1 #include <../src/mat/impls/aij/seq/aij.h>
2 #include <../src/mat/impls/baij/seq/baij.h>
3 #include <../src/mat/impls/sbaij/seq/sbaij.h>
4
MatConvert_SeqSBAIJ_SeqAIJ(Mat A,MatType newtype,MatReuse reuse,Mat * newmat)5 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
6 {
7 Mat B;
8 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
9 Mat_SeqAIJ *b;
10 PetscInt *ai = a->i, *aj = a->j, m = A->rmap->N, n = A->cmap->n, i, j, k, *bi, *bj, *rowlengths, nz, *rowstart, itmp;
11 PetscInt bs = A->rmap->bs, bs2 = bs * bs, mbs = A->rmap->N / bs, diagcnt = 0;
12 MatScalar *av, *bv;
13 const int aconj = PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0;
14
15 PetscFunctionBegin;
16 /* compute rowlengths of newmat */
17 PetscCall(PetscMalloc2(m, &rowlengths, m + 1, &rowstart));
18
19 for (i = 0; i < mbs; i++) rowlengths[i * bs] = 0;
20 k = 0;
21 for (i = 0; i < mbs; i++) {
22 nz = ai[i + 1] - ai[i];
23 if (nz) {
24 rowlengths[k] += nz; /* no. of upper triangular blocks */
25 if (*aj == i) {
26 aj++;
27 diagcnt++;
28 nz--;
29 } /* skip diagonal */
30 for (j = 0; j < nz; j++) { /* no. of lower triangular blocks */
31 rowlengths[(*aj) * bs]++;
32 aj++;
33 }
34 }
35 rowlengths[k] *= bs;
36 for (j = 1; j < bs; j++) rowlengths[k + j] = rowlengths[k];
37 k += bs;
38 }
39
40 if (reuse != MAT_REUSE_MATRIX) {
41 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
42 PetscCall(MatSetSizes(B, m, n, m, n));
43 PetscCall(MatSetType(B, MATSEQAIJ));
44 PetscCall(MatSeqAIJSetPreallocation(B, 0, rowlengths));
45 PetscCall(MatSetBlockSize(B, A->rmap->bs));
46 } else B = *newmat;
47
48 b = (Mat_SeqAIJ *)B->data;
49 bi = b->i;
50 bj = b->j;
51 bv = b->a;
52
53 /* set b->i */
54 bi[0] = 0;
55 rowstart[0] = 0;
56 for (i = 0; i < mbs; i++) {
57 for (j = 0; j < bs; j++) {
58 b->ilen[i * bs + j] = rowlengths[i * bs];
59 rowstart[i * bs + j + 1] = rowstart[i * bs + j] + rowlengths[i * bs];
60 }
61 bi[i + 1] = bi[i] + rowlengths[i * bs] / bs;
62 }
63 PetscCheck(bi[mbs] == 2 * a->nz - diagcnt, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bi[mbs]: %" PetscInt_FMT " != 2*a->nz-diagcnt: %" PetscInt_FMT, bi[mbs], 2 * a->nz - diagcnt);
64
65 /* set b->j and b->a */
66 aj = a->j;
67 av = a->a;
68 for (i = 0; i < mbs; i++) {
69 nz = ai[i + 1] - ai[i];
70 /* diagonal block */
71 if (nz && *aj == i) {
72 nz--;
73 for (j = 0; j < bs; j++) { /* row i*bs+j */
74 itmp = i * bs + j;
75 for (k = 0; k < bs; k++) { /* col i*bs+k */
76 *(bj + rowstart[itmp]) = (*aj) * bs + k;
77 *(bv + rowstart[itmp]) = *(av + k * bs + j);
78 rowstart[itmp]++;
79 }
80 }
81 aj++;
82 av += bs2;
83 }
84
85 while (nz--) {
86 /* lower triangular blocks */
87 for (j = 0; j < bs; j++) { /* row (*aj)*bs+j */
88 itmp = (*aj) * bs + j;
89 for (k = 0; k < bs; k++) { /* col i*bs+k */
90 *(bj + rowstart[itmp]) = i * bs + k;
91 *(bv + rowstart[itmp]) = aconj ? PetscConj(*(av + j * bs + k)) : *(av + j * bs + k);
92 rowstart[itmp]++;
93 }
94 }
95 /* upper triangular blocks */
96 for (j = 0; j < bs; j++) { /* row i*bs+j */
97 itmp = i * bs + j;
98 for (k = 0; k < bs; k++) { /* col (*aj)*bs+k */
99 *(bj + rowstart[itmp]) = (*aj) * bs + k;
100 *(bv + rowstart[itmp]) = *(av + k * bs + j);
101 rowstart[itmp]++;
102 }
103 }
104 aj++;
105 av += bs2;
106 }
107 }
108 PetscCall(PetscFree2(rowlengths, rowstart));
109 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
110 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
111
112 if (reuse == MAT_INPLACE_MATRIX) {
113 PetscCall(MatHeaderReplace(A, &B));
114 } else {
115 *newmat = B;
116 }
117 PetscFunctionReturn(PETSC_SUCCESS);
118 }
119
120 #include <../src/mat/impls/aij/seq/aij.h>
121
MatConvert_SeqAIJ_SeqSBAIJ_Preallocate(Mat A,PetscInt ** nnz)122 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ_Preallocate(Mat A, PetscInt **nnz)
123 {
124 Mat_SeqAIJ *Aa = (Mat_SeqAIJ *)A->data;
125 PetscInt m, n, bs = A->rmap->bs;
126 const PetscInt *ai = Aa->i, *aj = Aa->j;
127
128 PetscFunctionBegin;
129 PetscCall(MatGetSize(A, &m, &n));
130
131 if (bs == 1) {
132 const PetscInt *adiag;
133
134 PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
135 PetscCall(PetscMalloc1(m, nnz));
136 for (PetscInt i = 0; i < m; i++) {
137 if (adiag[i] == ai[i + 1]) {
138 (*nnz)[i] = 0;
139 for (PetscInt j = ai[i]; j < ai[i + 1]; j++) (*nnz)[i] += (aj[j] > i);
140 } else (*nnz)[i] = ai[i + 1] - adiag[i];
141 }
142 } else {
143 PetscHSetIJ ht;
144 PetscHashIJKey key;
145 PetscBool missing;
146
147 PetscCall(PetscHSetIJCreate(&ht));
148 PetscCall(PetscCalloc1(m / bs, nnz));
149 for (PetscInt i = 0; i < m; i++) {
150 key.i = i / bs;
151 for (PetscInt k = ai[i]; k < ai[i + 1]; k++) {
152 key.j = aj[k] / bs;
153 if (key.j < key.i) continue;
154 PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
155 if (missing) (*nnz)[key.i]++;
156 }
157 }
158 PetscCall(PetscHSetIJDestroy(&ht));
159 }
160 PetscFunctionReturn(PETSC_SUCCESS);
161 }
162
MatConvert_SeqAIJ_SeqSBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat * newmat)163 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
164 {
165 Mat B;
166 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
167 Mat_SeqSBAIJ *b;
168 PetscInt *ai = a->i, *aj, m = A->rmap->N, n = A->cmap->N, i, j, *bi, *bj, *rowlengths, bs = A->rmap->bs;
169 MatScalar *av, *bv;
170 PetscBool miss = PETSC_FALSE;
171 const PetscInt *adiag;
172
173 PetscFunctionBegin;
174 #if !defined(PETSC_USE_COMPLEX)
175 PetscCheck(A->symmetric == PETSC_BOOL3_TRUE, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
176 #else
177 PetscCheck(A->symmetric == PETSC_BOOL3_TRUE || A->hermitian == PETSC_BOOL3_TRUE, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Matrix must be either symmetric or Hermitian. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE) and/or MatSetOption(mat,MAT_HERMITIAN,PETSC_TRUE)");
178 #endif
179 PetscCheck(n == m, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix must be square");
180 PetscCall(MatGetDiagonalMarkers_SeqAIJ(A, &adiag, NULL));
181 if (bs == 1) {
182 PetscCall(PetscMalloc1(m, &rowlengths));
183 for (i = 0; i < m; i++) {
184 if (adiag[i] == ai[i + 1]) { /* missing diagonal */
185 rowlengths[i] = (ai[i + 1] - ai[i]) + 1; /* allocate some extra space */
186 miss = PETSC_TRUE;
187 } else {
188 rowlengths[i] = (ai[i + 1] - adiag[i]);
189 }
190 }
191 } else PetscCall(MatConvert_SeqAIJ_SeqSBAIJ_Preallocate(A, &rowlengths));
192 if (reuse != MAT_REUSE_MATRIX) {
193 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
194 PetscCall(MatSetSizes(B, m, n, m, n));
195 PetscCall(MatSetType(B, MATSEQSBAIJ));
196 PetscCall(MatSeqSBAIJSetPreallocation(B, bs, 0, rowlengths));
197 } else B = *newmat;
198
199 if (bs == 1 && !miss) {
200 b = (Mat_SeqSBAIJ *)B->data;
201 bi = b->i;
202 bj = b->j;
203 bv = b->a;
204
205 bi[0] = 0;
206 for (i = 0; i < m; i++) {
207 aj = a->j + adiag[i];
208 av = a->a + adiag[i];
209 for (j = 0; j < rowlengths[i]; j++) {
210 *bj = *aj;
211 bj++;
212 aj++;
213 *bv = *av;
214 bv++;
215 av++;
216 }
217 bi[i + 1] = bi[i] + rowlengths[i];
218 b->ilen[i] = rowlengths[i];
219 }
220 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
221 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
222 } else {
223 PetscCall(MatSetOption(B, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE));
224 /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
225 /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before */
226 /* MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below */
227 PetscCall(MatConvert_Basic(A, newtype, MAT_REUSE_MATRIX, &B));
228 }
229 PetscCall(PetscFree(rowlengths));
230 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
231 else *newmat = B;
232 PetscFunctionReturn(PETSC_SUCCESS);
233 }
234
MatConvert_SeqSBAIJ_SeqBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat * newmat)235 PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
236 {
237 Mat B;
238 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
239 Mat_SeqBAIJ *b;
240 PetscInt *ai = a->i, *aj = a->j, m = A->rmap->N, n = A->cmap->n, *bi, *bj, *browlengths, nz, *browstart, itmp;
241 PetscInt bs = A->rmap->bs, bs2 = bs * bs, mbs = m / bs;
242 MatScalar *av, *bv;
243 const int aconj = PetscDefined(USE_COMPLEX) && A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0;
244
245 PetscFunctionBegin;
246 /* compute browlengths of newmat */
247 PetscCall(PetscMalloc2(mbs, &browlengths, mbs, &browstart));
248 for (PetscInt i = 0; i < mbs; i++) browlengths[i] = 0;
249 for (PetscInt i = 0; i < mbs; i++) {
250 nz = ai[i + 1] - ai[i];
251 aj++; /* skip diagonal */
252 for (PetscInt k = 1; k < nz; k++) { /* no. of lower triangular blocks */
253 browlengths[*aj]++;
254 aj++;
255 }
256 browlengths[i] += nz; /* no. of upper triangular blocks */
257 }
258
259 if (reuse != MAT_REUSE_MATRIX) {
260 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
261 PetscCall(MatSetSizes(B, m, n, m, n));
262 PetscCall(MatSetType(B, MATSEQBAIJ));
263 PetscCall(MatSeqBAIJSetPreallocation(B, bs, 0, browlengths));
264 } else B = *newmat;
265
266 b = (Mat_SeqBAIJ *)B->data;
267 bi = b->i;
268 bj = b->j;
269 bv = b->a;
270
271 /* set b->i */
272 bi[0] = 0;
273 for (PetscInt i = 0; i < mbs; i++) {
274 b->ilen[i] = browlengths[i];
275 bi[i + 1] = bi[i] + browlengths[i];
276 browstart[i] = bi[i];
277 }
278 PetscCheck(bi[mbs] == 2 * a->nz - mbs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "bi[mbs]: %" PetscInt_FMT " != 2*a->nz - mbs: %" PetscInt_FMT, bi[mbs], 2 * a->nz - mbs);
279
280 /* set b->j and b->a */
281 aj = a->j;
282 av = a->a;
283 for (PetscInt i = 0; i < mbs; i++) {
284 /* diagonal block */
285 *(bj + browstart[i]) = *aj;
286 aj++;
287
288 itmp = bs2 * browstart[i];
289 for (PetscInt k = 0; k < bs2; k++) {
290 *(bv + itmp + k) = *av;
291 av++;
292 }
293 browstart[i]++;
294
295 nz = ai[i + 1] - ai[i] - 1;
296 while (nz--) {
297 /* lower triangular blocks - transpose blocks of A */
298 *(bj + browstart[*aj]) = i; /* block col index */
299
300 itmp = bs2 * browstart[*aj]; /* row index */
301 for (PetscInt col = 0; col < bs; col++) {
302 PetscInt k = col;
303 for (PetscInt row = 0; row < bs; row++) {
304 bv[itmp + col * bs + row] = aconj ? PetscConj(av[k]) : av[k];
305 k += bs;
306 }
307 }
308 browstart[*aj]++;
309
310 /* upper triangular blocks */
311 *(bj + browstart[i]) = *aj;
312 aj++;
313
314 itmp = bs2 * browstart[i];
315 for (PetscInt k = 0; k < bs2; k++) bv[itmp + k] = av[k];
316 av += bs2;
317 browstart[i]++;
318 }
319 }
320 PetscCall(PetscFree2(browlengths, browstart));
321 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
322 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
323
324 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
325 else *newmat = B;
326 PetscFunctionReturn(PETSC_SUCCESS);
327 }
328
MatConvert_SeqBAIJ_SeqSBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat * newmat)329 PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
330 {
331 Mat B;
332 Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
333 Mat_SeqSBAIJ *b;
334 PetscInt *ai = a->i, *aj, m = A->rmap->N, n = A->cmap->n, k, *bi, *bj, *browlengths;
335 PetscInt bs = A->rmap->bs, bs2 = bs * bs, mbs = m / bs;
336 MatScalar *av, *bv;
337 const PetscInt *adiag;
338
339 PetscFunctionBegin;
340 PetscCheck(A->symmetric == PETSC_BOOL3_TRUE, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
341 PetscCheck(n == m, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix must be square");
342 PetscCall(MatGetDiagonalMarkers_SeqBAIJ(A, &adiag, NULL));
343 PetscCall(PetscMalloc1(mbs, &browlengths));
344 for (PetscInt i = 0; i < mbs; i++) browlengths[i] = ai[i + 1] - adiag[i];
345
346 if (reuse != MAT_REUSE_MATRIX) {
347 PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
348 PetscCall(MatSetSizes(B, m, n, m, n));
349 PetscCall(MatSetType(B, MATSEQSBAIJ));
350 PetscCall(MatSeqSBAIJSetPreallocation(B, bs, 0, browlengths));
351 } else B = *newmat;
352
353 b = (Mat_SeqSBAIJ *)B->data;
354 bi = b->i;
355 bj = b->j;
356 bv = b->a;
357
358 bi[0] = 0;
359 for (PetscInt i = 0; i < mbs; i++) {
360 aj = a->j + adiag[i];
361 av = a->a + (adiag[i]) * bs2;
362 for (PetscInt j = 0; j < browlengths[i]; j++) {
363 *bj = *aj;
364 bj++;
365 aj++;
366 for (k = 0; k < bs2; k++) {
367 *bv = *av;
368 bv++;
369 av++;
370 }
371 }
372 bi[i + 1] = bi[i] + browlengths[i];
373 b->ilen[i] = browlengths[i];
374 }
375 PetscCall(PetscFree(browlengths));
376 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
377 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
378
379 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
380 else *newmat = B;
381 PetscFunctionReturn(PETSC_SUCCESS);
382 }
383