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