Actual source code: aijsbaij.c
petsc-3.9.4 2018-09-11
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;
12: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*rowlengths,nz,*rowstart,itmp;
13: PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=A->rmap->N/bs,diagcnt=0;
14: MatScalar *av,*bv;
17: /* compute rowlengths of newmat */
18: PetscMalloc2(m,&rowlengths,m+1,&rowstart);
20: for (i=0; i<mbs; i++) rowlengths[i*bs] = 0;
21: k = 0;
22: for (i=0; i<mbs; i++) {
23: nz = ai[i+1] - ai[i];
24: if (nz) {
25: rowlengths[k] += nz; /* no. of upper triangular blocks */
26: if (*aj == i) {aj++;diagcnt++;nz--;} /* skip diagonal */
27: for (j=0; j<nz; j++) { /* no. of lower triangular blocks */
28: rowlengths[(*aj)*bs]++; aj++;
29: }
30: }
31: rowlengths[k] *= bs;
32: for (j=1; j<bs; j++) {
33: rowlengths[k+j] = rowlengths[k];
34: }
35: k += bs;
36: /* printf(" rowlengths[%d]: %d\n",i, rowlengths[i]); */
37: }
39: MatCreate(PetscObjectComm((PetscObject)A),&B);
40: MatSetSizes(B,m,n,m,n);
41: MatSetType(B,MATSEQAIJ);
42: MatSeqAIJSetPreallocation(B,0,rowlengths);
43: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_FALSE);
45: B->rmap->bs = A->rmap->bs;
47: b = (Mat_SeqAIJ*)(B->data);
48: bi = b->i;
49: bj = b->j;
50: bv = b->a;
52: /* set b->i */
53: bi[0] = 0; rowstart[0] = 0;
54: for (i=0; i<mbs; i++) {
55: for (j=0; j<bs; j++) {
56: b->ilen[i*bs+j] = rowlengths[i*bs];
57: rowstart[i*bs+j+1] = rowstart[i*bs+j] + rowlengths[i*bs];
58: }
59: bi[i+1] = bi[i] + rowlengths[i*bs]/bs;
60: }
61: 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);
63: /* set b->j and b->a */
64: aj = a->j; av = a->a;
65: for (i=0; i<mbs; i++) {
66: nz = ai[i+1] - ai[i];
67: /* diagonal block */
68: if (nz && *aj == i) {
69: nz--;
70: for (j=0; j<bs; j++) { /* row i*bs+j */
71: itmp = i*bs+j;
72: for (k=0; k<bs; k++) { /* col i*bs+k */
73: *(bj + rowstart[itmp]) = (*aj)*bs+k;
74: *(bv + rowstart[itmp]) = *(av+k*bs+j);
75: rowstart[itmp]++;
76: }
77: }
78: aj++; av += bs2;
79: }
81: while (nz--) {
82: /* lower triangular blocks */
83: for (j=0; j<bs; j++) { /* row (*aj)*bs+j */
84: itmp = (*aj)*bs+j;
85: for (k=0; k<bs; k++) { /* col i*bs+k */
86: *(bj + rowstart[itmp]) = i*bs+k;
87: *(bv + rowstart[itmp]) = *(av+j*bs+k);
88: rowstart[itmp]++;
89: }
90: }
91: /* upper triangular blocks */
92: for (j=0; j<bs; j++) { /* row i*bs+j */
93: itmp = i*bs+j;
94: for (k=0; k<bs; k++) { /* col (*aj)*bs+k */
95: *(bj + rowstart[itmp]) = (*aj)*bs+k;
96: *(bv + rowstart[itmp]) = *(av+k*bs+j);
97: rowstart[itmp]++;
98: }
99: }
100: aj++; av += bs2;
101: }
102: }
103: PetscFree2(rowlengths,rowstart);
104: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
105: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
107: if (reuse == MAT_INPLACE_MATRIX) {
108: MatHeaderReplace(A,&B);
109: } else {
110: *newmat = B;
111: }
112: return(0);
113: }
115: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
116: {
117: Mat B;
118: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
119: Mat_SeqSBAIJ *b;
121: PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->N,i,j,*bi,*bj,*rowlengths;
122: MatScalar *av,*bv;
125: if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
126: if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
128: PetscMalloc1(m,&rowlengths);
129: for (i=0; i<m; i++) {
130: rowlengths[i] = ai[i+1] - a->diag[i];
131: }
132: MatCreate(PetscObjectComm((PetscObject)A),&B);
133: MatSetSizes(B,m,n,m,n);
134: MatSetType(B,MATSEQSBAIJ);
135: MatSeqSBAIJSetPreallocation(B,1,0,rowlengths);
137: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
139: b = (Mat_SeqSBAIJ*)(B->data);
140: bi = b->i;
141: bj = b->j;
142: bv = b->a;
144: bi[0] = 0;
145: for (i=0; i<m; i++) {
146: aj = a->j + a->diag[i];
147: av = a->a + a->diag[i];
148: for (j=0; j<rowlengths[i]; j++) {
149: *bj = *aj; bj++; aj++;
150: *bv = *av; bv++; av++;
151: }
152: bi[i+1] = bi[i] + rowlengths[i];
153: b->ilen[i] = rowlengths[i];
154: }
156: PetscFree(rowlengths);
157: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
158: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
160: if (reuse == MAT_INPLACE_MATRIX) {
161: MatHeaderReplace(A,&B);
162: } else {
163: *newmat = B;
164: }
165: return(0);
166: }
168: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
169: {
170: Mat B;
171: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
172: Mat_SeqBAIJ *b;
174: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp;
175: PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,col,row;
176: MatScalar *av,*bv;
179: /* compute browlengths of newmat */
180: PetscMalloc2(mbs,&browlengths,mbs,&browstart);
181: for (i=0; i<mbs; i++) browlengths[i] = 0;
182: for (i=0; i<mbs; i++) {
183: nz = ai[i+1] - ai[i];
184: aj++; /* skip diagonal */
185: for (k=1; k<nz; k++) { /* no. of lower triangular blocks */
186: browlengths[*aj]++; aj++;
187: }
188: browlengths[i] += nz; /* no. of upper triangular blocks */
189: }
191: MatCreate(PetscObjectComm((PetscObject)A),&B);
192: MatSetSizes(B,m,n,m,n);
193: MatSetType(B,MATSEQBAIJ);
194: MatSeqBAIJSetPreallocation(B,bs,0,browlengths);
195: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
197: b = (Mat_SeqBAIJ*)(B->data);
198: bi = b->i;
199: bj = b->j;
200: bv = b->a;
202: /* set b->i */
203: bi[0] = 0;
204: for (i=0; i<mbs; i++) {
205: b->ilen[i] = browlengths[i];
206: bi[i+1] = bi[i] + browlengths[i];
207: browstart[i] = bi[i];
208: }
209: 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);
211: /* set b->j and b->a */
212: aj = a->j; av = a->a;
213: for (i=0; i<mbs; i++) {
214: /* diagonal block */
215: *(bj + browstart[i]) = *aj; aj++;
217: itmp = bs2*browstart[i];
218: for (k=0; k<bs2; k++) {
219: *(bv + itmp + k) = *av; av++;
220: }
221: browstart[i]++;
223: nz = ai[i+1] - ai[i] -1;
224: while (nz--) {
225: /* lower triangular blocks - transpose blocks of A */
226: *(bj + browstart[*aj]) = i; /* block col index */
228: itmp = bs2*browstart[*aj]; /* row index */
229: for (col=0; col<bs; col++) {
230: k = col;
231: for (row=0; row<bs; row++) {
232: bv[itmp + col*bs+row] = av[k]; k+=bs;
233: }
234: }
235: browstart[*aj]++;
237: /* upper triangular blocks */
238: *(bj + browstart[i]) = *aj; aj++;
240: itmp = bs2*browstart[i];
241: for (k=0; k<bs2; k++) {
242: bv[itmp + k] = av[k];
243: }
244: av += bs2;
245: browstart[i]++;
246: }
247: }
248: PetscFree2(browlengths,browstart);
249: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
250: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
252: if (reuse == MAT_INPLACE_MATRIX) {
253: MatHeaderReplace(A,&B);
254: } else {
255: *newmat = B;
256: }
257: return(0);
258: }
260: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
261: {
262: Mat B;
263: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
264: Mat_SeqSBAIJ *b;
266: PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*browlengths;
267: PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,dd;
268: MatScalar *av,*bv;
269: PetscBool flg;
272: if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
273: if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
274: MatMissingDiagonal_SeqBAIJ(A,&flg,&dd); /* check for missing diagonals, then mark diag */
275: if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal %D",dd);
277: PetscMalloc1(mbs,&browlengths);
278: for (i=0; i<mbs; i++) {
279: browlengths[i] = ai[i+1] - a->diag[i];
280: }
282: MatCreate(PetscObjectComm((PetscObject)A),&B);
283: MatSetSizes(B,m,n,m,n);
284: MatSetType(B,MATSEQSBAIJ);
285: MatSeqSBAIJSetPreallocation(B,bs,0,browlengths);
286: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
288: b = (Mat_SeqSBAIJ*)(B->data);
289: bi = b->i;
290: bj = b->j;
291: bv = b->a;
293: bi[0] = 0;
294: for (i=0; i<mbs; i++) {
295: aj = a->j + a->diag[i];
296: av = a->a + (a->diag[i])*bs2;
297: for (j=0; j<browlengths[i]; j++) {
298: *bj = *aj; bj++; aj++;
299: for (k=0; k<bs2; k++) {
300: *bv = *av; bv++; av++;
301: }
302: }
303: bi[i+1] = bi[i] + browlengths[i];
304: b->ilen[i] = browlengths[i];
305: }
306: PetscFree(browlengths);
307: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
308: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
310: if (reuse == MAT_INPLACE_MATRIX) {
311: MatHeaderReplace(A,&B);
312: } else {
313: *newmat = B;
314: }
315: return(0);
316: }