Actual source code: aijsbaij.c
petsc-3.7.3 2016-08-01
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>
8: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
9: {
10: Mat B;
11: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
12: Mat_SeqAIJ *b;
14: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*rowlengths,nz,*rowstart,itmp;
15: PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=A->rmap->N/bs,diagcnt=0;
16: MatScalar *av,*bv;
19: /* compute rowlengths of newmat */
20: PetscMalloc2(m,&rowlengths,m+1,&rowstart);
22: for (i=0; i<mbs; i++) rowlengths[i*bs] = 0;
23: k = 0;
24: for (i=0; i<mbs; i++) {
25: nz = ai[i+1] - ai[i];
26: if (nz) {
27: rowlengths[k] += nz; /* no. of upper triangular blocks */
28: if (*aj == i) {aj++;diagcnt++;nz--;} /* skip diagonal */
29: for (j=0; j<nz; j++) { /* no. of lower triangular blocks */
30: rowlengths[(*aj)*bs]++; aj++;
31: }
32: }
33: rowlengths[k] *= bs;
34: for (j=1; j<bs; j++) {
35: rowlengths[k+j] = rowlengths[k];
36: }
37: k += bs;
38: /* printf(" rowlengths[%d]: %d\n",i, rowlengths[i]); */
39: }
41: MatCreate(PetscObjectComm((PetscObject)A),&B);
42: MatSetSizes(B,m,n,m,n);
43: MatSetType(B,MATSEQAIJ);
44: MatSeqAIJSetPreallocation(B,0,rowlengths);
45: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_FALSE);
47: B->rmap->bs = A->rmap->bs;
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: }
63: if (bi[mbs] != 2*a->nz - diagcnt) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bi[mbs]: %D != 2*a->nz-diagcnt: %D\n",bi[mbs],2*a->nz - diagcnt);
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]) = *(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: }
119: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
120: {
121: Mat B;
122: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
123: Mat_SeqSBAIJ *b;
125: PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->N,i,j,*bi,*bj,*rowlengths;
126: MatScalar *av,*bv;
129: if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
130: if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
132: PetscMalloc1(m,&rowlengths);
133: for (i=0; i<m; i++) {
134: rowlengths[i] = ai[i+1] - a->diag[i];
135: }
136: MatCreate(PetscObjectComm((PetscObject)A),&B);
137: MatSetSizes(B,m,n,m,n);
138: MatSetType(B,MATSEQSBAIJ);
139: MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,1,0,rowlengths);
141: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
143: b = (Mat_SeqSBAIJ*)(B->data);
144: bi = b->i;
145: bj = b->j;
146: bv = b->a;
148: bi[0] = 0;
149: for (i=0; i<m; i++) {
150: aj = a->j + a->diag[i];
151: av = a->a + a->diag[i];
152: for (j=0; j<rowlengths[i]; j++) {
153: *bj = *aj; bj++; aj++;
154: *bv = *av; bv++; av++;
155: }
156: bi[i+1] = bi[i] + rowlengths[i];
157: b->ilen[i] = rowlengths[i];
158: }
160: PetscFree(rowlengths);
161: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
162: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
164: if (reuse == MAT_INPLACE_MATRIX) {
165: MatHeaderReplace(A,&B);
166: } else {
167: *newmat = B;
168: }
169: return(0);
170: }
174: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
175: {
176: Mat B;
177: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
178: Mat_SeqBAIJ *b;
180: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp;
181: PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,col,row;
182: MatScalar *av,*bv;
185: /* compute browlengths of newmat */
186: PetscMalloc2(mbs,&browlengths,mbs,&browstart);
187: for (i=0; i<mbs; i++) browlengths[i] = 0;
188: for (i=0; i<mbs; i++) {
189: nz = ai[i+1] - ai[i];
190: aj++; /* skip diagonal */
191: for (k=1; k<nz; k++) { /* no. of lower triangular blocks */
192: browlengths[*aj]++; aj++;
193: }
194: browlengths[i] += nz; /* no. of upper triangular blocks */
195: }
197: MatCreate(PetscObjectComm((PetscObject)A),&B);
198: MatSetSizes(B,m,n,m,n);
199: MatSetType(B,MATSEQBAIJ);
200: MatSeqBAIJSetPreallocation(B,bs,0,browlengths);
201: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
203: b = (Mat_SeqBAIJ*)(B->data);
204: bi = b->i;
205: bj = b->j;
206: bv = b->a;
208: /* set b->i */
209: bi[0] = 0;
210: for (i=0; i<mbs; i++) {
211: b->ilen[i] = browlengths[i];
212: bi[i+1] = bi[i] + browlengths[i];
213: browstart[i] = bi[i];
214: }
215: if (bi[mbs] != 2*a->nz - mbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"bi[mbs]: %D != 2*a->nz - mbs: %D\n",bi[mbs],2*a->nz - mbs);
217: /* set b->j and b->a */
218: aj = a->j; av = a->a;
219: for (i=0; i<mbs; i++) {
220: /* diagonal block */
221: *(bj + browstart[i]) = *aj; aj++;
223: itmp = bs2*browstart[i];
224: for (k=0; k<bs2; k++) {
225: *(bv + itmp + k) = *av; av++;
226: }
227: browstart[i]++;
229: nz = ai[i+1] - ai[i] -1;
230: while (nz--) {
231: /* lower triangular blocks - transpose blocks of A */
232: *(bj + browstart[*aj]) = i; /* block col index */
234: itmp = bs2*browstart[*aj]; /* row index */
235: for (col=0; col<bs; col++) {
236: k = col;
237: for (row=0; row<bs; row++) {
238: bv[itmp + col*bs+row] = av[k]; k+=bs;
239: }
240: }
241: browstart[*aj]++;
243: /* upper triangular blocks */
244: *(bj + browstart[i]) = *aj; aj++;
246: itmp = bs2*browstart[i];
247: for (k=0; k<bs2; k++) {
248: bv[itmp + k] = av[k];
249: }
250: av += bs2;
251: browstart[i]++;
252: }
253: }
254: PetscFree2(browlengths,browstart);
255: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
256: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
258: if (reuse == MAT_INPLACE_MATRIX) {
259: MatHeaderReplace(A,&B);
260: } else {
261: *newmat = B;
262: }
263: return(0);
264: }
268: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
269: {
270: Mat B;
271: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
272: Mat_SeqSBAIJ *b;
274: PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*browlengths;
275: PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,dd;
276: MatScalar *av,*bv;
277: PetscBool flg;
280: if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
281: if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
282: MatMissingDiagonal_SeqBAIJ(A,&flg,&dd); /* check for missing diagonals, then mark diag */
283: if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal %D",dd);
285: PetscMalloc1(mbs,&browlengths);
286: for (i=0; i<mbs; i++) {
287: browlengths[i] = ai[i+1] - a->diag[i];
288: }
290: MatCreate(PetscObjectComm((PetscObject)A),&B);
291: MatSetSizes(B,m,n,m,n);
292: MatSetType(B,MATSEQSBAIJ);
293: MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,browlengths);
294: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
296: b = (Mat_SeqSBAIJ*)(B->data);
297: bi = b->i;
298: bj = b->j;
299: bv = b->a;
301: bi[0] = 0;
302: for (i=0; i<mbs; i++) {
303: aj = a->j + a->diag[i];
304: av = a->a + (a->diag[i])*bs2;
305: for (j=0; j<browlengths[i]; j++) {
306: *bj = *aj; bj++; aj++;
307: for (k=0; k<bs2; k++) {
308: *bv = *av; bv++; av++;
309: }
310: }
311: bi[i+1] = bi[i] + browlengths[i];
312: b->ilen[i] = browlengths[i];
313: }
314: PetscFree(browlengths);
315: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
316: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
318: if (reuse == MAT_INPLACE_MATRIX) {
319: MatHeaderReplace(A,&B);
320: } else {
321: *newmat = B;
322: }
323: return(0);
324: }