1 #define PETSCMAT_DLL 2 3 #include "../src/mat/impls/aij/seq/aij.h" 4 #include "../src/mat/impls/baij/seq/baij.h" 5 #include "../src/mat/impls/sbaij/seq/sbaij.h" 6 7 EXTERN_C_BEGIN 8 #undef __FUNCT__ 9 #define __FUNCT__ "MatConvert_SeqSBAIJ_SeqAIJ" 10 PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 11 { 12 Mat B; 13 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 14 Mat_SeqAIJ *b; 15 PetscErrorCode ierr; 16 PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*rowlengths,nz,*rowstart,itmp; 17 PetscInt bs=A->rmap->bs,bs2=bs*bs,mbs=A->rmap->N/bs; 18 MatScalar *av,*bv; 19 20 PetscFunctionBegin; 21 /* compute rowlengths of newmat */ 22 ierr = PetscMalloc2(m,PetscInt,&rowlengths,m+1,PetscInt,&rowstart);CHKERRQ(ierr); 23 24 for (i=0; i<mbs; i++) rowlengths[i*bs] = 0; 25 aj = a->j; 26 k = 0; 27 for (i=0; i<mbs; i++) { 28 nz = ai[i+1] - ai[i]; 29 aj++; /* skip diagonal */ 30 for (j=1; j<nz; j++) { /* no. of lower triangular blocks */ 31 rowlengths[(*aj)*bs]++; aj++; 32 } 33 rowlengths[k] += nz; /* no. of upper triangular blocks */ 34 rowlengths[k] *= bs; 35 for (j=1; j<bs; j++) { 36 rowlengths[k+j] = rowlengths[k]; 37 } 38 k += bs; 39 /* printf(" rowlengths[%d]: %d\n",i, rowlengths[i]); */ 40 } 41 42 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 43 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 44 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 45 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 46 ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_FALSE);CHKERRQ(ierr); 47 B->rmap->bs = A->rmap->bs; 48 49 b = (Mat_SeqAIJ*)(B->data); 50 bi = b->i; 51 bj = b->j; 52 bv = b->a; 53 54 /* set b->i */ 55 bi[0] = 0; rowstart[0] = 0; 56 for (i=0; i<mbs; i++){ 57 for (j=0; j<bs; j++){ 58 b->ilen[i*bs+j] = rowlengths[i*bs]; 59 rowstart[i*bs+j+1] = rowstart[i*bs+j] + rowlengths[i*bs]; 60 } 61 bi[i+1] = bi[i] + rowlengths[i*bs]/bs; 62 } 63 if (bi[mbs] != 2*a->nz - mbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bi[mbs]: %D != 2*a->nz-mbs: %D\n",bi[mbs],2*a->nz - mbs); 64 65 /* set b->j and b->a */ 66 aj = a->j; av = a->a; 67 for (i=0; i<mbs; i++) { 68 /* diagonal block */ 69 for (j=0; j<bs; j++){ /* row i*bs+j */ 70 itmp = i*bs+j; 71 for (k=0; k<bs; k++){ /* col i*bs+k */ 72 *(bj + rowstart[itmp]) = (*aj)*bs+k; 73 *(bv + rowstart[itmp]) = *(av+k*bs+j); 74 rowstart[itmp]++; 75 } 76 } 77 aj++; av += bs2; 78 79 nz = ai[i+1] - ai[i] -1; 80 while (nz--){ 81 /* lower triangular blocks */ 82 for (j=0; j<bs; j++){ /* row (*aj)*bs+j */ 83 itmp = (*aj)*bs+j; 84 for (k=0; k<bs; k++){ /* col i*bs+k */ 85 *(bj + rowstart[itmp]) = i*bs+k; 86 *(bv + rowstart[itmp]) = *(av+k*bs+j); 87 rowstart[itmp]++; 88 } 89 } 90 /* upper triangular blocks */ 91 for (j=0; j<bs; j++){ /* row i*bs+j */ 92 itmp = i*bs+j; 93 for (k=0; k<bs; k++){ /* col (*aj)*bs+k */ 94 *(bj + rowstart[itmp]) = (*aj)*bs+k; 95 *(bv + rowstart[itmp]) = *(av+k*bs+j); 96 rowstart[itmp]++; 97 } 98 } 99 aj++; av += bs2; 100 } 101 } 102 ierr = PetscFree2(rowlengths,rowstart);CHKERRQ(ierr); 103 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 104 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 105 106 if (reuse == MAT_REUSE_MATRIX) { 107 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 108 } else { 109 *newmat = B; 110 } 111 PetscFunctionReturn(0); 112 } 113 EXTERN_C_END 114 115 EXTERN_C_BEGIN 116 #undef __FUNCT__ 117 #define __FUNCT__ "MatConvert_SeqAIJ_SeqSBAIJ" 118 PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqSBAIJ(Mat A,const MatType newtype,MatReuse reuse,Mat *newmat) { 119 Mat B; 120 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 121 Mat_SeqSBAIJ *b; 122 PetscErrorCode ierr; 123 PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->N,i,j,*bi,*bj,*rowlengths; 124 MatScalar *av,*bv; 125 126 PetscFunctionBegin; 127 if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square"); 128 129 ierr = PetscMalloc(m*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 130 for (i=0; i<m; i++) { 131 rowlengths[i] = ai[i+1] - a->diag[i]; 132 } 133 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 134 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 135 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 136 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,1,0,rowlengths);CHKERRQ(ierr); 137 138 ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 139 140 b = (Mat_SeqSBAIJ*)(B->data); 141 bi = b->i; 142 bj = b->j; 143 bv = b->a; 144 145 bi[0] = 0; 146 for (i=0; i<m; i++) { 147 aj = a->j + a->diag[i]; 148 av = a->a + a->diag[i]; 149 for (j=0; j<rowlengths[i]; j++){ 150 *bj = *aj; bj++; aj++; 151 *bv = *av; bv++; av++; 152 } 153 bi[i+1] = bi[i] + rowlengths[i]; 154 b->ilen[i] = rowlengths[i]; 155 } 156 157 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 158 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 159 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 160 if (A->hermitian){ 161 ierr = MatSetOption(B,MAT_HERMITIAN,PETSC_TRUE);CHKERRQ(ierr); 162 } 163 164 if (reuse == MAT_REUSE_MATRIX) { 165 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 166 } else { 167 *newmat = B; 168 } 169 PetscFunctionReturn(0); 170 } 171 EXTERN_C_END 172 173 EXTERN_C_BEGIN 174 #undef __FUNCT__ 175 #define __FUNCT__ "MatConvert_SeqSBAI_SeqBAIJ" 176 PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 177 { 178 Mat B; 179 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 180 Mat_SeqBAIJ *b; 181 PetscErrorCode ierr; 182 PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp; 183 PetscInt bs=A->rmap->bs,bs2=bs*bs,mbs=m/bs; 184 MatScalar *av,*bv; 185 186 PetscFunctionBegin; 187 /* compute browlengths of newmat */ 188 ierr = PetscMalloc2(mbs,PetscInt,&browlengths,mbs,PetscInt,&browstart);CHKERRQ(ierr); 189 for (i=0; i<mbs; i++) browlengths[i] = 0; 190 aj = a->j; 191 for (i=0; i<mbs; i++) { 192 nz = ai[i+1] - ai[i]; 193 aj++; /* skip diagonal */ 194 for (k=1; k<nz; k++) { /* no. of lower triangular blocks */ 195 browlengths[*aj]++; aj++; 196 } 197 browlengths[i] += nz; /* no. of upper triangular blocks */ 198 } 199 200 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 201 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 202 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 203 ierr = MatSeqBAIJSetPreallocation(B,bs,0,browlengths);CHKERRQ(ierr); 204 ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 205 206 b = (Mat_SeqBAIJ*)(B->data); 207 bi = b->i; 208 bj = b->j; 209 bv = b->a; 210 211 /* set b->i */ 212 bi[0] = 0; 213 for (i=0; i<mbs; i++){ 214 b->ilen[i] = browlengths[i]; 215 bi[i+1] = bi[i] + browlengths[i]; 216 browstart[i] = bi[i]; 217 } 218 if (bi[mbs] != 2*a->nz - mbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bi[mbs]: %D != 2*a->nz - mbs: %D\n",bi[mbs],2*a->nz - mbs); 219 220 /* set b->j and b->a */ 221 aj = a->j; av = a->a; 222 for (i=0; i<mbs; i++) { 223 /* diagonal block */ 224 *(bj + browstart[i]) = *aj; aj++; 225 itmp = bs2*browstart[i]; 226 for (k=0; k<bs2; k++){ 227 *(bv + itmp + k) = *av; av++; 228 } 229 browstart[i]++; 230 231 nz = ai[i+1] - ai[i] -1; 232 while (nz--){ 233 /* lower triangular blocks */ 234 *(bj + browstart[*aj]) = i; 235 itmp = bs2*browstart[*aj]; 236 for (k=0; k<bs2; k++){ 237 *(bv + itmp + k) = *(av + k); 238 } 239 browstart[*aj]++; 240 241 /* upper triangular blocks */ 242 *(bj + browstart[i]) = *aj; aj++; 243 itmp = bs2*browstart[i]; 244 for (k=0; k<bs2; k++){ 245 *(bv + itmp + k) = *av; av++; 246 } 247 browstart[i]++; 248 } 249 } 250 ierr = PetscFree2(browlengths,browstart);CHKERRQ(ierr); 251 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 252 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 253 254 if (reuse == MAT_REUSE_MATRIX) { 255 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 256 } else { 257 *newmat = B; 258 } 259 PetscFunctionReturn(0); 260 } 261 EXTERN_C_END 262 263 EXTERN_C_BEGIN 264 #undef __FUNCT__ 265 #define __FUNCT__ "MatConvert_SeqBAIJ_SeqSBAIJ" 266 PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 267 { 268 Mat B; 269 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 270 Mat_SeqSBAIJ *b; 271 PetscErrorCode ierr; 272 PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*browlengths; 273 PetscInt bs=A->rmap->bs,bs2=bs*bs,mbs=m/bs,dd; 274 MatScalar *av,*bv; 275 PetscTruth flg; 276 277 PetscFunctionBegin; 278 if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square"); 279 ierr = MatMissingDiagonal_SeqBAIJ(A,&flg,&dd);CHKERRQ(ierr); /* check for missing diagonals, then mark diag */ 280 if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal %D",dd); 281 282 ierr = PetscMalloc(mbs*sizeof(PetscInt),&browlengths);CHKERRQ(ierr); 283 for (i=0; i<mbs; i++) { 284 browlengths[i] = ai[i+1] - a->diag[i]; 285 } 286 287 ierr = MatCreate(((PetscObject)A)->comm,&B);CHKERRQ(ierr); 288 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 289 ierr = MatSetType(B,newtype);CHKERRQ(ierr); 290 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,browlengths);CHKERRQ(ierr); 291 ierr = MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);CHKERRQ(ierr); 292 293 b = (Mat_SeqSBAIJ*)(B->data); 294 bi = b->i; 295 bj = b->j; 296 bv = b->a; 297 298 bi[0] = 0; 299 for (i=0; i<mbs; i++) { 300 aj = a->j + a->diag[i]; 301 av = a->a + (a->diag[i])*bs2; 302 for (j=0; j<browlengths[i]; j++){ 303 *bj = *aj; bj++; aj++; 304 for (k=0; k<bs2; k++){ 305 *bv = *av; bv++; av++; 306 } 307 } 308 bi[i+1] = bi[i] + browlengths[i]; 309 b->ilen[i] = browlengths[i]; 310 } 311 ierr = PetscFree(browlengths);CHKERRQ(ierr); 312 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 313 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 314 315 if (reuse == MAT_REUSE_MATRIX) { 316 ierr = MatHeaderReplace(A,B);CHKERRQ(ierr); 317 } else { 318 *newmat = B; 319 } 320 321 PetscFunctionReturn(0); 322 } 323 EXTERN_C_END 324