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: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
125: {
126: Mat B;
127: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
128: Mat_SeqSBAIJ *b;
129: PetscInt *ai = a->i, *aj, m = A->rmap->N, n = A->cmap->N, i, j, *bi, *bj, *rowlengths, bs = PetscAbs(A->rmap->bs);
130: MatScalar *av, *bv;
131: PetscBool miss = PETSC_FALSE;
133: #if !defined(PETSC_USE_COMPLEX)
135: #else
137: #endif
140: PetscMalloc1(m / bs, &rowlengths);
141: for (i = 0; i < m / bs; i++) {
142: if (a->diag[i * bs] == ai[i * bs + 1]) { /* missing diagonal */
143: rowlengths[i] = (ai[i * bs + 1] - ai[i * bs]) / bs; /* allocate some extra space */
144: miss = PETSC_TRUE;
145: } else {
146: rowlengths[i] = (ai[i * bs + 1] - a->diag[i * bs]) / bs;
147: }
148: }
149: if (reuse != MAT_REUSE_MATRIX) {
150: MatCreate(PetscObjectComm((PetscObject)A), &B);
151: MatSetSizes(B, m, n, m, n);
152: MatSetType(B, MATSEQSBAIJ);
153: MatSeqSBAIJSetPreallocation(B, bs, 0, rowlengths);
154: } else B = *newmat;
156: if (bs == 1 && !miss) {
157: b = (Mat_SeqSBAIJ *)(B->data);
158: bi = b->i;
159: bj = b->j;
160: bv = b->a;
162: bi[0] = 0;
163: for (i = 0; i < m; i++) {
164: aj = a->j + a->diag[i];
165: av = a->a + a->diag[i];
166: for (j = 0; j < rowlengths[i]; j++) {
167: *bj = *aj;
168: bj++;
169: aj++;
170: *bv = *av;
171: bv++;
172: av++;
173: }
174: bi[i + 1] = bi[i] + rowlengths[i];
175: b->ilen[i] = rowlengths[i];
176: }
177: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
178: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
179: } else {
180: MatSetOption(B, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);
181: /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
182: /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before */
183: /* MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below */
184: MatConvert_Basic(A, newtype, MAT_REUSE_MATRIX, &B);
185: }
186: PetscFree(rowlengths);
187: if (reuse == MAT_INPLACE_MATRIX) {
188: MatHeaderReplace(A, &B);
189: } else *newmat = B;
190: return 0;
191: }
193: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
194: {
195: Mat B;
196: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ *)A->data;
197: Mat_SeqBAIJ *b;
198: PetscInt *ai = a->i, *aj = a->j, m = A->rmap->N, n = A->cmap->n, i, k, *bi, *bj, *browlengths, nz, *browstart, itmp;
199: PetscInt bs = A->rmap->bs, bs2 = bs * bs, mbs = m / bs, col, row;
200: MatScalar *av, *bv;
201: #if defined(PETSC_USE_COMPLEX)
202: const int aconj = A->hermitian == PETSC_BOOL3_TRUE ? 1 : 0;
203: #else
204: const int aconj = 0;
205: #endif
207: /* compute browlengths of newmat */
208: PetscMalloc2(mbs, &browlengths, mbs, &browstart);
209: for (i = 0; i < mbs; i++) browlengths[i] = 0;
210: for (i = 0; i < mbs; i++) {
211: nz = ai[i + 1] - ai[i];
212: aj++; /* skip diagonal */
213: for (k = 1; k < nz; k++) { /* no. of lower triangular blocks */
214: browlengths[*aj]++;
215: aj++;
216: }
217: browlengths[i] += nz; /* no. of upper triangular blocks */
218: }
220: if (reuse != MAT_REUSE_MATRIX) {
221: MatCreate(PetscObjectComm((PetscObject)A), &B);
222: MatSetSizes(B, m, n, m, n);
223: MatSetType(B, MATSEQBAIJ);
224: MatSeqBAIJSetPreallocation(B, bs, 0, browlengths);
225: } else B = *newmat;
227: b = (Mat_SeqBAIJ *)(B->data);
228: bi = b->i;
229: bj = b->j;
230: bv = b->a;
232: /* set b->i */
233: bi[0] = 0;
234: for (i = 0; i < mbs; i++) {
235: b->ilen[i] = browlengths[i];
236: bi[i + 1] = bi[i] + browlengths[i];
237: browstart[i] = bi[i];
238: }
241: /* set b->j and b->a */
242: aj = a->j;
243: av = a->a;
244: for (i = 0; i < mbs; i++) {
245: /* diagonal block */
246: *(bj + browstart[i]) = *aj;
247: aj++;
249: itmp = bs2 * browstart[i];
250: for (k = 0; k < bs2; k++) {
251: *(bv + itmp + k) = *av;
252: av++;
253: }
254: browstart[i]++;
256: nz = ai[i + 1] - ai[i] - 1;
257: while (nz--) {
258: /* lower triangular blocks - transpose blocks of A */
259: *(bj + browstart[*aj]) = i; /* block col index */
261: itmp = bs2 * browstart[*aj]; /* row index */
262: for (col = 0; col < bs; col++) {
263: k = col;
264: for (row = 0; row < bs; row++) {
265: bv[itmp + col * bs + row] = aconj ? PetscConj(av[k]) : av[k];
266: k += bs;
267: }
268: }
269: browstart[*aj]++;
271: /* upper triangular blocks */
272: *(bj + browstart[i]) = *aj;
273: aj++;
275: itmp = bs2 * browstart[i];
276: for (k = 0; k < bs2; k++) bv[itmp + k] = av[k];
277: av += bs2;
278: browstart[i]++;
279: }
280: }
281: PetscFree2(browlengths, browstart);
282: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
283: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
285: if (reuse == MAT_INPLACE_MATRIX) {
286: MatHeaderReplace(A, &B);
287: } else *newmat = B;
288: return 0;
289: }
291: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
292: {
293: Mat B;
294: Mat_SeqBAIJ *a = (Mat_SeqBAIJ *)A->data;
295: Mat_SeqSBAIJ *b;
296: PetscInt *ai = a->i, *aj, m = A->rmap->N, n = A->cmap->n, i, j, k, *bi, *bj, *browlengths;
297: PetscInt bs = A->rmap->bs, bs2 = bs * bs, mbs = m / bs, dd;
298: MatScalar *av, *bv;
299: PetscBool flg;
303: MatMissingDiagonal_SeqBAIJ(A, &flg, &dd); /* check for missing diagonals, then mark diag */
306: PetscMalloc1(mbs, &browlengths);
307: for (i = 0; i < mbs; i++) browlengths[i] = ai[i + 1] - a->diag[i];
309: if (reuse != MAT_REUSE_MATRIX) {
310: MatCreate(PetscObjectComm((PetscObject)A), &B);
311: MatSetSizes(B, m, n, m, n);
312: MatSetType(B, MATSEQSBAIJ);
313: MatSeqSBAIJSetPreallocation(B, bs, 0, browlengths);
314: } else B = *newmat;
316: b = (Mat_SeqSBAIJ *)(B->data);
317: bi = b->i;
318: bj = b->j;
319: bv = b->a;
321: bi[0] = 0;
322: for (i = 0; i < mbs; i++) {
323: aj = a->j + a->diag[i];
324: av = a->a + (a->diag[i]) * bs2;
325: for (j = 0; j < browlengths[i]; j++) {
326: *bj = *aj;
327: bj++;
328: aj++;
329: for (k = 0; k < bs2; k++) {
330: *bv = *av;
331: bv++;
332: av++;
333: }
334: }
335: bi[i + 1] = bi[i] + browlengths[i];
336: b->ilen[i] = browlengths[i];
337: }
338: PetscFree(browlengths);
339: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
340: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
342: if (reuse == MAT_INPLACE_MATRIX) {
343: MatHeaderReplace(A, &B);
344: } else *newmat = B;
345: return 0;
346: }