Actual source code: aijsbaij.c


  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;
 11:   PetscInt     *ai = a->i, *aj = a->j, m = A->rmap->N, n = A->cmap->n, i, j, k, *bi, *bj, *rowlengths, nz, *rowstart, itmp;
 12:   PetscInt      bs = A->rmap->bs, bs2 = bs * bs, mbs = A->rmap->N / bs, diagcnt = 0;
 13:   MatScalar    *av, *bv;
 14: #if defined(PETSC_USE_COMPLEX)
 15:   const int aconj = A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0;
 16: #else
 17:   const int aconj = 0;
 18: #endif

 20:   PetscFunctionBegin;
 21:   /* compute rowlengths of newmat */
 22:   PetscCall(PetscMalloc2(m, &rowlengths, m + 1, &rowstart));

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

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

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

 58:   /* set b->i */
 59:   bi[0]       = 0;
 60:   rowstart[0] = 0;
 61:   for (i = 0; i < mbs; i++) {
 62:     for (j = 0; j < bs; j++) {
 63:       b->ilen[i * bs + j]      = rowlengths[i * bs];
 64:       rowstart[i * bs + j + 1] = rowstart[i * bs + j] + rowlengths[i * bs];
 65:     }
 66:     bi[i + 1] = bi[i] + rowlengths[i * bs] / bs;
 67:   }
 68:   PetscCheck(bi[mbs] == 2 * a->nz - diagcnt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bi[mbs]: %" PetscInt_FMT " != 2*a->nz-diagcnt: %" PetscInt_FMT, bi[mbs], 2 * a->nz - diagcnt);

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

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

117:   if (reuse == MAT_INPLACE_MATRIX) {
118:     PetscCall(MatHeaderReplace(A, &B));
119:   } else {
120:     *newmat = B;
121:   }
122:   PetscFunctionReturn(PETSC_SUCCESS);
123: }

125: #include <../src/mat/impls/aij/seq/aij.h>

127: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ_Preallocate(Mat A, PetscInt **nnz)
128: {
129:   Mat_SeqAIJ     *Aa = (Mat_SeqAIJ *)A->data;
130:   PetscInt        m, n, bs = PetscAbs(A->rmap->bs);
131:   const PetscInt *ai = Aa->i, *aj = Aa->j;

133:   PetscFunctionBegin;
134:   PetscCall(MatGetSize(A, &m, &n));

136:   if (bs == 1) {
137:     const PetscInt *adiag = Aa->diag;

139:     PetscCall(PetscMalloc1(m, nnz));
140:     for (PetscInt i = 0; i < m; i++) {
141:       if (adiag[i] == ai[i + 1]) {
142:         (*nnz)[i] = 0;
143:         for (PetscInt j = ai[i]; j < ai[i + 1]; j++) (*nnz)[i] += (aj[j] > i);
144:       } else (*nnz)[i] = ai[i + 1] - adiag[i];
145:     }
146:   } else {
147:     PetscHSetIJ    ht;
148:     PetscHashIJKey key;
149:     PetscBool      missing;

151:     PetscCall(PetscHSetIJCreate(&ht));
152:     PetscCall(PetscCalloc1(m / bs, nnz));
153:     for (PetscInt i = 0; i < m; i++) {
154:       key.i = i / bs;
155:       for (PetscInt k = ai[i]; k < ai[i + 1]; k++) {
156:         key.j = aj[k] / bs;
157:         if (key.j < key.i) continue;
158:         PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
159:         if (missing) (*nnz)[key.i]++;
160:       }
161:     }
162:     PetscCall(PetscHSetIJDestroy(&ht));
163:   }
164:   PetscFunctionReturn(PETSC_SUCCESS);
165: }

167: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
168: {
169:   Mat           B;
170:   Mat_SeqAIJ   *a = (Mat_SeqAIJ *)A->data;
171:   Mat_SeqSBAIJ *b;
172:   PetscInt     *ai = a->i, *aj, m = A->rmap->N, n = A->cmap->N, i, j, *bi, *bj, *rowlengths, bs = PetscAbs(A->rmap->bs);
173:   MatScalar    *av, *bv;
174:   PetscBool     miss = PETSC_FALSE;

176:   PetscFunctionBegin;
177: #if !defined(PETSC_USE_COMPLEX)
178:   PetscCheck(A->symmetric == PETSC_BOOL3_TRUE, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
179: #else
180:   PetscCheck(A->symmetric == PETSC_BOOL3_TRUE || A->hermitian == PETSC_BOOL3_TRUE, 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)");
181: #endif
182:   PetscCheck(n == m, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix must be square");

184:   if (bs == 1) {
185:     PetscCall(PetscMalloc1(m, &rowlengths));
186:     for (i = 0; i < m; i++) {
187:       if (a->diag[i] == ai[i + 1]) {             /* missing diagonal */
188:         rowlengths[i] = (ai[i + 1] - ai[i]) + 1; /* allocate some extra space */
189:         miss          = PETSC_TRUE;
190:       } else {
191:         rowlengths[i] = (ai[i + 1] - a->diag[i]);
192:       }
193:     }
194:   } else PetscCall(MatConvert_SeqAIJ_SeqSBAIJ_Preallocate(A, &rowlengths));
195:   if (reuse != MAT_REUSE_MATRIX) {
196:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
197:     PetscCall(MatSetSizes(B, m, n, m, n));
198:     PetscCall(MatSetType(B, MATSEQSBAIJ));
199:     PetscCall(MatSeqSBAIJSetPreallocation(B, bs, 0, rowlengths));
200:   } else B = *newmat;

202:   if (bs == 1 && !miss) {
203:     b  = (Mat_SeqSBAIJ *)(B->data);
204:     bi = b->i;
205:     bj = b->j;
206:     bv = b->a;

208:     bi[0] = 0;
209:     for (i = 0; i < m; i++) {
210:       aj = a->j + a->diag[i];
211:       av = a->a + a->diag[i];
212:       for (j = 0; j < rowlengths[i]; j++) {
213:         *bj = *aj;
214:         bj++;
215:         aj++;
216:         *bv = *av;
217:         bv++;
218:         av++;
219:       }
220:       bi[i + 1]  = bi[i] + rowlengths[i];
221:       b->ilen[i] = rowlengths[i];
222:     }
223:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
224:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
225:   } else {
226:     PetscCall(MatSetOption(B, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE));
227:     /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
228:     /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before                                          */
229:     /*                      MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below                                        */
230:     PetscCall(MatConvert_Basic(A, newtype, MAT_REUSE_MATRIX, &B));
231:   }
232:   PetscCall(PetscFree(rowlengths));
233:   if (reuse == MAT_INPLACE_MATRIX) {
234:     PetscCall(MatHeaderReplace(A, &B));
235:   } else *newmat = B;
236:   PetscFunctionReturn(PETSC_SUCCESS);
237: }

239: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
240: {
241:   Mat           B;
242:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
243:   Mat_SeqBAIJ  *b;
244:   PetscInt     *ai = a->i, *aj = a->j, m = A->rmap->N, n = A->cmap->n, i, k, *bi, *bj, *browlengths, nz, *browstart, itmp;
245:   PetscInt      bs = A->rmap->bs, bs2 = bs * bs, mbs = m / bs, col, row;
246:   MatScalar    *av, *bv;
247: #if defined(PETSC_USE_COMPLEX)
248:   const int aconj = A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0;
249: #else
250:   const int aconj = 0;
251: #endif

253:   PetscFunctionBegin;
254:   /* compute browlengths of newmat */
255:   PetscCall(PetscMalloc2(mbs, &browlengths, mbs, &browstart));
256:   for (i = 0; i < mbs; i++) browlengths[i] = 0;
257:   for (i = 0; i < mbs; i++) {
258:     nz = ai[i + 1] - ai[i];
259:     aj++;                      /* skip diagonal */
260:     for (k = 1; k < nz; k++) { /* no. of lower triangular blocks */
261:       browlengths[*aj]++;
262:       aj++;
263:     }
264:     browlengths[i] += nz; /* no. of upper triangular blocks */
265:   }

267:   if (reuse != MAT_REUSE_MATRIX) {
268:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
269:     PetscCall(MatSetSizes(B, m, n, m, n));
270:     PetscCall(MatSetType(B, MATSEQBAIJ));
271:     PetscCall(MatSeqBAIJSetPreallocation(B, bs, 0, browlengths));
272:   } else B = *newmat;

274:   b  = (Mat_SeqBAIJ *)(B->data);
275:   bi = b->i;
276:   bj = b->j;
277:   bv = b->a;

279:   /* set b->i */
280:   bi[0] = 0;
281:   for (i = 0; i < mbs; i++) {
282:     b->ilen[i]   = browlengths[i];
283:     bi[i + 1]    = bi[i] + browlengths[i];
284:     browstart[i] = bi[i];
285:   }
286:   PetscCheck(bi[mbs] == 2 * a->nz - mbs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bi[mbs]: %" PetscInt_FMT " != 2*a->nz - mbs: %" PetscInt_FMT, bi[mbs], 2 * a->nz - mbs);

288:   /* set b->j and b->a */
289:   aj = a->j;
290:   av = a->a;
291:   for (i = 0; i < mbs; i++) {
292:     /* diagonal block */
293:     *(bj + browstart[i]) = *aj;
294:     aj++;

296:     itmp = bs2 * browstart[i];
297:     for (k = 0; k < bs2; k++) {
298:       *(bv + itmp + k) = *av;
299:       av++;
300:     }
301:     browstart[i]++;

303:     nz = ai[i + 1] - ai[i] - 1;
304:     while (nz--) {
305:       /* lower triangular blocks - transpose blocks of A */
306:       *(bj + browstart[*aj]) = i; /* block col index */

308:       itmp = bs2 * browstart[*aj]; /* row index */
309:       for (col = 0; col < bs; col++) {
310:         k = col;
311:         for (row = 0; row < bs; row++) {
312:           bv[itmp + col * bs + row] = aconj ? PetscConj(av[k]) : av[k];
313:           k += bs;
314:         }
315:       }
316:       browstart[*aj]++;

318:       /* upper triangular blocks */
319:       *(bj + browstart[i]) = *aj;
320:       aj++;

322:       itmp = bs2 * browstart[i];
323:       for (k = 0; k < bs2; k++) bv[itmp + k] = av[k];
324:       av += bs2;
325:       browstart[i]++;
326:     }
327:   }
328:   PetscCall(PetscFree2(browlengths, browstart));
329:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
330:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

332:   if (reuse == MAT_INPLACE_MATRIX) {
333:     PetscCall(MatHeaderReplace(A, &B));
334:   } else *newmat = B;
335:   PetscFunctionReturn(PETSC_SUCCESS);
336: }

338: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
339: {
340:   Mat           B;
341:   Mat_SeqBAIJ  *a = (Mat_SeqBAIJ *)A->data;
342:   Mat_SeqSBAIJ *b;
343:   PetscInt     *ai = a->i, *aj, m = A->rmap->N, n = A->cmap->n, i, j, k, *bi, *bj, *browlengths;
344:   PetscInt      bs = A->rmap->bs, bs2 = bs * bs, mbs = m / bs, dd;
345:   MatScalar    *av, *bv;
346:   PetscBool     flg;

348:   PetscFunctionBegin;
349:   PetscCheck(A->symmetric, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
350:   PetscCheck(n == m, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix must be square");
351:   PetscCall(MatMissingDiagonal_SeqBAIJ(A, &flg, &dd)); /* check for missing diagonals, then mark diag */
352:   PetscCheck(!flg, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal %" PetscInt_FMT, dd);

354:   PetscCall(PetscMalloc1(mbs, &browlengths));
355:   for (i = 0; i < mbs; i++) browlengths[i] = ai[i + 1] - a->diag[i];

357:   if (reuse != MAT_REUSE_MATRIX) {
358:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
359:     PetscCall(MatSetSizes(B, m, n, m, n));
360:     PetscCall(MatSetType(B, MATSEQSBAIJ));
361:     PetscCall(MatSeqSBAIJSetPreallocation(B, bs, 0, browlengths));
362:   } else B = *newmat;

364:   b  = (Mat_SeqSBAIJ *)(B->data);
365:   bi = b->i;
366:   bj = b->j;
367:   bv = b->a;

369:   bi[0] = 0;
370:   for (i = 0; i < mbs; i++) {
371:     aj = a->j + a->diag[i];
372:     av = a->a + (a->diag[i]) * bs2;
373:     for (j = 0; j < browlengths[i]; j++) {
374:       *bj = *aj;
375:       bj++;
376:       aj++;
377:       for (k = 0; k < bs2; k++) {
378:         *bv = *av;
379:         bv++;
380:         av++;
381:       }
382:     }
383:     bi[i + 1]  = bi[i] + browlengths[i];
384:     b->ilen[i] = browlengths[i];
385:   }
386:   PetscCall(PetscFree(browlengths));
387:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
388:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

390:   if (reuse == MAT_INPLACE_MATRIX) {
391:     PetscCall(MatHeaderReplace(A, &B));
392:   } else *newmat = B;
393:   PetscFunctionReturn(PETSC_SUCCESS);
394: }