1 #include <../src/mat/impls/baij/seq/baij.h> 2 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; 12 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; 32 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)); 51 52 if (reuse == MAT_INPLACE_MATRIX) { 53 PetscCall(MatHeaderReplace(A, &B)); 54 } else *newmat = B; 55 PetscFunctionReturn(PETSC_SUCCESS); 56 } 57 58 #include <../src/mat/impls/aij/seq/aij.h> 59 60 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ_Preallocate(Mat A, PetscInt **nnz) 61 { 62 Mat_SeqAIJ *Aa = (Mat_SeqAIJ *)A->data; 63 PetscInt m, n, bs = PetscAbs(A->rmap->bs); 64 const PetscInt *ai = Aa->i, *aj = Aa->j; 65 66 PetscFunctionBegin; 67 PetscCall(MatGetSize(A, &m, &n)); 68 69 if (bs == 1) { 70 PetscCall(PetscMalloc1(m, nnz)); 71 for (PetscInt i = 0; i < m; i++) (*nnz)[i] = ai[i + 1] - ai[i]; 72 } else { 73 PetscHSetIJ ht; 74 PetscInt k; 75 PetscHashIJKey key; 76 PetscBool missing; 77 78 PetscCall(PetscHSetIJCreate(&ht)); 79 PetscCall(PetscCalloc1(m / bs, nnz)); 80 for (PetscInt i = 0; i < m; i++) { 81 key.i = i / bs; 82 for (k = ai[i]; k < ai[i + 1]; k++) { 83 key.j = aj[k] / bs; 84 PetscCall(PetscHSetIJQueryAdd(ht, key, &missing)); 85 if (missing) (*nnz)[key.i]++; 86 } 87 } 88 PetscCall(PetscHSetIJDestroy(&ht)); 89 } 90 PetscFunctionReturn(PETSC_SUCCESS); 91 } 92 93 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat) 94 { 95 Mat B; 96 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data; 97 Mat_SeqBAIJ *b; 98 PetscInt m = A->rmap->N, n = A->cmap->N, *rowlengths, bs = PetscAbs(A->rmap->bs); 99 100 PetscFunctionBegin; 101 if (reuse != MAT_REUSE_MATRIX) { 102 PetscCall(MatConvert_SeqAIJ_SeqBAIJ_Preallocate(A, &rowlengths)); 103 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B)); 104 PetscCall(MatSetSizes(B, m, n, m, n)); 105 PetscCall(MatSetType(B, MATSEQBAIJ)); 106 PetscCall(MatSeqBAIJSetPreallocation(B, bs, 0, rowlengths)); 107 PetscCall(PetscFree(rowlengths)); 108 } else B = *newmat; 109 110 if (bs == 1) { 111 b = (Mat_SeqBAIJ *)B->data; 112 113 PetscCall(PetscArraycpy(b->i, a->i, m + 1)); 114 PetscCall(PetscArraycpy(b->ilen, a->ilen, m)); 115 PetscCall(PetscArraycpy(b->j, a->j, a->nz)); 116 PetscCall(PetscArraycpy(b->a, a->a, a->nz)); 117 118 PetscCall(MatSetOption(B, MAT_ROW_ORIENTED, PETSC_TRUE)); 119 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 120 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 121 } else { 122 /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */ 123 /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before */ 124 /* MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below */ 125 PetscCall(MatConvert_Basic(A, newtype, MAT_REUSE_MATRIX, &B)); 126 } 127 128 if (reuse == MAT_INPLACE_MATRIX) { 129 PetscCall(MatHeaderReplace(A, &B)); 130 } else *newmat = B; 131 PetscFunctionReturn(PETSC_SUCCESS); 132 } 133