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:   /* compute rowlengths of newmat */
 21:   PetscMalloc2(m, &rowlengths, m + 1, &rowstart);

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

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

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

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

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

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

116:   if (reuse == MAT_INPLACE_MATRIX) {
117:     MatHeaderReplace(A, &B);
118:   } else {
119:     *newmat = B;
120:   }
121:   return 0;
122: }

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

126: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ_Preallocate(Mat A, PetscInt **nnz)
127: {
128:   PetscInt    m, n, bs = PetscAbs(A->rmap->bs);
129:   PetscInt   *ii;
130:   Mat_SeqAIJ *Aa = (Mat_SeqAIJ *)A->data;

132:   MatGetSize(A, &m, &n);
133:   PetscCalloc1(m / bs, nnz);
134:   PetscMalloc1(bs, &ii);

136:   /* loop over all potential blocks in the matrix to see if the potential block has a nonzero */
137:   for (PetscInt i = 0; i < m / bs; i++) {
138:     for (PetscInt k = 0; k < bs; k++) ii[k] = Aa->i[i * bs + k];
139:     for (PetscInt j = 0; j < n / bs; j++) {
140:       for (PetscInt k = 0; k < bs; k++) {
141:         for (; ii[k] < Aa->i[i * bs + k + 1] && Aa->j[ii[k]] < (j + 1) * bs; ii[k]++) {
142:           if (j >= i && Aa->j[ii[k]] >= j * bs) {
143:             /* found a nonzero in the potential block so allocate for that block and jump to check the next potential block */
144:             (*nnz)[i]++;
145:             goto theend;
146:           }
147:         }
148:       }
149:     theend:;
150:     }
151:   }
152:   PetscFree(ii);
153:   return 0;
154: }

156: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
157: {
158:   Mat           B;
159:   Mat_SeqAIJ   *a = (Mat_SeqAIJ *)A->data;
160:   Mat_SeqSBAIJ *b;
161:   PetscInt     *ai = a->i, *aj, m = A->rmap->N, n = A->cmap->N, i, j, *bi, *bj, *rowlengths, bs = PetscAbs(A->rmap->bs);
162:   MatScalar    *av, *bv;
163:   PetscBool     miss = PETSC_FALSE;

165: #if !defined(PETSC_USE_COMPLEX)
167: #else
169: #endif

172:   if (bs == 1) {
173:     PetscMalloc1(m, &rowlengths);
174:     for (i = 0; i < m; i++) {
175:       if (a->diag[i] == ai[i + 1]) {             /* missing diagonal */
176:         rowlengths[i] = (ai[i + 1] - ai[i]) + 1; /* allocate some extra space */
177:         miss          = PETSC_TRUE;
178:       } else {
179:         rowlengths[i] = (ai[i + 1] - a->diag[i]);
180:       }
181:     }
182:   } else MatConvert_SeqAIJ_SeqSBAIJ_Preallocate(A, &rowlengths);
183:   if (reuse != MAT_REUSE_MATRIX) {
184:     MatCreate(PetscObjectComm((PetscObject)A), &B);
185:     MatSetSizes(B, m, n, m, n);
186:     MatSetType(B, MATSEQSBAIJ);
187:     MatSeqSBAIJSetPreallocation(B, bs, 0, rowlengths);
188:   } else B = *newmat;

190:   if (bs == 1 && !miss) {
191:     b  = (Mat_SeqSBAIJ *)(B->data);
192:     bi = b->i;
193:     bj = b->j;
194:     bv = b->a;

196:     bi[0] = 0;
197:     for (i = 0; i < m; i++) {
198:       aj = a->j + a->diag[i];
199:       av = a->a + a->diag[i];
200:       for (j = 0; j < rowlengths[i]; j++) {
201:         *bj = *aj;
202:         bj++;
203:         aj++;
204:         *bv = *av;
205:         bv++;
206:         av++;
207:       }
208:       bi[i + 1]  = bi[i] + rowlengths[i];
209:       b->ilen[i] = rowlengths[i];
210:     }
211:     MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
212:     MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
213:   } else {
214:     MatSetOption(B, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
215:     /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
216:     /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before                                          */
217:     /*                      MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below                                        */
218:     MatConvert_Basic(A, newtype, MAT_REUSE_MATRIX, &B);
219:   }
220:   PetscFree(rowlengths);
221:   if (reuse == MAT_INPLACE_MATRIX) {
222:     MatHeaderReplace(A, &B);
223:   } else *newmat = B;
224:   return 0;
225: }

227: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
228: {
229:   Mat           B;
230:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
231:   Mat_SeqBAIJ  *b;
232:   PetscInt     *ai = a->i, *aj = a->j, m = A->rmap->N, n = A->cmap->n, i, k, *bi, *bj, *browlengths, nz, *browstart, itmp;
233:   PetscInt      bs = A->rmap->bs, bs2 = bs * bs, mbs = m / bs, col, row;
234:   MatScalar    *av, *bv;
235: #if defined(PETSC_USE_COMPLEX)
236:   const int aconj = A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0;
237: #else
238:   const int aconj = 0;
239: #endif

241:   /* compute browlengths of newmat */
242:   PetscMalloc2(mbs, &browlengths, mbs, &browstart);
243:   for (i = 0; i < mbs; i++) browlengths[i] = 0;
244:   for (i = 0; i < mbs; i++) {
245:     nz = ai[i + 1] - ai[i];
246:     aj++;                      /* skip diagonal */
247:     for (k = 1; k < nz; k++) { /* no. of lower triangular blocks */
248:       browlengths[*aj]++;
249:       aj++;
250:     }
251:     browlengths[i] += nz; /* no. of upper triangular blocks */
252:   }

254:   if (reuse != MAT_REUSE_MATRIX) {
255:     MatCreate(PetscObjectComm((PetscObject)A), &B);
256:     MatSetSizes(B, m, n, m, n);
257:     MatSetType(B, MATSEQBAIJ);
258:     MatSeqBAIJSetPreallocation(B, bs, 0, browlengths);
259:   } else B = *newmat;

261:   b  = (Mat_SeqBAIJ *)(B->data);
262:   bi = b->i;
263:   bj = b->j;
264:   bv = b->a;

266:   /* set b->i */
267:   bi[0] = 0;
268:   for (i = 0; i < mbs; i++) {
269:     b->ilen[i]   = browlengths[i];
270:     bi[i + 1]    = bi[i] + browlengths[i];
271:     browstart[i] = bi[i];
272:   }

275:   /* set b->j and b->a */
276:   aj = a->j;
277:   av = a->a;
278:   for (i = 0; i < mbs; i++) {
279:     /* diagonal block */
280:     *(bj + browstart[i]) = *aj;
281:     aj++;

283:     itmp = bs2 * browstart[i];
284:     for (k = 0; k < bs2; k++) {
285:       *(bv + itmp + k) = *av;
286:       av++;
287:     }
288:     browstart[i]++;

290:     nz = ai[i + 1] - ai[i] - 1;
291:     while (nz--) {
292:       /* lower triangular blocks - transpose blocks of A */
293:       *(bj + browstart[*aj]) = i; /* block col index */

295:       itmp = bs2 * browstart[*aj]; /* row index */
296:       for (col = 0; col < bs; col++) {
297:         k = col;
298:         for (row = 0; row < bs; row++) {
299:           bv[itmp + col * bs + row] = aconj ? PetscConj(av[k]) : av[k];
300:           k += bs;
301:         }
302:       }
303:       browstart[*aj]++;

305:       /* upper triangular blocks */
306:       *(bj + browstart[i]) = *aj;
307:       aj++;

309:       itmp = bs2 * browstart[i];
310:       for (k = 0; k < bs2; k++) bv[itmp + k] = av[k];
311:       av += bs2;
312:       browstart[i]++;
313:     }
314:   }
315:   PetscFree2(browlengths, browstart);
316:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
317:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);

319:   if (reuse == MAT_INPLACE_MATRIX) {
320:     MatHeaderReplace(A, &B);
321:   } else *newmat = B;
322:   return 0;
323: }

325: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
326: {
327:   Mat           B;
328:   Mat_SeqBAIJ  *a = (Mat_SeqBAIJ *)A->data;
329:   Mat_SeqSBAIJ *b;
330:   PetscInt     *ai = a->i, *aj, m = A->rmap->N, n = A->cmap->n, i, j, k, *bi, *bj, *browlengths;
331:   PetscInt      bs = A->rmap->bs, bs2 = bs * bs, mbs = m / bs, dd;
332:   MatScalar    *av, *bv;
333:   PetscBool     flg;

337:   MatMissingDiagonal_SeqBAIJ(A, &flg, &dd); /* check for missing diagonals, then mark diag */

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

343:   if (reuse != MAT_REUSE_MATRIX) {
344:     MatCreate(PetscObjectComm((PetscObject)A), &B);
345:     MatSetSizes(B, m, n, m, n);
346:     MatSetType(B, MATSEQSBAIJ);
347:     MatSeqSBAIJSetPreallocation(B, bs, 0, browlengths);
348:   } else B = *newmat;

350:   b  = (Mat_SeqSBAIJ *)(B->data);
351:   bi = b->i;
352:   bj = b->j;
353:   bv = b->a;

355:   bi[0] = 0;
356:   for (i = 0; i < mbs; i++) {
357:     aj = a->j + a->diag[i];
358:     av = a->a + (a->diag[i]) * bs2;
359:     for (j = 0; j < browlengths[i]; j++) {
360:       *bj = *aj;
361:       bj++;
362:       aj++;
363:       for (k = 0; k < bs2; k++) {
364:         *bv = *av;
365:         bv++;
366:         av++;
367:       }
368:     }
369:     bi[i + 1]  = bi[i] + browlengths[i];
370:     b->ilen[i] = browlengths[i];
371:   }
372:   PetscFree(browlengths);
373:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
374:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);

376:   if (reuse == MAT_INPLACE_MATRIX) {
377:     MatHeaderReplace(A, &B);
378:   } else *newmat = B;
379:   return 0;
380: }