Actual source code: matnest.c

petsc-master 2016-08-28
Report Typos and Errors
 2:  #include <../src/mat/impls/nest/matnestimpl.h>
 3:  #include <petscsf.h>

  5: static PetscErrorCode MatSetUp_NestIS_Private(Mat,PetscInt,const IS[],PetscInt,const IS[]);
  6: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left);

  8: /* private functions */
 11: static PetscErrorCode MatNestGetSizes_Private(Mat A,PetscInt *m,PetscInt *n,PetscInt *M,PetscInt *N)
 12: {
 13:   Mat_Nest       *bA = (Mat_Nest*)A->data;
 14:   PetscInt       i,j;

 18:   *m = *n = *M = *N = 0;
 19:   for (i=0; i<bA->nr; i++) {  /* rows */
 20:     PetscInt sm,sM;
 21:     ISGetLocalSize(bA->isglobal.row[i],&sm);
 22:     ISGetSize(bA->isglobal.row[i],&sM);
 23:     *m  += sm;
 24:     *M  += sM;
 25:   }
 26:   for (j=0; j<bA->nc; j++) {  /* cols */
 27:     PetscInt sn,sN;
 28:     ISGetLocalSize(bA->isglobal.col[j],&sn);
 29:     ISGetSize(bA->isglobal.col[j],&sN);
 30:     *n  += sn;
 31:     *N  += sN;
 32:   }
 33:   return(0);
 34: }

 36: /* operations */
 39: static PetscErrorCode MatMult_Nest(Mat A,Vec x,Vec y)
 40: {
 41:   Mat_Nest       *bA = (Mat_Nest*)A->data;
 42:   Vec            *bx = bA->right,*by = bA->left;
 43:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

 47:   for (i=0; i<nr; i++) {VecGetSubVector(y,bA->isglobal.row[i],&by[i]);}
 48:   for (i=0; i<nc; i++) {VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);}
 49:   for (i=0; i<nr; i++) {
 50:     VecZeroEntries(by[i]);
 51:     for (j=0; j<nc; j++) {
 52:       if (!bA->m[i][j]) continue;
 53:       /* y[i] <- y[i] + A[i][j] * x[j] */
 54:       MatMultAdd(bA->m[i][j],bx[j],by[i],by[i]);
 55:     }
 56:   }
 57:   for (i=0; i<nr; i++) {VecRestoreSubVector(y,bA->isglobal.row[i],&by[i]);}
 58:   for (i=0; i<nc; i++) {VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);}
 59:   return(0);
 60: }

 64: static PetscErrorCode MatMultAdd_Nest(Mat A,Vec x,Vec y,Vec z)
 65: {
 66:   Mat_Nest       *bA = (Mat_Nest*)A->data;
 67:   Vec            *bx = bA->right,*bz = bA->left;
 68:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

 72:   for (i=0; i<nr; i++) {VecGetSubVector(z,bA->isglobal.row[i],&bz[i]);}
 73:   for (i=0; i<nc; i++) {VecGetSubVector(x,bA->isglobal.col[i],&bx[i]);}
 74:   for (i=0; i<nr; i++) {
 75:     if (y != z) {
 76:       Vec by;
 77:       VecGetSubVector(y,bA->isglobal.row[i],&by);
 78:       VecCopy(by,bz[i]);
 79:       VecRestoreSubVector(y,bA->isglobal.row[i],&by);
 80:     }
 81:     for (j=0; j<nc; j++) {
 82:       if (!bA->m[i][j]) continue;
 83:       /* y[i] <- y[i] + A[i][j] * x[j] */
 84:       MatMultAdd(bA->m[i][j],bx[j],bz[i],bz[i]);
 85:     }
 86:   }
 87:   for (i=0; i<nr; i++) {VecRestoreSubVector(z,bA->isglobal.row[i],&bz[i]);}
 88:   for (i=0; i<nc; i++) {VecRestoreSubVector(x,bA->isglobal.col[i],&bx[i]);}
 89:   return(0);
 90: }

 94: static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y)
 95: {
 96:   Mat_Nest       *bA = (Mat_Nest*)A->data;
 97:   Vec            *bx = bA->left,*by = bA->right;
 98:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

102:   for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
103:   for (i=0; i<nc; i++) {VecGetSubVector(y,bA->isglobal.col[i],&by[i]);}
104:   for (j=0; j<nc; j++) {
105:     VecZeroEntries(by[j]);
106:     for (i=0; i<nr; i++) {
107:       if (!bA->m[i][j]) continue;
108:       /* y[j] <- y[j] + (A[i][j])^T * x[i] */
109:       MatMultTransposeAdd(bA->m[i][j],bx[i],by[j],by[j]);
110:     }
111:   }
112:   for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
113:   for (i=0; i<nc; i++) {VecRestoreSubVector(y,bA->isglobal.col[i],&by[i]);}
114:   return(0);
115: }

119: static PetscErrorCode MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z)
120: {
121:   Mat_Nest       *bA = (Mat_Nest*)A->data;
122:   Vec            *bx = bA->left,*bz = bA->right;
123:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

127:   for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
128:   for (i=0; i<nc; i++) {VecGetSubVector(z,bA->isglobal.col[i],&bz[i]);}
129:   for (j=0; j<nc; j++) {
130:     if (y != z) {
131:       Vec by;
132:       VecGetSubVector(y,bA->isglobal.col[j],&by);
133:       VecCopy(by,bz[j]);
134:       VecRestoreSubVector(y,bA->isglobal.col[j],&by);
135:     }
136:     for (i=0; i<nr; i++) {
137:       if (!bA->m[i][j]) continue;
138:       /* z[j] <- y[j] + (A[i][j])^T * x[i] */
139:       MatMultTransposeAdd(bA->m[i][j],bx[i],bz[j],bz[j]);
140:     }
141:   }
142:   for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
143:   for (i=0; i<nc; i++) {VecRestoreSubVector(z,bA->isglobal.col[i],&bz[i]);}
144:   return(0);
145: }

149: static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list)
150: {
152:   IS             *lst = *list;
153:   PetscInt       i;

156:   if (!lst) return(0);
157:   for (i=0; i<n; i++) if (lst[i]) {ISDestroy(&lst[i]);}
158:   PetscFree(lst);
159:   *list = NULL;
160:   return(0);
161: }

165: static PetscErrorCode MatDestroy_Nest(Mat A)
166: {
167:   Mat_Nest       *vs = (Mat_Nest*)A->data;
168:   PetscInt       i,j;

172:   /* release the matrices and the place holders */
173:   MatNestDestroyISList(vs->nr,&vs->isglobal.row);
174:   MatNestDestroyISList(vs->nc,&vs->isglobal.col);
175:   MatNestDestroyISList(vs->nr,&vs->islocal.row);
176:   MatNestDestroyISList(vs->nc,&vs->islocal.col);

178:   PetscFree(vs->row_len);
179:   PetscFree(vs->col_len);

181:   PetscFree2(vs->left,vs->right);

183:   /* release the matrices and the place holders */
184:   if (vs->m) {
185:     for (i=0; i<vs->nr; i++) {
186:       for (j=0; j<vs->nc; j++) {
187:         MatDestroy(&vs->m[i][j]);
188:       }
189:       PetscFree(vs->m[i]);
190:     }
191:     PetscFree(vs->m);
192:   }
193:   PetscFree(A->data);

195:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",0);
196:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",0);
197:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",0);
198:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",0);
199:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",0);
200:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",0);
201:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",0);
202:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",0);
203:   return(0);
204: }

208: static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)
209: {
210:   Mat_Nest       *vs = (Mat_Nest*)A->data;
211:   PetscInt       i,j;

215:   for (i=0; i<vs->nr; i++) {
216:     for (j=0; j<vs->nc; j++) {
217:       if (vs->m[i][j]) {
218:         MatAssemblyBegin(vs->m[i][j],type);
219:         if (!vs->splitassembly) {
220:           /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
221:            * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
222:            * already performing an assembly, but the result would by more complicated and appears to offer less
223:            * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
224:            * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
225:            */
226:           MatAssemblyEnd(vs->m[i][j],type);
227:         }
228:       }
229:     }
230:   }
231:   return(0);
232: }

236: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
237: {
238:   Mat_Nest       *vs = (Mat_Nest*)A->data;
239:   PetscInt       i,j;

243:   for (i=0; i<vs->nr; i++) {
244:     for (j=0; j<vs->nc; j++) {
245:       if (vs->m[i][j]) {
246:         if (vs->splitassembly) {
247:           MatAssemblyEnd(vs->m[i][j],type);
248:         }
249:       }
250:     }
251:   }
252:   return(0);
253: }

257: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B)
258: {
260:   Mat_Nest       *vs = (Mat_Nest*)A->data;
261:   PetscInt       j;
262:   Mat            sub;

265:   sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
266:   for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j];
267:   if (sub) {MatSetUp(sub);}       /* Ensure that the sizes are available */
268:   *B = sub;
269:   return(0);
270: }

274: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B)
275: {
277:   Mat_Nest       *vs = (Mat_Nest*)A->data;
278:   PetscInt       i;
279:   Mat            sub;

282:   sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
283:   for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col];
284:   if (sub) {MatSetUp(sub);}       /* Ensure that the sizes are available */
285:   *B = sub;
286:   return(0);
287: }

291: static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found)
292: {
294:   PetscInt       i;
295:   PetscBool      flg;

301:   *found = -1;
302:   for (i=0; i<n; i++) {
303:     if (!list[i]) continue;
304:     ISEqual(list[i],is,&flg);
305:     if (flg) {
306:       *found = i;
307:       return(0);
308:     }
309:   }
310:   SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set");
311:   return(0);
312: }

316: /* Get a block row as a new MatNest */
317: static PetscErrorCode MatNestGetRow(Mat A,PetscInt row,Mat *B)
318: {
319:   Mat_Nest       *vs = (Mat_Nest*)A->data;
320:   char           keyname[256];

324:   *B   = NULL;
325:   PetscSNPrintf(keyname,sizeof(keyname),"NestRow_%D",row);
326:   PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);
327:   if (*B) return(0);

329:   MatCreateNest(PetscObjectComm((PetscObject)A),1,NULL,vs->nc,vs->isglobal.col,vs->m[row],B);

331:   (*B)->assembled = A->assembled;

333:   PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);
334:   PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
335:   return(0);
336: }

340: static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
341: {
342:   Mat_Nest       *vs = (Mat_Nest*)A->data;
344:   PetscInt       row,col;
345:   PetscBool      same,isFullCol,isFullColGlobal;

348:   /* Check if full column space. This is a hack */
349:   isFullCol = PETSC_FALSE;
350:   PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&same);
351:   if (same) {
352:     PetscInt n,first,step,i,an,am,afirst,astep;
353:     ISStrideGetInfo(iscol,&first,&step);
354:     ISGetLocalSize(iscol,&n);
355:     isFullCol = PETSC_TRUE;
356:     for (i=0,an=A->cmap->rstart; i<vs->nc; i++) {
357:       ISStrideGetInfo(is->col[i],&afirst,&astep);
358:       ISGetLocalSize(is->col[i],&am);
359:       if (afirst != an || astep != step) isFullCol = PETSC_FALSE;
360:       an += am;
361:     }
362:     if (an != A->cmap->rstart+n) isFullCol = PETSC_FALSE;
363:   }
364:   MPIU_Allreduce(&isFullCol,&isFullColGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)iscol));

366:   if (isFullColGlobal) {
367:     PetscInt row;
368:     MatNestFindIS(A,vs->nr,is->row,isrow,&row);
369:     MatNestGetRow(A,row,B);
370:   } else {
371:     MatNestFindIS(A,vs->nr,is->row,isrow,&row);
372:     MatNestFindIS(A,vs->nc,is->col,iscol,&col);
373:     if (!vs->m[row][col]) {
374:       PetscInt lr,lc;

376:       MatCreate(PetscObjectComm((PetscObject)A),&vs->m[row][col]);
377:       ISGetLocalSize(vs->isglobal.row[row],&lr);
378:       ISGetLocalSize(vs->isglobal.col[col],&lc);
379:       MatSetSizes(vs->m[row][col],lr,lc,PETSC_DECIDE,PETSC_DECIDE);
380:       MatSetUp(vs->m[row][col]);
381:       MatAssemblyBegin(vs->m[row][col],MAT_FINAL_ASSEMBLY);
382:       MatAssemblyEnd(vs->m[row][col],MAT_FINAL_ASSEMBLY);
383:     }
384:     *B = vs->m[row][col];
385:   }
386:   return(0);
387: }

391: static PetscErrorCode MatGetSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
392: {
394:   Mat_Nest       *vs = (Mat_Nest*)A->data;
395:   Mat            sub;

398:   MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
399:   switch (reuse) {
400:   case MAT_INITIAL_MATRIX:
401:     if (sub) { PetscObjectReference((PetscObject)sub); }
402:     *B = sub;
403:     break;
404:   case MAT_REUSE_MATRIX:
405:     if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call");
406:     break;
407:   case MAT_IGNORE_MATRIX:       /* Nothing to do */
408:     break;
409:   case MAT_INPLACE_MATRIX:       /* Nothing to do */
410:     SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet");
411:     break;
412:   }
413:   return(0);
414: }

418: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
419: {
421:   Mat_Nest       *vs = (Mat_Nest*)A->data;
422:   Mat            sub;

425:   MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
426:   /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
427:   if (sub) {PetscObjectReference((PetscObject)sub);}
428:   *B = sub;
429:   return(0);
430: }

434: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
435: {
437:   Mat_Nest       *vs = (Mat_Nest*)A->data;
438:   Mat            sub;

441:   MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
442:   if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten");
443:   if (sub) {
444:     if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times");
445:     MatDestroy(B);
446:   }
447:   return(0);
448: }

452: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
453: {
454:   Mat_Nest       *bA = (Mat_Nest*)A->data;
455:   PetscInt       i;

459:   for (i=0; i<bA->nr; i++) {
460:     Vec bv;
461:     VecGetSubVector(v,bA->isglobal.row[i],&bv);
462:     if (bA->m[i][i]) {
463:       MatGetDiagonal(bA->m[i][i],bv);
464:     } else {
465:       VecSet(bv,0.0);
466:     }
467:     VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
468:   }
469:   return(0);
470: }

474: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
475: {
476:   Mat_Nest       *bA = (Mat_Nest*)A->data;
477:   Vec            bl,*br;
478:   PetscInt       i,j;

482:   PetscCalloc1(bA->nc,&br);
483:   if (r) {
484:     for (j=0; j<bA->nc; j++) {VecGetSubVector(r,bA->isglobal.col[j],&br[j]);}
485:   }
486:   bl = NULL;
487:   for (i=0; i<bA->nr; i++) {
488:     if (l) {
489:       VecGetSubVector(l,bA->isglobal.row[i],&bl);
490:     }
491:     for (j=0; j<bA->nc; j++) {
492:       if (bA->m[i][j]) {
493:         MatDiagonalScale(bA->m[i][j],bl,br[j]);
494:       }
495:     }
496:     if (l) {
497:       VecRestoreSubVector(l,bA->isglobal.row[i],&bl);
498:     }
499:   }
500:   if (r) {
501:     for (j=0; j<bA->nc; j++) {VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);}
502:   }
503:   PetscFree(br);
504:   return(0);
505: }

509: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
510: {
511:   Mat_Nest       *bA = (Mat_Nest*)A->data;
512:   PetscInt       i,j;

516:   for (i=0; i<bA->nr; i++) {
517:     for (j=0; j<bA->nc; j++) {
518:       if (bA->m[i][j]) {
519:         MatScale(bA->m[i][j],a);
520:       }
521:     }
522:   }
523:   return(0);
524: }

528: static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
529: {
530:   Mat_Nest       *bA = (Mat_Nest*)A->data;
531:   PetscInt       i;

535:   for (i=0; i<bA->nr; i++) {
536:     if (!bA->m[i][i]) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support for shifting an empty diagonal block, insert a matrix in block (%D,%D)",i,i);
537:     MatShift(bA->m[i][i],a);
538:   }
539:   return(0);
540: }

544: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left)
545: {
546:   Mat_Nest       *bA = (Mat_Nest*)A->data;
547:   Vec            *L,*R;
548:   MPI_Comm       comm;
549:   PetscInt       i,j;

553:   PetscObjectGetComm((PetscObject)A,&comm);
554:   if (right) {
555:     /* allocate R */
556:     PetscMalloc1(bA->nc, &R);
557:     /* Create the right vectors */
558:     for (j=0; j<bA->nc; j++) {
559:       for (i=0; i<bA->nr; i++) {
560:         if (bA->m[i][j]) {
561:           MatCreateVecs(bA->m[i][j],&R[j],NULL);
562:           break;
563:         }
564:       }
565:       if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
566:     }
567:     VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
568:     /* hand back control to the nest vector */
569:     for (j=0; j<bA->nc; j++) {
570:       VecDestroy(&R[j]);
571:     }
572:     PetscFree(R);
573:   }

575:   if (left) {
576:     /* allocate L */
577:     PetscMalloc1(bA->nr, &L);
578:     /* Create the left vectors */
579:     for (i=0; i<bA->nr; i++) {
580:       for (j=0; j<bA->nc; j++) {
581:         if (bA->m[i][j]) {
582:           MatCreateVecs(bA->m[i][j],NULL,&L[i]);
583:           break;
584:         }
585:       }
586:       if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
587:     }

589:     VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
590:     for (i=0; i<bA->nr; i++) {
591:       VecDestroy(&L[i]);
592:     }

594:     PetscFree(L);
595:   }
596:   return(0);
597: }

601: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
602: {
603:   Mat_Nest       *bA = (Mat_Nest*)A->data;
604:   PetscBool      isascii;
605:   PetscInt       i,j;

609:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
610:   if (isascii) {

612:     PetscViewerASCIIPrintf(viewer,"Matrix object: \n");
613:     PetscViewerASCIIPushTab(viewer);    /* push0 */
614:     PetscViewerASCIIPrintf(viewer, "type=nest, rows=%d, cols=%d \n",bA->nr,bA->nc);

616:     PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");
617:     for (i=0; i<bA->nr; i++) {
618:       for (j=0; j<bA->nc; j++) {
619:         MatType   type;
620:         char      name[256] = "",prefix[256] = "";
621:         PetscInt  NR,NC;
622:         PetscBool isNest = PETSC_FALSE;

624:         if (!bA->m[i][j]) {
625:           PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);
626:           continue;
627:         }
628:         MatGetSize(bA->m[i][j],&NR,&NC);
629:         MatGetType(bA->m[i][j], &type);
630:         if (((PetscObject)bA->m[i][j])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);}
631:         if (((PetscObject)bA->m[i][j])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);}
632:         PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);

634:         PetscViewerASCIIPrintf(viewer,"(%D,%D) : %s%stype=%s, rows=%D, cols=%D \n",i,j,name,prefix,type,NR,NC);

636:         if (isNest) {
637:           PetscViewerASCIIPushTab(viewer);  /* push1 */
638:           MatView(bA->m[i][j],viewer);
639:           PetscViewerASCIIPopTab(viewer);    /* pop1 */
640:         }
641:       }
642:     }
643:     PetscViewerASCIIPopTab(viewer);    /* pop0 */
644:   }
645:   return(0);
646: }

650: static PetscErrorCode MatZeroEntries_Nest(Mat A)
651: {
652:   Mat_Nest       *bA = (Mat_Nest*)A->data;
653:   PetscInt       i,j;

657:   for (i=0; i<bA->nr; i++) {
658:     for (j=0; j<bA->nc; j++) {
659:       if (!bA->m[i][j]) continue;
660:       MatZeroEntries(bA->m[i][j]);
661:     }
662:   }
663:   return(0);
664: }

668: static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str)
669: {
670:   Mat_Nest       *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data;
671:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

675:   if (nr != bB->nr || nc != bB->nc) SETERRQ4(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Cannot copy a Mat_Nest of block size (%D,%D) to a Mat_Nest of block size (%D,%D)",bB->nr,bB->nc,nr,nc);
676:   for (i=0; i<nr; i++) {
677:     for (j=0; j<nc; j++) {
678:       if (bA->m[i][j] && bB->m[i][j]) {
679:         MatCopy(bA->m[i][j],bB->m[i][j],str);
680:       } else if (bA->m[i][j] || bB->m[i][j]) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Matrix block does not exist at %D,%D",i,j);
681:     }
682:   }
683:   return(0);
684: }

688: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
689: {
690:   Mat_Nest       *bA = (Mat_Nest*)A->data;
691:   Mat            *b;
692:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

696:   PetscMalloc1(nr*nc,&b);
697:   for (i=0; i<nr; i++) {
698:     for (j=0; j<nc; j++) {
699:       if (bA->m[i][j]) {
700:         MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
701:       } else {
702:         b[i*nc+j] = NULL;
703:       }
704:     }
705:   }
706:   MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
707:   /* Give the new MatNest exclusive ownership */
708:   for (i=0; i<nr*nc; i++) {
709:     MatDestroy(&b[i]);
710:   }
711:   PetscFree(b);

713:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
714:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
715:   return(0);
716: }

718: /* nest api */
721: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
722: {
723:   Mat_Nest *bA = (Mat_Nest*)A->data;

726:   if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
727:   if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
728:   *mat = bA->m[idxm][jdxm];
729:   return(0);
730: }

734: /*@
735:  MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.

737:  Not collective

739:  Input Parameters:
740: +   A  - nest matrix
741: .   idxm - index of the matrix within the nest matrix
742: -   jdxm - index of the matrix within the nest matrix

744:  Output Parameter:
745: .   sub - matrix at index idxm,jdxm within the nest matrix

747:  Level: developer

749: .seealso: MatNestGetSize(), MatNestGetSubMats()
750: @*/
751: PetscErrorCode  MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
752: {

756:   PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
757:   return(0);
758: }

762: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
763: {
764:   Mat_Nest       *bA = (Mat_Nest*)A->data;
765:   PetscInt       m,n,M,N,mi,ni,Mi,Ni;

769:   if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
770:   if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
771:   MatGetLocalSize(mat,&m,&n);
772:   MatGetSize(mat,&M,&N);
773:   ISGetLocalSize(bA->isglobal.row[idxm],&mi);
774:   ISGetSize(bA->isglobal.row[idxm],&Mi);
775:   ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
776:   ISGetSize(bA->isglobal.col[jdxm],&Ni);
777:   if (M != Mi || N != Ni) SETERRQ4(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_INCOMP,"Submatrix dimension (%D,%D) incompatible with nest block (%D,%D)",M,N,Mi,Ni);
778:   if (m != mi || n != ni) SETERRQ4(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_INCOMP,"Submatrix local dimension (%D,%D) incompatible with nest block (%D,%D)",m,n,mi,ni);

780:   PetscObjectReference((PetscObject)mat);
781:   MatDestroy(&bA->m[idxm][jdxm]);
782:   bA->m[idxm][jdxm] = mat;
783:   return(0);
784: }

788: /*@
789:  MatNestSetSubMat - Set a single submatrix in the nest matrix.

791:  Logically collective on the submatrix communicator

793:  Input Parameters:
794: +   A  - nest matrix
795: .   idxm - index of the matrix within the nest matrix
796: .   jdxm - index of the matrix within the nest matrix
797: -   sub - matrix at index idxm,jdxm within the nest matrix

799:  Notes:
800:  The new submatrix must have the same size and communicator as that block of the nest.

802:  This increments the reference count of the submatrix.

804:  Level: developer

806: .seealso: MatNestSetSubMats(), MatNestGetSubMat()
807: @*/
808: PetscErrorCode  MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
809: {

813:   PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
814:   return(0);
815: }

819: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
820: {
821:   Mat_Nest *bA = (Mat_Nest*)A->data;

824:   if (M)   *M   = bA->nr;
825:   if (N)   *N   = bA->nc;
826:   if (mat) *mat = bA->m;
827:   return(0);
828: }

832: /*@C
833:  MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a nest matrix.

835:  Not collective

837:  Input Parameters:
838: .   A  - nest matrix

840:  Output Parameter:
841: +   M - number of rows in the nest matrix
842: .   N - number of cols in the nest matrix
843: -   mat - 2d array of matrices

845:  Notes:

847:  The user should not free the array mat.

849:  In Fortran, this routine has a calling sequence
850: $   call MatNestGetSubMats(A, M, N, mat, ierr)
851:  where the space allocated for the optional argument mat is assumed large enough (if provided).

853:  Level: developer

855: .seealso: MatNestGetSize(), MatNestGetSubMat()
856: @*/
857: PetscErrorCode  MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
858: {

862:   PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
863:   return(0);
864: }

868: PetscErrorCode  MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
869: {
870:   Mat_Nest *bA = (Mat_Nest*)A->data;

873:   if (M) *M = bA->nr;
874:   if (N) *N = bA->nc;
875:   return(0);
876: }

880: /*@
881:  MatNestGetSize - Returns the size of the nest matrix.

883:  Not collective

885:  Input Parameters:
886: .   A  - nest matrix

888:  Output Parameter:
889: +   M - number of rows in the nested mat
890: -   N - number of cols in the nested mat

892:  Notes:

894:  Level: developer

896: .seealso: MatNestGetSubMat(), MatNestGetSubMats()
897: @*/
898: PetscErrorCode  MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
899: {

903:   PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
904:   return(0);
905: }

909: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
910: {
911:   Mat_Nest *vs = (Mat_Nest*)A->data;
912:   PetscInt i;

915:   if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
916:   if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
917:   return(0);
918: }

922: /*@C
923:  MatNestGetISs - Returns the index sets partitioning the row and column spaces

925:  Not collective

927:  Input Parameters:
928: .   A  - nest matrix

930:  Output Parameter:
931: +   rows - array of row index sets
932: -   cols - array of column index sets

934:  Level: advanced

936:  Notes:
937:  The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.

939: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs()
940: @*/
941: PetscErrorCode  MatNestGetISs(Mat A,IS rows[],IS cols[])
942: {

947:   PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
948:   return(0);
949: }

953: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
954: {
955:   Mat_Nest *vs = (Mat_Nest*)A->data;
956:   PetscInt i;

959:   if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
960:   if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
961:   return(0);
962: }

966: /*@C
967:  MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces

969:  Not collective

971:  Input Parameters:
972: .   A  - nest matrix

974:  Output Parameter:
975: +   rows - array of row index sets (or NULL to ignore)
976: -   cols - array of column index sets (or NULL to ignore)

978:  Level: advanced

980:  Notes:
981:  The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.

983: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs()
984: @*/
985: PetscErrorCode  MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
986: {

991:   PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
992:   return(0);
993: }

997: PetscErrorCode  MatNestSetVecType_Nest(Mat A,VecType vtype)
998: {
1000:   PetscBool      flg;

1003:   PetscStrcmp(vtype,VECNEST,&flg);
1004:   /* In reality, this only distinguishes VECNEST and "other" */
1005:   if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1006:   else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
1007:   return(0);
1008: }

1012: /*@C
1013:  MatNestSetVecType - Sets the type of Vec returned by MatCreateVecs()

1015:  Not collective

1017:  Input Parameters:
1018: +  A  - nest matrix
1019: -  vtype - type to use for creating vectors

1021:  Notes:

1023:  Level: developer

1025: .seealso: MatCreateVecs()
1026: @*/
1027: PetscErrorCode  MatNestSetVecType(Mat A,VecType vtype)
1028: {

1032:   PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
1033:   return(0);
1034: }

1038: PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1039: {
1040:   Mat_Nest       *s = (Mat_Nest*)A->data;
1041:   PetscInt       i,j,m,n,M,N;

1045:   s->nr = nr;
1046:   s->nc = nc;

1048:   /* Create space for submatrices */
1049:   PetscMalloc1(nr,&s->m);
1050:   for (i=0; i<nr; i++) {
1051:     PetscMalloc1(nc,&s->m[i]);
1052:   }
1053:   for (i=0; i<nr; i++) {
1054:     for (j=0; j<nc; j++) {
1055:       s->m[i][j] = a[i*nc+j];
1056:       if (a[i*nc+j]) {
1057:         PetscObjectReference((PetscObject)a[i*nc+j]);
1058:       }
1059:     }
1060:   }

1062:   MatSetUp_NestIS_Private(A,nr,is_row,nc,is_col);

1064:   PetscMalloc1(nr,&s->row_len);
1065:   PetscMalloc1(nc,&s->col_len);
1066:   for (i=0; i<nr; i++) s->row_len[i]=-1;
1067:   for (j=0; j<nc; j++) s->col_len[j]=-1;

1069:   MatNestGetSizes_Private(A,&m,&n,&M,&N);

1071:   PetscLayoutSetSize(A->rmap,M);
1072:   PetscLayoutSetLocalSize(A->rmap,m);
1073:   PetscLayoutSetSize(A->cmap,N);
1074:   PetscLayoutSetLocalSize(A->cmap,n);

1076:   PetscLayoutSetUp(A->rmap);
1077:   PetscLayoutSetUp(A->cmap);

1079:   PetscCalloc2(nr,&s->left,nc,&s->right);
1080:   return(0);
1081: }

1085: /*@
1086:    MatNestSetSubMats - Sets the nested submatrices

1088:    Collective on Mat

1090:    Input Parameter:
1091: +  N - nested matrix
1092: .  nr - number of nested row blocks
1093: .  is_row - index sets for each nested row block, or NULL to make contiguous
1094: .  nc - number of nested column blocks
1095: .  is_col - index sets for each nested column block, or NULL to make contiguous
1096: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

1098:    Level: advanced

1100: .seealso: MatCreateNest(), MATNEST
1101: @*/
1102: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1103: {
1105:   PetscInt       i;

1109:   if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1110:   if (nr && is_row) {
1113:   }
1114:   if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1115:   if (nc && is_col) {
1118:   }
1120:   PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1121:   return(0);
1122: }

1126: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1127: {
1129:   PetscBool      flg;
1130:   PetscInt       i,j,m,mi,*ix;

1133:   for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1134:     if (islocal[i]) {
1135:       ISGetSize(islocal[i],&mi);
1136:       flg  = PETSC_TRUE;      /* We found a non-trivial entry */
1137:     } else {
1138:       ISGetSize(isglobal[i],&mi);
1139:     }
1140:     m += mi;
1141:   }
1142:   if (flg) {
1143:     PetscMalloc1(m,&ix);
1144:     for (i=0,n=0; i<n; i++) {
1145:       ISLocalToGlobalMapping smap = NULL;
1146:       VecScatter             scat;
1147:       IS                     isreq;
1148:       Vec                    lvec,gvec;
1149:       union {char padding[sizeof(PetscScalar)]; PetscInt integer;} *x;
1150:       Mat sub;

1152:       if (sizeof(*x) != sizeof(PetscScalar)) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"No support when scalars smaller than integers");
1153:       if (colflg) {
1154:         MatNestFindNonzeroSubMatRow(A,i,&sub);
1155:       } else {
1156:         MatNestFindNonzeroSubMatCol(A,i,&sub);
1157:       }
1158:       if (sub) {MatGetLocalToGlobalMapping(sub,&smap,NULL);}
1159:       if (islocal[i]) {
1160:         ISGetSize(islocal[i],&mi);
1161:       } else {
1162:         ISGetSize(isglobal[i],&mi);
1163:       }
1164:       for (j=0; j<mi; j++) ix[m+j] = j;
1165:       if (smap) {ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);}
1166:       /*
1167:         Now we need to extract the monolithic global indices that correspond to the given split global indices.
1168:         In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1169:         The approach here is ugly because it uses VecScatter to move indices.
1170:        */
1171:       VecCreateSeq(PETSC_COMM_SELF,mi,&lvec);
1172:       VecCreateMPI(((PetscObject)isglobal[i])->comm,mi,PETSC_DECIDE,&gvec);
1173:       ISCreateGeneral(((PetscObject)isglobal[i])->comm,mi,ix+m,PETSC_COPY_VALUES,&isreq);
1174:       VecScatterCreate(gvec,isreq,lvec,NULL,&scat);
1175:       VecGetArray(gvec,(PetscScalar**)&x);
1176:       for (j=0; j<mi; j++) x[j].integer = ix[m+j];
1177:       VecRestoreArray(gvec,(PetscScalar**)&x);
1178:       VecScatterBegin(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);
1179:       VecScatterEnd(scat,gvec,lvec,INSERT_VALUES,SCATTER_FORWARD);
1180:       VecGetArray(lvec,(PetscScalar**)&x);
1181:       for (j=0; j<mi; j++) ix[m+j] = x[j].integer;
1182:       VecRestoreArray(lvec,(PetscScalar**)&x);
1183:       VecDestroy(&lvec);
1184:       VecDestroy(&gvec);
1185:       ISDestroy(&isreq);
1186:       VecScatterDestroy(&scat);
1187:       m   += mi;
1188:     }
1189:     ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);
1190:   } else {
1191:     *ltog  = NULL;
1192:   }
1193:   return(0);
1194: }


1197: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1198: /*
1199:   nprocessors = NP
1200:   Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1201:        proc 0: => (g_0,h_0,)
1202:        proc 1: => (g_1,h_1,)
1203:        ...
1204:        proc nprocs-1: => (g_NP-1,h_NP-1,)

1206:             proc 0:                      proc 1:                    proc nprocs-1:
1207:     is[0] = (0,1,2,...,nlocal(g_0)-1)  (0,1,...,nlocal(g_1)-1)  (0,1,...,nlocal(g_NP-1))

1209:             proc 0:
1210:     is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1211:             proc 1:
1212:     is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)

1214:             proc NP-1:
1215:     is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1216: */
1219: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1220: {
1221:   Mat_Nest       *vs = (Mat_Nest*)A->data;
1222:   PetscInt       i,j,offset,n,nsum,bs;
1224:   Mat            sub = NULL;

1227:   PetscMalloc1(nr,&vs->isglobal.row);
1228:   PetscMalloc1(nc,&vs->isglobal.col);
1229:   if (is_row) { /* valid IS is passed in */
1230:     /* refs on is[] are incremeneted */
1231:     for (i=0; i<vs->nr; i++) {
1232:       PetscObjectReference((PetscObject)is_row[i]);

1234:       vs->isglobal.row[i] = is_row[i];
1235:     }
1236:   } else {                      /* Create the ISs by inspecting sizes of a submatrix in each row */
1237:     nsum = 0;
1238:     for (i=0; i<vs->nr; i++) {  /* Add up the local sizes to compute the aggregate offset */
1239:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1240:       if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i);
1241:       MatGetLocalSize(sub,&n,NULL);
1242:       if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1243:       nsum += n;
1244:     }
1245:     MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1246:     offset -= nsum;
1247:     for (i=0; i<vs->nr; i++) {
1248:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1249:       MatGetLocalSize(sub,&n,NULL);
1250:       MatGetBlockSize(sub,&bs);
1251:       ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);
1252:       ISSetBlockSize(vs->isglobal.row[i],bs);
1253:       offset += n;
1254:     }
1255:   }

1257:   if (is_col) { /* valid IS is passed in */
1258:     /* refs on is[] are incremeneted */
1259:     for (j=0; j<vs->nc; j++) {
1260:       PetscObjectReference((PetscObject)is_col[j]);

1262:       vs->isglobal.col[j] = is_col[j];
1263:     }
1264:   } else {                      /* Create the ISs by inspecting sizes of a submatrix in each column */
1265:     offset = A->cmap->rstart;
1266:     nsum   = 0;
1267:     for (j=0; j<vs->nc; j++) {
1268:       MatNestFindNonzeroSubMatCol(A,j,&sub);
1269:       if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i);
1270:       MatGetLocalSize(sub,NULL,&n);
1271:       if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1272:       nsum += n;
1273:     }
1274:     MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1275:     offset -= nsum;
1276:     for (j=0; j<vs->nc; j++) {
1277:       MatNestFindNonzeroSubMatCol(A,j,&sub);
1278:       MatGetLocalSize(sub,NULL,&n);
1279:       MatGetBlockSize(sub,&bs);
1280:       ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);
1281:       ISSetBlockSize(vs->isglobal.col[j],bs);
1282:       offset += n;
1283:     }
1284:   }

1286:   /* Set up the local ISs */
1287:   PetscMalloc1(vs->nr,&vs->islocal.row);
1288:   PetscMalloc1(vs->nc,&vs->islocal.col);
1289:   for (i=0,offset=0; i<vs->nr; i++) {
1290:     IS                     isloc;
1291:     ISLocalToGlobalMapping rmap = NULL;
1292:     PetscInt               nlocal,bs;
1293:     MatNestFindNonzeroSubMatRow(A,i,&sub);
1294:     if (sub) {MatGetLocalToGlobalMapping(sub,&rmap,NULL);}
1295:     if (rmap) {
1296:       MatGetBlockSize(sub,&bs);
1297:       ISLocalToGlobalMappingGetSize(rmap,&nlocal);
1298:       ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1299:       ISSetBlockSize(isloc,bs);
1300:     } else {
1301:       nlocal = 0;
1302:       isloc  = NULL;
1303:     }
1304:     vs->islocal.row[i] = isloc;
1305:     offset            += nlocal;
1306:   }
1307:   for (i=0,offset=0; i<vs->nc; i++) {
1308:     IS                     isloc;
1309:     ISLocalToGlobalMapping cmap = NULL;
1310:     PetscInt               nlocal,bs;
1311:     MatNestFindNonzeroSubMatCol(A,i,&sub);
1312:     if (sub) {MatGetLocalToGlobalMapping(sub,NULL,&cmap);}
1313:     if (cmap) {
1314:       MatGetBlockSize(sub,&bs);
1315:       ISLocalToGlobalMappingGetSize(cmap,&nlocal);
1316:       ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1317:       ISSetBlockSize(isloc,bs);
1318:     } else {
1319:       nlocal = 0;
1320:       isloc  = NULL;
1321:     }
1322:     vs->islocal.col[i] = isloc;
1323:     offset            += nlocal;
1324:   }

1326:   /* Set up the aggregate ISLocalToGlobalMapping */
1327:   {
1328:     ISLocalToGlobalMapping rmap,cmap;
1329:     MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);
1330:     MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);
1331:     if (rmap && cmap) {MatSetLocalToGlobalMapping(A,rmap,cmap);}
1332:     ISLocalToGlobalMappingDestroy(&rmap);
1333:     ISLocalToGlobalMappingDestroy(&cmap);
1334:   }

1336: #if defined(PETSC_USE_DEBUG)
1337:   for (i=0; i<vs->nr; i++) {
1338:     for (j=0; j<vs->nc; j++) {
1339:       PetscInt m,n,M,N,mi,ni,Mi,Ni;
1340:       Mat      B = vs->m[i][j];
1341:       if (!B) continue;
1342:       MatGetSize(B,&M,&N);
1343:       MatGetLocalSize(B,&m,&n);
1344:       ISGetSize(vs->isglobal.row[i],&Mi);
1345:       ISGetSize(vs->isglobal.col[j],&Ni);
1346:       ISGetLocalSize(vs->isglobal.row[i],&mi);
1347:       ISGetLocalSize(vs->isglobal.col[j],&ni);
1348:       if (M != Mi || N != Ni) SETERRQ6(PetscObjectComm((PetscObject)sub),PETSC_ERR_ARG_INCOMP,"Global sizes (%D,%D) of nested submatrix (%D,%D) do not agree with space defined by index sets (%D,%D)",M,N,i,j,Mi,Ni);
1349:       if (m != mi || n != ni) SETERRQ6(PetscObjectComm((PetscObject)sub),PETSC_ERR_ARG_INCOMP,"Local sizes (%D,%D) of nested submatrix (%D,%D) do not agree with space defined by index sets (%D,%D)",m,n,i,j,mi,ni);
1350:     }
1351:   }
1352: #endif

1354:   /* Set A->assembled if all non-null blocks are currently assembled */
1355:   for (i=0; i<vs->nr; i++) {
1356:     for (j=0; j<vs->nc; j++) {
1357:       if (vs->m[i][j] && !vs->m[i][j]->assembled) return(0);
1358:     }
1359:   }
1360:   A->assembled = PETSC_TRUE;
1361:   return(0);
1362: }

1366: /*@C
1367:    MatCreateNest - Creates a new matrix containing several nested submatrices, each stored separately

1369:    Collective on Mat

1371:    Input Parameter:
1372: +  comm - Communicator for the new Mat
1373: .  nr - number of nested row blocks
1374: .  is_row - index sets for each nested row block, or NULL to make contiguous
1375: .  nc - number of nested column blocks
1376: .  is_col - index sets for each nested column block, or NULL to make contiguous
1377: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

1379:    Output Parameter:
1380: .  B - new matrix

1382:    Level: advanced

1384: .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST
1385: @*/
1386: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1387: {
1388:   Mat            A;

1392:   *B   = 0;
1393:   MatCreate(comm,&A);
1394:   MatSetType(A,MATNEST);
1395:   A->preallocated = PETSC_TRUE;
1396:   MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1397:   *B   = A;
1398:   return(0);
1399: }

1403: PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1404: {
1406:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1407:   PetscInt       m,n,M,N,i,j,k,*dnnz,*onnz,rstart;
1408:   PetscInt       cstart,cend;
1409:   Mat            C;

1412:   MatGetSize(A,&M,&N);
1413:   MatGetLocalSize(A,&m,&n);
1414:   MatGetOwnershipRangeColumn(A,&cstart,&cend);
1415:   switch (reuse) {
1416:   case MAT_INITIAL_MATRIX:
1417:     MatCreate(PetscObjectComm((PetscObject)A),&C);
1418:     MatSetType(C,newtype);
1419:     MatSetSizes(C,m,n,M,N);
1420:     *newmat = C;
1421:     break;
1422:   case MAT_REUSE_MATRIX:
1423:     C = *newmat;
1424:     break;
1425:   default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse");
1426:   }
1427:   PetscMalloc1(2*m,&dnnz);
1428:   onnz = dnnz + m;
1429:   for (k=0; k<m; k++) {
1430:     dnnz[k] = 0;
1431:     onnz[k] = 0;
1432:   }
1433:   for (j=0; j<nest->nc; ++j) {
1434:     IS             bNis;
1435:     PetscInt       bN;
1436:     const PetscInt *bNindices;
1437:     /* Using global column indices and ISAllGather() is not scalable. */
1438:     ISAllGather(nest->isglobal.col[j], &bNis);
1439:     ISGetSize(bNis, &bN);
1440:     ISGetIndices(bNis,&bNindices);
1441:     for (i=0; i<nest->nr; ++i) {
1442:       PetscSF        bmsf;
1443:       PetscSFNode    *iremote;
1444:       Mat            B;
1445:       PetscInt       bm, *sub_dnnz,*sub_onnz, br;
1446:       const PetscInt *bmindices;
1447:       B = nest->m[i][j];
1448:       if (!B) continue;
1449:       ISGetLocalSize(nest->isglobal.row[i],&bm);
1450:       ISGetIndices(nest->isglobal.row[i],&bmindices);
1451:       PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);
1452:       PetscMalloc1(bm,&iremote);
1453:       PetscMalloc1(bm,&sub_dnnz);
1454:       PetscMalloc1(bm,&sub_onnz);
1455:       for (k = 0; k < bm; ++k){
1456:             sub_dnnz[k] = 0;
1457:             sub_onnz[k] = 0;
1458:       }
1459:       /*
1460:        Locate the owners for all of the locally-owned global row indices for this row block.
1461:        These determine the roots of PetscSF used to communicate preallocation data to row owners.
1462:        The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
1463:        */
1464:       MatGetOwnershipRange(B,&rstart,NULL);
1465:       for (br = 0; br < bm; ++br) {
1466:         PetscInt       row = bmindices[br], rowowner = 0, brncols, col;
1467:         const PetscInt *brcols;
1468:         PetscInt       rowrel = 0; /* row's relative index on its owner rank */
1469:         PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);
1470:         /* how many roots  */
1471:         iremote[br].rank = rowowner; iremote[br].index = rowrel;           /* edge from bmdnnz to dnnz */
1472:         /* get nonzero pattern */
1473:         MatGetRow(B,br+rstart,&brncols,&brcols,NULL);
1474:         for (k=0; k<brncols; k++) {
1475:           col  = bNindices[brcols[k]];
1476:           if(col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]){
1477:                 sub_dnnz[br]++;
1478:           }else{
1479:                 sub_onnz[br]++;
1480:           }
1481:         }
1482:         MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);
1483:       }
1484:       ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1485:       /* bsf will have to take care of disposing of bedges. */
1486:       PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
1487:       PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
1488:       PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
1489:       PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
1490:       PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
1491:       PetscFree(sub_dnnz);
1492:       PetscFree(sub_onnz);
1493:       PetscSFDestroy(&bmsf);
1494:     }
1495:     ISRestoreIndices(bNis,&bNindices);
1496:     ISDestroy(&bNis);
1497:   }
1498:   MatSeqAIJSetPreallocation(C,0,dnnz);
1499:   MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);
1500:   PetscFree(dnnz);

1502:   /* Fill by row */
1503:   for (j=0; j<nest->nc; ++j) {
1504:     /* Using global column indices and ISAllGather() is not scalable. */
1505:     IS             bNis;
1506:     PetscInt       bN;
1507:     const PetscInt *bNindices;
1508:     ISAllGather(nest->isglobal.col[j], &bNis);
1509:     ISGetSize(bNis,&bN);
1510:     ISGetIndices(bNis,&bNindices);
1511:     for (i=0; i<nest->nr; ++i) {
1512:       Mat            B;
1513:       PetscInt       bm, br;
1514:       const PetscInt *bmindices;
1515:       B = nest->m[i][j];
1516:       if (!B) continue;
1517:       ISGetLocalSize(nest->isglobal.row[i],&bm);
1518:       ISGetIndices(nest->isglobal.row[i],&bmindices);
1519:       MatGetOwnershipRange(B,&rstart,NULL);
1520:       for (br = 0; br < bm; ++br) {
1521:         PetscInt          row = bmindices[br], brncols,  *cols;
1522:         const PetscInt    *brcols;
1523:         const PetscScalar *brcoldata;
1524:         MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1525:         PetscMalloc1(brncols,&cols);
1526:         for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
1527:         /*
1528:           Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
1529:           Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
1530:          */
1531:         MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);
1532:         MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1533:         PetscFree(cols);
1534:       }
1535:       ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1536:     }
1537:     ISRestoreIndices(bNis,&bNindices);
1538:     ISDestroy(&bNis);
1539:   }
1540:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1541:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1542:   return(0);
1543: }

1545: /*MC
1546:   MATNEST - MATNEST = "nest" - Matrix type consisting of nested submatrices, each stored separately.

1548:   Level: intermediate

1550:   Notes:
1551:   This matrix type permits scalable use of PCFieldSplit and avoids the large memory costs of extracting submatrices.
1552:   It allows the use of symmetric and block formats for parts of multi-physics simulations.
1553:   It is usually used with DMComposite and DMCreateMatrix()

1555: .seealso: MatCreate(), MatType, MatCreateNest()
1556: M*/
1559: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
1560: {
1561:   Mat_Nest       *s;

1565:   PetscNewLog(A,&s);
1566:   A->data = (void*)s;

1568:   s->nr            = -1;
1569:   s->nc            = -1;
1570:   s->m             = NULL;
1571:   s->splitassembly = PETSC_FALSE;

1573:   PetscMemzero(A->ops,sizeof(*A->ops));

1575:   A->ops->mult                  = MatMult_Nest;
1576:   A->ops->multadd               = MatMultAdd_Nest;
1577:   A->ops->multtranspose         = MatMultTranspose_Nest;
1578:   A->ops->multtransposeadd      = MatMultTransposeAdd_Nest;
1579:   A->ops->assemblybegin         = MatAssemblyBegin_Nest;
1580:   A->ops->assemblyend           = MatAssemblyEnd_Nest;
1581:   A->ops->zeroentries           = MatZeroEntries_Nest;
1582:   A->ops->copy                  = MatCopy_Nest;
1583:   A->ops->duplicate             = MatDuplicate_Nest;
1584:   A->ops->getsubmatrix          = MatGetSubMatrix_Nest;
1585:   A->ops->destroy               = MatDestroy_Nest;
1586:   A->ops->view                  = MatView_Nest;
1587:   A->ops->getvecs               = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
1588:   A->ops->getlocalsubmatrix     = MatGetLocalSubMatrix_Nest;
1589:   A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
1590:   A->ops->getdiagonal           = MatGetDiagonal_Nest;
1591:   A->ops->diagonalscale         = MatDiagonalScale_Nest;
1592:   A->ops->scale                 = MatScale_Nest;
1593:   A->ops->shift                 = MatShift_Nest;

1595:   A->spptr        = 0;
1596:   A->assembled    = PETSC_FALSE;

1598:   /* expose Nest api's */
1599:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",   MatNestGetSubMat_Nest);
1600:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",   MatNestSetSubMat_Nest);
1601:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",  MatNestGetSubMats_Nest);
1602:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",     MatNestGetSize_Nest);
1603:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",      MatNestGetISs_Nest);
1604:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);
1605:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",  MatNestSetVecType_Nest);
1606:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",  MatNestSetSubMats_Nest);
1607:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",MatConvert_Nest_AIJ);

1609:   PetscObjectChangeTypeName((PetscObject)A,MATNEST);
1610:   return(0);
1611: }