Actual source code: aijsbaij.c
petsc-3.4.5 2014-06-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>
8: PETSC_EXTERN 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,PetscInt,&rowlengths,m+1,PetscInt,&rowstart);
22: for (i=0; i<mbs; i++) rowlengths[i*bs] = 0;
23: aj = a->j;
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: /* printf(" rowlengths[%d]: %d\n",i, rowlengths[i]); */
40: }
42: MatCreate(PetscObjectComm((PetscObject)A),&B);
43: MatSetSizes(B,m,n,m,n);
44: MatSetType(B,MATSEQAIJ);
45: MatSeqAIJSetPreallocation(B,0,rowlengths);
46: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_FALSE);
48: B->rmap->bs = A->rmap->bs;
50: b = (Mat_SeqAIJ*)(B->data);
51: bi = b->i;
52: bj = b->j;
53: bv = b->a;
55: /* set b->i */
56: bi[0] = 0; rowstart[0] = 0;
57: for (i=0; i<mbs; i++) {
58: for (j=0; j<bs; j++) {
59: b->ilen[i*bs+j] = rowlengths[i*bs];
60: rowstart[i*bs+j+1] = rowstart[i*bs+j] + rowlengths[i*bs];
61: }
62: bi[i+1] = bi[i] + rowlengths[i*bs]/bs;
63: }
64: 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);
66: /* set b->j and b->a */
67: aj = a->j; av = a->a;
68: for (i=0; i<mbs; i++) {
69: nz = ai[i+1] - ai[i];
70: /* diagonal block */
71: if (nz && *aj == i) {
72: nz--;
73: for (j=0; j<bs; j++) { /* row i*bs+j */
74: itmp = i*bs+j;
75: for (k=0; k<bs; k++) { /* col i*bs+k */
76: *(bj + rowstart[itmp]) = (*aj)*bs+k;
77: *(bv + rowstart[itmp]) = *(av+k*bs+j);
78: rowstart[itmp]++;
79: }
80: }
81: aj++; av += bs2;
82: }
84: while (nz--) {
85: /* lower triangular blocks */
86: for (j=0; j<bs; j++) { /* row (*aj)*bs+j */
87: itmp = (*aj)*bs+j;
88: for (k=0; k<bs; k++) { /* col i*bs+k */
89: *(bj + rowstart[itmp]) = i*bs+k;
90: *(bv + rowstart[itmp]) = *(av+j*bs+k);
91: rowstart[itmp]++;
92: }
93: }
94: /* upper triangular blocks */
95: for (j=0; j<bs; j++) { /* row i*bs+j */
96: itmp = i*bs+j;
97: for (k=0; k<bs; k++) { /* col (*aj)*bs+k */
98: *(bj + rowstart[itmp]) = (*aj)*bs+k;
99: *(bv + rowstart[itmp]) = *(av+k*bs+j);
100: rowstart[itmp]++;
101: }
102: }
103: aj++; av += bs2;
104: }
105: }
106: PetscFree2(rowlengths,rowstart);
107: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
108: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
110: if (reuse == MAT_REUSE_MATRIX) {
111: MatHeaderReplace(A,B);
112: } else {
113: *newmat = B;
114: }
115: return(0);
116: }
120: PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
121: {
122: Mat B;
123: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
124: Mat_SeqSBAIJ *b;
126: PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->N,i,j,*bi,*bj,*rowlengths;
127: MatScalar *av,*bv;
130: if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
131: if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
133: PetscMalloc(m*sizeof(PetscInt),&rowlengths);
134: for (i=0; i<m; i++) {
135: rowlengths[i] = ai[i+1] - a->diag[i];
136: }
137: MatCreate(PetscObjectComm((PetscObject)A),&B);
138: MatSetSizes(B,m,n,m,n);
139: MatSetType(B,MATSEQSBAIJ);
140: MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,1,0,rowlengths);
142: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
144: b = (Mat_SeqSBAIJ*)(B->data);
145: bi = b->i;
146: bj = b->j;
147: bv = b->a;
149: bi[0] = 0;
150: for (i=0; i<m; i++) {
151: aj = a->j + a->diag[i];
152: av = a->a + a->diag[i];
153: for (j=0; j<rowlengths[i]; j++) {
154: *bj = *aj; bj++; aj++;
155: *bv = *av; bv++; av++;
156: }
157: bi[i+1] = bi[i] + rowlengths[i];
158: b->ilen[i] = rowlengths[i];
159: }
161: PetscFree(rowlengths);
162: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
163: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
165: if (reuse == MAT_REUSE_MATRIX) {
166: MatHeaderReplace(A,B);
167: } else {
168: *newmat = B;
169: }
170: return(0);
171: }
175: PETSC_EXTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
176: {
177: Mat B;
178: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
179: Mat_SeqBAIJ *b;
181: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp;
182: PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,col,row;
183: MatScalar *av,*bv;
186: /* compute browlengths of newmat */
187: PetscMalloc2(mbs,PetscInt,&browlengths,mbs,PetscInt,&browstart);
188: for (i=0; i<mbs; i++) browlengths[i] = 0;
189: aj = a->j;
190: for (i=0; i<mbs; i++) {
191: nz = ai[i+1] - ai[i];
192: aj++; /* skip diagonal */
193: for (k=1; k<nz; k++) { /* no. of lower triangular blocks */
194: browlengths[*aj]++; aj++;
195: }
196: browlengths[i] += nz; /* no. of upper triangular blocks */
197: }
199: MatCreate(PetscObjectComm((PetscObject)A),&B);
200: MatSetSizes(B,m,n,m,n);
201: MatSetType(B,MATSEQBAIJ);
202: MatSeqBAIJSetPreallocation(B,bs,0,browlengths);
203: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
205: b = (Mat_SeqBAIJ*)(B->data);
206: bi = b->i;
207: bj = b->j;
208: bv = b->a;
210: /* set b->i */
211: bi[0] = 0;
212: for (i=0; i<mbs; i++) {
213: b->ilen[i] = browlengths[i];
214: bi[i+1] = bi[i] + browlengths[i];
215: browstart[i] = bi[i];
216: }
217: 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);
219: /* set b->j and b->a */
220: aj = a->j; av = a->a;
221: for (i=0; i<mbs; i++) {
222: /* diagonal block */
223: *(bj + browstart[i]) = *aj; aj++;
225: itmp = bs2*browstart[i];
226: for (k=0; k<bs2; k++) {
227: *(bv + itmp + k) = *av; av++;
228: }
229: browstart[i]++;
231: nz = ai[i+1] - ai[i] -1;
232: while (nz--) {
233: /* lower triangular blocks - transpose blocks of A */
234: *(bj + browstart[*aj]) = i; /* block col index */
236: itmp = bs2*browstart[*aj]; /* row index */
237: for (col=0; col<bs; col++) {
238: k = col;
239: for (row=0; row<bs; row++) {
240: bv[itmp + col*bs+row] = av[k]; k+=bs;
241: }
242: }
243: browstart[*aj]++;
245: /* upper triangular blocks */
246: *(bj + browstart[i]) = *aj; aj++;
248: itmp = bs2*browstart[i];
249: for (k=0; k<bs2; k++) {
250: bv[itmp + k] = av[k];
251: }
252: av += bs2;
253: browstart[i]++;
254: }
255: }
256: PetscFree2(browlengths,browstart);
257: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
258: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
260: if (reuse == MAT_REUSE_MATRIX) {
261: MatHeaderReplace(A,B);
262: } else {
263: *newmat = B;
264: }
265: return(0);
266: }
270: PETSC_EXTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
271: {
272: Mat B;
273: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
274: Mat_SeqSBAIJ *b;
276: PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*browlengths;
277: PetscInt bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,dd;
278: MatScalar *av,*bv;
279: PetscBool flg;
282: if (!A->symmetric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
283: if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
284: MatMissingDiagonal_SeqBAIJ(A,&flg,&dd); /* check for missing diagonals, then mark diag */
285: if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal %D",dd);
287: PetscMalloc(mbs*sizeof(PetscInt),&browlengths);
288: for (i=0; i<mbs; i++) {
289: browlengths[i] = ai[i+1] - a->diag[i];
290: }
292: MatCreate(PetscObjectComm((PetscObject)A),&B);
293: MatSetSizes(B,m,n,m,n);
294: MatSetType(B,MATSEQSBAIJ);
295: MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,browlengths);
296: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
298: b = (Mat_SeqSBAIJ*)(B->data);
299: bi = b->i;
300: bj = b->j;
301: bv = b->a;
303: bi[0] = 0;
304: for (i=0; i<mbs; i++) {
305: aj = a->j + a->diag[i];
306: av = a->a + (a->diag[i])*bs2;
307: for (j=0; j<browlengths[i]; j++) {
308: *bj = *aj; bj++; aj++;
309: for (k=0; k<bs2; k++) {
310: *bv = *av; bv++; av++;
311: }
312: }
313: bi[i+1] = bi[i] + browlengths[i];
314: b->ilen[i] = browlengths[i];
315: }
316: PetscFree(browlengths);
317: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
318: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
320: if (reuse == MAT_REUSE_MATRIX) {
321: MatHeaderReplace(A,B);
322: } else {
323: *newmat = B;
324: }
325: return(0);
326: }