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: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
125: {
126:   Mat           B;
127:   Mat_SeqAIJ   *a = (Mat_SeqAIJ *)A->data;
128:   Mat_SeqSBAIJ *b;
129:   PetscInt     *ai = a->i, *aj, m = A->rmap->N, n = A->cmap->N, i, j, *bi, *bj, *rowlengths, bs = PetscAbs(A->rmap->bs);
130:   MatScalar    *av, *bv;
131:   PetscBool     miss = PETSC_FALSE;

133: #if !defined(PETSC_USE_COMPLEX)
135: #else
137: #endif

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

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

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

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

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

220:   if (reuse != MAT_REUSE_MATRIX) {
221:     MatCreate(PetscObjectComm((PetscObject)A), &B);
222:     MatSetSizes(B, m, n, m, n);
223:     MatSetType(B, MATSEQBAIJ);
224:     MatSeqBAIJSetPreallocation(B, bs, 0, browlengths);
225:   } else B = *newmat;

227:   b  = (Mat_SeqBAIJ *)(B->data);
228:   bi = b->i;
229:   bj = b->j;
230:   bv = b->a;

232:   /* set b->i */
233:   bi[0] = 0;
234:   for (i = 0; i < mbs; i++) {
235:     b->ilen[i]   = browlengths[i];
236:     bi[i + 1]    = bi[i] + browlengths[i];
237:     browstart[i] = bi[i];
238:   }

241:   /* set b->j and b->a */
242:   aj = a->j;
243:   av = a->a;
244:   for (i = 0; i < mbs; i++) {
245:     /* diagonal block */
246:     *(bj + browstart[i]) = *aj;
247:     aj++;

249:     itmp = bs2 * browstart[i];
250:     for (k = 0; k < bs2; k++) {
251:       *(bv + itmp + k) = *av;
252:       av++;
253:     }
254:     browstart[i]++;

256:     nz = ai[i + 1] - ai[i] - 1;
257:     while (nz--) {
258:       /* lower triangular blocks - transpose blocks of A */
259:       *(bj + browstart[*aj]) = i; /* block col index */

261:       itmp = bs2 * browstart[*aj]; /* row index */
262:       for (col = 0; col < bs; col++) {
263:         k = col;
264:         for (row = 0; row < bs; row++) {
265:           bv[itmp + col * bs + row] = aconj ? PetscConj(av[k]) : av[k];
266:           k += bs;
267:         }
268:       }
269:       browstart[*aj]++;

271:       /* upper triangular blocks */
272:       *(bj + browstart[i]) = *aj;
273:       aj++;

275:       itmp = bs2 * browstart[i];
276:       for (k = 0; k < bs2; k++) bv[itmp + k] = av[k];
277:       av += bs2;
278:       browstart[i]++;
279:     }
280:   }
281:   PetscFree2(browlengths, browstart);
282:   MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
283:   MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);

285:   if (reuse == MAT_INPLACE_MATRIX) {
286:     MatHeaderReplace(A, &B);
287:   } else *newmat = B;
288:   return 0;
289: }

291: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
292: {
293:   Mat           B;
294:   Mat_SeqBAIJ  *a = (Mat_SeqBAIJ *)A->data;
295:   Mat_SeqSBAIJ *b;
296:   PetscInt     *ai = a->i, *aj, m = A->rmap->N, n = A->cmap->n, i, j, k, *bi, *bj, *browlengths;
297:   PetscInt      bs = A->rmap->bs, bs2 = bs * bs, mbs = m / bs, dd;
298:   MatScalar    *av, *bv;
299:   PetscBool     flg;

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

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

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

316:   b  = (Mat_SeqSBAIJ *)(B->data);
317:   bi = b->i;
318:   bj = b->j;
319:   bv = b->a;

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

342:   if (reuse == MAT_INPLACE_MATRIX) {
343:     MatHeaderReplace(A, &B);
344:   } else *newmat = B;
345:   return 0;
346: }