Actual source code: aijbaij.c


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

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

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

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

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

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

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

 67:   PetscFunctionBegin;
 68:   PetscCall(MatGetSize(A, &m, &n));

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

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

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

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

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

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

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

129:   if (reuse == MAT_INPLACE_MATRIX) {
130:     PetscCall(MatHeaderReplace(A, &B));
131:   } else *newmat = B;
132:   PetscFunctionReturn(PETSC_SUCCESS);
133: }