Actual source code: aijsbaij.c
petsc-3.12.5 2020-03-29
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: }
38: if (reuse != MAT_REUSE_MATRIX) {
39: MatCreate(PetscObjectComm((PetscObject)A),&B);
40: MatSetSizes(B,m,n,m,n);
41: MatSetType(B,MATSEQAIJ);
42: MatSeqAIJSetPreallocation(B,0,rowlengths);
43: MatSetBlockSize(B,A->rmap->bs);
44: } else B = *newmat;
46: b = (Mat_SeqAIJ*)(B->data);
47: bi = b->i;
48: bj = b->j;
49: bv = b->a;
51: /* set b->i */
52: bi[0] = 0; rowstart[0] = 0;
53: for (i=0; i<mbs; i++) {
54: for (j=0; j<bs; j++) {
55: b->ilen[i*bs+j] = rowlengths[i*bs];
56: rowstart[i*bs+j+1] = rowstart[i*bs+j] + rowlengths[i*bs];
57: }
58: bi[i+1] = bi[i] + rowlengths[i*bs]/bs;
59: }
60: 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);
62: /* set b->j and b->a */
63: aj = a->j; av = a->a;
64: for (i=0; i<mbs; i++) {
65: nz = ai[i+1] - ai[i];
66: /* diagonal block */
67: if (nz && *aj == i) {
68: nz--;
69: for (j=0; j<bs; j++) { /* row i*bs+j */
70: itmp = i*bs+j;
71: for (k=0; k<bs; k++) { /* col i*bs+k */
72: *(bj + rowstart[itmp]) = (*aj)*bs+k;
73: *(bv + rowstart[itmp]) = *(av+k*bs+j);
74: rowstart[itmp]++;
75: }
76: }
77: aj++; av += bs2;
78: }
80: while (nz--) {
81: /* lower triangular blocks */
82: for (j=0; j<bs; j++) { /* row (*aj)*bs+j */
83: itmp = (*aj)*bs+j;
84: for (k=0; k<bs; k++) { /* col i*bs+k */
85: *(bj + rowstart[itmp]) = i*bs+k;
86: *(bv + rowstart[itmp]) = *(av+j*bs+k);
87: rowstart[itmp]++;
88: }
89: }
90: /* upper triangular blocks */
91: for (j=0; j<bs; j++) { /* row i*bs+j */
92: itmp = i*bs+j;
93: for (k=0; k<bs; k++) { /* col (*aj)*bs+k */
94: *(bj + rowstart[itmp]) = (*aj)*bs+k;
95: *(bv + rowstart[itmp]) = *(av+k*bs+j);
96: rowstart[itmp]++;
97: }
98: }
99: aj++; av += bs2;
100: }
101: }
102: PetscFree2(rowlengths,rowstart);
103: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
104: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
106: if (reuse == MAT_INPLACE_MATRIX) {
107: MatHeaderReplace(A,&B);
108: } else {
109: *newmat = B;
110: }
111: return(0);
112: }
114: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
115: {
116: Mat B;
117: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
118: Mat_SeqSBAIJ *b;
120: PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->N,i,j,*bi,*bj,*rowlengths;
121: MatScalar *av,*bv;
124: if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
125: if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
127: PetscMalloc1(m,&rowlengths);
128: for (i=0; i<m; i++) {
129: rowlengths[i] = ai[i+1] - a->diag[i];
130: }
131: if (reuse != MAT_REUSE_MATRIX) {
132: MatCreate(PetscObjectComm((PetscObject)A),&B);
133: MatSetSizes(B,m,n,m,n);
134: MatSetType(B,MATSEQSBAIJ);
135: MatSeqSBAIJSetPreallocation(B,1,0,rowlengths);
136: } else B = *newmat;
138: b = (Mat_SeqSBAIJ*)(B->data);
139: bi = b->i;
140: bj = b->j;
141: bv = b->a;
143: bi[0] = 0;
144: for (i=0; i<m; i++) {
145: aj = a->j + a->diag[i];
146: av = a->a + a->diag[i];
147: for (j=0; j<rowlengths[i]; j++) {
148: *bj = *aj; bj++; aj++;
149: *bv = *av; bv++; av++;
150: }
151: bi[i+1] = bi[i] + rowlengths[i];
152: b->ilen[i] = rowlengths[i];
153: }
155: PetscFree(rowlengths);
156: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
157: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
159: if (reuse == MAT_INPLACE_MATRIX) {
160: MatHeaderReplace(A,&B);
161: } else {
162: *newmat = B;
163: }
164: return(0);
165: }
167: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
168: {
169: Mat B;
170: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
171: Mat_SeqBAIJ *b;
173: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp;
174: PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,col,row;
175: MatScalar *av,*bv;
178: /* compute browlengths of newmat */
179: PetscMalloc2(mbs,&browlengths,mbs,&browstart);
180: for (i=0; i<mbs; i++) browlengths[i] = 0;
181: for (i=0; i<mbs; i++) {
182: nz = ai[i+1] - ai[i];
183: aj++; /* skip diagonal */
184: for (k=1; k<nz; k++) { /* no. of lower triangular blocks */
185: browlengths[*aj]++; aj++;
186: }
187: browlengths[i] += nz; /* no. of upper triangular blocks */
188: }
190: if (reuse != MAT_REUSE_MATRIX) {
191: MatCreate(PetscObjectComm((PetscObject)A),&B);
192: MatSetSizes(B,m,n,m,n);
193: MatSetType(B,MATSEQBAIJ);
194: MatSeqBAIJSetPreallocation(B,bs,0,browlengths);
195: } else B = *newmat;
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: if (reuse != MAT_REUSE_MATRIX) {
283: MatCreate(PetscObjectComm((PetscObject)A),&B);
284: MatSetSizes(B,m,n,m,n);
285: MatSetType(B,MATSEQSBAIJ);
286: MatSeqSBAIJSetPreallocation(B,bs,0,browlengths);
287: } else B = *newmat;
289: b = (Mat_SeqSBAIJ*)(B->data);
290: bi = b->i;
291: bj = b->j;
292: bv = b->a;
294: bi[0] = 0;
295: for (i=0; i<mbs; i++) {
296: aj = a->j + a->diag[i];
297: av = a->a + (a->diag[i])*bs2;
298: for (j=0; j<browlengths[i]; j++) {
299: *bj = *aj; bj++; aj++;
300: for (k=0; k<bs2; k++) {
301: *bv = *av; bv++; av++;
302: }
303: }
304: bi[i+1] = bi[i] + browlengths[i];
305: b->ilen[i] = browlengths[i];
306: }
307: PetscFree(browlengths);
308: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
309: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
311: if (reuse == MAT_INPLACE_MATRIX) {
312: MatHeaderReplace(A,&B);
313: } else {
314: *newmat = B;
315: }
316: return(0);
317: }