Actual source code: aijbaij.c

      1: #include <../src/mat/impls/baij/seq/baij.h>

      3: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
      4: {
      5:   Mat          B;
      6:   Mat_SeqAIJ  *b;
      7:   PetscBool    roworiented;
      8:   Mat_SeqBAIJ *a  = (Mat_SeqBAIJ *)A->data;
      9:   PetscInt     bs = A->rmap->bs, *ai = a->i, *aj = a->j, n = A->rmap->N / bs, i, j, k;
     10:   PetscInt    *rowlengths, *rows, *cols, maxlen            = 0, ncols;
     11:   MatScalar   *aa = a->a;

     13:   PetscFunctionBegin;
     14:   if (reuse == MAT_REUSE_MATRIX) {
     15:     B = *newmat;
     16:     for (i = 0; i < n; i++) maxlen = PetscMax(maxlen, ai[i + 1] - ai[i]);
     17:   } else {
     18:     PetscCall(PetscMalloc1(n * bs, &rowlengths));
     19:     for (i = 0; i < n; i++) {
     20:       maxlen = PetscMax(maxlen, ai[i + 1] - ai[i]);
     21:       for (j = 0; j < bs; j++) rowlengths[i * bs + j] = bs * (ai[i + 1] - ai[i]);
     22:     }
     23:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
     24:     PetscCall(MatSetType(B, MATSEQAIJ));
     25:     PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
     26:     PetscCall(MatSetBlockSizes(B, A->rmap->bs, A->cmap->bs));
     27:     PetscCall(MatSeqAIJSetPreallocation(B, 0, rowlengths));
     28:     PetscCall(PetscFree(rowlengths));
     29:   }
     30:   b           = (Mat_SeqAIJ *)B->data;
     31:   roworiented = b->roworiented;

     33:   PetscCall(MatSetOption(B, MAT_ROW_ORIENTED, PETSC_FALSE));
     34:   PetscCall(PetscMalloc1(bs, &rows));
     35:   PetscCall(PetscMalloc1(bs * maxlen, &cols));
     36:   for (i = 0; i < n; i++) {
     37:     for (j = 0; j < bs; j++) rows[j] = i * bs + j;
     38:     ncols = ai[i + 1] - ai[i];
     39:     for (k = 0; k < ncols; k++) {
     40:       for (j = 0; j < bs; j++) cols[k * bs + j] = bs * (*aj) + j;
     41:       aj++;
     42:     }
     43:     PetscCall(MatSetValues(B, bs, rows, bs * ncols, cols, aa, INSERT_VALUES));
     44:     aa = PetscSafePointerPlusOffset(aa, ncols * bs * bs);
     45:   }
     46:   PetscCall(PetscFree(cols));
     47:   PetscCall(PetscFree(rows));
     48:   PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
     49:   PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
     50:   PetscCall(MatSetOption(B, MAT_ROW_ORIENTED, roworiented));

     52:   if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
     53:   else *newmat = B;
     54:   PetscFunctionReturn(PETSC_SUCCESS);
     55: }

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

     59: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ_Preallocate(Mat A, PetscInt **nnz)
     60: {
     61:   Mat_SeqAIJ     *Aa = (Mat_SeqAIJ *)A->data;
     62:   PetscInt        m, n, bs = A->rmap->bs;
     63:   const PetscInt *ai = Aa->i, *aj = Aa->j;

     65:   PetscFunctionBegin;
     66:   PetscCall(MatGetSize(A, &m, &n));

     68:   if (bs == 1) {
     69:     PetscCall(PetscMalloc1(m, nnz));
     70:     for (PetscInt i = 0; i < m; i++) (*nnz)[i] = ai[i + 1] - ai[i];
     71:   } else {
     72:     PetscHSetIJ    ht;
     73:     PetscInt       k;
     74:     PetscHashIJKey key;
     75:     PetscBool      missing;

     77:     PetscCall(PetscHSetIJCreate(&ht));
     78:     PetscCall(PetscCalloc1(m / bs, nnz));
     79:     for (PetscInt i = 0; i < m; i++) {
     80:       key.i = i / bs;
     81:       for (k = ai[i]; k < ai[i + 1]; k++) {
     82:         key.j = aj[k] / bs;
     83:         PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
     84:         if (missing) (*nnz)[key.i]++;
     85:       }
     86:     }
     87:     PetscCall(PetscHSetIJDestroy(&ht));
     88:   }
     89:   PetscFunctionReturn(PETSC_SUCCESS);
     90: }

     92: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
     93: {
     94:   Mat          B;
     95:   Mat_SeqAIJ  *a = (Mat_SeqAIJ *)A->data;
     96:   Mat_SeqBAIJ *b;
     97:   PetscInt     m = A->rmap->N, n = A->cmap->N, *rowlengths, bs = A->rmap->bs;

     99:   PetscFunctionBegin;
    100:   if (reuse != MAT_REUSE_MATRIX) {
    101:     PetscCall(MatConvert_SeqAIJ_SeqBAIJ_Preallocate(A, &rowlengths));
    102:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
    103:     PetscCall(MatSetSizes(B, m, n, m, n));
    104:     PetscCall(MatSetType(B, MATSEQBAIJ));
    105:     PetscCall(MatSeqBAIJSetPreallocation(B, bs, 0, rowlengths));
    106:     PetscCall(PetscFree(rowlengths));
    107:   } else B = *newmat;

    109:   if (bs == 1) {
    110:     b = (Mat_SeqBAIJ *)B->data;

    112:     PetscCall(PetscArraycpy(b->i, a->i, m + 1));
    113:     PetscCall(PetscArraycpy(b->ilen, a->ilen, m));
    114:     PetscCall(PetscArraycpy(b->j, a->j, a->nz));
    115:     PetscCall(PetscArraycpy(b->a, a->a, a->nz));

    117:     PetscCall(MatSetOption(B, MAT_ROW_ORIENTED, PETSC_TRUE));
    118:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
    119:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
    120:   } else {
    121:     /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
    122:     /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before                                          */
    123:     /*                      MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below                                        */
    124:     PetscCall(MatConvert_Basic(A, newtype, MAT_REUSE_MATRIX, &B));
    125:   }

    127:   if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
    128:   else *newmat = B;
    129:   PetscFunctionReturn(PETSC_SUCCESS);
    130: }