Actual source code: aijsbaij.c
petsc-3.3-p7 2013-05-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: EXTERN_C_BEGIN
9: PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
10: {
11: Mat B;
12: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
13: Mat_SeqAIJ *b;
15: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*rowlengths,nz,*rowstart,itmp;
16: PetscInt bs=A->rmap->bs,bs2=bs*bs,mbs=A->rmap->N/bs,diagcnt=0;
17: MatScalar *av,*bv;
20: /* compute rowlengths of newmat */
21: PetscMalloc2(m,PetscInt,&rowlengths,m+1,PetscInt,&rowstart);
22:
23: for (i=0; i<mbs; i++) rowlengths[i*bs] = 0;
24: aj = a->j;
25: k = 0;
26: for (i=0; i<mbs; i++) {
27: nz = ai[i+1] - ai[i];
28: if (nz) {
29: rowlengths[k] += nz; /* no. of upper triangular blocks */
30: if (*aj == i) {aj++;diagcnt++;nz--;} /* skip diagonal */
31: for (j=0; j<nz; j++) { /* no. of lower triangular blocks */
32: rowlengths[(*aj)*bs]++; aj++;
33: }
34: }
35: rowlengths[k] *= bs;
36: for (j=1; j<bs; j++) {
37: rowlengths[k+j] = rowlengths[k];
38: }
39: k += bs;
40: /* printf(" rowlengths[%d]: %d\n",i, rowlengths[i]); */
41: }
42:
43: MatCreate(((PetscObject)A)->comm,&B);
44: MatSetSizes(B,m,n,m,n);
45: MatSetType(B,MATSEQAIJ);
46: MatSeqAIJSetPreallocation(B,0,rowlengths);
47: 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: }
83:
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: }
117: EXTERN_C_END
119: EXTERN_C_BEGIN
122: PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A,const MatType newtype,MatReuse reuse,Mat *newmat) {
123: Mat B;
124: Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data;
125: Mat_SeqSBAIJ *b;
127: PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->N,i,j,*bi,*bj,*rowlengths;
128: MatScalar *av,*bv;
131: if (!A->symmetric) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
132: if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
134: PetscMalloc(m*sizeof(PetscInt),&rowlengths);
135: for (i=0; i<m; i++) {
136: rowlengths[i] = ai[i+1] - a->diag[i];
137: }
138: MatCreate(((PetscObject)A)->comm,&B);
139: MatSetSizes(B,m,n,m,n);
140: MatSetType(B,MATSEQSBAIJ);
141: MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,1,0,rowlengths);
143: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
144:
145: b = (Mat_SeqSBAIJ*)(B->data);
146: bi = b->i;
147: bj = b->j;
148: bv = b->a;
150: bi[0] = 0;
151: for (i=0; i<m; i++) {
152: aj = a->j + a->diag[i];
153: av = a->a + a->diag[i];
154: for (j=0; j<rowlengths[i]; j++){
155: *bj = *aj; bj++; aj++;
156: *bv = *av; bv++; av++;
157: }
158: bi[i+1] = bi[i] + rowlengths[i];
159: b->ilen[i] = rowlengths[i];
160: }
161:
162: PetscFree(rowlengths);
163: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
164: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
166: if (reuse == MAT_REUSE_MATRIX) {
167: MatHeaderReplace(A,B);
168: } else {
169: *newmat = B;
170: }
171: return(0);
172: }
173: EXTERN_C_END
175: EXTERN_C_BEGIN
178: PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
179: {
180: Mat B;
181: Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data;
182: Mat_SeqBAIJ *b;
184: PetscInt *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp;
185: PetscInt bs=A->rmap->bs,bs2=bs*bs,mbs=m/bs,col,row;
186: MatScalar *av,*bv;
189: /* compute browlengths of newmat */
190: PetscMalloc2(mbs,PetscInt,&browlengths,mbs,PetscInt,&browstart);
191: for (i=0; i<mbs; i++) browlengths[i] = 0;
192: aj = a->j;
193: for (i=0; i<mbs; i++) {
194: nz = ai[i+1] - ai[i];
195: aj++; /* skip diagonal */
196: for (k=1; k<nz; k++) { /* no. of lower triangular blocks */
197: browlengths[*aj]++; aj++;
198: }
199: browlengths[i] += nz; /* no. of upper triangular blocks */
200: }
201:
202: MatCreate(((PetscObject)A)->comm,&B);
203: MatSetSizes(B,m,n,m,n);
204: MatSetType(B,MATSEQBAIJ);
205: MatSeqBAIJSetPreallocation(B,bs,0,browlengths);
206: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
207:
208: b = (Mat_SeqBAIJ*)(B->data);
209: bi = b->i;
210: bj = b->j;
211: bv = b->a;
213: /* set b->i */
214: bi[0] = 0;
215: for (i=0; i<mbs; i++){
216: b->ilen[i] = browlengths[i];
217: bi[i+1] = bi[i] + browlengths[i];
218: browstart[i] = bi[i];
219: }
220: 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);
221:
222: /* set b->j and b->a */
223: aj = a->j; av = a->a;
224: for (i=0; i<mbs; i++) {
225: /* diagonal block */
226: *(bj + browstart[i]) = *aj; aj++;
227: itmp = bs2*browstart[i];
228: for (k=0; k<bs2; k++){
229: *(bv + itmp + k) = *av; av++;
230: }
231: browstart[i]++;
232:
233: nz = ai[i+1] - ai[i] -1;
234: while (nz--){
235: /* lower triangular blocks - transpose blocks of A */
236: *(bj + browstart[*aj]) = i; /* block col index */
237: itmp = bs2*browstart[*aj]; /* row index */
238: for (col=0; col<bs; col++){
239: k = col;
240: for (row=0; row<bs; row++){
241: bv[itmp + col*bs+row] = av[k]; k+=bs;
242: }
243: }
244: browstart[*aj]++;
246: /* upper triangular blocks */
247: *(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: }
267: EXTERN_C_END
269: EXTERN_C_BEGIN
272: PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
273: {
274: Mat B;
275: Mat_SeqBAIJ *a = (Mat_SeqBAIJ*)A->data;
276: Mat_SeqSBAIJ *b;
278: PetscInt *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*browlengths;
279: PetscInt bs=A->rmap->bs,bs2=bs*bs,mbs=m/bs,dd;
280: MatScalar *av,*bv;
281: PetscBool flg;
284: if (!A->symmetric) SETERRQ(((PetscObject)A)->comm,PETSC_ERR_USER,"Matrix must be symmetric. Call MatSetOption(mat,MAT_SYMMETRIC,PETSC_TRUE)");
285: if (n != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Matrix must be square");
286: MatMissingDiagonal_SeqBAIJ(A,&flg,&dd); /* check for missing diagonals, then mark diag */
287: if (flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal %D",dd);
289: PetscMalloc(mbs*sizeof(PetscInt),&browlengths);
290: for (i=0; i<mbs; i++) {
291: browlengths[i] = ai[i+1] - a->diag[i];
292: }
294: MatCreate(((PetscObject)A)->comm,&B);
295: MatSetSizes(B,m,n,m,n);
296: MatSetType(B,MATSEQSBAIJ);
297: MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,browlengths);
298: MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
300: b = (Mat_SeqSBAIJ*)(B->data);
301: bi = b->i;
302: bj = b->j;
303: bv = b->a;
305: bi[0] = 0;
306: for (i=0; i<mbs; i++) {
307: aj = a->j + a->diag[i];
308: av = a->a + (a->diag[i])*bs2;
309: for (j=0; j<browlengths[i]; j++){
310: *bj = *aj; bj++; aj++;
311: for (k=0; k<bs2; k++){
312: *bv = *av; bv++; av++;
313: }
314: }
315: bi[i+1] = bi[i] + browlengths[i];
316: b->ilen[i] = browlengths[i];
317: }
318: PetscFree(browlengths);
319: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
320: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
322: if (reuse == MAT_REUSE_MATRIX) {
323: MatHeaderReplace(A,B);
324: } else {
325: *newmat = B;
326: }
328: return(0);
329: }
330: EXTERN_C_END