Actual source code: matnest.c

petsc-3.12.0 2019-09-29
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_mpiaij_C",0);
243:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",0);
244:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",0);
245:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",0);
246:   return(0);
247: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

443: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
444: {
446:   Mat_Nest       *vs = (Mat_Nest*)A->data;
447:   Mat            sub;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

642:     PetscFree(L);
643:   }
644:   return(0);
645: }

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

655:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
656:   if (isascii) {

658:     PetscOptionsGetBool(((PetscObject)A)->options,((PetscObject)A)->prefix,"-mat_view_nest_sub",&viewSub,NULL);
659:     PetscViewerASCIIPrintf(viewer,"Matrix object: \n");
660:     PetscViewerASCIIPushTab(viewer);
661:     PetscViewerASCIIPrintf(viewer, "type=nest, rows=%D, cols=%D \n",bA->nr,bA->nc);

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

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

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

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

695: static PetscErrorCode MatZeroEntries_Nest(Mat A)
696: {
697:   Mat_Nest       *bA = (Mat_Nest*)A->data;
698:   PetscInt       i,j;

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

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

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

730: static PetscErrorCode MatAXPY_Nest(Mat Y,PetscScalar a,Mat X,MatStructure str)
731: {
732:   Mat_Nest       *bY = (Mat_Nest*)Y->data,*bX = (Mat_Nest*)X->data;
733:   PetscInt       i,j,nr = bY->nr,nc = bY->nc;

737:   if (nr != bX->nr || nc != bX->nc) SETERRQ4(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Cannot AXPY a MatNest of block size (%D,%D) with a MatNest of block size (%D,%D)",bX->nr,bX->nc,nr,nc);
738:   for (i=0; i<nr; i++) {
739:     for (j=0; j<nc; j++) {
740:       if (bY->m[i][j] && bX->m[i][j]) {
741:         MatAXPY(bY->m[i][j],a,bX->m[i][j],str);
742:       } else if (bX->m[i][j]) {
743:         Mat M;

745:         if (str != DIFFERENT_NONZERO_PATTERN) SETERRQ2(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Matrix block does not exist at %D,%D",i,j);
746:         MatDuplicate(bX->m[i][j],MAT_COPY_VALUES,&M);
747:         MatNestSetSubMat(Y,i,j,M);
748:         MatDestroy(&M);
749:       }
750:     }
751:   }
752:   return(0);
753: }

755: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
756: {
757:   Mat_Nest       *bA = (Mat_Nest*)A->data;
758:   Mat            *b;
759:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

763:   PetscMalloc1(nr*nc,&b);
764:   for (i=0; i<nr; i++) {
765:     for (j=0; j<nc; j++) {
766:       if (bA->m[i][j]) {
767:         MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
768:       } else {
769:         b[i*nc+j] = NULL;
770:       }
771:     }
772:   }
773:   MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
774:   /* Give the new MatNest exclusive ownership */
775:   for (i=0; i<nr*nc; i++) {
776:     MatDestroy(&b[i]);
777:   }
778:   PetscFree(b);

780:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
781:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
782:   return(0);
783: }

785: /* nest api */
786: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
787: {
788:   Mat_Nest *bA = (Mat_Nest*)A->data;

791:   if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
792:   if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
793:   *mat = bA->m[idxm][jdxm];
794:   return(0);
795: }

797: /*@
798:  MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.

800:  Not collective

802:  Input Parameters:
803: +   A  - nest matrix
804: .   idxm - index of the matrix within the nest matrix
805: -   jdxm - index of the matrix within the nest matrix

807:  Output Parameter:
808: .   sub - matrix at index idxm,jdxm within the nest matrix

810:  Level: developer

812: .seealso: MatNestGetSize(), MatNestGetSubMats()
813: @*/
814: PetscErrorCode  MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
815: {

819:   PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
820:   return(0);
821: }

823: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
824: {
825:   Mat_Nest       *bA = (Mat_Nest*)A->data;
826:   PetscInt       m,n,M,N,mi,ni,Mi,Ni;

830:   if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
831:   if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
832:   MatGetLocalSize(mat,&m,&n);
833:   MatGetSize(mat,&M,&N);
834:   ISGetLocalSize(bA->isglobal.row[idxm],&mi);
835:   ISGetSize(bA->isglobal.row[idxm],&Mi);
836:   ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
837:   ISGetSize(bA->isglobal.col[jdxm],&Ni);
838:   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);
839:   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);

841:   PetscObjectReference((PetscObject)mat);
842:   MatDestroy(&bA->m[idxm][jdxm]);
843:   bA->m[idxm][jdxm] = mat;
844:   return(0);
845: }

847: /*@
848:  MatNestSetSubMat - Set a single submatrix in the nest matrix.

850:  Logically collective on the submatrix communicator

852:  Input Parameters:
853: +   A  - nest matrix
854: .   idxm - index of the matrix within the nest matrix
855: .   jdxm - index of the matrix within the nest matrix
856: -   sub - matrix at index idxm,jdxm within the nest matrix

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

861:  This increments the reference count of the submatrix.

863:  Level: developer

865: .seealso: MatNestSetSubMats(), MatNestGetSubMats()
866: @*/
867: PetscErrorCode  MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
868: {

872:   PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
873:   return(0);
874: }

876: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
877: {
878:   Mat_Nest *bA = (Mat_Nest*)A->data;

881:   if (M)   *M   = bA->nr;
882:   if (N)   *N   = bA->nc;
883:   if (mat) *mat = bA->m;
884:   return(0);
885: }

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

890:  Not collective

892:  Input Parameters:
893: .   A  - nest matrix

895:  Output Parameter:
896: +   M - number of rows in the nest matrix
897: .   N - number of cols in the nest matrix
898: -   mat - 2d array of matrices

900:  Notes:

902:  The user should not free the array mat.

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

908:  Level: developer

910: .seealso: MatNestGetSize(), MatNestGetSubMat()
911: @*/
912: PetscErrorCode  MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
913: {

917:   PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
918:   return(0);
919: }

921: PetscErrorCode  MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
922: {
923:   Mat_Nest *bA = (Mat_Nest*)A->data;

926:   if (M) *M = bA->nr;
927:   if (N) *N = bA->nc;
928:   return(0);
929: }

931: /*@
932:  MatNestGetSize - Returns the size of the nest matrix.

934:  Not collective

936:  Input Parameters:
937: .   A  - nest matrix

939:  Output Parameter:
940: +   M - number of rows in the nested mat
941: -   N - number of cols in the nested mat

943:  Notes:

945:  Level: developer

947: .seealso: MatNestGetSubMat(), MatNestGetSubMats()
948: @*/
949: PetscErrorCode  MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
950: {

954:   PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
955:   return(0);
956: }

958: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
959: {
960:   Mat_Nest *vs = (Mat_Nest*)A->data;
961:   PetscInt i;

964:   if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
965:   if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
966:   return(0);
967: }

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

972:  Not collective

974:  Input Parameters:
975: .   A  - nest matrix

977:  Output Parameter:
978: +   rows - array of row index sets
979: -   cols - array of column index sets

981:  Level: advanced

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

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

994:   PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
995:   return(0);
996: }

998: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
999: {
1000:   Mat_Nest *vs = (Mat_Nest*)A->data;
1001:   PetscInt i;

1004:   if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
1005:   if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
1006:   return(0);
1007: }

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

1012:  Not collective

1014:  Input Parameters:
1015: .   A  - nest matrix

1017:  Output Parameter:
1018: +   rows - array of row index sets (or NULL to ignore)
1019: -   cols - array of column index sets (or NULL to ignore)

1021:  Level: advanced

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

1026: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs()
1027: @*/
1028: PetscErrorCode  MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
1029: {

1034:   PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1035:   return(0);
1036: }

1038: PetscErrorCode  MatNestSetVecType_Nest(Mat A,VecType vtype)
1039: {
1041:   PetscBool      flg;

1044:   PetscStrcmp(vtype,VECNEST,&flg);
1045:   /* In reality, this only distinguishes VECNEST and "other" */
1046:   if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1047:   else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
1048:   return(0);
1049: }

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

1054:  Not collective

1056:  Input Parameters:
1057: +  A  - nest matrix
1058: -  vtype - type to use for creating vectors

1060:  Notes:

1062:  Level: developer

1064: .seealso: MatCreateVecs()
1065: @*/
1066: PetscErrorCode  MatNestSetVecType(Mat A,VecType vtype)
1067: {

1071:   PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
1072:   return(0);
1073: }

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

1082:   s->nr = nr;
1083:   s->nc = nc;

1085:   /* Create space for submatrices */
1086:   PetscMalloc1(nr,&s->m);
1087:   for (i=0; i<nr; i++) {
1088:     PetscMalloc1(nc,&s->m[i]);
1089:   }
1090:   for (i=0; i<nr; i++) {
1091:     for (j=0; j<nc; j++) {
1092:       s->m[i][j] = a[i*nc+j];
1093:       if (a[i*nc+j]) {
1094:         PetscObjectReference((PetscObject)a[i*nc+j]);
1095:       }
1096:     }
1097:   }

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

1101:   PetscMalloc1(nr,&s->row_len);
1102:   PetscMalloc1(nc,&s->col_len);
1103:   for (i=0; i<nr; i++) s->row_len[i]=-1;
1104:   for (j=0; j<nc; j++) s->col_len[j]=-1;

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

1108:   PetscLayoutSetSize(A->rmap,M);
1109:   PetscLayoutSetLocalSize(A->rmap,m);
1110:   PetscLayoutSetSize(A->cmap,N);
1111:   PetscLayoutSetLocalSize(A->cmap,n);

1113:   PetscLayoutSetUp(A->rmap);
1114:   PetscLayoutSetUp(A->cmap);

1116:   PetscCalloc2(nr,&s->left,nc,&s->right);
1117:   return(0);
1118: }

1120: /*@
1121:    MatNestSetSubMats - Sets the nested submatrices

1123:    Collective on Mat

1125:    Input Parameter:
1126: +  A - nested matrix
1127: .  nr - number of nested row blocks
1128: .  is_row - index sets for each nested row block, or NULL to make contiguous
1129: .  nc - number of nested column blocks
1130: .  is_col - index sets for each nested column block, or NULL to make contiguous
1131: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

1133:    Level: advanced

1135: .seealso: MatCreateNest(), MATNEST
1136: @*/
1137: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1138: {
1140:   PetscInt       i,nr_nc;

1144:   if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1145:   if (nr && is_row) {
1148:   }
1149:   if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1150:   if (nc && is_col) {
1153:   }
1154:   nr_nc=nr*nc;
1156:   PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1157:   return(0);
1158: }

1160: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1161: {
1163:   PetscBool      flg;
1164:   PetscInt       i,j,m,mi,*ix;

1167:   for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1168:     if (islocal[i]) {
1169:       ISGetSize(islocal[i],&mi);
1170:       flg  = PETSC_TRUE;      /* We found a non-trivial entry */
1171:     } else {
1172:       ISGetSize(isglobal[i],&mi);
1173:     }
1174:     m += mi;
1175:   }
1176:   if (flg) {
1177:     PetscMalloc1(m,&ix);
1178:     for (i=0,m=0; i<n; i++) {
1179:       ISLocalToGlobalMapping smap = NULL;
1180:       Mat                    sub = NULL;
1181:       PetscSF                sf;
1182:       PetscLayout            map;
1183:       PetscInt               *ix2;

1185:       if (!colflg) {
1186:         MatNestFindNonzeroSubMatRow(A,i,&sub);
1187:       } else {
1188:         MatNestFindNonzeroSubMatCol(A,i,&sub);
1189:       }
1190:       if (sub) {
1191:         if (!colflg) {
1192:           MatGetLocalToGlobalMapping(sub,&smap,NULL);
1193:         } else {
1194:           MatGetLocalToGlobalMapping(sub,NULL,&smap);
1195:         }
1196:       }
1197:       if (islocal[i]) {
1198:         ISGetSize(islocal[i],&mi);
1199:       } else {
1200:         ISGetSize(isglobal[i],&mi);
1201:       }
1202:       for (j=0; j<mi; j++) ix[m+j] = j;
1203:       if (smap) {ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);}

1205:       /*
1206:         Now we need to extract the monolithic global indices that correspond to the given split global indices.
1207:         In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1208:        */
1209:       PetscMalloc1(mi,&ix2);
1210:       PetscSFCreate(((PetscObject)isglobal[i])->comm,&sf);
1211:       PetscLayoutCreate(((PetscObject)isglobal[i])->comm,&map);
1212:       PetscLayoutSetLocalSize(map,mi);
1213:       PetscLayoutSetUp(map);
1214:       PetscSFSetGraphLayout(sf,map,mi,NULL,PETSC_USE_POINTER,ix+m);
1215:       PetscLayoutDestroy(&map);
1216:       for (j=0; j<mi; j++) ix2[j] = ix[m+j];
1217:       PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m);
1218:       PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m);
1219:       PetscSFDestroy(&sf);
1220:       PetscFree(ix2);
1221:       m   += mi;
1222:     }
1223:     ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);
1224:   } else {
1225:     *ltog = NULL;
1226:   }
1227:   return(0);
1228: }


1231: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1232: /*
1233:   nprocessors = NP
1234:   Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1235:        proc 0: => (g_0,h_0,)
1236:        proc 1: => (g_1,h_1,)
1237:        ...
1238:        proc nprocs-1: => (g_NP-1,h_NP-1,)

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

1243:             proc 0:
1244:     is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1245:             proc 1:
1246:     is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)

1248:             proc NP-1:
1249:     is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1250: */
1251: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1252: {
1253:   Mat_Nest       *vs = (Mat_Nest*)A->data;
1254:   PetscInt       i,j,offset,n,nsum,bs;
1256:   Mat            sub = NULL;

1259:   PetscMalloc1(nr,&vs->isglobal.row);
1260:   PetscMalloc1(nc,&vs->isglobal.col);
1261:   if (is_row) { /* valid IS is passed in */
1262:     /* refs on is[] are incremeneted */
1263:     for (i=0; i<vs->nr; i++) {
1264:       PetscObjectReference((PetscObject)is_row[i]);

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

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

1294:       vs->isglobal.col[j] = is_col[j];
1295:     }
1296:   } else {                      /* Create the ISs by inspecting sizes of a submatrix in each column */
1297:     offset = A->cmap->rstart;
1298:     nsum   = 0;
1299:     for (j=0; j<vs->nc; j++) {
1300:       MatNestFindNonzeroSubMatCol(A,j,&sub);
1301:       if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i);
1302:       MatGetLocalSize(sub,NULL,&n);
1303:       if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1304:       nsum += n;
1305:     }
1306:     MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1307:     offset -= nsum;
1308:     for (j=0; j<vs->nc; j++) {
1309:       MatNestFindNonzeroSubMatCol(A,j,&sub);
1310:       MatGetLocalSize(sub,NULL,&n);
1311:       MatGetBlockSize(sub,&bs);
1312:       ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);
1313:       ISSetBlockSize(vs->isglobal.col[j],bs);
1314:       offset += n;
1315:     }
1316:   }

1318:   /* Set up the local ISs */
1319:   PetscMalloc1(vs->nr,&vs->islocal.row);
1320:   PetscMalloc1(vs->nc,&vs->islocal.col);
1321:   for (i=0,offset=0; i<vs->nr; i++) {
1322:     IS                     isloc;
1323:     ISLocalToGlobalMapping rmap = NULL;
1324:     PetscInt               nlocal,bs;
1325:     MatNestFindNonzeroSubMatRow(A,i,&sub);
1326:     if (sub) {MatGetLocalToGlobalMapping(sub,&rmap,NULL);}
1327:     if (rmap) {
1328:       MatGetBlockSize(sub,&bs);
1329:       ISLocalToGlobalMappingGetSize(rmap,&nlocal);
1330:       ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1331:       ISSetBlockSize(isloc,bs);
1332:     } else {
1333:       nlocal = 0;
1334:       isloc  = NULL;
1335:     }
1336:     vs->islocal.row[i] = isloc;
1337:     offset            += nlocal;
1338:   }
1339:   for (i=0,offset=0; i<vs->nc; i++) {
1340:     IS                     isloc;
1341:     ISLocalToGlobalMapping cmap = NULL;
1342:     PetscInt               nlocal,bs;
1343:     MatNestFindNonzeroSubMatCol(A,i,&sub);
1344:     if (sub) {MatGetLocalToGlobalMapping(sub,NULL,&cmap);}
1345:     if (cmap) {
1346:       MatGetBlockSize(sub,&bs);
1347:       ISLocalToGlobalMappingGetSize(cmap,&nlocal);
1348:       ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1349:       ISSetBlockSize(isloc,bs);
1350:     } else {
1351:       nlocal = 0;
1352:       isloc  = NULL;
1353:     }
1354:     vs->islocal.col[i] = isloc;
1355:     offset            += nlocal;
1356:   }

1358:   /* Set up the aggregate ISLocalToGlobalMapping */
1359:   {
1360:     ISLocalToGlobalMapping rmap,cmap;
1361:     MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);
1362:     MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);
1363:     if (rmap && cmap) {MatSetLocalToGlobalMapping(A,rmap,cmap);}
1364:     ISLocalToGlobalMappingDestroy(&rmap);
1365:     ISLocalToGlobalMappingDestroy(&cmap);
1366:   }

1368: #if defined(PETSC_USE_DEBUG)
1369:   for (i=0; i<vs->nr; i++) {
1370:     for (j=0; j<vs->nc; j++) {
1371:       PetscInt m,n,M,N,mi,ni,Mi,Ni;
1372:       Mat      B = vs->m[i][j];
1373:       if (!B) continue;
1374:       MatGetSize(B,&M,&N);
1375:       MatGetLocalSize(B,&m,&n);
1376:       ISGetSize(vs->isglobal.row[i],&Mi);
1377:       ISGetSize(vs->isglobal.col[j],&Ni);
1378:       ISGetLocalSize(vs->isglobal.row[i],&mi);
1379:       ISGetLocalSize(vs->isglobal.col[j],&ni);
1380:       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);
1381:       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);
1382:     }
1383:   }
1384: #endif

1386:   /* Set A->assembled if all non-null blocks are currently assembled */
1387:   for (i=0; i<vs->nr; i++) {
1388:     for (j=0; j<vs->nc; j++) {
1389:       if (vs->m[i][j] && !vs->m[i][j]->assembled) return(0);
1390:     }
1391:   }
1392:   A->assembled = PETSC_TRUE;
1393:   return(0);
1394: }

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

1399:    Collective on Mat

1401:    Input Parameter:
1402: +  comm - Communicator for the new Mat
1403: .  nr - number of nested row blocks
1404: .  is_row - index sets for each nested row block, or NULL to make contiguous
1405: .  nc - number of nested column blocks
1406: .  is_col - index sets for each nested column block, or NULL to make contiguous
1407: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

1409:    Output Parameter:
1410: .  B - new matrix

1412:    Level: advanced

1414: .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST
1415: @*/
1416: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1417: {
1418:   Mat            A;

1422:   *B   = 0;
1423:   MatCreate(comm,&A);
1424:   MatSetType(A,MATNEST);
1425:   A->preallocated = PETSC_TRUE;
1426:   MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1427:   *B   = A;
1428:   return(0);
1429: }

1431: static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1432: {
1433:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1434:   Mat            *trans;
1435:   PetscScalar    **avv;
1436:   PetscScalar    *vv;
1437:   PetscInt       **aii,**ajj;
1438:   PetscInt       *ii,*jj,*ci;
1439:   PetscInt       nr,nc,nnz,i,j;
1440:   PetscBool      done;

1444:   MatGetSize(A,&nr,&nc);
1445:   if (reuse == MAT_REUSE_MATRIX) {
1446:     PetscInt rnr;

1448:     MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);
1449:     if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ");
1450:     if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows");
1451:     MatSeqAIJGetArray(*newmat,&vv);
1452:   }
1453:   /* extract CSR for nested SeqAIJ matrices */
1454:   nnz  = 0;
1455:   PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);
1456:   for (i=0; i<nest->nr; ++i) {
1457:     for (j=0; j<nest->nc; ++j) {
1458:       Mat B = nest->m[i][j];
1459:       if (B) {
1460:         PetscScalar *naa;
1461:         PetscInt    *nii,*njj,nnr;
1462:         PetscBool   istrans;

1464:         PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1465:         if (istrans) {
1466:           Mat Bt;

1468:           MatTransposeGetMat(B,&Bt);
1469:           MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);
1470:           B    = trans[i*nest->nc+j];
1471:         }
1472:         MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);
1473:         if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ");
1474:         MatSeqAIJGetArray(B,&naa);
1475:         nnz += nii[nnr];

1477:         aii[i*nest->nc+j] = nii;
1478:         ajj[i*nest->nc+j] = njj;
1479:         avv[i*nest->nc+j] = naa;
1480:       }
1481:     }
1482:   }
1483:   if (reuse != MAT_REUSE_MATRIX) {
1484:     PetscMalloc1(nr+1,&ii);
1485:     PetscMalloc1(nnz,&jj);
1486:     PetscMalloc1(nnz,&vv);
1487:   } else {
1488:     if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros");
1489:   }

1491:   /* new row pointer */
1492:   PetscArrayzero(ii,nr+1);
1493:   for (i=0; i<nest->nr; ++i) {
1494:     PetscInt       ncr,rst;

1496:     ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1497:     ISGetLocalSize(nest->isglobal.row[i],&ncr);
1498:     for (j=0; j<nest->nc; ++j) {
1499:       if (aii[i*nest->nc+j]) {
1500:         PetscInt    *nii = aii[i*nest->nc+j];
1501:         PetscInt    ir;

1503:         for (ir=rst; ir<ncr+rst; ++ir) {
1504:           ii[ir+1] += nii[1]-nii[0];
1505:           nii++;
1506:         }
1507:       }
1508:     }
1509:   }
1510:   for (i=0; i<nr; i++) ii[i+1] += ii[i];

1512:   /* construct CSR for the new matrix */
1513:   PetscCalloc1(nr,&ci);
1514:   for (i=0; i<nest->nr; ++i) {
1515:     PetscInt       ncr,rst;

1517:     ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1518:     ISGetLocalSize(nest->isglobal.row[i],&ncr);
1519:     for (j=0; j<nest->nc; ++j) {
1520:       if (aii[i*nest->nc+j]) {
1521:         PetscScalar *nvv = avv[i*nest->nc+j];
1522:         PetscInt    *nii = aii[i*nest->nc+j];
1523:         PetscInt    *njj = ajj[i*nest->nc+j];
1524:         PetscInt    ir,cst;

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

1530:           for (ij=0;ij<rsize;ij++) {
1531:             jj[ist+ij] = *njj+cst;
1532:             vv[ist+ij] = *nvv;
1533:             njj++;
1534:             nvv++;
1535:           }
1536:           ci[ir] += rsize;
1537:           nii++;
1538:         }
1539:       }
1540:     }
1541:   }
1542:   PetscFree(ci);

1544:   /* restore info */
1545:   for (i=0; i<nest->nr; ++i) {
1546:     for (j=0; j<nest->nc; ++j) {
1547:       Mat B = nest->m[i][j];
1548:       if (B) {
1549:         PetscInt nnr = 0, k = i*nest->nc+j;

1551:         B    = (trans[k] ? trans[k] : B);
1552:         MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);
1553:         if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ");
1554:         MatSeqAIJRestoreArray(B,&avv[k]);
1555:         MatDestroy(&trans[k]);
1556:       }
1557:     }
1558:   }
1559:   PetscFree4(aii,ajj,avv,trans);

1561:   /* finalize newmat */
1562:   if (reuse == MAT_INITIAL_MATRIX) {
1563:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);
1564:   } else if (reuse == MAT_INPLACE_MATRIX) {
1565:     Mat B;

1567:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);
1568:     MatHeaderReplace(A,&B);
1569:   }
1570:   MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1571:   MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1572:   {
1573:     Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data);
1574:     a->free_a     = PETSC_TRUE;
1575:     a->free_ij    = PETSC_TRUE;
1576:   }
1577:   return(0);
1578: }

1580: PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1581: {
1583:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1584:   PetscInt       m,n,M,N,i,j,k,*dnnz,*onnz,rstart;
1585:   PetscInt       cstart,cend;
1586:   PetscMPIInt    size;
1587:   Mat            C;

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

1595:     PetscStrcmp(newtype,MATAIJ,&fast);
1596:     if (!fast) {
1597:       PetscStrcmp(newtype,MATSEQAIJ,&fast);
1598:     }
1599:     for (i=0; i<nest->nr && fast; ++i) {
1600:       for (j=0; j<nest->nc && fast; ++j) {
1601:         Mat B = nest->m[i][j];
1602:         if (B) {
1603:           PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);
1604:           if (!fast) {
1605:             PetscBool istrans;

1607:             PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1608:             if (istrans) {
1609:               Mat Bt;

1611:               MatTransposeGetMat(B,&Bt);
1612:               PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);
1613:             }
1614:           }
1615:         }
1616:       }
1617:     }
1618:     for (i=0, nf=0; i<nest->nr && fast; ++i) {
1619:       PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);
1620:       if (fast) {
1621:         PetscInt f,s;

1623:         ISStrideGetInfo(nest->isglobal.row[i],&f,&s);
1624:         if (f != nf || s != 1) { fast = PETSC_FALSE; }
1625:         else {
1626:           ISGetSize(nest->isglobal.row[i],&f);
1627:           nf  += f;
1628:         }
1629:       }
1630:     }
1631:     for (i=0, nf=0; i<nest->nc && fast; ++i) {
1632:       PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);
1633:       if (fast) {
1634:         PetscInt f,s;

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

1744:   /* Fill by row */
1745:   for (j=0; j<nest->nc; ++j) {
1746:     /* Using global column indices and ISAllGather() is not scalable. */
1747:     IS             bNis;
1748:     PetscInt       bN;
1749:     const PetscInt *bNindices;
1750:     ISAllGather(nest->isglobal.col[j], &bNis);
1751:     ISGetSize(bNis,&bN);
1752:     ISGetIndices(bNis,&bNindices);
1753:     for (i=0; i<nest->nr; ++i) {
1754:       Mat            B;
1755:       PetscInt       bm, br;
1756:       const PetscInt *bmindices;
1757:       B = nest->m[i][j];
1758:       if (!B) continue;
1759:       ISGetLocalSize(nest->isglobal.row[i],&bm);
1760:       ISGetIndices(nest->isglobal.row[i],&bmindices);
1761:       MatGetOwnershipRange(B,&rstart,NULL);
1762:       for (br = 0; br < bm; ++br) {
1763:         PetscInt          row = bmindices[br], brncols,  *cols;
1764:         const PetscInt    *brcols;
1765:         const PetscScalar *brcoldata;
1766:         MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1767:         PetscMalloc1(brncols,&cols);
1768:         for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
1769:         /*
1770:           Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
1771:           Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
1772:          */
1773:         MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);
1774:         MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1775:         PetscFree(cols);
1776:       }
1777:       ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1778:     }
1779:     ISRestoreIndices(bNis,&bNindices);
1780:     ISDestroy(&bNis);
1781:   }
1782:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1783:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1784:   return(0);
1785: }

1787: PetscErrorCode  MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool  *has)
1788: {
1790:   *has = PETSC_FALSE;
1791:   if (op == MATOP_MULT_TRANSPOSE) {
1792:     Mat_Nest       *bA = (Mat_Nest*)mat->data;
1793:     PetscInt       i,j,nr = bA->nr,nc = bA->nc;
1795:     PetscBool      flg;

1797:     for (j=0; j<nc; j++) {
1798:       for (i=0; i<nr; i++) {
1799:         if (!bA->m[i][j]) continue;
1800:         MatHasOperation(bA->m[i][j],MATOP_MULT_TRANSPOSE_ADD,&flg);
1801:         if (!flg) return(0);
1802:       }
1803:     }
1804:   }
1805:   if (((void**)mat->ops)[op]) *has =  PETSC_TRUE;
1806:   return(0);
1807: }

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

1812:   Level: intermediate

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

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

1823: .seealso: MatCreate(), MatType, MatCreateNest()
1824: M*/
1825: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
1826: {
1827:   Mat_Nest       *s;

1831:   PetscNewLog(A,&s);
1832:   A->data = (void*)s;

1834:   s->nr            = -1;
1835:   s->nc            = -1;
1836:   s->m             = NULL;
1837:   s->splitassembly = PETSC_FALSE;

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

1841:   A->ops->mult                  = MatMult_Nest;
1842:   A->ops->multadd               = MatMultAdd_Nest;
1843:   A->ops->multtranspose         = MatMultTranspose_Nest;
1844:   A->ops->multtransposeadd      = MatMultTransposeAdd_Nest;
1845:   A->ops->transpose             = MatTranspose_Nest;
1846:   A->ops->assemblybegin         = MatAssemblyBegin_Nest;
1847:   A->ops->assemblyend           = MatAssemblyEnd_Nest;
1848:   A->ops->zeroentries           = MatZeroEntries_Nest;
1849:   A->ops->copy                  = MatCopy_Nest;
1850:   A->ops->axpy                  = MatAXPY_Nest;
1851:   A->ops->duplicate             = MatDuplicate_Nest;
1852:   A->ops->createsubmatrix       = MatCreateSubMatrix_Nest;
1853:   A->ops->destroy               = MatDestroy_Nest;
1854:   A->ops->view                  = MatView_Nest;
1855:   A->ops->getvecs               = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
1856:   A->ops->getlocalsubmatrix     = MatGetLocalSubMatrix_Nest;
1857:   A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
1858:   A->ops->getdiagonal           = MatGetDiagonal_Nest;
1859:   A->ops->diagonalscale         = MatDiagonalScale_Nest;
1860:   A->ops->scale                 = MatScale_Nest;
1861:   A->ops->shift                 = MatShift_Nest;
1862:   A->ops->diagonalset           = MatDiagonalSet_Nest;
1863:   A->ops->setrandom             = MatSetRandom_Nest;
1864:   A->ops->hasoperation          = MatHasOperation_Nest;

1866:   A->spptr        = 0;
1867:   A->assembled    = PETSC_FALSE;

1869:   /* expose Nest api's */
1870:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",      MatNestGetSubMat_Nest);
1871:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",      MatNestSetSubMat_Nest);
1872:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",     MatNestGetSubMats_Nest);
1873:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",        MatNestGetSize_Nest);
1874:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",         MatNestGetISs_Nest);
1875:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",    MatNestGetLocalISs_Nest);
1876:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",     MatNestSetVecType_Nest);
1877:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",     MatNestSetSubMats_Nest);
1878:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",MatConvert_Nest_AIJ);
1879:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",MatConvert_Nest_AIJ);
1880:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",   MatConvert_Nest_AIJ);
1881:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",    MatConvert_Nest_IS);

1883:   PetscObjectChangeTypeName((PetscObject)A,MATNEST);
1884:   return(0);
1885: }