Actual source code: aijsbaij.c
2: #include <../src/mat/impls/aij/seq/aij.h>
3: #include <../src/mat/impls/baij/seq/baij.h>
4: #include <../src/mat/impls/sbaij/seq/sbaij.h>
6: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
7: {
8: Mat B;
9: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
10: Mat_SeqAIJ *b;
11: PetscInt *ai = a->i, *aj = a->j, m = A->rmap->N, n = A->cmap->n, i, j, k, *bi, *bj, *rowlengths, nz, *rowstart, itmp;
12: PetscInt bs = A->rmap->bs, bs2 = bs * bs, mbs = A->rmap->N / bs, diagcnt = 0;
13: MatScalar *av, *bv;
14: #if defined(PETSC_USE_COMPLEX)
15: const int aconj = A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0;
16: #else
17: const int aconj = 0;
18: #endif
20: /* compute rowlengths of newmat */
21: 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: MatCreate(PetscObjectComm((PetscObject)A), &B);
46: MatSetSizes(B, m, n, m, n);
47: MatSetType(B, MATSEQAIJ);
48: MatSeqAIJSetPreallocation(B, 0, rowlengths);
49: 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: }
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: PetscFree2(rowlengths, rowstart);
113: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
114: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
116: if (reuse == MAT_INPLACE_MATRIX) {
117: MatHeaderReplace(A, &B);
118: } else {
119: *newmat = B;
120: }
121: return 0;
122: }
124: #include <../src/mat/impls/aij/seq/aij.h>
126: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ_Preallocate(Mat A, PetscInt **nnz)
127: {
128: PetscInt m, n, bs = PetscAbs(A->rmap->bs);
129: PetscInt *ii;
130: Mat_SeqAIJ *Aa = (Mat_SeqAIJ *)A->data;
132: MatGetSize(A, &m, &n);
133: PetscCalloc1(m / bs, nnz);
134: PetscMalloc1(bs, &ii);
136: /* loop over all potential blocks in the matrix to see if the potential block has a nonzero */
137: for (PetscInt i = 0; i < m / bs; i++) {
138: for (PetscInt k = 0; k < bs; k++) ii[k] = Aa->i[i * bs + k];
139: for (PetscInt j = 0; j < n / bs; j++) {
140: for (PetscInt k = 0; k < bs; k++) {
141: for (; ii[k] < Aa->i[i * bs + k + 1] && Aa->j[ii[k]] < (j + 1) * bs; ii[k]++) {
142: if (j >= i && Aa->j[ii[k]] >= j * bs) {
143: /* found a nonzero in the potential block so allocate for that block and jump to check the next potential block */
144: (*nnz)[i]++;
145: goto theend;
146: }
147: }
148: }
149: theend:;
150: }
151: }
152: PetscFree(ii);
153: return 0;
154: }
156: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
157: {
158: Mat B;
159: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
160: Mat_SeqSBAIJ *b;
161: PetscInt *ai = a->i, *aj, m = A->rmap->N, n = A->cmap->N, i, j, *bi, *bj, *rowlengths, bs = PetscAbs(A->rmap->bs);
162: MatScalar *av, *bv;
163: PetscBool miss = PETSC_FALSE;
165: #if !defined(PETSC_USE_COMPLEX)
167: #else
169: #endif
172: if (bs == 1) {
173: PetscMalloc1(m, &rowlengths);
174: for (i = 0; i < m; i++) {
175: if (a->diag[i] == ai[i + 1]) { /* missing diagonal */
176: rowlengths[i] = (ai[i + 1] - ai[i]) + 1; /* allocate some extra space */
177: miss = PETSC_TRUE;
178: } else {
179: rowlengths[i] = (ai[i + 1] - a->diag[i]);
180: }
181: }
182: } else MatConvert_SeqAIJ_SeqSBAIJ_Preallocate(A, &rowlengths);
183: if (reuse != MAT_REUSE_MATRIX) {
184: MatCreate(PetscObjectComm((PetscObject)A), &B);
185: MatSetSizes(B, m, n, m, n);
186: MatSetType(B, MATSEQSBAIJ);
187: MatSeqSBAIJSetPreallocation(B, bs, 0, rowlengths);
188: } else B = *newmat;
190: if (bs == 1 && !miss) {
191: b = (Mat_SeqSBAIJ *)(B->data);
192: bi = b->i;
193: bj = b->j;
194: bv = b->a;
196: bi[0] = 0;
197: for (i = 0; i < m; i++) {
198: aj = a->j + a->diag[i];
199: av = a->a + a->diag[i];
200: for (j = 0; j < rowlengths[i]; j++) {
201: *bj = *aj;
202: bj++;
203: aj++;
204: *bv = *av;
205: bv++;
206: av++;
207: }
208: bi[i + 1] = bi[i] + rowlengths[i];
209: b->ilen[i] = rowlengths[i];
210: }
211: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
212: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
213: } else {
214: MatSetOption(B, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
215: /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
216: /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before */
217: /* MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below */
218: MatConvert_Basic(A, newtype, MAT_REUSE_MATRIX, &B);
219: }
220: PetscFree(rowlengths);
221: if (reuse == MAT_INPLACE_MATRIX) {
222: MatHeaderReplace(A, &B);
223: } else *newmat = B;
224: return 0;
225: }
227: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
228: {
229: Mat B;
230: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
231: Mat_SeqBAIJ *b;
232: PetscInt *ai = a->i, *aj = a->j, m = A->rmap->N, n = A->cmap->n, i, k, *bi, *bj, *browlengths, nz, *browstart, itmp;
233: PetscInt bs = A->rmap->bs, bs2 = bs * bs, mbs = m / bs, col, row;
234: MatScalar *av, *bv;
235: #if defined(PETSC_USE_COMPLEX)
236: const int aconj = A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0;
237: #else
238: const int aconj = 0;
239: #endif
241: /* compute browlengths of newmat */
242: PetscMalloc2(mbs, &browlengths, mbs, &browstart);
243: for (i = 0; i < mbs; i++) browlengths[i] = 0;
244: for (i = 0; i < mbs; i++) {
245: nz = ai[i + 1] - ai[i];
246: aj++; /* skip diagonal */
247: for (k = 1; k < nz; k++) { /* no. of lower triangular blocks */
248: browlengths[*aj]++;
249: aj++;
250: }
251: browlengths[i] += nz; /* no. of upper triangular blocks */
252: }
254: if (reuse != MAT_REUSE_MATRIX) {
255: MatCreate(PetscObjectComm((PetscObject)A), &B);
256: MatSetSizes(B, m, n, m, n);
257: MatSetType(B, MATSEQBAIJ);
258: MatSeqBAIJSetPreallocation(B, bs, 0, browlengths);
259: } else B = *newmat;
261: b = (Mat_SeqBAIJ *)(B->data);
262: bi = b->i;
263: bj = b->j;
264: bv = b->a;
266: /* set b->i */
267: bi[0] = 0;
268: for (i = 0; i < mbs; i++) {
269: b->ilen[i] = browlengths[i];
270: bi[i + 1] = bi[i] + browlengths[i];
271: browstart[i] = bi[i];
272: }
275: /* set b->j and b->a */
276: aj = a->j;
277: av = a->a;
278: for (i = 0; i < mbs; i++) {
279: /* diagonal block */
280: *(bj + browstart[i]) = *aj;
281: aj++;
283: itmp = bs2 * browstart[i];
284: for (k = 0; k < bs2; k++) {
285: *(bv + itmp + k) = *av;
286: av++;
287: }
288: browstart[i]++;
290: nz = ai[i + 1] - ai[i] - 1;
291: while (nz--) {
292: /* lower triangular blocks - transpose blocks of A */
293: *(bj + browstart[*aj]) = i; /* block col index */
295: itmp = bs2 * browstart[*aj]; /* row index */
296: for (col = 0; col < bs; col++) {
297: k = col;
298: for (row = 0; row < bs; row++) {
299: bv[itmp + col * bs + row] = aconj ? PetscConj(av[k]) : av[k];
300: k += bs;
301: }
302: }
303: browstart[*aj]++;
305: /* upper triangular blocks */
306: *(bj + browstart[i]) = *aj;
307: aj++;
309: itmp = bs2 * browstart[i];
310: for (k = 0; k < bs2; k++) bv[itmp + k] = av[k];
311: av += bs2;
312: browstart[i]++;
313: }
314: }
315: PetscFree2(browlengths, browstart);
316: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
317: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
319: if (reuse == MAT_INPLACE_MATRIX) {
320: MatHeaderReplace(A, &B);
321: } else *newmat = B;
322: return 0;
323: }
325: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
326: {
327: Mat B;
328: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
329: Mat_SeqSBAIJ *b;
330: PetscInt *ai = a->i, *aj, m = A->rmap->N, n = A->cmap->n, i, j, k, *bi, *bj, *browlengths;
331: PetscInt bs = A->rmap->bs, bs2 = bs * bs, mbs = m / bs, dd;
332: MatScalar *av, *bv;
333: PetscBool flg;
337: MatMissingDiagonal_SeqBAIJ(A, &flg, &dd); /* check for missing diagonals, then mark diag */
340: PetscMalloc1(mbs, &browlengths);
341: for (i = 0; i < mbs; i++) browlengths[i] = ai[i + 1] - a->diag[i];
343: if (reuse != MAT_REUSE_MATRIX) {
344: MatCreate(PetscObjectComm((PetscObject)A), &B);
345: MatSetSizes(B, m, n, m, n);
346: MatSetType(B, MATSEQSBAIJ);
347: MatSeqSBAIJSetPreallocation(B, bs, 0, browlengths);
348: } else B = *newmat;
350: b = (Mat_SeqSBAIJ *)(B->data);
351: bi = b->i;
352: bj = b->j;
353: bv = b->a;
355: bi[0] = 0;
356: for (i = 0; i < mbs; i++) {
357: aj = a->j + a->diag[i];
358: av = a->a + (a->diag[i]) * bs2;
359: for (j = 0; j < browlengths[i]; j++) {
360: *bj = *aj;
361: bj++;
362: aj++;
363: for (k = 0; k < bs2; k++) {
364: *bv = *av;
365: bv++;
366: av++;
367: }
368: }
369: bi[i + 1] = bi[i] + browlengths[i];
370: b->ilen[i] = browlengths[i];
371: }
372: PetscFree(browlengths);
373: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
374: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
376: if (reuse == MAT_INPLACE_MATRIX) {
377: MatHeaderReplace(A, &B);
378: } else *newmat = B;
379: return 0;
380: }