Actual source code: aijsbaij.c


  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;
 11:   PetscInt       *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*rowlengths,nz,*rowstart,itmp;
 12:   PetscInt       bs =A->rmap->bs,bs2=bs*bs,mbs=A->rmap->N/bs,diagcnt=0;
 13:   MatScalar      *av,*bv;
 14: #if defined(PETSC_USE_COMPLEX)
 15:   const int      aconj = A->hermitian ? 1 : 0;
 16: #else
 17:   const int      aconj = 0;
 18: #endif

 20:   /* compute rowlengths of newmat */
 21:   PetscMalloc2(m,&rowlengths,m+1,&rowstart);

 23:   for (i=0; i<mbs; i++) rowlengths[i*bs] = 0;
 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:   }

 41:   if (reuse != MAT_REUSE_MATRIX) {
 42:     MatCreate(PetscObjectComm((PetscObject)A),&B);
 43:     MatSetSizes(B,m,n,m,n);
 44:     MatSetType(B,MATSEQAIJ);
 45:     MatSeqAIJSetPreallocation(B,0,rowlengths);
 46:     MatSetBlockSize(B,A->rmap->bs);
 47:   } else B = *newmat;

 49:   b  = (Mat_SeqAIJ*)(B->data);
 50:   bi = b->i;
 51:   bj = b->j;
 52:   bv = b->a;

 54:   /* set b->i */
 55:   bi[0] = 0; rowstart[0] = 0;
 56:   for (i=0; i<mbs; i++) {
 57:     for (j=0; j<bs; j++) {
 58:       b->ilen[i*bs+j]    = rowlengths[i*bs];
 59:       rowstart[i*bs+j+1] = rowstart[i*bs+j] + rowlengths[i*bs];
 60:     }
 61:     bi[i+1] = bi[i] + rowlengths[i*bs]/bs;
 62:   }

 65:   /* set b->j and b->a */
 66:   aj = a->j; av = a->a;
 67:   for (i=0; i<mbs; i++) {
 68:     nz = ai[i+1] - ai[i];
 69:     /* diagonal block */
 70:     if (nz && *aj == i) {
 71:       nz--;
 72:       for (j=0; j<bs; j++) {   /* row i*bs+j */
 73:         itmp = i*bs+j;
 74:         for (k=0; k<bs; k++) { /* col i*bs+k */
 75:           *(bj + rowstart[itmp]) = (*aj)*bs+k;
 76:           *(bv + rowstart[itmp]) = *(av+k*bs+j);
 77:           rowstart[itmp]++;
 78:         }
 79:       }
 80:       aj++; av += bs2;
 81:     }

 83:     while (nz--) {
 84:       /* lower triangular blocks */
 85:       for (j=0; j<bs; j++) {   /* row (*aj)*bs+j */
 86:         itmp = (*aj)*bs+j;
 87:         for (k=0; k<bs; k++) { /* col i*bs+k */
 88:           *(bj + rowstart[itmp]) = i*bs+k;
 89:           *(bv + rowstart[itmp]) = aconj ? PetscConj(*(av+j*bs+k)) : *(av+j*bs+k);
 90:           rowstart[itmp]++;
 91:         }
 92:       }
 93:       /* upper triangular blocks */
 94:       for (j=0; j<bs; j++) {   /* row i*bs+j */
 95:         itmp = i*bs+j;
 96:         for (k=0; k<bs; k++) { /* col (*aj)*bs+k */
 97:           *(bj + rowstart[itmp]) = (*aj)*bs+k;
 98:           *(bv + rowstart[itmp]) = *(av+k*bs+j);
 99:           rowstart[itmp]++;
100:         }
101:       }
102:       aj++; av += bs2;
103:     }
104:   }
105:   PetscFree2(rowlengths,rowstart);
106:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
107:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

109:   if (reuse == MAT_INPLACE_MATRIX) {
110:     MatHeaderReplace(A,&B);
111:   } else {
112:     *newmat = B;
113:   }
114:   return 0;
115: }

117: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
118: {
119:   Mat            B;
120:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
121:   Mat_SeqSBAIJ   *b;
122:   PetscInt       *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->N,i,j,*bi,*bj,*rowlengths,bs=PetscAbs(A->rmap->bs);
123:   MatScalar      *av,*bv;
124:   PetscBool      miss = PETSC_FALSE;

126: #if !defined(PETSC_USE_COMPLEX)
128: #else
130: #endif

133:   PetscMalloc1(m/bs,&rowlengths);
134:   for (i=0; i<m/bs; i++) {
135:     if (a->diag[i*bs] == ai[i*bs+1]) { /* missing diagonal */
136:       rowlengths[i] = (ai[i*bs+1] - ai[i*bs])/bs; /* allocate some extra space */
137:       miss = PETSC_TRUE;
138:     } else {
139:       rowlengths[i] = (ai[i*bs+1] - a->diag[i*bs])/bs;
140:     }
141:   }
142:   if (reuse != MAT_REUSE_MATRIX) {
143:     MatCreate(PetscObjectComm((PetscObject)A),&B);
144:     MatSetSizes(B,m,n,m,n);
145:     MatSetType(B,MATSEQSBAIJ);
146:     MatSeqSBAIJSetPreallocation(B,bs,0,rowlengths);
147:   } else B = *newmat;

149:   if (bs == 1 && !miss) {
150:     b  = (Mat_SeqSBAIJ*)(B->data);
151:     bi = b->i;
152:     bj = b->j;
153:     bv = b->a;

155:     bi[0] = 0;
156:     for (i=0; i<m; i++) {
157:       aj = a->j + a->diag[i];
158:       av = a->a + a->diag[i];
159:       for (j=0; j<rowlengths[i]; j++) {
160:         *bj = *aj; bj++; aj++;
161:         *bv = *av; bv++; av++;
162:       }
163:       bi[i+1]    = bi[i] + rowlengths[i];
164:       b->ilen[i] = rowlengths[i];
165:     }
166:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
167:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
168:   } else {
169:     MatSetOption(B,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE);
170:     /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
171:     /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before                                          */
172:     /*                      MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below                                        */
173:     MatConvert_Basic(A,newtype,MAT_REUSE_MATRIX,&B);
174:   }
175:   PetscFree(rowlengths);
176:   if (reuse == MAT_INPLACE_MATRIX) {
177:     MatHeaderReplace(A,&B);
178:   } else *newmat = B;
179:   return 0;
180: }

182: PETSC_INTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
183: {
184:   Mat            B;
185:   Mat_SeqSBAIJ   *a = (Mat_SeqSBAIJ*)A->data;
186:   Mat_SeqBAIJ    *b;
187:   PetscInt       *ai=a->i,*aj=a->j,m=A->rmap->N,n=A->cmap->n,i,k,*bi,*bj,*browlengths,nz,*browstart,itmp;
188:   PetscInt       bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,col,row;
189:   MatScalar      *av,*bv;
190: #if defined(PETSC_USE_COMPLEX)
191:   const int      aconj = A->hermitian ? 1 : 0;
192: #else
193:   const int      aconj = 0;
194: #endif

196:   /* compute browlengths of newmat */
197:   PetscMalloc2(mbs,&browlengths,mbs,&browstart);
198:   for (i=0; i<mbs; i++) browlengths[i] = 0;
199:   for (i=0; i<mbs; i++) {
200:     nz = ai[i+1] - ai[i];
201:     aj++; /* skip diagonal */
202:     for (k=1; k<nz; k++) { /* no. of lower triangular blocks */
203:       browlengths[*aj]++; aj++;
204:     }
205:     browlengths[i] += nz;   /* no. of upper triangular blocks */
206:   }

208:   if (reuse != MAT_REUSE_MATRIX) {
209:     MatCreate(PetscObjectComm((PetscObject)A),&B);
210:     MatSetSizes(B,m,n,m,n);
211:     MatSetType(B,MATSEQBAIJ);
212:     MatSeqBAIJSetPreallocation(B,bs,0,browlengths);
213:   } else B = *newmat;

215:   b  = (Mat_SeqBAIJ*)(B->data);
216:   bi = b->i;
217:   bj = b->j;
218:   bv = b->a;

220:   /* set b->i */
221:   bi[0] = 0;
222:   for (i=0; i<mbs; i++) {
223:     b->ilen[i]   = browlengths[i];
224:     bi[i+1]      = bi[i] + browlengths[i];
225:     browstart[i] = bi[i];
226:   }

229:   /* set b->j and b->a */
230:   aj = a->j; av = a->a;
231:   for (i=0; i<mbs; i++) {
232:     /* diagonal block */
233:     *(bj + browstart[i]) = *aj; aj++;

235:     itmp = bs2*browstart[i];
236:     for (k=0; k<bs2; k++) {
237:       *(bv + itmp + k) = *av; av++;
238:     }
239:     browstart[i]++;

241:     nz = ai[i+1] - ai[i] -1;
242:     while (nz--) {
243:       /* lower triangular blocks - transpose blocks of A */
244:       *(bj + browstart[*aj]) = i; /* block col index */

246:       itmp = bs2*browstart[*aj];  /* row index */
247:       for (col=0; col<bs; col++) {
248:         k = col;
249:         for (row=0; row<bs; row++) {
250:           bv[itmp + col*bs+row] = aconj ? PetscConj(av[k]) : av[k];
251:           k+=bs;
252:         }
253:       }
254:       browstart[*aj]++;

256:       /* upper triangular blocks */
257:       *(bj + browstart[i]) = *aj; aj++;

259:       itmp = bs2*browstart[i];
260:       for (k=0; k<bs2; k++) {
261:         bv[itmp + k] = av[k];
262:       }
263:       av += bs2;
264:       browstart[i]++;
265:     }
266:   }
267:   PetscFree2(browlengths,browstart);
268:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
269:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

271:   if (reuse == MAT_INPLACE_MATRIX) {
272:     MatHeaderReplace(A,&B);
273:   } else *newmat = B;
274:   return 0;
275: }

277: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqSBAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
278: {
279:   Mat            B;
280:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
281:   Mat_SeqSBAIJ   *b;
282:   PetscInt       *ai=a->i,*aj,m=A->rmap->N,n=A->cmap->n,i,j,k,*bi,*bj,*browlengths;
283:   PetscInt       bs =A->rmap->bs,bs2=bs*bs,mbs=m/bs,dd;
284:   MatScalar      *av,*bv;
285:   PetscBool      flg;

289:   MatMissingDiagonal_SeqBAIJ(A,&flg,&dd); /* check for missing diagonals, then mark diag */

292:   PetscMalloc1(mbs,&browlengths);
293:   for (i=0; i<mbs; i++) {
294:     browlengths[i] = ai[i+1] - a->diag[i];
295:   }

297:   if (reuse != MAT_REUSE_MATRIX) {
298:     MatCreate(PetscObjectComm((PetscObject)A),&B);
299:     MatSetSizes(B,m,n,m,n);
300:     MatSetType(B,MATSEQSBAIJ);
301:     MatSeqSBAIJSetPreallocation(B,bs,0,browlengths);
302:   } else B = *newmat;

304:   b  = (Mat_SeqSBAIJ*)(B->data);
305:   bi = b->i;
306:   bj = b->j;
307:   bv = b->a;

309:   bi[0] = 0;
310:   for (i=0; i<mbs; i++) {
311:     aj = a->j + a->diag[i];
312:     av = a->a + (a->diag[i])*bs2;
313:     for (j=0; j<browlengths[i]; j++) {
314:       *bj = *aj; bj++; aj++;
315:       for (k=0; k<bs2; k++) {
316:         *bv = *av; bv++; av++;
317:       }
318:     }
319:     bi[i+1]    = bi[i] + browlengths[i];
320:     b->ilen[i] = browlengths[i];
321:   }
322:   PetscFree(browlengths);
323:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
324:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

326:   if (reuse == MAT_INPLACE_MATRIX) {
327:     MatHeaderReplace(A,&B);
328:   } else *newmat = B;
329:   return 0;
330: }