Actual source code: blockmat.c

petsc-3.6.3 2015-12-03
Report Typos and Errors
  2: /*
  3:    This provides a matrix that consists of Mats
  4: */

  6: #include <petsc/private/matimpl.h>              /*I "petscmat.h" I*/
  7: #include <../src/mat/impls/baij/seq/baij.h>    /* use the common AIJ data-structure */

  9: typedef struct {
 10:   SEQAIJHEADER(Mat);
 11:   SEQBAIJHEADER;
 12:   Mat *diags;

 14:   Vec left,right,middle,workb;                 /* dummy vectors to perform local parts of product */
 15: } Mat_BlockMat;

 17: extern PetscErrorCode  MatBlockMatSetPreallocation(Mat,PetscInt,PetscInt,const PetscInt*);

 21: PetscErrorCode MatSOR_BlockMat_Symmetric(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
 22: {
 23:   Mat_BlockMat      *a = (Mat_BlockMat*)A->data;
 24:   PetscScalar       *x;
 25:   const Mat         *v;
 26:   const PetscScalar *b;
 27:   PetscErrorCode    ierr;
 28:   PetscInt          n = A->cmap->n,i,mbs = n/A->rmap->bs,j,bs = A->rmap->bs;
 29:   const PetscInt    *idx;
 30:   IS                row,col;
 31:   MatFactorInfo     info;
 32:   Vec               left = a->left,right = a->right, middle = a->middle;
 33:   Mat               *diag;

 36:   its = its*lits;
 37:   if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
 38:   if (flag & SOR_EISENSTAT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for Eisenstat");
 39:   if (omega != 1.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for omega not equal to 1.0");
 40:   if (fshift) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for fshift");
 41:   if ((flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) && !(flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP)) {
 42:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot do backward sweep without forward sweep");
 43:   }

 45:   if (!a->diags) {
 46:     PetscMalloc1(mbs,&a->diags);
 47:     MatFactorInfoInitialize(&info);
 48:     for (i=0; i<mbs; i++) {
 49:       MatGetOrdering(a->a[a->diag[i]], MATORDERINGND,&row,&col);
 50:       MatCholeskyFactorSymbolic(a->diags[i],a->a[a->diag[i]],row,&info);
 51:       MatCholeskyFactorNumeric(a->diags[i],a->a[a->diag[i]],&info);
 52:       ISDestroy(&row);
 53:       ISDestroy(&col);
 54:     }
 55:     VecDuplicate(bb,&a->workb);
 56:   }
 57:   diag = a->diags;

 59:   VecSet(xx,0.0);
 60:   VecGetArray(xx,&x);
 61:   /* copy right hand side because it must be modified during iteration */
 62:   VecCopy(bb,a->workb);
 63:   VecGetArrayRead(a->workb,&b);

 65:   /* need to add code for when initial guess is zero, see MatSOR_SeqAIJ */
 66:   while (its--) {
 67:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {

 69:       for (i=0; i<mbs; i++) {
 70:         n   = a->i[i+1] - a->i[i] - 1;
 71:         idx = a->j + a->i[i] + 1;
 72:         v   = a->a + a->i[i] + 1;

 74:         VecSet(left,0.0);
 75:         for (j=0; j<n; j++) {
 76:           VecPlaceArray(right,x + idx[j]*bs);
 77:           MatMultAdd(v[j],right,left,left);
 78:           VecResetArray(right);
 79:         }
 80:         VecPlaceArray(right,b + i*bs);
 81:         VecAYPX(left,-1.0,right);
 82:         VecResetArray(right);

 84:         VecPlaceArray(right,x + i*bs);
 85:         MatSolve(diag[i],left,right);

 87:         /* now adjust right hand side, see MatSOR_SeqSBAIJ */
 88:         for (j=0; j<n; j++) {
 89:           MatMultTranspose(v[j],right,left);
 90:           VecPlaceArray(middle,b + idx[j]*bs);
 91:           VecAXPY(middle,-1.0,left);
 92:           VecResetArray(middle);
 93:         }
 94:         VecResetArray(right);

 96:       }
 97:     }
 98:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {

100:       for (i=mbs-1; i>=0; i--) {
101:         n   = a->i[i+1] - a->i[i] - 1;
102:         idx = a->j + a->i[i] + 1;
103:         v   = a->a + a->i[i] + 1;

105:         VecSet(left,0.0);
106:         for (j=0; j<n; j++) {
107:           VecPlaceArray(right,x + idx[j]*bs);
108:           MatMultAdd(v[j],right,left,left);
109:           VecResetArray(right);
110:         }
111:         VecPlaceArray(right,b + i*bs);
112:         VecAYPX(left,-1.0,right);
113:         VecResetArray(right);

115:         VecPlaceArray(right,x + i*bs);
116:         MatSolve(diag[i],left,right);
117:         VecResetArray(right);

119:       }
120:     }
121:   }
122:   VecRestoreArray(xx,&x);
123:   VecRestoreArrayRead(a->workb,&b);
124:   return(0);
125: }

129: PetscErrorCode MatSOR_BlockMat(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx)
130: {
131:   Mat_BlockMat      *a = (Mat_BlockMat*)A->data;
132:   PetscScalar       *x;
133:   const Mat         *v;
134:   const PetscScalar *b;
135:   PetscErrorCode    ierr;
136:   PetscInt          n = A->cmap->n,i,mbs = n/A->rmap->bs,j,bs = A->rmap->bs;
137:   const PetscInt    *idx;
138:   IS                row,col;
139:   MatFactorInfo     info;
140:   Vec               left = a->left,right = a->right;
141:   Mat               *diag;

144:   its = its*lits;
145:   if (its <= 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits);
146:   if (flag & SOR_EISENSTAT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for Eisenstat");
147:   if (omega != 1.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for omega not equal to 1.0");
148:   if (fshift) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for fshift");

150:   if (!a->diags) {
151:     PetscMalloc1(mbs,&a->diags);
152:     MatFactorInfoInitialize(&info);
153:     for (i=0; i<mbs; i++) {
154:       MatGetOrdering(a->a[a->diag[i]], MATORDERINGND,&row,&col);
155:       MatLUFactorSymbolic(a->diags[i],a->a[a->diag[i]],row,col,&info);
156:       MatLUFactorNumeric(a->diags[i],a->a[a->diag[i]],&info);
157:       ISDestroy(&row);
158:       ISDestroy(&col);
159:     }
160:   }
161:   diag = a->diags;

163:   VecSet(xx,0.0);
164:   VecGetArray(xx,&x);
165:   VecGetArrayRead(bb,&b);

167:   /* need to add code for when initial guess is zero, see MatSOR_SeqAIJ */
168:   while (its--) {
169:     if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) {

171:       for (i=0; i<mbs; i++) {
172:         n   = a->i[i+1] - a->i[i];
173:         idx = a->j + a->i[i];
174:         v   = a->a + a->i[i];

176:         VecSet(left,0.0);
177:         for (j=0; j<n; j++) {
178:           if (idx[j] != i) {
179:             VecPlaceArray(right,x + idx[j]*bs);
180:             MatMultAdd(v[j],right,left,left);
181:             VecResetArray(right);
182:           }
183:         }
184:         VecPlaceArray(right,b + i*bs);
185:         VecAYPX(left,-1.0,right);
186:         VecResetArray(right);

188:         VecPlaceArray(right,x + i*bs);
189:         MatSolve(diag[i],left,right);
190:         VecResetArray(right);
191:       }
192:     }
193:     if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) {

195:       for (i=mbs-1; i>=0; i--) {
196:         n   = a->i[i+1] - a->i[i];
197:         idx = a->j + a->i[i];
198:         v   = a->a + a->i[i];

200:         VecSet(left,0.0);
201:         for (j=0; j<n; j++) {
202:           if (idx[j] != i) {
203:             VecPlaceArray(right,x + idx[j]*bs);
204:             MatMultAdd(v[j],right,left,left);
205:             VecResetArray(right);
206:           }
207:         }
208:         VecPlaceArray(right,b + i*bs);
209:         VecAYPX(left,-1.0,right);
210:         VecResetArray(right);

212:         VecPlaceArray(right,x + i*bs);
213:         MatSolve(diag[i],left,right);
214:         VecResetArray(right);

216:       }
217:     }
218:   }
219:   VecRestoreArray(xx,&x);
220:   VecRestoreArrayRead(bb,&b);
221:   return(0);
222: }

226: PetscErrorCode MatSetValues_BlockMat(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is)
227: {
228:   Mat_BlockMat   *a = (Mat_BlockMat*)A->data;
229:   PetscInt       *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1;
230:   PetscInt       *imax=a->imax,*ai=a->i,*ailen=a->ilen;
231:   PetscInt       *aj  =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol;
233:   PetscInt       ridx,cidx;
234:   PetscBool      roworiented=a->roworiented;
235:   MatScalar      value;
236:   Mat            *ap,*aa = a->a;

239:   for (k=0; k<m; k++) { /* loop over added rows */
240:     row  = im[k];
241:     brow = row/bs;
242:     if (row < 0) continue;
243: #if defined(PETSC_USE_DEBUG)
244:     if (row >= A->rmap->N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->N-1);
245: #endif
246:     rp   = aj + ai[brow];
247:     ap   = aa + ai[brow];
248:     rmax = imax[brow];
249:     nrow = ailen[brow];
250:     low  = 0;
251:     high = nrow;
252:     for (l=0; l<n; l++) { /* loop over added columns */
253:       if (in[l] < 0) continue;
254: #if defined(PETSC_USE_DEBUG)
255:       if (in[l] >= A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1);
256: #endif
257:       col = in[l]; bcol = col/bs;
258:       if (A->symmetric && brow > bcol) continue;
259:       ridx = row % bs; cidx = col % bs;
260:       if (roworiented) value = v[l + k*n];
261:       else value = v[k + l*m];

263:       if (col <= lastcol) low = 0;
264:       else high = nrow;
265:       lastcol = col;
266:       while (high-low > 7) {
267:         t = (low+high)/2;
268:         if (rp[t] > bcol) high = t;
269:         else              low  = t;
270:       }
271:       for (i=low; i<high; i++) {
272:         if (rp[i] > bcol) break;
273:         if (rp[i] == bcol) goto noinsert1;
274:       }
275:       if (nonew == 1) goto noinsert1;
276:       if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col);
277:       MatSeqXAIJReallocateAIJ(A,a->mbs,1,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,Mat);
278:       N = nrow++ - 1; high++;
279:       /* shift up all the later entries in this row */
280:       for (ii=N; ii>=i; ii--) {
281:         rp[ii+1] = rp[ii];
282:         ap[ii+1] = ap[ii];
283:       }
284:       if (N>=i) ap[i] = 0;
285:       rp[i] = bcol;
286:       a->nz++;
287:       A->nonzerostate++;
288: noinsert1:;
289:       if (!*(ap+i)) {
290:         MatCreateSeqAIJ(PETSC_COMM_SELF,bs,bs,0,0,ap+i);
291:       }
292:       MatSetValues(ap[i],1,&ridx,1,&cidx,&value,is);
293:       low  = i;
294:     }
295:     ailen[brow] = nrow;
296:   }
297:   return(0);
298: }

302: PetscErrorCode MatLoad_BlockMat(Mat newmat, PetscViewer viewer)
303: {
304:   PetscErrorCode    ierr;
305:   Mat               tmpA;
306:   PetscInt          i,j,m,n,bs = 1,ncols,*lens,currentcol,mbs,**ii,*ilens,nextcol,*llens,cnt = 0;
307:   const PetscInt    *cols;
308:   const PetscScalar *values;
309:   PetscBool         flg = PETSC_FALSE,notdone;
310:   Mat_SeqAIJ        *a;
311:   Mat_BlockMat      *amat;

314:   /* force binary viewer to load .info file if it has not yet done so */
315:   PetscViewerSetUp(viewer);
316:   MatCreate(PETSC_COMM_SELF,&tmpA);
317:   MatSetType(tmpA,MATSEQAIJ);
318:   MatLoad_SeqAIJ(tmpA,viewer);

320:   MatGetLocalSize(tmpA,&m,&n);
321:   PetscOptionsBegin(PETSC_COMM_SELF,NULL,"Options for loading BlockMat matrix 1","Mat");
322:   PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);
323:   PetscOptionsBool("-matload_symmetric","Store the matrix as symmetric","MatLoad",flg,&flg,NULL);
324:   PetscOptionsEnd();

326:   /* Determine number of nonzero blocks for each block row */
327:   a    = (Mat_SeqAIJ*) tmpA->data;
328:   mbs  = m/bs;
329:   PetscMalloc3(mbs,&lens,bs,&ii,bs,&ilens);
330:   PetscMemzero(lens,mbs*sizeof(PetscInt));

332:   for (i=0; i<mbs; i++) {
333:     for (j=0; j<bs; j++) {
334:       ii[j]    = a->j + a->i[i*bs + j];
335:       ilens[j] = a->i[i*bs + j + 1] - a->i[i*bs + j];
336:     }

338:     currentcol = -1;
339:     notdone    = PETSC_TRUE;
340:     while (PETSC_TRUE) {
341:       notdone = PETSC_FALSE;
342:       nextcol = 1000000000;
343:       for (j=0; j<bs; j++) {
344:         while ((ilens[j] > 0 && ii[j][0]/bs <= currentcol)) {
345:           ii[j]++;
346:           ilens[j]--;
347:         }
348:         if (ilens[j] > 0) {
349:           notdone = PETSC_TRUE;
350:           nextcol = PetscMin(nextcol,ii[j][0]/bs);
351:         }
352:       }
353:       if (!notdone) break;
354:       if (!flg || (nextcol >= i)) lens[i]++;
355:       currentcol = nextcol;
356:     }
357:   }

359:   if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) {
360:     MatSetSizes(newmat,m,n,PETSC_DETERMINE,PETSC_DETERMINE);
361:   }
362:   MatBlockMatSetPreallocation(newmat,bs,0,lens);
363:   if (flg) {
364:     MatSetOption(newmat,MAT_SYMMETRIC,PETSC_TRUE);
365:   }
366:   amat = (Mat_BlockMat*)(newmat)->data;

368:   /* preallocate the submatrices */
369:   PetscMalloc1(bs,&llens);
370:   for (i=0; i<mbs; i++) { /* loops for block rows */
371:     for (j=0; j<bs; j++) {
372:       ii[j]    = a->j + a->i[i*bs + j];
373:       ilens[j] = a->i[i*bs + j + 1] - a->i[i*bs + j];
374:     }

376:     currentcol = 1000000000;
377:     for (j=0; j<bs; j++) { /* loop over rows in block finding first nonzero block */
378:       if (ilens[j] > 0) {
379:         currentcol = PetscMin(currentcol,ii[j][0]/bs);
380:       }
381:     }

383:     notdone = PETSC_TRUE;
384:     while (PETSC_TRUE) {  /* loops over blocks in block row */

386:       notdone = PETSC_FALSE;
387:       nextcol = 1000000000;
388:       PetscMemzero(llens,bs*sizeof(PetscInt));
389:       for (j=0; j<bs; j++) { /* loop over rows in block */
390:         while ((ilens[j] > 0 && ii[j][0]/bs <= currentcol)) { /* loop over columns in row */
391:           ii[j]++;
392:           ilens[j]--;
393:           llens[j]++;
394:         }
395:         if (ilens[j] > 0) {
396:           notdone = PETSC_TRUE;
397:           nextcol = PetscMin(nextcol,ii[j][0]/bs);
398:         }
399:       }
400:       if (cnt >= amat->maxnz) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Number of blocks found greater than expected %D",cnt);
401:       if (!flg || currentcol >= i) {
402:         amat->j[cnt] = currentcol;
403:         MatCreateSeqAIJ(PETSC_COMM_SELF,bs,bs,0,llens,amat->a+cnt++);
404:       }

406:       if (!notdone) break;
407:       currentcol = nextcol;
408:     }
409:     amat->ilen[i] = lens[i];
410:   }

412:   PetscFree3(lens,ii,ilens);
413:   PetscFree(llens);

415:   /* copy over the matrix, one row at a time */
416:   for (i=0; i<m; i++) {
417:     MatGetRow(tmpA,i,&ncols,&cols,&values);
418:     MatSetValues(newmat,1,&i,ncols,cols,values,INSERT_VALUES);
419:     MatRestoreRow(tmpA,i,&ncols,&cols,&values);
420:   }
421:   MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);
422:   MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);
423:   return(0);
424: }

428: PetscErrorCode MatView_BlockMat(Mat A,PetscViewer viewer)
429: {
430:   Mat_BlockMat      *a = (Mat_BlockMat*)A->data;
431:   PetscErrorCode    ierr;
432:   const char        *name;
433:   PetscViewerFormat format;

436:   PetscObjectGetName((PetscObject)A,&name);
437:   PetscViewerGetFormat(viewer,&format);
438:   if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) {
439:     PetscViewerASCIIPrintf(viewer,"Nonzero block matrices = %D \n",a->nz);
440:     if (A->symmetric) {
441:       PetscViewerASCIIPrintf(viewer,"Only upper triangular part of symmetric matrix is stored\n");
442:     }
443:   }
444:   return(0);
445: }

449: PetscErrorCode MatDestroy_BlockMat(Mat mat)
450: {
452:   Mat_BlockMat   *bmat = (Mat_BlockMat*)mat->data;
453:   PetscInt       i;

456:   VecDestroy(&bmat->right);
457:   VecDestroy(&bmat->left);
458:   VecDestroy(&bmat->middle);
459:   VecDestroy(&bmat->workb);
460:   if (bmat->diags) {
461:     for (i=0; i<mat->rmap->n/mat->rmap->bs; i++) {
462:       MatDestroy(&bmat->diags[i]);
463:     }
464:   }
465:   if (bmat->a) {
466:     for (i=0; i<bmat->nz; i++) {
467:       MatDestroy(&bmat->a[i]);
468:     }
469:   }
470:   MatSeqXAIJFreeAIJ(mat,(PetscScalar**)&bmat->a,&bmat->j,&bmat->i);
471:   PetscFree(mat->data);
472:   return(0);
473: }

477: PetscErrorCode MatMult_BlockMat(Mat A,Vec x,Vec y)
478: {
479:   Mat_BlockMat   *bmat = (Mat_BlockMat*)A->data;
481:   PetscScalar    *xx,*yy;
482:   PetscInt       *aj,i,*ii,jrow,m = A->rmap->n/A->rmap->bs,bs = A->rmap->bs,n,j;
483:   Mat            *aa;

486:   /*
487:      Standard CSR multiply except each entry is a Mat
488:   */
489:   VecGetArray(x,&xx);

491:   VecSet(y,0.0);
492:   VecGetArray(y,&yy);
493:   aj   = bmat->j;
494:   aa   = bmat->a;
495:   ii   = bmat->i;
496:   for (i=0; i<m; i++) {
497:     jrow = ii[i];
498:     VecPlaceArray(bmat->left,yy + bs*i);
499:     n    = ii[i+1] - jrow;
500:     for (j=0; j<n; j++) {
501:       VecPlaceArray(bmat->right,xx + bs*aj[jrow]);
502:       MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);
503:       VecResetArray(bmat->right);
504:       jrow++;
505:     }
506:     VecResetArray(bmat->left);
507:   }
508:   VecRestoreArray(x,&xx);
509:   VecRestoreArray(y,&yy);
510:   return(0);
511: }

515: PetscErrorCode MatMult_BlockMat_Symmetric(Mat A,Vec x,Vec y)
516: {
517:   Mat_BlockMat   *bmat = (Mat_BlockMat*)A->data;
519:   PetscScalar    *xx,*yy;
520:   PetscInt       *aj,i,*ii,jrow,m = A->rmap->n/A->rmap->bs,bs = A->rmap->bs,n,j;
521:   Mat            *aa;

524:   /*
525:      Standard CSR multiply except each entry is a Mat
526:   */
527:   VecGetArray(x,&xx);

529:   VecSet(y,0.0);
530:   VecGetArray(y,&yy);
531:   aj   = bmat->j;
532:   aa   = bmat->a;
533:   ii   = bmat->i;
534:   for (i=0; i<m; i++) {
535:     jrow = ii[i];
536:     n    = ii[i+1] - jrow;
537:     VecPlaceArray(bmat->left,yy + bs*i);
538:     VecPlaceArray(bmat->middle,xx + bs*i);
539:     /* if we ALWAYS required a diagonal entry then could remove this if test */
540:     if (aj[jrow] == i) {
541:       VecPlaceArray(bmat->right,xx + bs*aj[jrow]);
542:       MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);
543:       VecResetArray(bmat->right);
544:       jrow++;
545:       n--;
546:     }
547:     for (j=0; j<n; j++) {
548:       VecPlaceArray(bmat->right,xx + bs*aj[jrow]);            /* upper triangular part */
549:       MatMultAdd(aa[jrow],bmat->right,bmat->left,bmat->left);
550:       VecResetArray(bmat->right);

552:       VecPlaceArray(bmat->right,yy + bs*aj[jrow]);            /* lower triangular part */
553:       MatMultTransposeAdd(aa[jrow],bmat->middle,bmat->right,bmat->right);
554:       VecResetArray(bmat->right);
555:       jrow++;
556:     }
557:     VecResetArray(bmat->left);
558:     VecResetArray(bmat->middle);
559:   }
560:   VecRestoreArray(x,&xx);
561:   VecRestoreArray(y,&yy);
562:   return(0);
563: }

567: PetscErrorCode MatMultAdd_BlockMat(Mat A,Vec x,Vec y,Vec z)
568: {
570:   return(0);
571: }

575: PetscErrorCode MatMultTranspose_BlockMat(Mat A,Vec x,Vec y)
576: {
578:   return(0);
579: }

583: PetscErrorCode MatMultTransposeAdd_BlockMat(Mat A,Vec x,Vec y,Vec z)
584: {
586:   return(0);
587: }

589: /*
590:      Adds diagonal pointers to sparse matrix structure.
591: */
594: PetscErrorCode MatMarkDiagonal_BlockMat(Mat A)
595: {
596:   Mat_BlockMat   *a = (Mat_BlockMat*)A->data;
598:   PetscInt       i,j,mbs = A->rmap->n/A->rmap->bs;

601:   if (!a->diag) {
602:     PetscMalloc1(mbs,&a->diag);
603:   }
604:   for (i=0; i<mbs; i++) {
605:     a->diag[i] = a->i[i+1];
606:     for (j=a->i[i]; j<a->i[i+1]; j++) {
607:       if (a->j[j] == i) {
608:         a->diag[i] = j;
609:         break;
610:       }
611:     }
612:   }
613:   return(0);
614: }

618: PetscErrorCode MatGetSubMatrix_BlockMat(Mat A,IS isrow,IS iscol,MatReuse scall,Mat *B)
619: {
620:   Mat_BlockMat   *a = (Mat_BlockMat*)A->data;
621:   Mat_SeqAIJ     *c;
623:   PetscInt       i,k,first,step,lensi,nrows,ncols;
624:   PetscInt       *j_new,*i_new,*aj = a->j,*ailen = a->ilen;
625:   PetscScalar    *a_new;
626:   Mat            C,*aa = a->a;
627:   PetscBool      stride,equal;

630:   ISEqual(isrow,iscol,&equal);
631:   if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only for idential column and row indices");
632:   PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);
633:   if (!stride) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only for stride indices");
634:   ISStrideGetInfo(iscol,&first,&step);
635:   if (step != A->rmap->bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only select one entry from each block");

637:   ISGetLocalSize(isrow,&nrows);
638:   ncols = nrows;

640:   /* create submatrix */
641:   if (scall == MAT_REUSE_MATRIX) {
642:     PetscInt n_cols,n_rows;
643:     C    = *B;
644:     MatGetSize(C,&n_rows,&n_cols);
645:     if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size");
646:     MatZeroEntries(C);
647:   } else {
648:     MatCreate(PetscObjectComm((PetscObject)A),&C);
649:     MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);
650:     if (A->symmetric) {
651:       MatSetType(C,MATSEQSBAIJ);
652:     } else {
653:       MatSetType(C,MATSEQAIJ);
654:     }
655:     MatSeqAIJSetPreallocation(C,0,ailen);
656:     MatSeqSBAIJSetPreallocation(C,1,0,ailen);
657:   }
658:   c = (Mat_SeqAIJ*)C->data;

660:   /* loop over rows inserting into submatrix */
661:   a_new = c->a;
662:   j_new = c->j;
663:   i_new = c->i;

665:   for (i=0; i<nrows; i++) {
666:     lensi = ailen[i];
667:     for (k=0; k<lensi; k++) {
668:       *j_new++ = *aj++;
669:       MatGetValue(*aa++,first,first,a_new++);
670:     }
671:     i_new[i+1] = i_new[i] + lensi;
672:     c->ilen[i] = lensi;
673:   }

675:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
676:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
677:   *B   = C;
678:   return(0);
679: }

683: PetscErrorCode MatAssemblyEnd_BlockMat(Mat A,MatAssemblyType mode)
684: {
685:   Mat_BlockMat   *a = (Mat_BlockMat*)A->data;
687:   PetscInt       fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax;
688:   PetscInt       m      = a->mbs,*ip,N,*ailen = a->ilen,rmax = 0;
689:   Mat            *aa    = a->a,*ap;

692:   if (mode == MAT_FLUSH_ASSEMBLY) return(0);

694:   if (m) rmax = ailen[0]; /* determine row with most nonzeros */
695:   for (i=1; i<m; i++) {
696:     /* move each row back by the amount of empty slots (fshift) before it*/
697:     fshift += imax[i-1] - ailen[i-1];
698:     rmax    = PetscMax(rmax,ailen[i]);
699:     if (fshift) {
700:       ip = aj + ai[i];
701:       ap = aa + ai[i];
702:       N  = ailen[i];
703:       for (j=0; j<N; j++) {
704:         ip[j-fshift] = ip[j];
705:         ap[j-fshift] = ap[j];
706:       }
707:     }
708:     ai[i] = ai[i-1] + ailen[i-1];
709:   }
710:   if (m) {
711:     fshift += imax[m-1] - ailen[m-1];
712:     ai[m]   = ai[m-1] + ailen[m-1];
713:   }
714:   /* reset ilen and imax for each row */
715:   for (i=0; i<m; i++) {
716:     ailen[i] = imax[i] = ai[i+1] - ai[i];
717:   }
718:   a->nz = ai[m];
719:   for (i=0; i<a->nz; i++) {
720: #if defined(PETSC_USE_DEBUG)
721:     if (!aa[i]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Null matrix at location %D column %D nz %D",i,aj[i],a->nz);
722: #endif
723:     MatAssemblyBegin(aa[i],MAT_FINAL_ASSEMBLY);
724:     MatAssemblyEnd(aa[i],MAT_FINAL_ASSEMBLY);
725:   }
726:   PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap->n/A->cmap->bs,fshift,a->nz);
727:   PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);
728:   PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);

730:   A->info.mallocs    += a->reallocs;
731:   a->reallocs         = 0;
732:   A->info.nz_unneeded = (double)fshift;
733:   a->rmax             = rmax;
734:   MatMarkDiagonal_BlockMat(A);
735:   return(0);
736: }

740: PetscErrorCode MatSetOption_BlockMat(Mat A,MatOption opt,PetscBool flg)
741: {
743:   if (opt == MAT_SYMMETRIC && flg) {
744:     A->ops->sor  = MatSOR_BlockMat_Symmetric;
745:     A->ops->mult = MatMult_BlockMat_Symmetric;
746:   } else {
747:     PetscInfo1(A,"Unused matrix option %s\n",MatOptions[opt]);
748:   }
749:   return(0);
750: }


753: static struct _MatOps MatOps_Values = {MatSetValues_BlockMat,
754:                                        0,
755:                                        0,
756:                                        MatMult_BlockMat,
757:                                /*  4*/ MatMultAdd_BlockMat,
758:                                        MatMultTranspose_BlockMat,
759:                                        MatMultTransposeAdd_BlockMat,
760:                                        0,
761:                                        0,
762:                                        0,
763:                                /* 10*/ 0,
764:                                        0,
765:                                        0,
766:                                        MatSOR_BlockMat,
767:                                        0,
768:                                /* 15*/ 0,
769:                                        0,
770:                                        0,
771:                                        0,
772:                                        0,
773:                                /* 20*/ 0,
774:                                        MatAssemblyEnd_BlockMat,
775:                                        MatSetOption_BlockMat,
776:                                        0,
777:                                /* 24*/ 0,
778:                                        0,
779:                                        0,
780:                                        0,
781:                                        0,
782:                                /* 29*/ 0,
783:                                        0,
784:                                        0,
785:                                        0,
786:                                        0,
787:                                /* 34*/ 0,
788:                                        0,
789:                                        0,
790:                                        0,
791:                                        0,
792:                                /* 39*/ 0,
793:                                        0,
794:                                        0,
795:                                        0,
796:                                        0,
797:                                /* 44*/ 0,
798:                                        0,
799:                                        MatShift_Basic,
800:                                        0,
801:                                        0,
802:                                /* 49*/ 0,
803:                                        0,
804:                                        0,
805:                                        0,
806:                                        0,
807:                                /* 54*/ 0,
808:                                        0,
809:                                        0,
810:                                        0,
811:                                        0,
812:                                /* 59*/ MatGetSubMatrix_BlockMat,
813:                                        MatDestroy_BlockMat,
814:                                        MatView_BlockMat,
815:                                        0,
816:                                        0,
817:                                /* 64*/ 0,
818:                                        0,
819:                                        0,
820:                                        0,
821:                                        0,
822:                                /* 69*/ 0,
823:                                        0,
824:                                        0,
825:                                        0,
826:                                        0,
827:                                /* 74*/ 0,
828:                                        0,
829:                                        0,
830:                                        0,
831:                                        0,
832:                                /* 79*/ 0,
833:                                        0,
834:                                        0,
835:                                        0,
836:                                        MatLoad_BlockMat,
837:                                /* 84*/ 0,
838:                                        0,
839:                                        0,
840:                                        0,
841:                                        0,
842:                                /* 89*/ 0,
843:                                        0,
844:                                        0,
845:                                        0,
846:                                        0,
847:                                /* 94*/ 0,
848:                                        0,
849:                                        0,
850:                                        0,
851:                                        0,
852:                                /* 99*/ 0,
853:                                        0,
854:                                        0,
855:                                        0,
856:                                        0,
857:                                /*104*/ 0,
858:                                        0,
859:                                        0,
860:                                        0,
861:                                        0,
862:                                /*109*/ 0,
863:                                        0,
864:                                        0,
865:                                        0,
866:                                        0,
867:                                /*114*/ 0,
868:                                        0,
869:                                        0,
870:                                        0,
871:                                        0,
872:                                /*119*/ 0,
873:                                        0,
874:                                        0,
875:                                        0,
876:                                        0,
877:                                /*124*/ 0,
878:                                        0,
879:                                        0,
880:                                        0,
881:                                        0,
882:                                /*129*/ 0,
883:                                        0,
884:                                        0,
885:                                        0,
886:                                        0,
887:                                /*134*/ 0,
888:                                        0,
889:                                        0,
890:                                        0,
891:                                        0,
892:                                /*139*/ 0,
893:                                        0,
894:                                        0
895: };

899: /*@C
900:    MatBlockMatSetPreallocation - For good matrix assembly performance
901:    the user should preallocate the matrix storage by setting the parameter nz
902:    (or the array nnz).  By setting these parameters accurately, performance
903:    during matrix assembly can be increased by more than a factor of 50.

905:    Collective on MPI_Comm

907:    Input Parameters:
908: +  B - The matrix
909: .  bs - size of each block in matrix
910: .  nz - number of nonzeros per block row (same for all rows)
911: -  nnz - array containing the number of nonzeros in the various block rows
912:          (possibly different for each row) or NULL

914:    Notes:
915:      If nnz is given then nz is ignored

917:    Specify the preallocated storage with either nz or nnz (not both).
918:    Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory
919:    allocation.  For large problems you MUST preallocate memory or you
920:    will get TERRIBLE performance, see the users' manual chapter on matrices.

922:    Level: intermediate

924: .seealso: MatCreate(), MatCreateBlockMat(), MatSetValues()

926: @*/
927: PetscErrorCode  MatBlockMatSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[])
928: {

932:   PetscTryMethod(B,"MatBlockMatSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));
933:   return(0);
934: }

938: PetscErrorCode  MatBlockMatSetPreallocation_BlockMat(Mat A,PetscInt bs,PetscInt nz,PetscInt *nnz)
939: {
940:   Mat_BlockMat   *bmat = (Mat_BlockMat*)A->data;
942:   PetscInt       i;

945:   PetscLayoutSetBlockSize(A->rmap,bs);
946:   PetscLayoutSetBlockSize(A->cmap,bs);
947:   PetscLayoutSetUp(A->rmap);
948:   PetscLayoutSetUp(A->cmap);
949:   PetscLayoutGetBlockSize(A->rmap,&bs);

951:   if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5;
952:   if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz);
953:   if (nnz) {
954:     for (i=0; i<A->rmap->n/bs; i++) {
955:       if (nnz[i] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]);
956:       if (nnz[i] > A->cmap->n/bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],A->cmap->n/bs);
957:     }
958:   }
959:   bmat->mbs = A->rmap->n/bs;

961:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,bs,NULL,&bmat->right);
962:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,bs,NULL,&bmat->middle);
963:   VecCreateSeq(PETSC_COMM_SELF,bs,&bmat->left);

965:   if (!bmat->imax) {
966:     PetscMalloc2(A->rmap->n,&bmat->imax,A->rmap->n,&bmat->ilen);
967:     PetscLogObjectMemory((PetscObject)A,2*A->rmap->n*sizeof(PetscInt));
968:   }
969:   if (nnz) {
970:     nz = 0;
971:     for (i=0; i<A->rmap->n/A->rmap->bs; i++) {
972:       bmat->imax[i] = nnz[i];
973:       nz           += nnz[i];
974:     }
975:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently requires block row by row preallocation");

977:   /* bmat->ilen will count nonzeros in each row so far. */
978:   for (i=0; i<bmat->mbs; i++) bmat->ilen[i] = 0;

980:   /* allocate the matrix space */
981:   MatSeqXAIJFreeAIJ(A,(PetscScalar**)&bmat->a,&bmat->j,&bmat->i);
982:   PetscMalloc3(nz,&bmat->a,nz,&bmat->j,A->rmap->n+1,&bmat->i);
983:   PetscLogObjectMemory((PetscObject)A,(A->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));
984:   bmat->i[0] = 0;
985:   for (i=1; i<bmat->mbs+1; i++) {
986:     bmat->i[i] = bmat->i[i-1] + bmat->imax[i-1];
987:   }
988:   bmat->singlemalloc = PETSC_TRUE;
989:   bmat->free_a       = PETSC_TRUE;
990:   bmat->free_ij      = PETSC_TRUE;

992:   bmat->nz            = 0;
993:   bmat->maxnz         = nz;
994:   A->info.nz_unneeded = (double)bmat->maxnz;
995:   MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
996:   return(0);
997: }

999: /*MC
1000:    MATBLOCKMAT - A matrix that is defined by a set of Mat's that represents a sparse block matrix
1001:                  consisting of (usually) sparse blocks.

1003:   Level: advanced

1005: .seealso: MatCreateBlockMat()

1007: M*/

1011: PETSC_EXTERN PetscErrorCode MatCreate_BlockMat(Mat A)
1012: {
1013:   Mat_BlockMat   *b;

1017:   PetscNewLog(A,&b);
1018:   A->data = (void*)b;
1019:   PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));

1021:   A->assembled    = PETSC_TRUE;
1022:   A->preallocated = PETSC_FALSE;
1023:   PetscObjectChangeTypeName((PetscObject)A,MATBLOCKMAT);

1025:   PetscObjectComposeFunction((PetscObject)A,"MatBlockMatSetPreallocation_C",MatBlockMatSetPreallocation_BlockMat);
1026:   return(0);
1027: }

1031: /*@C
1032:    MatCreateBlockMat - Creates a new matrix in which each block contains a uniform-size sequential Mat object

1034:   Collective on MPI_Comm

1036:    Input Parameters:
1037: +  comm - MPI communicator
1038: .  m - number of rows
1039: .  n  - number of columns
1040: .  bs - size of each submatrix
1041: .  nz  - expected maximum number of nonzero blocks in row (use PETSC_DEFAULT if not known)
1042: -  nnz - expected number of nonzers per block row if known (use NULL otherwise)


1045:    Output Parameter:
1046: .  A - the matrix

1048:    Level: intermediate

1050:    Notes: Matrices of this type are nominally-sparse matrices in which each "entry" is a Mat object.  Each Mat must
1051:    have the same size and be sequential.  The local and global sizes must be compatible with this decomposition.

1053:    For matrices containing parallel submatrices and variable block sizes, see MATNEST.

1055: .keywords: matrix, bmat, create

1057: .seealso: MATBLOCKMAT, MatCreateNest()
1058: @*/
1059: PetscErrorCode  MatCreateBlockMat(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt bs,PetscInt nz,PetscInt *nnz, Mat *A)
1060: {

1064:   MatCreate(comm,A);
1065:   MatSetSizes(*A,m,n,PETSC_DETERMINE,PETSC_DETERMINE);
1066:   MatSetType(*A,MATBLOCKMAT);
1067:   MatBlockMatSetPreallocation(*A,bs,nz,nnz);
1068:   return(0);
1069: }