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 ? 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) {aj++;diagcnt++;nz--;} /* skip diagonal */
30: for (j=0; j<nz; j++) { /* no. of lower triangular blocks */
31: rowlengths[(*aj)*bs]++; aj++;
32: }
33: }
34: rowlengths[k] *= bs;
35: for (j=1; j<bs; j++) {
36: rowlengths[k+j] = rowlengths[k];
37: }
38: k += bs;
39: }
41: if (reuse != MAT_REUSE_MATRIX) {
42: MatCreate(PetscObjectComm((PetscObject)A),&B);
43: MatSetSizes(B,m,n,m,n);
44: MatSetType(B,MATSEQAIJ);
45: MatSeqAIJSetPreallocation(B,0,rowlengths);
46: MatSetBlockSize(B,A->rmap->bs);
47: } else B = *newmat;
49: b = (Mat_SeqAIJ*)(B->data);
50: bi = b->i;
51: bj = b->j;
52: bv = b->a;
54: /* set b->i */
55: bi[0] = 0; 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: }
65: /* set b->j and b->a */
66: aj = a->j; av = a->a;
67: for (i=0; i<mbs; i++) {
68: nz = ai[i+1] - ai[i];
69: /* diagonal block */
70: if (nz && *aj == i) {
71: nz--;
72: for (j=0; j<bs; j++) { /* row i*bs+j */
73: itmp = i*bs+j;
74: for (k=0; k<bs; k++) { /* col i*bs+k */
75: *(bj + rowstart[itmp]) = (*aj)*bs+k;
76: *(bv + rowstart[itmp]) = *(av+k*bs+j);
77: rowstart[itmp]++;
78: }
79: }
80: aj++; av += bs2;
81: }
83: while (nz--) {
84: /* lower triangular blocks */
85: for (j=0; j<bs; j++) { /* row (*aj)*bs+j */
86: itmp = (*aj)*bs+j;
87: for (k=0; k<bs; k++) { /* col i*bs+k */
88: *(bj + rowstart[itmp]) = i*bs+k;
89: *(bv + rowstart[itmp]) = aconj ? PetscConj(*(av+j*bs+k)) : *(av+j*bs+k);
90: rowstart[itmp]++;
91: }
92: }
93: /* upper triangular blocks */
94: for (j=0; j<bs; j++) { /* row i*bs+j */
95: itmp = i*bs+j;
96: for (k=0; k<bs; k++) { /* col (*aj)*bs+k */
97: *(bj + rowstart[itmp]) = (*aj)*bs+k;
98: *(bv + rowstart[itmp]) = *(av+k*bs+j);
99: rowstart[itmp]++;
100: }
101: }
102: aj++; av += bs2;
103: }
104: }
105: PetscFree2(rowlengths,rowstart);
106: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
107: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
109: if (reuse == MAT_INPLACE_MATRIX) {
110: MatHeaderReplace(A,&B);
111: } else {
112: *newmat = B;
113: }
114: return 0;
115: }
117: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
118: {
119: Mat B;
120: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
121: Mat_SeqSBAIJ *b;
122: PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->N,i,j,*bi,*bj,*rowlengths,bs=PetscAbs(A->rmap->bs);
123: MatScalar *av,*bv;
124: PetscBool miss = PETSC_FALSE;
126: #if !defined(PETSC_USE_COMPLEX)
128: #else
130: #endif
133: PetscMalloc1(m/bs,&rowlengths);
134: for (i=0; i<m/bs; i++) {
135: if (a->diag[i*bs] == ai[i*bs+1]) { /* missing diagonal */
136: rowlengths[i] = (ai[i*bs+1] - ai[i*bs])/bs; /* allocate some extra space */
137: miss = PETSC_TRUE;
138: } else {
139: rowlengths[i] = (ai[i*bs+1] - a->diag[i*bs])/bs;
140: }
141: }
142: if (reuse != MAT_REUSE_MATRIX) {
143: MatCreate(PetscObjectComm((PetscObject)A),&B);
144: MatSetSizes(B,m,n,m,n);
145: MatSetType(B,MATSEQSBAIJ);
146: MatSeqSBAIJSetPreallocation(B,bs,0,rowlengths);
147: } else B = *newmat;
149: if (bs == 1 && !miss) {
150: b = (Mat_SeqSBAIJ*)(B->data);
151: bi = b->i;
152: bj = b->j;
153: bv = b->a;
155: bi[0] = 0;
156: for (i=0; i<m; i++) {
157: aj = a->j + a->diag[i];
158: av = a->a + a->diag[i];
159: for (j=0; j<rowlengths[i]; j++) {
160: *bj = *aj; bj++; aj++;
161: *bv = *av; bv++; av++;
162: }
163: bi[i+1] = bi[i] + rowlengths[i];
164: b->ilen[i] = rowlengths[i];
165: }
166: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
167: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
168: } else {
169: MatSetOption(B,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);
170: /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
171: /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before */
172: /* MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below */
173: MatConvert_Basic(A,newtype,MAT_REUSE_MATRIX,&B);
174: }
175: PetscFree(rowlengths);
176: if (reuse == MAT_INPLACE_MATRIX) {
177: MatHeaderReplace(A,&B);
178: } else *newmat = B;
179: return 0;
180: }
182: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
183: {
184: Mat B;
185: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
186: Mat_SeqBAIJ *b;
187: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp;
188: PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,col,row;
189: MatScalar *av,*bv;
190: #if defined(PETSC_USE_COMPLEX)
191: const int aconj = A->hermitian ? 1 : 0;
192: #else
193: const int aconj = 0;
194: #endif
196: /* compute browlengths of newmat */
197: PetscMalloc2(mbs,&browlengths,mbs,&browstart);
198: for (i=0; i<mbs; i++) browlengths[i] = 0;
199: for (i=0; i<mbs; i++) {
200: nz = ai[i+1] - ai[i];
201: aj++; /* skip diagonal */
202: for (k=1; k<nz; k++) { /* no. of lower triangular blocks */
203: browlengths[*aj]++; aj++;
204: }
205: browlengths[i] += nz; /* no. of upper triangular blocks */
206: }
208: if (reuse != MAT_REUSE_MATRIX) {
209: MatCreate(PetscObjectComm((PetscObject)A),&B);
210: MatSetSizes(B,m,n,m,n);
211: MatSetType(B,MATSEQBAIJ);
212: MatSeqBAIJSetPreallocation(B,bs,0,browlengths);
213: } else B = *newmat;
215: b = (Mat_SeqBAIJ*)(B->data);
216: bi = b->i;
217: bj = b->j;
218: bv = b->a;
220: /* set b->i */
221: bi[0] = 0;
222: for (i=0; i<mbs; i++) {
223: b->ilen[i] = browlengths[i];
224: bi[i+1] = bi[i] + browlengths[i];
225: browstart[i] = bi[i];
226: }
229: /* set b->j and b->a */
230: aj = a->j; av = a->a;
231: for (i=0; i<mbs; i++) {
232: /* diagonal block */
233: *(bj + browstart[i]) = *aj; aj++;
235: itmp = bs2*browstart[i];
236: for (k=0; k<bs2; k++) {
237: *(bv + itmp + k) = *av; av++;
238: }
239: browstart[i]++;
241: nz = ai[i+1] - ai[i] -1;
242: while (nz--) {
243: /* lower triangular blocks - transpose blocks of A */
244: *(bj + browstart[*aj]) = i; /* block col index */
246: itmp = bs2*browstart[*aj]; /* row index */
247: for (col=0; col<bs; col++) {
248: k = col;
249: for (row=0; row<bs; row++) {
250: bv[itmp + col*bs+row] = aconj ? PetscConj(av[k]) : av[k];
251: k+=bs;
252: }
253: }
254: browstart[*aj]++;
256: /* upper triangular blocks */
257: *(bj + browstart[i]) = *aj; aj++;
259: itmp = bs2*browstart[i];
260: for (k=0; k<bs2; k++) {
261: bv[itmp + k] = av[k];
262: }
263: av += bs2;
264: browstart[i]++;
265: }
266: }
267: PetscFree2(browlengths,browstart);
268: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
269: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
271: if (reuse == MAT_INPLACE_MATRIX) {
272: MatHeaderReplace(A,&B);
273: } else *newmat = B;
274: return 0;
275: }
277: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
278: {
279: Mat B;
280: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
281: Mat_SeqSBAIJ *b;
282: PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*browlengths;
283: PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,dd;
284: MatScalar *av,*bv;
285: PetscBool flg;
289: MatMissingDiagonal_SeqBAIJ(A,&flg,&dd); /* check for missing diagonals, then mark diag */
292: PetscMalloc1(mbs,&browlengths);
293: for (i=0; i<mbs; i++) {
294: browlengths[i] = ai[i+1] - a->diag[i];
295: }
297: if (reuse != MAT_REUSE_MATRIX) {
298: MatCreate(PetscObjectComm((PetscObject)A),&B);
299: MatSetSizes(B,m,n,m,n);
300: MatSetType(B,MATSEQSBAIJ);
301: MatSeqSBAIJSetPreallocation(B,bs,0,browlengths);
302: } else B = *newmat;
304: b = (Mat_SeqSBAIJ*)(B->data);
305: bi = b->i;
306: bj = b->j;
307: bv = b->a;
309: bi[0] = 0;
310: for (i=0; i<mbs; i++) {
311: aj = a->j + a->diag[i];
312: av = a->a + (a->diag[i])*bs2;
313: for (j=0; j<browlengths[i]; j++) {
314: *bj = *aj; bj++; aj++;
315: for (k=0; k<bs2; k++) {
316: *bv = *av; bv++; av++;
317: }
318: }
319: bi[i+1] = bi[i] + browlengths[i];
320: b->ilen[i] = browlengths[i];
321: }
322: PetscFree(browlengths);
323: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
324: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
326: if (reuse == MAT_INPLACE_MATRIX) {
327: MatHeaderReplace(A,&B);
328: } else *newmat = B;
329: return 0;
330: }