Actual source code: matnest.c

petsc-master 2017-09-23
Report Typos and Errors

  2:  #include <../src/mat/impls/nest/matnestimpl.h>
  3:  #include <../src/mat/impls/aij/seq/aij.h>
  4:  #include <petscsf.h>

  6: static PetscErrorCode MatSetUp_NestIS_Private(Mat,PetscInt,const IS[],PetscInt,const IS[]);
  7: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left);
  8: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat,MatType,MatReuse,Mat*);

 10: /* 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 */
 37: static PetscErrorCode MatMult_Nest(Mat A,Vec x,Vec y)
 38: {
 39:   Mat_Nest       *bA = (Mat_Nest*)A->data;
 40:   Vec            *bx = bA->right,*by = bA->left;
 41:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

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

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

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

 88: static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y)
 89: {
 90:   Mat_Nest       *bA = (Mat_Nest*)A->data;
 91:   Vec            *bx = bA->left,*by = bA->right;
 92:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

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

111: static PetscErrorCode MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z)
112: {
113:   Mat_Nest       *bA = (Mat_Nest*)A->data;
114:   Vec            *bx = bA->left,*bz = bA->right;
115:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

119:   for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
120:   for (i=0; i<nc; i++) {VecGetSubVector(z,bA->isglobal.col[i],&bz[i]);}
121:   for (j=0; j<nc; j++) {
122:     if (y != z) {
123:       Vec by;
124:       VecGetSubVector(y,bA->isglobal.col[j],&by);
125:       VecCopy(by,bz[j]);
126:       VecRestoreSubVector(y,bA->isglobal.col[j],&by);
127:     }
128:     for (i=0; i<nr; i++) {
129:       if (!bA->m[i][j]) continue;
130:       /* z[j] <- y[j] + (A[i][j])^T * x[i] */
131:       MatMultTransposeAdd(bA->m[i][j],bx[i],bz[j],bz[j]);
132:     }
133:   }
134:   for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
135:   for (i=0; i<nc; i++) {VecRestoreSubVector(z,bA->isglobal.col[i],&bz[i]);}
136:   return(0);
137: }

139: static PetscErrorCode MatTranspose_Nest(Mat A,MatReuse reuse,Mat *B)
140: {
141:   Mat_Nest       *bA = (Mat_Nest*)A->data, *bC;
142:   Mat            C;
143:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

147:   if (reuse == MAT_INPLACE_MATRIX && nr != nc) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_SIZ,"Square nested matrix only for in-place");

149:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
150:     Mat *subs;
151:     IS  *is_row,*is_col;

153:     PetscCalloc1(nr * nc,&subs);
154:     PetscMalloc2(nr,&is_row,nc,&is_col);
155:     MatNestGetISs(A,is_row,is_col);
156:     if (reuse == MAT_INPLACE_MATRIX) {
157:       for (i=0; i<nr; i++) {
158:         for (j=0; j<nc; j++) {
159:           subs[i + nr * j] = bA->m[i][j];
160:         }
161:       }
162:     }

164:     MatCreateNest(PetscObjectComm((PetscObject)A),nc,is_col,nr,is_row,subs,&C);
165:     PetscFree(subs);
166:     PetscFree2(is_row,is_col);
167:   } else {
168:     C = *B;
169:   }

171:   bC = (Mat_Nest*)C->data;
172:   for (i=0; i<nr; i++) {
173:     for (j=0; j<nc; j++) {
174:       if (bA->m[i][j]) {
175:         MatTranspose(bA->m[i][j], reuse, &(bC->m[j][i]));
176:       } else {
177:         bC->m[j][i] = NULL;
178:       }
179:     }
180:   }

182:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
183:     *B = C;
184:   } else {
185:     MatHeaderMerge(A, &C);
186:   }
187:   return(0);
188: }

190: static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list)
191: {
193:   IS             *lst = *list;
194:   PetscInt       i;

197:   if (!lst) return(0);
198:   for (i=0; i<n; i++) if (lst[i]) {ISDestroy(&lst[i]);}
199:   PetscFree(lst);
200:   *list = NULL;
201:   return(0);
202: }

204: static PetscErrorCode MatDestroy_Nest(Mat A)
205: {
206:   Mat_Nest       *vs = (Mat_Nest*)A->data;
207:   PetscInt       i,j;

211:   /* release the matrices and the place holders */
212:   MatNestDestroyISList(vs->nr,&vs->isglobal.row);
213:   MatNestDestroyISList(vs->nc,&vs->isglobal.col);
214:   MatNestDestroyISList(vs->nr,&vs->islocal.row);
215:   MatNestDestroyISList(vs->nc,&vs->islocal.col);

217:   PetscFree(vs->row_len);
218:   PetscFree(vs->col_len);

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

222:   /* release the matrices and the place holders */
223:   if (vs->m) {
224:     for (i=0; i<vs->nr; i++) {
225:       for (j=0; j<vs->nc; j++) {
226:         MatDestroy(&vs->m[i][j]);
227:       }
228:       PetscFree(vs->m[i]);
229:     }
230:     PetscFree(vs->m);
231:   }
232:   PetscFree(A->data);

234:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",0);
235:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",0);
236:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",0);
237:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",0);
238:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",0);
239:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",0);
240:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",0);
241:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",0);
242:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",0);
243:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",0);
244:   return(0);
245: }

247: static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)
248: {
249:   Mat_Nest       *vs = (Mat_Nest*)A->data;
250:   PetscInt       i,j;

254:   for (i=0; i<vs->nr; i++) {
255:     for (j=0; j<vs->nc; j++) {
256:       if (vs->m[i][j]) {
257:         MatAssemblyBegin(vs->m[i][j],type);
258:         if (!vs->splitassembly) {
259:           /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
260:            * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
261:            * already performing an assembly, but the result would by more complicated and appears to offer less
262:            * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
263:            * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
264:            */
265:           MatAssemblyEnd(vs->m[i][j],type);
266:         }
267:       }
268:     }
269:   }
270:   return(0);
271: }

273: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
274: {
275:   Mat_Nest       *vs = (Mat_Nest*)A->data;
276:   PetscInt       i,j;

280:   for (i=0; i<vs->nr; i++) {
281:     for (j=0; j<vs->nc; j++) {
282:       if (vs->m[i][j]) {
283:         if (vs->splitassembly) {
284:           MatAssemblyEnd(vs->m[i][j],type);
285:         }
286:       }
287:     }
288:   }
289:   return(0);
290: }

292: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B)
293: {
295:   Mat_Nest       *vs = (Mat_Nest*)A->data;
296:   PetscInt       j;
297:   Mat            sub;

300:   sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
301:   for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j];
302:   if (sub) {MatSetUp(sub);}       /* Ensure that the sizes are available */
303:   *B = sub;
304:   return(0);
305: }

307: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B)
308: {
310:   Mat_Nest       *vs = (Mat_Nest*)A->data;
311:   PetscInt       i;
312:   Mat            sub;

315:   sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
316:   for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col];
317:   if (sub) {MatSetUp(sub);}       /* Ensure that the sizes are available */
318:   *B = sub;
319:   return(0);
320: }

322: static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found)
323: {
325:   PetscInt       i;
326:   PetscBool      flg;

332:   *found = -1;
333:   for (i=0; i<n; i++) {
334:     if (!list[i]) continue;
335:     ISEqual(list[i],is,&flg);
336:     if (flg) {
337:       *found = i;
338:       return(0);
339:     }
340:   }
341:   SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set");
342:   return(0);
343: }

345: /* Get a block row as a new MatNest */
346: static PetscErrorCode MatNestGetRow(Mat A,PetscInt row,Mat *B)
347: {
348:   Mat_Nest       *vs = (Mat_Nest*)A->data;
349:   char           keyname[256];

353:   *B   = NULL;
354:   PetscSNPrintf(keyname,sizeof(keyname),"NestRow_%D",row);
355:   PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);
356:   if (*B) return(0);

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

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

362:   PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);
363:   PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
364:   return(0);
365: }

367: static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
368: {
369:   Mat_Nest       *vs = (Mat_Nest*)A->data;
371:   PetscInt       row,col;
372:   PetscBool      same,isFullCol,isFullColGlobal;

375:   /* Check if full column space. This is a hack */
376:   isFullCol = PETSC_FALSE;
377:   PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&same);
378:   if (same) {
379:     PetscInt n,first,step,i,an,am,afirst,astep;
380:     ISStrideGetInfo(iscol,&first,&step);
381:     ISGetLocalSize(iscol,&n);
382:     isFullCol = PETSC_TRUE;
383:     for (i=0,an=A->cmap->rstart; i<vs->nc; i++) {
384:       ISStrideGetInfo(is->col[i],&afirst,&astep);
385:       ISGetLocalSize(is->col[i],&am);
386:       if (afirst != an || astep != step) isFullCol = PETSC_FALSE;
387:       an += am;
388:     }
389:     if (an != A->cmap->rstart+n) isFullCol = PETSC_FALSE;
390:   }
391:   MPIU_Allreduce(&isFullCol,&isFullColGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)iscol));

393:   if (isFullColGlobal && vs->nc > 1) {
394:     PetscInt row;
395:     MatNestFindIS(A,vs->nr,is->row,isrow,&row);
396:     MatNestGetRow(A,row,B);
397:   } else {
398:     MatNestFindIS(A,vs->nr,is->row,isrow,&row);
399:     MatNestFindIS(A,vs->nc,is->col,iscol,&col);
400:     if (!vs->m[row][col]) {
401:       PetscInt lr,lc;

403:       MatCreate(PetscObjectComm((PetscObject)A),&vs->m[row][col]);
404:       ISGetLocalSize(vs->isglobal.row[row],&lr);
405:       ISGetLocalSize(vs->isglobal.col[col],&lc);
406:       MatSetSizes(vs->m[row][col],lr,lc,PETSC_DECIDE,PETSC_DECIDE);
407:       MatSetUp(vs->m[row][col]);
408:       MatAssemblyBegin(vs->m[row][col],MAT_FINAL_ASSEMBLY);
409:       MatAssemblyEnd(vs->m[row][col],MAT_FINAL_ASSEMBLY);
410:     }
411:     *B = vs->m[row][col];
412:   }
413:   return(0);
414: }

416: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
417: {
419:   Mat_Nest       *vs = (Mat_Nest*)A->data;
420:   Mat            sub;

423:   MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
424:   switch (reuse) {
425:   case MAT_INITIAL_MATRIX:
426:     if (sub) { PetscObjectReference((PetscObject)sub); }
427:     *B = sub;
428:     break;
429:   case MAT_REUSE_MATRIX:
430:     if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call");
431:     break;
432:   case MAT_IGNORE_MATRIX:       /* Nothing to do */
433:     break;
434:   case MAT_INPLACE_MATRIX:       /* Nothing to do */
435:     SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet");
436:     break;
437:   }
438:   return(0);
439: }

441: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
442: {
444:   Mat_Nest       *vs = (Mat_Nest*)A->data;
445:   Mat            sub;

448:   MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
449:   /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
450:   if (sub) {PetscObjectReference((PetscObject)sub);}
451:   *B = sub;
452:   return(0);
453: }

455: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
456: {
458:   Mat_Nest       *vs = (Mat_Nest*)A->data;
459:   Mat            sub;

462:   MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
463:   if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten");
464:   if (sub) {
465:     if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times");
466:     MatDestroy(B);
467:   }
468:   return(0);
469: }

471: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
472: {
473:   Mat_Nest       *bA = (Mat_Nest*)A->data;
474:   PetscInt       i;

478:   for (i=0; i<bA->nr; i++) {
479:     Vec bv;
480:     VecGetSubVector(v,bA->isglobal.row[i],&bv);
481:     if (bA->m[i][i]) {
482:       MatGetDiagonal(bA->m[i][i],bv);
483:     } else {
484:       VecSet(bv,0.0);
485:     }
486:     VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
487:   }
488:   return(0);
489: }

491: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
492: {
493:   Mat_Nest       *bA = (Mat_Nest*)A->data;
494:   Vec            bl,*br;
495:   PetscInt       i,j;

499:   PetscCalloc1(bA->nc,&br);
500:   if (r) {
501:     for (j=0; j<bA->nc; j++) {VecGetSubVector(r,bA->isglobal.col[j],&br[j]);}
502:   }
503:   bl = NULL;
504:   for (i=0; i<bA->nr; i++) {
505:     if (l) {
506:       VecGetSubVector(l,bA->isglobal.row[i],&bl);
507:     }
508:     for (j=0; j<bA->nc; j++) {
509:       if (bA->m[i][j]) {
510:         MatDiagonalScale(bA->m[i][j],bl,br[j]);
511:       }
512:     }
513:     if (l) {
514:       VecRestoreSubVector(l,bA->isglobal.row[i],&bl);
515:     }
516:   }
517:   if (r) {
518:     for (j=0; j<bA->nc; j++) {VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);}
519:   }
520:   PetscFree(br);
521:   return(0);
522: }

524: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
525: {
526:   Mat_Nest       *bA = (Mat_Nest*)A->data;
527:   PetscInt       i,j;

531:   for (i=0; i<bA->nr; i++) {
532:     for (j=0; j<bA->nc; j++) {
533:       if (bA->m[i][j]) {
534:         MatScale(bA->m[i][j],a);
535:       }
536:     }
537:   }
538:   return(0);
539: }

541: static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
542: {
543:   Mat_Nest       *bA = (Mat_Nest*)A->data;
544:   PetscInt       i;

548:   for (i=0; i<bA->nr; i++) {
549:     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);
550:     MatShift(bA->m[i][i],a);
551:   }
552:   return(0);
553: }

555: static PetscErrorCode MatDiagonalSet_Nest(Mat A,Vec D,InsertMode is)
556: {
557:   Mat_Nest       *bA = (Mat_Nest*)A->data;
558:   PetscInt       i;

562:   for (i=0; i<bA->nr; i++) {
563:     Vec bv;
564:     VecGetSubVector(D,bA->isglobal.row[i],&bv);
565:     if (bA->m[i][i]) {
566:       MatDiagonalSet(bA->m[i][i],bv,is);
567:     }
568:     VecRestoreSubVector(D,bA->isglobal.row[i],&bv);
569:   }
570:   return(0);
571: }

573: static PetscErrorCode MatSetRandom_Nest(Mat A,PetscRandom rctx)
574: {
575:   Mat_Nest       *bA = (Mat_Nest*)A->data;
576:   PetscInt       i,j;

580:   for (i=0; i<bA->nr; i++) {
581:     for (j=0; j<bA->nc; j++) {
582:       if (bA->m[i][j]) {
583:         MatSetRandom(bA->m[i][j],rctx);
584:       }
585:     }
586:   }
587:   return(0);
588: }

590: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left)
591: {
592:   Mat_Nest       *bA = (Mat_Nest*)A->data;
593:   Vec            *L,*R;
594:   MPI_Comm       comm;
595:   PetscInt       i,j;

599:   PetscObjectGetComm((PetscObject)A,&comm);
600:   if (right) {
601:     /* allocate R */
602:     PetscMalloc1(bA->nc, &R);
603:     /* Create the right vectors */
604:     for (j=0; j<bA->nc; j++) {
605:       for (i=0; i<bA->nr; i++) {
606:         if (bA->m[i][j]) {
607:           MatCreateVecs(bA->m[i][j],&R[j],NULL);
608:           break;
609:         }
610:       }
611:       if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
612:     }
613:     VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
614:     /* hand back control to the nest vector */
615:     for (j=0; j<bA->nc; j++) {
616:       VecDestroy(&R[j]);
617:     }
618:     PetscFree(R);
619:   }

621:   if (left) {
622:     /* allocate L */
623:     PetscMalloc1(bA->nr, &L);
624:     /* Create the left vectors */
625:     for (i=0; i<bA->nr; i++) {
626:       for (j=0; j<bA->nc; j++) {
627:         if (bA->m[i][j]) {
628:           MatCreateVecs(bA->m[i][j],NULL,&L[i]);
629:           break;
630:         }
631:       }
632:       if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
633:     }

635:     VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
636:     for (i=0; i<bA->nr; i++) {
637:       VecDestroy(&L[i]);
638:     }

640:     PetscFree(L);
641:   }
642:   return(0);
643: }

645: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
646: {
647:   Mat_Nest       *bA = (Mat_Nest*)A->data;
648:   PetscBool      isascii;
649:   PetscInt       i,j;

653:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
654:   if (isascii) {

656:     PetscViewerASCIIPrintf(viewer,"Matrix object: \n");
657:     PetscViewerASCIIPushTab(viewer);
658:     PetscViewerASCIIPrintf(viewer, "type=nest, rows=%D, cols=%D \n",bA->nr,bA->nc);

660:     PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");
661:     for (i=0; i<bA->nr; i++) {
662:       for (j=0; j<bA->nc; j++) {
663:         MatType   type;
664:         char      name[256] = "",prefix[256] = "";
665:         PetscInt  NR,NC;
666:         PetscBool isNest = PETSC_FALSE;

668:         if (!bA->m[i][j]) {
669:           PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);
670:           continue;
671:         }
672:         MatGetSize(bA->m[i][j],&NR,&NC);
673:         MatGetType(bA->m[i][j], &type);
674:         if (((PetscObject)bA->m[i][j])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);}
675:         if (((PetscObject)bA->m[i][j])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);}
676:         PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);

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

680:         if (isNest) {
681:           PetscViewerASCIIPushTab(viewer);  /* push1 */
682:           MatView(bA->m[i][j],viewer);
683:           PetscViewerASCIIPopTab(viewer);    /* pop1 */
684:         }
685:       }
686:     }
687:     PetscViewerASCIIPopTab(viewer);    /* pop0 */
688:   }
689:   return(0);
690: }

692: static PetscErrorCode MatZeroEntries_Nest(Mat A)
693: {
694:   Mat_Nest       *bA = (Mat_Nest*)A->data;
695:   PetscInt       i,j;

699:   for (i=0; i<bA->nr; i++) {
700:     for (j=0; j<bA->nc; j++) {
701:       if (!bA->m[i][j]) continue;
702:       MatZeroEntries(bA->m[i][j]);
703:     }
704:   }
705:   return(0);
706: }

708: static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str)
709: {
710:   Mat_Nest       *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data;
711:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

715:   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);
716:   for (i=0; i<nr; i++) {
717:     for (j=0; j<nc; j++) {
718:       if (bA->m[i][j] && bB->m[i][j]) {
719:         MatCopy(bA->m[i][j],bB->m[i][j],str);
720:       } 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);
721:     }
722:   }
723:   PetscObjectStateIncrease((PetscObject)B);
724:   return(0);
725: }

727: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
728: {
729:   Mat_Nest       *bA = (Mat_Nest*)A->data;
730:   Mat            *b;
731:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

735:   PetscMalloc1(nr*nc,&b);
736:   for (i=0; i<nr; i++) {
737:     for (j=0; j<nc; j++) {
738:       if (bA->m[i][j]) {
739:         MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
740:       } else {
741:         b[i*nc+j] = NULL;
742:       }
743:     }
744:   }
745:   MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
746:   /* Give the new MatNest exclusive ownership */
747:   for (i=0; i<nr*nc; i++) {
748:     MatDestroy(&b[i]);
749:   }
750:   PetscFree(b);

752:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
753:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
754:   return(0);
755: }

757: /* nest api */
758: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
759: {
760:   Mat_Nest *bA = (Mat_Nest*)A->data;

763:   if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
764:   if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
765:   *mat = bA->m[idxm][jdxm];
766:   return(0);
767: }

769: /*@
770:  MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.

772:  Not collective

774:  Input Parameters:
775: +   A  - nest matrix
776: .   idxm - index of the matrix within the nest matrix
777: -   jdxm - index of the matrix within the nest matrix

779:  Output Parameter:
780: .   sub - matrix at index idxm,jdxm within the nest matrix

782:  Level: developer

784: .seealso: MatNestGetSize(), MatNestGetSubMats()
785: @*/
786: PetscErrorCode  MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
787: {

791:   PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
792:   return(0);
793: }

795: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
796: {
797:   Mat_Nest       *bA = (Mat_Nest*)A->data;
798:   PetscInt       m,n,M,N,mi,ni,Mi,Ni;

802:   if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
803:   if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
804:   MatGetLocalSize(mat,&m,&n);
805:   MatGetSize(mat,&M,&N);
806:   ISGetLocalSize(bA->isglobal.row[idxm],&mi);
807:   ISGetSize(bA->isglobal.row[idxm],&Mi);
808:   ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
809:   ISGetSize(bA->isglobal.col[jdxm],&Ni);
810:   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);
811:   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);

813:   PetscObjectReference((PetscObject)mat);
814:   MatDestroy(&bA->m[idxm][jdxm]);
815:   bA->m[idxm][jdxm] = mat;
816:   return(0);
817: }

819: /*@
820:  MatNestSetSubMat - Set a single submatrix in the nest matrix.

822:  Logically collective on the submatrix communicator

824:  Input Parameters:
825: +   A  - nest matrix
826: .   idxm - index of the matrix within the nest matrix
827: .   jdxm - index of the matrix within the nest matrix
828: -   sub - matrix at index idxm,jdxm within the nest matrix

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

833:  This increments the reference count of the submatrix.

835:  Level: developer

837: .seealso: MatNestSetSubMats(), MatNestGetSubMats()
838: @*/
839: PetscErrorCode  MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
840: {

844:   PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
845:   return(0);
846: }

848: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
849: {
850:   Mat_Nest *bA = (Mat_Nest*)A->data;

853:   if (M)   *M   = bA->nr;
854:   if (N)   *N   = bA->nc;
855:   if (mat) *mat = bA->m;
856:   return(0);
857: }

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

862:  Not collective

864:  Input Parameters:
865: .   A  - nest matrix

867:  Output Parameter:
868: +   M - number of rows in the nest matrix
869: .   N - number of cols in the nest matrix
870: -   mat - 2d array of matrices

872:  Notes:

874:  The user should not free the array mat.

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

880:  Level: developer

882: .seealso: MatNestGetSize(), MatNestGetSubMat()
883: @*/
884: PetscErrorCode  MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
885: {

889:   PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
890:   return(0);
891: }

893: PetscErrorCode  MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
894: {
895:   Mat_Nest *bA = (Mat_Nest*)A->data;

898:   if (M) *M = bA->nr;
899:   if (N) *N = bA->nc;
900:   return(0);
901: }

903: /*@
904:  MatNestGetSize - Returns the size of the nest matrix.

906:  Not collective

908:  Input Parameters:
909: .   A  - nest matrix

911:  Output Parameter:
912: +   M - number of rows in the nested mat
913: -   N - number of cols in the nested mat

915:  Notes:

917:  Level: developer

919: .seealso: MatNestGetSubMat(), MatNestGetSubMats()
920: @*/
921: PetscErrorCode  MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
922: {

926:   PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
927:   return(0);
928: }

930: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
931: {
932:   Mat_Nest *vs = (Mat_Nest*)A->data;
933:   PetscInt i;

936:   if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
937:   if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
938:   return(0);
939: }

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

944:  Not collective

946:  Input Parameters:
947: .   A  - nest matrix

949:  Output Parameter:
950: +   rows - array of row index sets
951: -   cols - array of column index sets

953:  Level: advanced

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

958: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs()
959: @*/
960: PetscErrorCode  MatNestGetISs(Mat A,IS rows[],IS cols[])
961: {

966:   PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
967:   return(0);
968: }

970: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
971: {
972:   Mat_Nest *vs = (Mat_Nest*)A->data;
973:   PetscInt i;

976:   if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
977:   if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
978:   return(0);
979: }

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

984:  Not collective

986:  Input Parameters:
987: .   A  - nest matrix

989:  Output Parameter:
990: +   rows - array of row index sets (or NULL to ignore)
991: -   cols - array of column index sets (or NULL to ignore)

993:  Level: advanced

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

998: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs()
999: @*/
1000: PetscErrorCode  MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
1001: {

1006:   PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1007:   return(0);
1008: }

1010: PetscErrorCode  MatNestSetVecType_Nest(Mat A,VecType vtype)
1011: {
1013:   PetscBool      flg;

1016:   PetscStrcmp(vtype,VECNEST,&flg);
1017:   /* In reality, this only distinguishes VECNEST and "other" */
1018:   if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1019:   else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
1020:   return(0);
1021: }

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

1026:  Not collective

1028:  Input Parameters:
1029: +  A  - nest matrix
1030: -  vtype - type to use for creating vectors

1032:  Notes:

1034:  Level: developer

1036: .seealso: MatCreateVecs()
1037: @*/
1038: PetscErrorCode  MatNestSetVecType(Mat A,VecType vtype)
1039: {

1043:   PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
1044:   return(0);
1045: }

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

1054:   s->nr = nr;
1055:   s->nc = nc;

1057:   /* Create space for submatrices */
1058:   PetscMalloc1(nr,&s->m);
1059:   for (i=0; i<nr; i++) {
1060:     PetscMalloc1(nc,&s->m[i]);
1061:   }
1062:   for (i=0; i<nr; i++) {
1063:     for (j=0; j<nc; j++) {
1064:       s->m[i][j] = a[i*nc+j];
1065:       if (a[i*nc+j]) {
1066:         PetscObjectReference((PetscObject)a[i*nc+j]);
1067:       }
1068:     }
1069:   }

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

1073:   PetscMalloc1(nr,&s->row_len);
1074:   PetscMalloc1(nc,&s->col_len);
1075:   for (i=0; i<nr; i++) s->row_len[i]=-1;
1076:   for (j=0; j<nc; j++) s->col_len[j]=-1;

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

1080:   PetscLayoutSetSize(A->rmap,M);
1081:   PetscLayoutSetLocalSize(A->rmap,m);
1082:   PetscLayoutSetSize(A->cmap,N);
1083:   PetscLayoutSetLocalSize(A->cmap,n);

1085:   PetscLayoutSetUp(A->rmap);
1086:   PetscLayoutSetUp(A->cmap);

1088:   PetscCalloc2(nr,&s->left,nc,&s->right);
1089:   return(0);
1090: }

1092: /*@
1093:    MatNestSetSubMats - Sets the nested submatrices

1095:    Collective on Mat

1097:    Input Parameter:
1098: +  N - nested matrix
1099: .  nr - number of nested row blocks
1100: .  is_row - index sets for each nested row block, or NULL to make contiguous
1101: .  nc - number of nested column blocks
1102: .  is_col - index sets for each nested column block, or NULL to make contiguous
1103: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

1105:    Level: advanced

1107: .seealso: MatCreateNest(), MATNEST
1108: @*/
1109: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1110: {
1112:   PetscInt       i,nr_nc;

1116:   if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1117:   if (nr && is_row) {
1120:   }
1121:   if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1122:   if (nc && is_col) {
1125:   }
1126:   nr_nc=nr*nc;
1128:   PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1129:   return(0);
1130: }

1132: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1133: {
1135:   PetscBool      flg;
1136:   PetscInt       i,j,m,mi,*ix;

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

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


1203: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1204: /*
1205:   nprocessors = NP
1206:   Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1207:        proc 0: => (g_0,h_0,)
1208:        proc 1: => (g_1,h_1,)
1209:        ...
1210:        proc nprocs-1: => (g_NP-1,h_NP-1,)

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

1215:             proc 0:
1216:     is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1217:             proc 1:
1218:     is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)

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

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

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

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

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

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

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

1340: #if defined(PETSC_USE_DEBUG)
1341:   for (i=0; i<vs->nr; i++) {
1342:     for (j=0; j<vs->nc; j++) {
1343:       PetscInt m,n,M,N,mi,ni,Mi,Ni;
1344:       Mat      B = vs->m[i][j];
1345:       if (!B) continue;
1346:       MatGetSize(B,&M,&N);
1347:       MatGetLocalSize(B,&m,&n);
1348:       ISGetSize(vs->isglobal.row[i],&Mi);
1349:       ISGetSize(vs->isglobal.col[j],&Ni);
1350:       ISGetLocalSize(vs->isglobal.row[i],&mi);
1351:       ISGetLocalSize(vs->isglobal.col[j],&ni);
1352:       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);
1353:       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);
1354:     }
1355:   }
1356: #endif

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

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

1371:    Collective on Mat

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

1381:    Output Parameter:
1382: .  B - new matrix

1384:    Level: advanced

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

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

1403: static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1404: {
1405:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1406:   Mat            *trans;
1407:   PetscScalar    **avv;
1408:   PetscScalar    *vv;
1409:   PetscInt       **aii,**ajj;
1410:   PetscInt       *ii,*jj,*ci;
1411:   PetscInt       nr,nc,nnz,i,j;
1412:   PetscBool      done;

1416:   MatGetSize(A,&nr,&nc);
1417:   if (reuse == MAT_REUSE_MATRIX) {
1418:     PetscInt rnr;

1420:     MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);
1421:     if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ");
1422:     if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows");
1423:     MatSeqAIJGetArray(*newmat,&vv);
1424:   }
1425:   /* extract CSR for nested SeqAIJ matrices */
1426:   nnz  = 0;
1427:   PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);
1428:   for (i=0; i<nest->nr; ++i) {
1429:     for (j=0; j<nest->nc; ++j) {
1430:       Mat B = nest->m[i][j];
1431:       if (B) {
1432:         PetscScalar *naa;
1433:         PetscInt    *nii,*njj,nnr;
1434:         PetscBool   istrans;

1436:         PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1437:         if (istrans) {
1438:           Mat Bt;

1440:           MatTransposeGetMat(B,&Bt);
1441:           MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);
1442:           B    = trans[i*nest->nc+j];
1443:         }
1444:         MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);
1445:         if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ");
1446:         MatSeqAIJGetArray(B,&naa);
1447:         nnz += nii[nnr];

1449:         aii[i*nest->nc+j] = nii;
1450:         ajj[i*nest->nc+j] = njj;
1451:         avv[i*nest->nc+j] = naa;
1452:       }
1453:     }
1454:   }
1455:   if (reuse != MAT_REUSE_MATRIX) {
1456:     PetscMalloc1(nr+1,&ii);
1457:     PetscMalloc1(nnz,&jj);
1458:     PetscMalloc1(nnz,&vv);
1459:   } else {
1460:     if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros");
1461:   }

1463:   /* new row pointer */
1464:   PetscMemzero(ii,(nr+1)*sizeof(PetscInt));
1465:   for (i=0; i<nest->nr; ++i) {
1466:     PetscInt       ncr,rst;

1468:     ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1469:     ISGetLocalSize(nest->isglobal.row[i],&ncr);
1470:     for (j=0; j<nest->nc; ++j) {
1471:       if (aii[i*nest->nc+j]) {
1472:         PetscInt    *nii = aii[i*nest->nc+j];
1473:         PetscInt    ir;

1475:         for (ir=rst; ir<ncr+rst; ++ir) {
1476:           ii[ir+1] += nii[1]-nii[0];
1477:           nii++;
1478:         }
1479:       }
1480:     }
1481:   }
1482:   for (i=0; i<nr; i++) ii[i+1] += ii[i];

1484:   /* construct CSR for the new matrix */
1485:   PetscCalloc1(nr,&ci);
1486:   for (i=0; i<nest->nr; ++i) {
1487:     PetscInt       ncr,rst;

1489:     ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1490:     ISGetLocalSize(nest->isglobal.row[i],&ncr);
1491:     for (j=0; j<nest->nc; ++j) {
1492:       if (aii[i*nest->nc+j]) {
1493:         PetscScalar *nvv = avv[i*nest->nc+j];
1494:         PetscInt    *nii = aii[i*nest->nc+j];
1495:         PetscInt    *njj = ajj[i*nest->nc+j];
1496:         PetscInt    ir,cst;

1498:         ISStrideGetInfo(nest->isglobal.col[j],&cst,NULL);
1499:         for (ir=rst; ir<ncr+rst; ++ir) {
1500:           PetscInt ij,rsize = nii[1]-nii[0],ist = ii[ir]+ci[ir];

1502:           for (ij=0;ij<rsize;ij++) {
1503:             jj[ist+ij] = *njj+cst;
1504:             vv[ist+ij] = *nvv;
1505:             njj++;
1506:             nvv++;
1507:           }
1508:           ci[ir] += rsize;
1509:           nii++;
1510:         }
1511:       }
1512:     }
1513:   }
1514:   PetscFree(ci);

1516:   /* restore info */
1517:   for (i=0; i<nest->nr; ++i) {
1518:     for (j=0; j<nest->nc; ++j) {
1519:       Mat B = nest->m[i][j];
1520:       if (B) {
1521:         PetscInt nnr = 0, k = i*nest->nc+j;

1523:         B    = (trans[k] ? trans[k] : B);
1524:         MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);
1525:         if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ");
1526:         MatSeqAIJRestoreArray(B,&avv[k]);
1527:         MatDestroy(&trans[k]);
1528:       }
1529:     }
1530:   }
1531:   PetscFree4(aii,ajj,avv,trans);

1533:   /* finalize newmat */
1534:   if (reuse == MAT_INITIAL_MATRIX) {
1535:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);
1536:   } else if (reuse == MAT_INPLACE_MATRIX) {
1537:     Mat B;

1539:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);
1540:     MatHeaderReplace(A,&B);
1541:   }
1542:   MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1543:   MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1544:   {
1545:     Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data);
1546:     a->free_a     = PETSC_TRUE;
1547:     a->free_ij    = PETSC_TRUE;
1548:   }
1549:   return(0);
1550: }

1552: PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1553: {
1555:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1556:   PetscInt       m,n,M,N,i,j,k,*dnnz,*onnz,rstart;
1557:   PetscInt       cstart,cend;
1558:   PetscMPIInt    size;
1559:   Mat            C;

1562:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1563:   if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */
1564:     PetscInt  nf;
1565:     PetscBool fast;

1567:     PetscStrcmp(newtype,MATAIJ,&fast);
1568:     if (!fast) {
1569:       PetscStrcmp(newtype,MATSEQAIJ,&fast);
1570:     }
1571:     for (i=0; i<nest->nr && fast; ++i) {
1572:       for (j=0; j<nest->nc && fast; ++j) {
1573:         Mat B = nest->m[i][j];
1574:         if (B) {
1575:           PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);
1576:           if (!fast) {
1577:             PetscBool istrans;

1579:             PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1580:             if (istrans) {
1581:               Mat Bt;

1583:               MatTransposeGetMat(B,&Bt);
1584:               PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);
1585:             }
1586:           }
1587:         }
1588:       }
1589:     }
1590:     for (i=0, nf=0; i<nest->nr && fast; ++i) {
1591:       PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);
1592:       if (fast) {
1593:         PetscInt f,s;

1595:         ISStrideGetInfo(nest->isglobal.row[i],&f,&s);
1596:         if (f != nf || s != 1) { fast = PETSC_FALSE; }
1597:         else {
1598:           ISGetSize(nest->isglobal.row[i],&f);
1599:           nf  += f;
1600:         }
1601:       }
1602:     }
1603:     for (i=0, nf=0; i<nest->nc && fast; ++i) {
1604:       PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);
1605:       if (fast) {
1606:         PetscInt f,s;

1608:         ISStrideGetInfo(nest->isglobal.col[i],&f,&s);
1609:         if (f != nf || s != 1) { fast = PETSC_FALSE; }
1610:         else {
1611:           ISGetSize(nest->isglobal.col[i],&f);
1612:           nf  += f;
1613:         }
1614:       }
1615:     }
1616:     if (fast) {
1617:       MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);
1618:       return(0);
1619:     }
1620:   }
1621:   MatGetSize(A,&M,&N);
1622:   MatGetLocalSize(A,&m,&n);
1623:   MatGetOwnershipRangeColumn(A,&cstart,&cend);
1624:   switch (reuse) {
1625:   case MAT_INITIAL_MATRIX:
1626:     MatCreate(PetscObjectComm((PetscObject)A),&C);
1627:     MatSetType(C,newtype);
1628:     MatSetSizes(C,m,n,M,N);
1629:     *newmat = C;
1630:     break;
1631:   case MAT_REUSE_MATRIX:
1632:     C = *newmat;
1633:     break;
1634:   default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse");
1635:   }
1636:   PetscMalloc1(2*m,&dnnz);
1637:   onnz = dnnz + m;
1638:   for (k=0; k<m; k++) {
1639:     dnnz[k] = 0;
1640:     onnz[k] = 0;
1641:   }
1642:   for (j=0; j<nest->nc; ++j) {
1643:     IS             bNis;
1644:     PetscInt       bN;
1645:     const PetscInt *bNindices;
1646:     /* Using global column indices and ISAllGather() is not scalable. */
1647:     ISAllGather(nest->isglobal.col[j], &bNis);
1648:     ISGetSize(bNis, &bN);
1649:     ISGetIndices(bNis,&bNindices);
1650:     for (i=0; i<nest->nr; ++i) {
1651:       PetscSF        bmsf;
1652:       PetscSFNode    *iremote;
1653:       Mat            B;
1654:       PetscInt       bm, *sub_dnnz,*sub_onnz, br;
1655:       const PetscInt *bmindices;
1656:       B = nest->m[i][j];
1657:       if (!B) continue;
1658:       ISGetLocalSize(nest->isglobal.row[i],&bm);
1659:       ISGetIndices(nest->isglobal.row[i],&bmindices);
1660:       PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);
1661:       PetscMalloc1(bm,&iremote);
1662:       PetscMalloc1(bm,&sub_dnnz);
1663:       PetscMalloc1(bm,&sub_onnz);
1664:       for (k = 0; k < bm; ++k){
1665:             sub_dnnz[k] = 0;
1666:             sub_onnz[k] = 0;
1667:       }
1668:       /*
1669:        Locate the owners for all of the locally-owned global row indices for this row block.
1670:        These determine the roots of PetscSF used to communicate preallocation data to row owners.
1671:        The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
1672:        */
1673:       MatGetOwnershipRange(B,&rstart,NULL);
1674:       for (br = 0; br < bm; ++br) {
1675:         PetscInt       row = bmindices[br], rowowner = 0, brncols, col;
1676:         const PetscInt *brcols;
1677:         PetscInt       rowrel = 0; /* row's relative index on its owner rank */
1678:         PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);
1679:         /* how many roots  */
1680:         iremote[br].rank = rowowner; iremote[br].index = rowrel;           /* edge from bmdnnz to dnnz */
1681:         /* get nonzero pattern */
1682:         MatGetRow(B,br+rstart,&brncols,&brcols,NULL);
1683:         for (k=0; k<brncols; k++) {
1684:           col  = bNindices[brcols[k]];
1685:           if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) {
1686:             sub_dnnz[br]++;
1687:           } else {
1688:             sub_onnz[br]++;
1689:           }
1690:         }
1691:         MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);
1692:       }
1693:       ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1694:       /* bsf will have to take care of disposing of bedges. */
1695:       PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
1696:       PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
1697:       PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
1698:       PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
1699:       PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
1700:       PetscFree(sub_dnnz);
1701:       PetscFree(sub_onnz);
1702:       PetscSFDestroy(&bmsf);
1703:     }
1704:     ISRestoreIndices(bNis,&bNindices);
1705:     ISDestroy(&bNis);
1706:   }
1707:   /* Resize preallocation if overestimated */
1708:   for (i=0;i<m;i++) {
1709:     dnnz[i] = PetscMin(dnnz[i],A->cmap->n);
1710:     onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n);
1711:   }
1712:   MatSeqAIJSetPreallocation(C,0,dnnz);
1713:   MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);
1714:   PetscFree(dnnz);

1716:   /* Fill by row */
1717:   for (j=0; j<nest->nc; ++j) {
1718:     /* Using global column indices and ISAllGather() is not scalable. */
1719:     IS             bNis;
1720:     PetscInt       bN;
1721:     const PetscInt *bNindices;
1722:     ISAllGather(nest->isglobal.col[j], &bNis);
1723:     ISGetSize(bNis,&bN);
1724:     ISGetIndices(bNis,&bNindices);
1725:     for (i=0; i<nest->nr; ++i) {
1726:       Mat            B;
1727:       PetscInt       bm, br;
1728:       const PetscInt *bmindices;
1729:       B = nest->m[i][j];
1730:       if (!B) continue;
1731:       ISGetLocalSize(nest->isglobal.row[i],&bm);
1732:       ISGetIndices(nest->isglobal.row[i],&bmindices);
1733:       MatGetOwnershipRange(B,&rstart,NULL);
1734:       for (br = 0; br < bm; ++br) {
1735:         PetscInt          row = bmindices[br], brncols,  *cols;
1736:         const PetscInt    *brcols;
1737:         const PetscScalar *brcoldata;
1738:         MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1739:         PetscMalloc1(brncols,&cols);
1740:         for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
1741:         /*
1742:           Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
1743:           Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
1744:          */
1745:         MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);
1746:         MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1747:         PetscFree(cols);
1748:       }
1749:       ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1750:     }
1751:     ISRestoreIndices(bNis,&bNindices);
1752:     ISDestroy(&bNis);
1753:   }
1754:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1755:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1756:   return(0);
1757: }

1759: PetscErrorCode  MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool  *has)
1760: {
1762:   *has = PETSC_FALSE;
1763:   if (op == MATOP_MULT_TRANSPOSE) {
1764:     Mat_Nest       *bA = (Mat_Nest*)mat->data;
1765:     PetscInt       i,j,nr = bA->nr,nc = bA->nc;
1767:     PetscBool      flg;

1769:     for (j=0; j<nc; j++) {
1770:       for (i=0; i<nr; i++) {
1771:         if (!bA->m[i][j]) continue;
1772:         MatHasOperation(bA->m[i][j],MATOP_MULT_TRANSPOSE_ADD,&flg);
1773:         if (!flg) return(0);
1774:       }
1775:     }
1776:   }
1777:   if (((void**)mat->ops)[op]) *has =  PETSC_TRUE;
1778:   return(0);
1779: }

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

1784:   Level: intermediate

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

1791:   Each of the submatrices lives on the same MPI communicator as the original nest matrix (though they can have zero
1792:   rows/columns on some processes.) Thus this is not meant for cases where the submatrices live on far fewer processes
1793:   than the nest matrix.

1795: .seealso: MatCreate(), MatType, MatCreateNest()
1796: M*/
1797: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
1798: {
1799:   Mat_Nest       *s;

1803:   PetscNewLog(A,&s);
1804:   A->data = (void*)s;

1806:   s->nr            = -1;
1807:   s->nc            = -1;
1808:   s->m             = NULL;
1809:   s->splitassembly = PETSC_FALSE;

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

1813:   A->ops->mult                  = MatMult_Nest;
1814:   A->ops->multadd               = MatMultAdd_Nest;
1815:   A->ops->multtranspose         = MatMultTranspose_Nest;
1816:   A->ops->multtransposeadd      = MatMultTransposeAdd_Nest;
1817:   A->ops->transpose             = MatTranspose_Nest;
1818:   A->ops->assemblybegin         = MatAssemblyBegin_Nest;
1819:   A->ops->assemblyend           = MatAssemblyEnd_Nest;
1820:   A->ops->zeroentries           = MatZeroEntries_Nest;
1821:   A->ops->copy                  = MatCopy_Nest;
1822:   A->ops->duplicate             = MatDuplicate_Nest;
1823:   A->ops->createsubmatrix       = MatCreateSubMatrix_Nest;
1824:   A->ops->destroy               = MatDestroy_Nest;
1825:   A->ops->view                  = MatView_Nest;
1826:   A->ops->getvecs               = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
1827:   A->ops->getlocalsubmatrix     = MatGetLocalSubMatrix_Nest;
1828:   A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
1829:   A->ops->getdiagonal           = MatGetDiagonal_Nest;
1830:   A->ops->diagonalscale         = MatDiagonalScale_Nest;
1831:   A->ops->scale                 = MatScale_Nest;
1832:   A->ops->shift                 = MatShift_Nest;
1833:   A->ops->diagonalset           = MatDiagonalSet_Nest;
1834:   A->ops->setrandom             = MatSetRandom_Nest;
1835:   A->ops->hasoperation          = MatHasOperation_Nest;

1837:   A->spptr        = 0;
1838:   A->assembled    = PETSC_FALSE;

1840:   /* expose Nest api's */
1841:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",   MatNestGetSubMat_Nest);
1842:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",   MatNestSetSubMat_Nest);
1843:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",  MatNestGetSubMats_Nest);
1844:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",     MatNestGetSize_Nest);
1845:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",      MatNestGetISs_Nest);
1846:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);
1847:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",  MatNestSetVecType_Nest);
1848:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",  MatNestSetSubMats_Nest);
1849:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",MatConvert_Nest_AIJ);
1850:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C", MatConvert_Nest_IS);

1852:   PetscObjectChangeTypeName((PetscObject)A,MATNEST);
1853:   return(0);
1854: }