Actual source code: aijsbaij.c
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>
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;
15: PetscFunctionBegin;
16: /* compute rowlengths of newmat */
17: PetscCall(PetscMalloc2(m, &rowlengths, m + 1, &rowstart));
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: }
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;
48: b = (Mat_SeqAIJ *)B->data;
49: bi = b->i;
50: bj = b->j;
51: bv = b->a;
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);
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: }
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));
112: if (reuse == MAT_INPLACE_MATRIX) {
113: PetscCall(MatHeaderReplace(A, &B));
114: } else {
115: *newmat = B;
116: }
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: #include <../src/mat/impls/aij/seq/aij.h>
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;
128: PetscFunctionBegin;
129: PetscCall(MatGetSize(A, &m, &n));
131: if (bs == 1) {
132: const PetscInt *adiag;
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;
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: }
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;
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;
199: if (bs == 1 && !miss) {
200: b = (Mat_SeqSBAIJ *)B->data;
201: bi = b->i;
202: bj = b->j;
203: bv = b->a;
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: }
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;
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: }
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;
266: b = (Mat_SeqBAIJ *)B->data;
267: bi = b->i;
268: bj = b->j;
269: bv = b->a;
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);
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++;
288: itmp = bs2 * browstart[i];
289: for (PetscInt k = 0; k < bs2; k++) {
290: *(bv + itmp + k) = *av;
291: av++;
292: }
293: browstart[i]++;
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 */
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]++;
310: /* upper triangular blocks */
311: *(bj + browstart[i]) = *aj;
312: aj++;
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));
324: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
325: else *newmat = B;
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
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;
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];
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;
353: b = (Mat_SeqSBAIJ *)B->data;
354: bi = b->i;
355: bj = b->j;
356: bv = b->a;
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));
379: if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &B));
380: else *newmat = B;
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }