Actual source code: aijsbaij.c

petsc-3.13.6 2020-09-29
Report Typos and Errors

  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>

  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;
 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: #if defined(PETSC_USE_COMPLEX)
 16:   const int      aconj = A->hermitian ? 1 : 0;
 17: #else
 18:   const int      aconj = 0;
 19: #endif

 22:   /* compute rowlengths of newmat */
 23:   PetscMalloc2(m,&rowlengths,m+1,&rowstart);

 25:   for (i=0; i<mbs; i++) rowlengths[i*bs] = 0;
 26:   k  = 0;
 27:   for (i=0; i<mbs; i++) {
 28:     nz = ai[i+1] - ai[i];
 29:     if (nz) {
 30:       rowlengths[k] += nz;   /* no. of upper triangular blocks */
 31:       if (*aj == i) {aj++;diagcnt++;nz--;} /* skip diagonal */
 32:       for (j=0; j<nz; j++) { /* no. of lower triangular blocks */
 33:         rowlengths[(*aj)*bs]++; aj++;
 34:       }
 35:     }
 36:     rowlengths[k] *= bs;
 37:     for (j=1; j<bs; j++) {
 38:       rowlengths[k+j] = rowlengths[k];
 39:     }
 40:     k += bs;
 41:   }

 43:   if (reuse != MAT_REUSE_MATRIX) {
 44:     MatCreate(PetscObjectComm((PetscObject)A),&B);
 45:     MatSetSizes(B,m,n,m,n);
 46:     MatSetType(B,MATSEQAIJ);
 47:     MatSeqAIJSetPreallocation(B,0,rowlengths);
 48:     MatSetBlockSize(B,A->rmap->bs);
 49:   } else B = *newmat;

 51:   b  = (Mat_SeqAIJ*)(B->data);
 52:   bi = b->i;
 53:   bj = b->j;
 54:   bv = b->a;

 56:   /* set b->i */
 57:   bi[0] = 0; rowstart[0] = 0;
 58:   for (i=0; i<mbs; i++) {
 59:     for (j=0; j<bs; j++) {
 60:       b->ilen[i*bs+j]    = rowlengths[i*bs];
 61:       rowstart[i*bs+j+1] = rowstart[i*bs+j] + rowlengths[i*bs];
 62:     }
 63:     bi[i+1] = bi[i] + rowlengths[i*bs]/bs;
 64:   }
 65:   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);

 67:   /* set b->j and b->a */
 68:   aj = a->j; av = a->a;
 69:   for (i=0; i<mbs; i++) {
 70:     nz = ai[i+1] - ai[i];
 71:     /* diagonal block */
 72:     if (nz && *aj == i) {
 73:       nz--;
 74:       for (j=0; j<bs; j++) {   /* row i*bs+j */
 75:         itmp = i*bs+j;
 76:         for (k=0; k<bs; k++) { /* col i*bs+k */
 77:           *(bj + rowstart[itmp]) = (*aj)*bs+k;
 78:           *(bv + rowstart[itmp]) = *(av+k*bs+j);
 79:           rowstart[itmp]++;
 80:         }
 81:       }
 82:       aj++; av += bs2;
 83:     }

 85:     while (nz--) {
 86:       /* lower triangular blocks */
 87:       for (j=0; j<bs; j++) {   /* row (*aj)*bs+j */
 88:         itmp = (*aj)*bs+j;
 89:         for (k=0; k<bs; k++) { /* col i*bs+k */
 90:           *(bj + rowstart[itmp]) = i*bs+k;
 91:           *(bv + rowstart[itmp]) = aconj ? PetscConj(*(av+j*bs+k)) : *(av+j*bs+k);
 92:           rowstart[itmp]++;
 93:         }
 94:       }
 95:       /* upper triangular blocks */
 96:       for (j=0; j<bs; j++) {   /* row i*bs+j */
 97:         itmp = i*bs+j;
 98:         for (k=0; k<bs; k++) { /* col (*aj)*bs+k */
 99:           *(bj + rowstart[itmp]) = (*aj)*bs+k;
100:           *(bv + rowstart[itmp]) = *(av+k*bs+j);
101:           rowstart[itmp]++;
102:         }
103:       }
104:       aj++; av += bs2;
105:     }
106:   }
107:   PetscFree2(rowlengths,rowstart);
108:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
109:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

111:   if (reuse == MAT_INPLACE_MATRIX) {
112:     MatHeaderReplace(A,&B);
113:   } else {
114:     *newmat = B;
115:   }
116:   return(0);
117: }

119: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
120: {
121:   Mat            B;
122:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
123:   Mat_SeqSBAIJ   *b;
125:   PetscInt       *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->N,i,j,*bi,*bj,*rowlengths,bs=PetscAbs(A->rmap->bs);
126:   MatScalar      *av,*bv;
127:   PetscBool      miss = PETSC_FALSE;

130: #if !defined(PETSC_USE_COMPLEX)
131:   if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
132: #else
133:   if (!A->symmetric && !A->hermitian) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be either symmetric or hermitian. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE) and/or MatSetOption(mat,MAT_HERMITIAN,PETSC_TRUE)");
134: #endif
135:   if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");

137:   PetscMalloc1(m/bs,&rowlengths);
138:   for (i=0; i<m/bs; i++) {
139:     if (a->diag[i*bs] == ai[i*bs+1]) { /* missing diagonal */
140:       rowlengths[i] = (ai[i*bs+1] - ai[i*bs])/bs; /* allocate some extra space */
141:       miss = PETSC_TRUE;
142:     } else {
143:       rowlengths[i] = (ai[i*bs+1] - a->diag[i*bs])/bs;
144:     }
145:   }
146:   if (reuse != MAT_REUSE_MATRIX) {
147:     MatCreate(PetscObjectComm((PetscObject)A),&B);
148:     MatSetSizes(B,m,n,m,n);
149:     MatSetType(B,MATSEQSBAIJ);
150:     MatSeqSBAIJSetPreallocation(B,bs,0,rowlengths);
151:   } else B = *newmat;

153:   if (bs == 1 && !miss) {
154:     b  = (Mat_SeqSBAIJ*)(B->data);
155:     bi = b->i;
156:     bj = b->j;
157:     bv = b->a;

159:     bi[0] = 0;
160:     for (i=0; i<m; i++) {
161:       aj = a->j + a->diag[i];
162:       av = a->a + a->diag[i];
163:       for (j=0; j<rowlengths[i]; j++) {
164:         *bj = *aj; bj++; aj++;
165:         *bv = *av; bv++; av++;
166:       }
167:       bi[i+1]    = bi[i] + rowlengths[i];
168:       b->ilen[i] = rowlengths[i];
169:     }
170:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
171:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
172:   } else {
173:     MatSetOption(B,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);
174:     /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
175:     /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before                                          */
176:     /*                      MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below                                        */
177:     MatConvert_Basic(A,newtype,MAT_REUSE_MATRIX,&B);
178:   }
179:   PetscFree(rowlengths);
180:   if (reuse == MAT_INPLACE_MATRIX) {
181:     MatHeaderReplace(A,&B);
182:   } else *newmat = B;
183:   return(0);
184: }

186: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
187: {
188:   Mat            B;
189:   Mat_SeqSBAIJ   *a = (Mat_SeqSBAIJ*)A->data;
190:   Mat_SeqBAIJ    *b;
192:   PetscInt       *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp;
193:   PetscInt       bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,col,row;
194:   MatScalar      *av,*bv;
195: #if defined(PETSC_USE_COMPLEX)
196:   const int      aconj = A->hermitian ? 1 : 0;
197: #else
198:   const int      aconj = 0;
199: #endif

202:   /* compute browlengths of newmat */
203:   PetscMalloc2(mbs,&browlengths,mbs,&browstart);
204:   for (i=0; i<mbs; i++) browlengths[i] = 0;
205:   for (i=0; i<mbs; i++) {
206:     nz = ai[i+1] - ai[i];
207:     aj++; /* skip diagonal */
208:     for (k=1; k<nz; k++) { /* no. of lower triangular blocks */
209:       browlengths[*aj]++; aj++;
210:     }
211:     browlengths[i] += nz;   /* no. of upper triangular blocks */
212:   }

214:   if (reuse != MAT_REUSE_MATRIX) {
215:     MatCreate(PetscObjectComm((PetscObject)A),&B);
216:     MatSetSizes(B,m,n,m,n);
217:     MatSetType(B,MATSEQBAIJ);
218:     MatSeqBAIJSetPreallocation(B,bs,0,browlengths);
219:   } else B = *newmat;

221:   b  = (Mat_SeqBAIJ*)(B->data);
222:   bi = b->i;
223:   bj = b->j;
224:   bv = b->a;

226:   /* set b->i */
227:   bi[0] = 0;
228:   for (i=0; i<mbs; i++) {
229:     b->ilen[i]   = browlengths[i];
230:     bi[i+1]      = bi[i] + browlengths[i];
231:     browstart[i] = bi[i];
232:   }
233:   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);

235:   /* set b->j and b->a */
236:   aj = a->j; av = a->a;
237:   for (i=0; i<mbs; i++) {
238:     /* diagonal block */
239:     *(bj + browstart[i]) = *aj; aj++;

241:     itmp = bs2*browstart[i];
242:     for (k=0; k<bs2; k++) {
243:       *(bv + itmp + k) = *av; av++;
244:     }
245:     browstart[i]++;

247:     nz = ai[i+1] - ai[i] -1;
248:     while (nz--) {
249:       /* lower triangular blocks - transpose blocks of A */
250:       *(bj + browstart[*aj]) = i; /* block col index */

252:       itmp = bs2*browstart[*aj];  /* row index */
253:       for (col=0; col<bs; col++) {
254:         k = col;
255:         for (row=0; row<bs; row++) {
256:           bv[itmp + col*bs+row] = aconj ? PetscConj(av[k]) : av[k];
257:           k+=bs;
258:         }
259:       }
260:       browstart[*aj]++;

262:       /* upper triangular blocks */
263:       *(bj + browstart[i]) = *aj; aj++;

265:       itmp = bs2*browstart[i];
266:       for (k=0; k<bs2; k++) {
267:         bv[itmp + k] = av[k];
268:       }
269:       av += bs2;
270:       browstart[i]++;
271:     }
272:   }
273:   PetscFree2(browlengths,browstart);
274:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
275:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

277:   if (reuse == MAT_INPLACE_MATRIX) {
278:     MatHeaderReplace(A,&B);
279:   } else *newmat = B;
280:   return(0);
281: }

283: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
284: {
285:   Mat            B;
286:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
287:   Mat_SeqSBAIJ   *b;
289:   PetscInt       *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*browlengths;
290:   PetscInt       bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,dd;
291:   MatScalar      *av,*bv;
292:   PetscBool      flg;

295:   if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
296:   if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
297:   MatMissingDiagonal_SeqBAIJ(A,&flg,&dd); /* check for missing diagonals, then mark diag */
298:   if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal %D",dd);

300:   PetscMalloc1(mbs,&browlengths);
301:   for (i=0; i<mbs; i++) {
302:     browlengths[i] = ai[i+1] - a->diag[i];
303:   }

305:   if (reuse != MAT_REUSE_MATRIX) {
306:     MatCreate(PetscObjectComm((PetscObject)A),&B);
307:     MatSetSizes(B,m,n,m,n);
308:     MatSetType(B,MATSEQSBAIJ);
309:     MatSeqSBAIJSetPreallocation(B,bs,0,browlengths);
310:   } else B = *newmat;

312:   b  = (Mat_SeqSBAIJ*)(B->data);
313:   bi = b->i;
314:   bj = b->j;
315:   bv = b->a;

317:   bi[0] = 0;
318:   for (i=0; i<mbs; i++) {
319:     aj = a->j + a->diag[i];
320:     av = a->a + (a->diag[i])*bs2;
321:     for (j=0; j<browlengths[i]; j++) {
322:       *bj = *aj; bj++; aj++;
323:       for (k=0; k<bs2; k++) {
324:         *bv = *av; bv++; av++;
325:       }
326:     }
327:     bi[i+1]    = bi[i] + browlengths[i];
328:     b->ilen[i] = browlengths[i];
329:   }
330:   PetscFree(browlengths);
331:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
332:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

334:   if (reuse == MAT_INPLACE_MATRIX) {
335:     MatHeaderReplace(A,&B);
336:   } else *newmat = B;
337:   return(0);
338: }