1 2 #include <../src/mat/impls/baij/seq/baij.h> 3 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; 13 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; 33 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)); 52 53 if (reuse == MAT_INPLACE_MATRIX) { 54 PetscCall(MatHeaderReplace(A, &B)); 55 } else *newmat = B; 56 PetscFunctionReturn(PETSC_SUCCESS); 57 } 58 59 #include <../src/mat/impls/aij/seq/aij.h> 60 61 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ_Preallocate(Mat A, PetscInt **nnz) 62 { 63 PetscInt m, n, bs = PetscAbs(A->rmap->bs); 64 PetscInt *ii; 65 Mat_SeqAIJ *Aa = (Mat_SeqAIJ *)A->data; 66 67 PetscFunctionBegin; 68 PetscCall(MatGetSize(A, &m, &n)); 69 PetscCall(PetscCalloc1(m / bs, nnz)); 70 PetscCall(PetscMalloc1(bs, &ii)); 71 72 /* loop over all potential blocks in the matrix to see if the potential block has a nonzero */ 73 for (PetscInt i = 0; i < m / bs; i++) { 74 for (PetscInt k = 0; k < bs; k++) ii[k] = Aa->i[i * bs + k]; 75 for (PetscInt j = 0; j < n / bs; j++) { 76 for (PetscInt k = 0; k < bs; k++) { 77 for (; ii[k] < Aa->i[i * bs + k + 1] && Aa->j[ii[k]] < (j + 1) * bs; ii[k]++) { 78 if (Aa->j[ii[k]] >= j * bs) { 79 /* found a nonzero in the potential block so allocate for that block and jump to check the next potential block */ 80 (*nnz)[i]++; 81 goto theend; 82 } 83 } 84 } 85 theend:; 86 } 87 } 88 PetscCall(PetscFree(ii)); 89 PetscFunctionReturn(PETSC_SUCCESS); 90 } 91 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 = PetscAbs(A->rmap->bs); 98 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; 108 109 if (bs == 1) { 110 b = (Mat_SeqBAIJ *)(B->data); 111 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)); 116 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 } 126 127 if (reuse == MAT_INPLACE_MATRIX) { 128 PetscCall(MatHeaderReplace(A, &B)); 129 } else *newmat = B; 130 PetscFunctionReturn(PETSC_SUCCESS); 131 } 132