1 2 #include <../src/mat/impls/aij/seq/aij.h> 3 #include <../src/mat/impls/baij/seq/baij.h> 4 #include <../src/mat/impls/sbaij/seq/sbaij.h> 5 6 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 7 { 8 Mat B; 9 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 10 Mat_SeqAIJ *b; 11 PetscErrorCode ierr; 12 PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*rowlengths,nz,*rowstart,itmp; 13 PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=A->rmap->N/bs,diagcnt=0; 14 MatScalar *av,*bv; 15 16 PetscFunctionBegin; 17 /* compute rowlengths of newmat */ 18 ierr = PetscMalloc2(m,&rowlengths,m+1,&rowstart);CHKERRQ(ierr); 19 20 for (i=0; i<mbs; i++) rowlengths[i*bs] = 0; 21 k = 0; 22 for (i=0; i<mbs; i++) { 23 nz = ai[i+1] - ai[i]; 24 if (nz) { 25 rowlengths[k] += nz; /* no. of upper triangular blocks */ 26 if (*aj == i) {aj++;diagcnt++;nz--;} /* skip diagonal */ 27 for (j=0; j<nz; j++) { /* no. of lower triangular blocks */ 28 rowlengths[(*aj)*bs]++; aj++; 29 } 30 } 31 rowlengths[k] *= bs; 32 for (j=1; j<bs; j++) { 33 rowlengths[k+j] = rowlengths[k]; 34 } 35 k += bs; 36 } 37 38 if (reuse != MAT_REUSE_MATRIX) { 39 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 40 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 41 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 42 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 43 ierr = MatSetBlockSize(B,A->rmap->bs);CHKERRQ(ierr); 44 } else B = *newmat; 45 46 b = (Mat_SeqAIJ*)(B->data); 47 bi = b->i; 48 bj = b->j; 49 bv = b->a; 50 51 /* set b->i */ 52 bi[0] = 0; rowstart[0] = 0; 53 for (i=0; i<mbs; i++) { 54 for (j=0; j<bs; j++) { 55 b->ilen[i*bs+j] = rowlengths[i*bs]; 56 rowstart[i*bs+j+1] = rowstart[i*bs+j] + rowlengths[i*bs]; 57 } 58 bi[i+1] = bi[i] + rowlengths[i*bs]/bs; 59 } 60 if (bi[mbs] != 2*a->nz - diagcnt) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bi[mbs]: %D != 2*a->nz-diagcnt: %D\n",bi[mbs],2*a->nz - diagcnt); 61 62 /* set b->j and b->a */ 63 aj = a->j; av = a->a; 64 for (i=0; i<mbs; i++) { 65 nz = ai[i+1] - ai[i]; 66 /* diagonal block */ 67 if (nz && *aj == i) { 68 nz--; 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 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+j*bs+k); 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_INPLACE_MATRIX) { 107 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 108 } else { 109 *newmat = B; 110 } 111 PetscFunctionReturn(0); 112 } 113 114 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat) 115 { 116 Mat B; 117 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 118 Mat_SeqSBAIJ *b; 119 PetscErrorCode ierr; 120 PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->N,i,j,*bi,*bj,*rowlengths,bs=PetscAbs(A->rmap->bs); 121 MatScalar *av,*bv; 122 123 PetscFunctionBegin; 124 if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)"); 125 if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square"); 126 127 ierr = PetscMalloc1(m/bs,&rowlengths);CHKERRQ(ierr); 128 for (i=0; i<m/bs; i++) { 129 rowlengths[i] = (ai[i*bs+1] - a->diag[i*bs])/bs; 130 } 131 if (reuse != MAT_REUSE_MATRIX) { 132 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 133 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 134 ierr = MatSetType(B,MATSEQSBAIJ);CHKERRQ(ierr); 135 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,rowlengths);CHKERRQ(ierr); 136 } else B = *newmat; 137 138 if (bs == 1) { 139 b = (Mat_SeqSBAIJ*)(B->data); 140 bi = b->i; 141 bj = b->j; 142 bv = b->a; 143 144 bi[0] = 0; 145 for (i=0; i<m; i++) { 146 aj = a->j + a->diag[i]; 147 av = a->a + a->diag[i]; 148 for (j=0; j<rowlengths[i]; j++) { 149 *bj = *aj; bj++; aj++; 150 *bv = *av; bv++; av++; 151 } 152 bi[i+1] = bi[i] + rowlengths[i]; 153 b->ilen[i] = rowlengths[i]; 154 } 155 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 156 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 157 } else { 158 ierr = MatSetOption(B,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);CHKERRQ(ierr); 159 /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */ 160 /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before */ 161 /* MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below */ 162 ierr = MatConvert_Basic(A,newtype,MAT_REUSE_MATRIX,&B);CHKERRQ(ierr); 163 } 164 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 165 if (reuse == MAT_INPLACE_MATRIX) { 166 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 167 } else *newmat = B; 168 PetscFunctionReturn(0); 169 } 170 171 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 172 { 173 Mat B; 174 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 175 Mat_SeqBAIJ *b; 176 PetscErrorCode ierr; 177 PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp; 178 PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,col,row; 179 MatScalar *av,*bv; 180 181 PetscFunctionBegin; 182 /* compute browlengths of newmat */ 183 ierr = PetscMalloc2(mbs,&browlengths,mbs,&browstart);CHKERRQ(ierr); 184 for (i=0; i<mbs; i++) browlengths[i] = 0; 185 for (i=0; i<mbs; i++) { 186 nz = ai[i+1] - ai[i]; 187 aj++; /* skip diagonal */ 188 for (k=1; k<nz; k++) { /* no. of lower triangular blocks */ 189 browlengths[*aj]++; aj++; 190 } 191 browlengths[i] += nz; /* no. of upper triangular blocks */ 192 } 193 194 if (reuse != MAT_REUSE_MATRIX) { 195 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 196 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 197 ierr = MatSetType(B,MATSEQBAIJ);CHKERRQ(ierr); 198 ierr = MatSeqBAIJSetPreallocation(B,bs,0,browlengths);CHKERRQ(ierr); 199 } else B = *newmat; 200 201 b = (Mat_SeqBAIJ*)(B->data); 202 bi = b->i; 203 bj = b->j; 204 bv = b->a; 205 206 /* set b->i */ 207 bi[0] = 0; 208 for (i=0; i<mbs; i++) { 209 b->ilen[i] = browlengths[i]; 210 bi[i+1] = bi[i] + browlengths[i]; 211 browstart[i] = bi[i]; 212 } 213 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); 214 215 /* set b->j and b->a */ 216 aj = a->j; av = a->a; 217 for (i=0; i<mbs; i++) { 218 /* diagonal block */ 219 *(bj + browstart[i]) = *aj; aj++; 220 221 itmp = bs2*browstart[i]; 222 for (k=0; k<bs2; k++) { 223 *(bv + itmp + k) = *av; av++; 224 } 225 browstart[i]++; 226 227 nz = ai[i+1] - ai[i] -1; 228 while (nz--) { 229 /* lower triangular blocks - transpose blocks of A */ 230 *(bj + browstart[*aj]) = i; /* block col index */ 231 232 itmp = bs2*browstart[*aj]; /* row index */ 233 for (col=0; col<bs; col++) { 234 k = col; 235 for (row=0; row<bs; row++) { 236 bv[itmp + col*bs+row] = av[k]; k+=bs; 237 } 238 } 239 browstart[*aj]++; 240 241 /* upper triangular blocks */ 242 *(bj + browstart[i]) = *aj; aj++; 243 244 itmp = bs2*browstart[i]; 245 for (k=0; k<bs2; k++) { 246 bv[itmp + k] = av[k]; 247 } 248 av += bs2; 249 browstart[i]++; 250 } 251 } 252 ierr = PetscFree2(browlengths,browstart);CHKERRQ(ierr); 253 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 254 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 255 256 if (reuse == MAT_INPLACE_MATRIX) { 257 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 258 } else *newmat = B; 259 PetscFunctionReturn(0); 260 } 261 262 PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat) 263 { 264 Mat B; 265 Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data; 266 Mat_SeqSBAIJ *b; 267 PetscErrorCode ierr; 268 PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*browlengths; 269 PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,dd; 270 MatScalar *av,*bv; 271 PetscBool flg; 272 273 PetscFunctionBegin; 274 if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)"); 275 if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square"); 276 ierr = MatMissingDiagonal_SeqBAIJ(A,&flg,&dd);CHKERRQ(ierr); /* check for missing diagonals, then mark diag */ 277 if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal %D",dd); 278 279 ierr = PetscMalloc1(mbs,&browlengths);CHKERRQ(ierr); 280 for (i=0; i<mbs; i++) { 281 browlengths[i] = ai[i+1] - a->diag[i]; 282 } 283 284 if (reuse != MAT_REUSE_MATRIX) { 285 ierr = MatCreate(PetscObjectComm((PetscObject)A),&B);CHKERRQ(ierr); 286 ierr = MatSetSizes(B,m,n,m,n);CHKERRQ(ierr); 287 ierr = MatSetType(B,MATSEQSBAIJ);CHKERRQ(ierr); 288 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,browlengths);CHKERRQ(ierr); 289 } else B = *newmat; 290 291 b = (Mat_SeqSBAIJ*)(B->data); 292 bi = b->i; 293 bj = b->j; 294 bv = b->a; 295 296 bi[0] = 0; 297 for (i=0; i<mbs; i++) { 298 aj = a->j + a->diag[i]; 299 av = a->a + (a->diag[i])*bs2; 300 for (j=0; j<browlengths[i]; j++) { 301 *bj = *aj; bj++; aj++; 302 for (k=0; k<bs2; k++) { 303 *bv = *av; bv++; av++; 304 } 305 } 306 bi[i+1] = bi[i] + browlengths[i]; 307 b->ilen[i] = browlengths[i]; 308 } 309 ierr = PetscFree(browlengths);CHKERRQ(ierr); 310 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 311 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 312 313 if (reuse == MAT_INPLACE_MATRIX) { 314 ierr = MatHeaderReplace(A,&B);CHKERRQ(ierr); 315 } else *newmat = B; 316 PetscFunctionReturn(0); 317 } 318