1 #include <../src/mat/impls/baij/seq/baij.h>
2
MatConvert_SeqBAIJ_SeqAIJ(Mat A,MatType newtype,MatReuse reuse,Mat * newmat)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) PetscCall(MatHeaderReplace(A, &B));
53 else *newmat = B;
54 PetscFunctionReturn(PETSC_SUCCESS);
55 }
56
57 #include <../src/mat/impls/aij/seq/aij.h>
58
MatConvert_SeqAIJ_SeqBAIJ_Preallocate(Mat A,PetscInt ** nnz)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;
64
65 PetscFunctionBegin;
66 PetscCall(MatGetSize(A, &m, &n));
67
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;
76
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 }
91
MatConvert_SeqAIJ_SeqBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat * newmat)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;
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) PetscCall(MatHeaderReplace(A, &B));
128 else *newmat = B;
129 PetscFunctionReturn(PETSC_SUCCESS);
130 }
131