Actual source code: aijsbaij.c

      1: #include <../src/mat/impls/aij/seq/aij.h>
      2: #include <../src/mat/impls/baij/seq/baij.h>
      3: #include <../src/mat/impls/sbaij/seq/sbaij.h>

      5: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
      6: {
      7:   Mat           B;
      8:   Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
      9:   Mat_SeqAIJ   *b;
     10:   PetscInt     *ai = a->i, *aj = a->j, m = A->rmap->N, n = A->cmap->n, i, j, k, *bi, *bj, *rowlengths, nz, *rowstart, itmp;
     11:   PetscInt      bs = A->rmap->bs, bs2 = bs * bs, mbs = A->rmap->N / bs, diagcnt = 0;
     12:   MatScalar    *av, *bv;
     13: #if defined(PETSC_USE_COMPLEX)
     14:   const int aconj = A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0;
     15: #else
     16:   const int aconj = 0;
     17: #endif

     19:   PetscFunctionBegin;
     20:   /* compute rowlengths of newmat */
     21:   PetscCall(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:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
     46:     PetscCall(MatSetSizes(B, m, n, m, n));
     47:     PetscCall(MatSetType(B, MATSEQAIJ));
     48:     PetscCall(MatSeqAIJSetPreallocation(B, 0, rowlengths));
     49:     PetscCall(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:   }
     67:   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);

     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:   PetscCall(PetscFree2(rowlengths, rowstart));
    113:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    114:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    278:   /* set b->i */
    279:   bi[0] = 0;
    280:   for (i = 0; i < mbs; i++) {
    281:     b->ilen[i]   = browlengths[i];
    282:     bi[i + 1]    = bi[i] + browlengths[i];
    283:     browstart[i] = bi[i];
    284:   }
    285:   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);

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

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

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

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

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

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

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

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

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

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

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

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

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

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