Actual source code: aijsbaij.c
petsc-3.14.6 2021-03-30
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;
15: #if defined(PETSC_USE_COMPLEX)
16: const int aconj = A->hermitian ? 1 : 0;
17: #else
18: const int aconj = 0;
19: #endif
22: /* compute rowlengths of newmat */
23: PetscMalloc2(m,&rowlengths,m+1,&rowstart);
25: for (i=0; i<mbs; i++) rowlengths[i*bs] = 0;
26: k = 0;
27: for (i=0; i<mbs; i++) {
28: nz = ai[i+1] - ai[i];
29: if (nz) {
30: rowlengths[k] += nz; /* no. of upper triangular blocks */
31: if (*aj == i) {aj++;diagcnt++;nz--;} /* skip diagonal */
32: for (j=0; j<nz; j++) { /* no. of lower triangular blocks */
33: rowlengths[(*aj)*bs]++; aj++;
34: }
35: }
36: rowlengths[k] *= bs;
37: for (j=1; j<bs; j++) {
38: rowlengths[k+j] = rowlengths[k];
39: }
40: k += bs;
41: }
43: if (reuse != MAT_REUSE_MATRIX) {
44: MatCreate(PetscObjectComm((PetscObject)A),&B);
45: MatSetSizes(B,m,n,m,n);
46: MatSetType(B,MATSEQAIJ);
47: MatSeqAIJSetPreallocation(B,0,rowlengths);
48: MatSetBlockSize(B,A->rmap->bs);
49: } else B = *newmat;
51: b = (Mat_SeqAIJ*)(B->data);
52: bi = b->i;
53: bj = b->j;
54: bv = b->a;
56: /* set b->i */
57: bi[0] = 0; rowstart[0] = 0;
58: for (i=0; i<mbs; i++) {
59: for (j=0; j<bs; j++) {
60: b->ilen[i*bs+j] = rowlengths[i*bs];
61: rowstart[i*bs+j+1] = rowstart[i*bs+j] + rowlengths[i*bs];
62: }
63: bi[i+1] = bi[i] + rowlengths[i*bs]/bs;
64: }
65: 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);
67: /* set b->j and b->a */
68: aj = a->j; av = a->a;
69: for (i=0; i<mbs; i++) {
70: nz = ai[i+1] - ai[i];
71: /* diagonal block */
72: if (nz && *aj == i) {
73: nz--;
74: for (j=0; j<bs; j++) { /* row i*bs+j */
75: itmp = i*bs+j;
76: for (k=0; k<bs; k++) { /* col i*bs+k */
77: *(bj + rowstart[itmp]) = (*aj)*bs+k;
78: *(bv + rowstart[itmp]) = *(av+k*bs+j);
79: rowstart[itmp]++;
80: }
81: }
82: aj++; av += bs2;
83: }
85: while (nz--) {
86: /* lower triangular blocks */
87: for (j=0; j<bs; j++) { /* row (*aj)*bs+j */
88: itmp = (*aj)*bs+j;
89: for (k=0; k<bs; k++) { /* col i*bs+k */
90: *(bj + rowstart[itmp]) = i*bs+k;
91: *(bv + rowstart[itmp]) = aconj ? PetscConj(*(av+j*bs+k)) : *(av+j*bs+k);
92: rowstart[itmp]++;
93: }
94: }
95: /* upper triangular blocks */
96: for (j=0; j<bs; j++) { /* row i*bs+j */
97: itmp = i*bs+j;
98: for (k=0; k<bs; k++) { /* col (*aj)*bs+k */
99: *(bj + rowstart[itmp]) = (*aj)*bs+k;
100: *(bv + rowstart[itmp]) = *(av+k*bs+j);
101: rowstart[itmp]++;
102: }
103: }
104: aj++; av += bs2;
105: }
106: }
107: PetscFree2(rowlengths,rowstart);
108: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
109: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
111: if (reuse == MAT_INPLACE_MATRIX) {
112: MatHeaderReplace(A,&B);
113: } else {
114: *newmat = B;
115: }
116: return(0);
117: }
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,bs=PetscAbs(A->rmap->bs);
126: MatScalar *av,*bv;
127: PetscBool miss = PETSC_FALSE;
130: #if !defined(PETSC_USE_COMPLEX)
131: if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
132: #else
133: if (!A->symmetric && !A->hermitian) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be either symmetric or hermitian. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE) and/or MatSetOption(mat,MAT_HERMITIAN,PETSC_TRUE)");
134: #endif
135: if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
137: PetscMalloc1(m/bs,&rowlengths);
138: for (i=0; i<m/bs; i++) {
139: if (a->diag[i*bs] == ai[i*bs+1]) { /* missing diagonal */
140: rowlengths[i] = (ai[i*bs+1] - ai[i*bs])/bs; /* allocate some extra space */
141: miss = PETSC_TRUE;
142: } else {
143: rowlengths[i] = (ai[i*bs+1] - a->diag[i*bs])/bs;
144: }
145: }
146: if (reuse != MAT_REUSE_MATRIX) {
147: MatCreate(PetscObjectComm((PetscObject)A),&B);
148: MatSetSizes(B,m,n,m,n);
149: MatSetType(B,MATSEQSBAIJ);
150: MatSeqSBAIJSetPreallocation(B,bs,0,rowlengths);
151: } else B = *newmat;
153: if (bs == 1 && !miss) {
154: b = (Mat_SeqSBAIJ*)(B->data);
155: bi = b->i;
156: bj = b->j;
157: bv = b->a;
159: bi[0] = 0;
160: for (i=0; i<m; i++) {
161: aj = a->j + a->diag[i];
162: av = a->a + a->diag[i];
163: for (j=0; j<rowlengths[i]; j++) {
164: *bj = *aj; bj++; aj++;
165: *bv = *av; bv++; av++;
166: }
167: bi[i+1] = bi[i] + rowlengths[i];
168: b->ilen[i] = rowlengths[i];
169: }
170: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
171: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
172: } else {
173: MatSetOption(B,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);
174: /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
175: /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before */
176: /* MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below */
177: MatConvert_Basic(A,newtype,MAT_REUSE_MATRIX,&B);
178: }
179: PetscFree(rowlengths);
180: if (reuse == MAT_INPLACE_MATRIX) {
181: MatHeaderReplace(A,&B);
182: } else *newmat = B;
183: return(0);
184: }
186: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
187: {
188: Mat B;
189: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
190: Mat_SeqBAIJ *b;
192: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp;
193: PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,col,row;
194: MatScalar *av,*bv;
195: #if defined(PETSC_USE_COMPLEX)
196: const int aconj = A->hermitian ? 1 : 0;
197: #else
198: const int aconj = 0;
199: #endif
202: /* compute browlengths of newmat */
203: PetscMalloc2(mbs,&browlengths,mbs,&browstart);
204: for (i=0; i<mbs; i++) browlengths[i] = 0;
205: for (i=0; i<mbs; i++) {
206: nz = ai[i+1] - ai[i];
207: aj++; /* skip diagonal */
208: for (k=1; k<nz; k++) { /* no. of lower triangular blocks */
209: browlengths[*aj]++; aj++;
210: }
211: browlengths[i] += nz; /* no. of upper triangular blocks */
212: }
214: if (reuse != MAT_REUSE_MATRIX) {
215: MatCreate(PetscObjectComm((PetscObject)A),&B);
216: MatSetSizes(B,m,n,m,n);
217: MatSetType(B,MATSEQBAIJ);
218: MatSeqBAIJSetPreallocation(B,bs,0,browlengths);
219: } else B = *newmat;
221: b = (Mat_SeqBAIJ*)(B->data);
222: bi = b->i;
223: bj = b->j;
224: bv = b->a;
226: /* set b->i */
227: bi[0] = 0;
228: for (i=0; i<mbs; i++) {
229: b->ilen[i] = browlengths[i];
230: bi[i+1] = bi[i] + browlengths[i];
231: browstart[i] = bi[i];
232: }
233: 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);
235: /* set b->j and b->a */
236: aj = a->j; av = a->a;
237: for (i=0; i<mbs; i++) {
238: /* diagonal block */
239: *(bj + browstart[i]) = *aj; aj++;
241: itmp = bs2*browstart[i];
242: for (k=0; k<bs2; k++) {
243: *(bv + itmp + k) = *av; av++;
244: }
245: browstart[i]++;
247: nz = ai[i+1] - ai[i] -1;
248: while (nz--) {
249: /* lower triangular blocks - transpose blocks of A */
250: *(bj + browstart[*aj]) = i; /* block col index */
252: itmp = bs2*browstart[*aj]; /* row index */
253: for (col=0; col<bs; col++) {
254: k = col;
255: for (row=0; row<bs; row++) {
256: bv[itmp + col*bs+row] = aconj ? PetscConj(av[k]) : av[k];
257: k+=bs;
258: }
259: }
260: browstart[*aj]++;
262: /* upper triangular blocks */
263: *(bj + browstart[i]) = *aj; aj++;
265: itmp = bs2*browstart[i];
266: for (k=0; k<bs2; k++) {
267: bv[itmp + k] = av[k];
268: }
269: av += bs2;
270: browstart[i]++;
271: }
272: }
273: PetscFree2(browlengths,browstart);
274: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
275: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
277: if (reuse == MAT_INPLACE_MATRIX) {
278: MatHeaderReplace(A,&B);
279: } else *newmat = B;
280: return(0);
281: }
283: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
284: {
285: Mat B;
286: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
287: Mat_SeqSBAIJ *b;
289: PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*browlengths;
290: PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,dd;
291: MatScalar *av,*bv;
292: PetscBool flg;
295: if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
296: if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
297: MatMissingDiagonal_SeqBAIJ(A,&flg,&dd); /* check for missing diagonals, then mark diag */
298: if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal %D",dd);
300: PetscMalloc1(mbs,&browlengths);
301: for (i=0; i<mbs; i++) {
302: browlengths[i] = ai[i+1] - a->diag[i];
303: }
305: if (reuse != MAT_REUSE_MATRIX) {
306: MatCreate(PetscObjectComm((PetscObject)A),&B);
307: MatSetSizes(B,m,n,m,n);
308: MatSetType(B,MATSEQSBAIJ);
309: MatSeqSBAIJSetPreallocation(B,bs,0,browlengths);
310: } else B = *newmat;
312: b = (Mat_SeqSBAIJ*)(B->data);
313: bi = b->i;
314: bj = b->j;
315: bv = b->a;
317: bi[0] = 0;
318: for (i=0; i<mbs; i++) {
319: aj = a->j + a->diag[i];
320: av = a->a + (a->diag[i])*bs2;
321: for (j=0; j<browlengths[i]; j++) {
322: *bj = *aj; bj++; aj++;
323: for (k=0; k<bs2; k++) {
324: *bv = *av; bv++; av++;
325: }
326: }
327: bi[i+1] = bi[i] + browlengths[i];
328: b->ilen[i] = browlengths[i];
329: }
330: PetscFree(browlengths);
331: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
332: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
334: if (reuse == MAT_INPLACE_MATRIX) {
335: MatHeaderReplace(A,&B);
336: } else *newmat = B;
337: return(0);
338: }