Actual source code: matnest.c

petsc-3.12.4 2020-02-04
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,Vec*,Vec*);
  8: static PetscErrorCode MatReset_Nest(Mat);

 10: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat,MatType,MatReuse,Mat*);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

192: static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list)
193: {
195:   IS             *lst = *list;
196:   PetscInt       i;

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

206: static PetscErrorCode MatReset_Nest(Mat A)
207: {
208:   Mat_Nest       *vs = (Mat_Nest*)A->data;
209:   PetscInt       i,j;

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

219:   PetscFree(vs->row_len);
220:   PetscFree(vs->col_len);
221:   PetscFree(vs->nnzstate);

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

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

236:   /* restore defaults */
237:   vs->nr = 0;
238:   vs->nc = 0;
239:   vs->splitassembly = PETSC_FALSE;
240:   return(0);
241: }

243: static PetscErrorCode MatDestroy_Nest(Mat A)
244: {

247:   MatReset_Nest(A);
248:   PetscFree(A->data);

250:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",0);
251:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",0);
252:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",0);
253:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",0);
254:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",0);
255:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",0);
256:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",0);
257:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",0);
258:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",0);
259:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",0);
260:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",0);
261:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",0);
262:   return(0);
263: }

265: static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)
266: {
267:   Mat_Nest       *vs = (Mat_Nest*)A->data;
268:   PetscInt       i,j;
270:   PetscBool      nnzstate = PETSC_FALSE;

273:   for (i=0; i<vs->nr; i++) {
274:     for (j=0; j<vs->nc; j++) {
275:       PetscObjectState subnnzstate = 0;
276:       if (vs->m[i][j]) {
277:         MatAssemblyBegin(vs->m[i][j],type);
278:         if (!vs->splitassembly) {
279:           /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
280:            * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
281:            * already performing an assembly, but the result would by more complicated and appears to offer less
282:            * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
283:            * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
284:            */
285:           MatAssemblyEnd(vs->m[i][j],type);
286:           MatGetNonzeroState(vs->m[i][j],&subnnzstate);
287:         }
288:       }
289:       nnzstate = (PetscBool)(nnzstate || vs->nnzstate[i*vs->nc+j] != subnnzstate);
290:       vs->nnzstate[i*vs->nc+j] = subnnzstate;
291:     }
292:   }
293:   if (nnzstate) A->nonzerostate++;
294:   return(0);
295: }

297: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
298: {
299:   Mat_Nest       *vs = (Mat_Nest*)A->data;
300:   PetscInt       i,j;

304:   for (i=0; i<vs->nr; i++) {
305:     for (j=0; j<vs->nc; j++) {
306:       if (vs->m[i][j]) {
307:         if (vs->splitassembly) {
308:           MatAssemblyEnd(vs->m[i][j],type);
309:         }
310:       }
311:     }
312:   }
313:   return(0);
314: }

316: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B)
317: {
319:   Mat_Nest       *vs = (Mat_Nest*)A->data;
320:   PetscInt       j;
321:   Mat            sub;

324:   sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
325:   for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j];
326:   if (sub) {MatSetUp(sub);}       /* Ensure that the sizes are available */
327:   *B = sub;
328:   return(0);
329: }

331: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B)
332: {
334:   Mat_Nest       *vs = (Mat_Nest*)A->data;
335:   PetscInt       i;
336:   Mat            sub;

339:   sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
340:   for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col];
341:   if (sub) {MatSetUp(sub);}       /* Ensure that the sizes are available */
342:   *B = sub;
343:   return(0);
344: }

346: static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found)
347: {
349:   PetscInt       i;
350:   PetscBool      flg;

356:   *found = -1;
357:   for (i=0; i<n; i++) {
358:     if (!list[i]) continue;
359:     ISEqualUnsorted(list[i],is,&flg);
360:     if (flg) {
361:       *found = i;
362:       return(0);
363:     }
364:   }
365:   SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set");
366:   return(0);
367: }

369: /* Get a block row as a new MatNest */
370: static PetscErrorCode MatNestGetRow(Mat A,PetscInt row,Mat *B)
371: {
372:   Mat_Nest       *vs = (Mat_Nest*)A->data;
373:   char           keyname[256];

377:   *B   = NULL;
378:   PetscSNPrintf(keyname,sizeof(keyname),"NestRow_%D",row);
379:   PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);
380:   if (*B) return(0);

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

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

386:   PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);
387:   PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
388:   return(0);
389: }

391: static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
392: {
393:   Mat_Nest       *vs = (Mat_Nest*)A->data;
395:   PetscInt       row,col;
396:   PetscBool      same,isFullCol,isFullColGlobal;

399:   /* Check if full column space. This is a hack */
400:   isFullCol = PETSC_FALSE;
401:   PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&same);
402:   if (same) {
403:     PetscInt n,first,step,i,an,am,afirst,astep;
404:     ISStrideGetInfo(iscol,&first,&step);
405:     ISGetLocalSize(iscol,&n);
406:     isFullCol = PETSC_TRUE;
407:     for (i=0,an=A->cmap->rstart; i<vs->nc; i++) {
408:       PetscObjectTypeCompare((PetscObject)is->col[i],ISSTRIDE,&same);
409:       ISGetLocalSize(is->col[i],&am);
410:       if (same) {
411:         ISStrideGetInfo(is->col[i],&afirst,&astep);
412:         if (afirst != an || astep != step) isFullCol = PETSC_FALSE;
413:       } else isFullCol = PETSC_FALSE;
414:       an += am;
415:     }
416:     if (an != A->cmap->rstart+n) isFullCol = PETSC_FALSE;
417:   }
418:   MPIU_Allreduce(&isFullCol,&isFullColGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)iscol));

420:   if (isFullColGlobal && vs->nc > 1) {
421:     PetscInt row;
422:     MatNestFindIS(A,vs->nr,is->row,isrow,&row);
423:     MatNestGetRow(A,row,B);
424:   } else {
425:     MatNestFindIS(A,vs->nr,is->row,isrow,&row);
426:     MatNestFindIS(A,vs->nc,is->col,iscol,&col);
427:     if (!vs->m[row][col]) {
428:       PetscInt lr,lc;

430:       MatCreate(PetscObjectComm((PetscObject)A),&vs->m[row][col]);
431:       ISGetLocalSize(vs->isglobal.row[row],&lr);
432:       ISGetLocalSize(vs->isglobal.col[col],&lc);
433:       MatSetSizes(vs->m[row][col],lr,lc,PETSC_DECIDE,PETSC_DECIDE);
434:       MatSetType(vs->m[row][col],MATAIJ);
435:       MatSeqAIJSetPreallocation(vs->m[row][col],0,NULL);
436:       MatMPIAIJSetPreallocation(vs->m[row][col],0,NULL,0,NULL);
437:       MatSetUp(vs->m[row][col]);
438:       MatAssemblyBegin(vs->m[row][col],MAT_FINAL_ASSEMBLY);
439:       MatAssemblyEnd(vs->m[row][col],MAT_FINAL_ASSEMBLY);
440:     }
441:     *B = vs->m[row][col];
442:   }
443:   return(0);
444: }

446: /*
447:    TODO: This does not actually returns a submatrix we can modify
448: */
449: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
450: {
452:   Mat_Nest       *vs = (Mat_Nest*)A->data;
453:   Mat            sub;

456:   MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
457:   switch (reuse) {
458:   case MAT_INITIAL_MATRIX:
459:     if (sub) { PetscObjectReference((PetscObject)sub); }
460:     *B = sub;
461:     break;
462:   case MAT_REUSE_MATRIX:
463:     if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call");
464:     break;
465:   case MAT_IGNORE_MATRIX:       /* Nothing to do */
466:     break;
467:   case MAT_INPLACE_MATRIX:       /* Nothing to do */
468:     SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet");
469:     break;
470:   }
471:   return(0);
472: }

474: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
475: {
477:   Mat_Nest       *vs = (Mat_Nest*)A->data;
478:   Mat            sub;

481:   MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
482:   /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
483:   if (sub) {PetscObjectReference((PetscObject)sub);}
484:   *B = sub;
485:   return(0);
486: }

488: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
489: {
491:   Mat_Nest       *vs = (Mat_Nest*)A->data;
492:   Mat            sub;

495:   MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
496:   if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten");
497:   if (sub) {
498:     if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times");
499:     MatDestroy(B);
500:   }
501:   return(0);
502: }

504: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
505: {
506:   Mat_Nest       *bA = (Mat_Nest*)A->data;
507:   PetscInt       i;

511:   for (i=0; i<bA->nr; i++) {
512:     Vec bv;
513:     VecGetSubVector(v,bA->isglobal.row[i],&bv);
514:     if (bA->m[i][i]) {
515:       MatGetDiagonal(bA->m[i][i],bv);
516:     } else {
517:       VecSet(bv,0.0);
518:     }
519:     VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
520:   }
521:   return(0);
522: }

524: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
525: {
526:   Mat_Nest       *bA = (Mat_Nest*)A->data;
527:   Vec            bl,*br;
528:   PetscInt       i,j;

532:   PetscCalloc1(bA->nc,&br);
533:   if (r) {
534:     for (j=0; j<bA->nc; j++) {VecGetSubVector(r,bA->isglobal.col[j],&br[j]);}
535:   }
536:   bl = NULL;
537:   for (i=0; i<bA->nr; i++) {
538:     if (l) {
539:       VecGetSubVector(l,bA->isglobal.row[i],&bl);
540:     }
541:     for (j=0; j<bA->nc; j++) {
542:       if (bA->m[i][j]) {
543:         MatDiagonalScale(bA->m[i][j],bl,br[j]);
544:       }
545:     }
546:     if (l) {
547:       VecRestoreSubVector(l,bA->isglobal.row[i],&bl);
548:     }
549:   }
550:   if (r) {
551:     for (j=0; j<bA->nc; j++) {VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);}
552:   }
553:   PetscFree(br);
554:   return(0);
555: }

557: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
558: {
559:   Mat_Nest       *bA = (Mat_Nest*)A->data;
560:   PetscInt       i,j;

564:   for (i=0; i<bA->nr; i++) {
565:     for (j=0; j<bA->nc; j++) {
566:       if (bA->m[i][j]) {
567:         MatScale(bA->m[i][j],a);
568:       }
569:     }
570:   }
571:   return(0);
572: }

574: static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
575: {
576:   Mat_Nest       *bA = (Mat_Nest*)A->data;
577:   PetscInt       i;
579:   PetscBool      nnzstate = PETSC_FALSE;

582:   for (i=0; i<bA->nr; i++) {
583:     PetscObjectState subnnzstate = 0;
584:     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);
585:     MatShift(bA->m[i][i],a);
586:     MatGetNonzeroState(bA->m[i][i],&subnnzstate);
587:     nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
588:     bA->nnzstate[i*bA->nc+i] = subnnzstate;
589:   }
590:   if (nnzstate) A->nonzerostate++;
591:   return(0);
592: }

594: static PetscErrorCode MatDiagonalSet_Nest(Mat A,Vec D,InsertMode is)
595: {
596:   Mat_Nest       *bA = (Mat_Nest*)A->data;
597:   PetscInt       i;
599:   PetscBool      nnzstate = PETSC_FALSE;

602:   for (i=0; i<bA->nr; i++) {
603:     PetscObjectState subnnzstate = 0;
604:     Vec              bv;
605:     VecGetSubVector(D,bA->isglobal.row[i],&bv);
606:     if (bA->m[i][i]) {
607:       MatDiagonalSet(bA->m[i][i],bv,is);
608:       MatGetNonzeroState(bA->m[i][i],&subnnzstate);
609:     }
610:     VecRestoreSubVector(D,bA->isglobal.row[i],&bv);
611:     nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
612:     bA->nnzstate[i*bA->nc+i] = subnnzstate;
613:   }
614:   if (nnzstate) A->nonzerostate++;
615:   return(0);
616: }

618: static PetscErrorCode MatSetRandom_Nest(Mat A,PetscRandom rctx)
619: {
620:   Mat_Nest       *bA = (Mat_Nest*)A->data;
621:   PetscInt       i,j;

625:   for (i=0; i<bA->nr; i++) {
626:     for (j=0; j<bA->nc; j++) {
627:       if (bA->m[i][j]) {
628:         MatSetRandom(bA->m[i][j],rctx);
629:       }
630:     }
631:   }
632:   return(0);
633: }

635: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left)
636: {
637:   Mat_Nest       *bA = (Mat_Nest*)A->data;
638:   Vec            *L,*R;
639:   MPI_Comm       comm;
640:   PetscInt       i,j;

644:   PetscObjectGetComm((PetscObject)A,&comm);
645:   if (right) {
646:     /* allocate R */
647:     PetscMalloc1(bA->nc, &R);
648:     /* Create the right vectors */
649:     for (j=0; j<bA->nc; j++) {
650:       for (i=0; i<bA->nr; i++) {
651:         if (bA->m[i][j]) {
652:           MatCreateVecs(bA->m[i][j],&R[j],NULL);
653:           break;
654:         }
655:       }
656:       if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
657:     }
658:     VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
659:     /* hand back control to the nest vector */
660:     for (j=0; j<bA->nc; j++) {
661:       VecDestroy(&R[j]);
662:     }
663:     PetscFree(R);
664:   }

666:   if (left) {
667:     /* allocate L */
668:     PetscMalloc1(bA->nr, &L);
669:     /* Create the left vectors */
670:     for (i=0; i<bA->nr; i++) {
671:       for (j=0; j<bA->nc; j++) {
672:         if (bA->m[i][j]) {
673:           MatCreateVecs(bA->m[i][j],NULL,&L[i]);
674:           break;
675:         }
676:       }
677:       if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
678:     }

680:     VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
681:     for (i=0; i<bA->nr; i++) {
682:       VecDestroy(&L[i]);
683:     }

685:     PetscFree(L);
686:   }
687:   return(0);
688: }

690: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
691: {
692:   Mat_Nest       *bA = (Mat_Nest*)A->data;
693:   PetscBool      isascii,viewSub = PETSC_FALSE;
694:   PetscInt       i,j;

698:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
699:   if (isascii) {

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

706:     PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");
707:     for (i=0; i<bA->nr; i++) {
708:       for (j=0; j<bA->nc; j++) {
709:         MatType   type;
710:         char      name[256] = "",prefix[256] = "";
711:         PetscInt  NR,NC;
712:         PetscBool isNest = PETSC_FALSE;

714:         if (!bA->m[i][j]) {
715:           PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);
716:           continue;
717:         }
718:         MatGetSize(bA->m[i][j],&NR,&NC);
719:         MatGetType(bA->m[i][j], &type);
720:         if (((PetscObject)bA->m[i][j])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);}
721:         if (((PetscObject)bA->m[i][j])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);}
722:         PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);

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

726:         if (isNest || viewSub) {
727:           PetscViewerASCIIPushTab(viewer);  /* push1 */
728:           MatView(bA->m[i][j],viewer);
729:           PetscViewerASCIIPopTab(viewer);    /* pop1 */
730:         }
731:       }
732:     }
733:     PetscViewerASCIIPopTab(viewer);    /* pop0 */
734:   }
735:   return(0);
736: }

738: static PetscErrorCode MatZeroEntries_Nest(Mat A)
739: {
740:   Mat_Nest       *bA = (Mat_Nest*)A->data;
741:   PetscInt       i,j;

745:   for (i=0; i<bA->nr; i++) {
746:     for (j=0; j<bA->nc; j++) {
747:       if (!bA->m[i][j]) continue;
748:       MatZeroEntries(bA->m[i][j]);
749:     }
750:   }
751:   return(0);
752: }

754: static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str)
755: {
756:   Mat_Nest       *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data;
757:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;
759:   PetscBool      nnzstate = PETSC_FALSE;

762:   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);
763:   for (i=0; i<nr; i++) {
764:     for (j=0; j<nc; j++) {
765:       PetscObjectState subnnzstate = 0;
766:       if (bA->m[i][j] && bB->m[i][j]) {
767:         MatCopy(bA->m[i][j],bB->m[i][j],str);
768:       } 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);
769:       MatGetNonzeroState(bB->m[i][j],&subnnzstate);
770:       nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i*nc+j] != subnnzstate);
771:       bB->nnzstate[i*nc+j] = subnnzstate;
772:     }
773:   }
774:   if (nnzstate) B->nonzerostate++;
775:   return(0);
776: }

778: static PetscErrorCode MatAXPY_Nest(Mat Y,PetscScalar a,Mat X,MatStructure str)
779: {
780:   Mat_Nest       *bY = (Mat_Nest*)Y->data,*bX = (Mat_Nest*)X->data;
781:   PetscInt       i,j,nr = bY->nr,nc = bY->nc;
783:   PetscBool      nnzstate = PETSC_FALSE;

786:   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);
787:   for (i=0; i<nr; i++) {
788:     for (j=0; j<nc; j++) {
789:       PetscObjectState subnnzstate = 0;
790:       if (bY->m[i][j] && bX->m[i][j]) {
791:         MatAXPY(bY->m[i][j],a,bX->m[i][j],str);
792:       } else if (bX->m[i][j]) {
793:         Mat M;

795:         if (str != DIFFERENT_NONZERO_PATTERN) SETERRQ2(PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Matrix block does not exist at %D,%D. Use DIFFERENT_NONZERO_PATTERN",i,j);
796:         MatDuplicate(bX->m[i][j],MAT_COPY_VALUES,&M);
797:         MatNestSetSubMat(Y,i,j,M);
798:         MatDestroy(&M);
799:       }
800:       if (bY->m[i][j]) { MatGetNonzeroState(bY->m[i][j],&subnnzstate); }
801:       nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i*nc+j] != subnnzstate);
802:       bY->nnzstate[i*nc+j] = subnnzstate;
803:     }
804:   }
805:   if (nnzstate) Y->nonzerostate++;
806:   return(0);
807: }

809: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
810: {
811:   Mat_Nest       *bA = (Mat_Nest*)A->data;
812:   Mat            *b;
813:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

817:   PetscMalloc1(nr*nc,&b);
818:   for (i=0; i<nr; i++) {
819:     for (j=0; j<nc; j++) {
820:       if (bA->m[i][j]) {
821:         MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
822:       } else {
823:         b[i*nc+j] = NULL;
824:       }
825:     }
826:   }
827:   MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
828:   /* Give the new MatNest exclusive ownership */
829:   for (i=0; i<nr*nc; i++) {
830:     MatDestroy(&b[i]);
831:   }
832:   PetscFree(b);

834:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
835:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
836:   return(0);
837: }

839: /* nest api */
840: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
841: {
842:   Mat_Nest *bA = (Mat_Nest*)A->data;

845:   if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
846:   if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
847:   *mat = bA->m[idxm][jdxm];
848:   return(0);
849: }

851: /*@
852:  MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.

854:  Not collective

856:  Input Parameters:
857: +   A  - nest matrix
858: .   idxm - index of the matrix within the nest matrix
859: -   jdxm - index of the matrix within the nest matrix

861:  Output Parameter:
862: .   sub - matrix at index idxm,jdxm within the nest matrix

864:  Level: developer

866: .seealso: MatNestGetSize(), MatNestGetSubMats()
867: @*/
868: PetscErrorCode  MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
869: {

873:   PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
874:   return(0);
875: }

877: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
878: {
879:   Mat_Nest       *bA = (Mat_Nest*)A->data;
880:   PetscInt       m,n,M,N,mi,ni,Mi,Ni;

884:   if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
885:   if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
886:   MatGetLocalSize(mat,&m,&n);
887:   MatGetSize(mat,&M,&N);
888:   ISGetLocalSize(bA->isglobal.row[idxm],&mi);
889:   ISGetSize(bA->isglobal.row[idxm],&Mi);
890:   ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
891:   ISGetSize(bA->isglobal.col[jdxm],&Ni);
892:   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);
893:   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);

895:   /* do not increase object state */
896:   if (mat == bA->m[idxm][jdxm]) return(0);

898:   PetscObjectReference((PetscObject)mat);
899:   MatDestroy(&bA->m[idxm][jdxm]);
900:   bA->m[idxm][jdxm] = mat;
901:   PetscObjectStateIncrease((PetscObject)A);
902:   MatGetNonzeroState(mat,&bA->nnzstate[idxm*bA->nc+jdxm]);
903:   A->nonzerostate++;
904:   return(0);
905: }

907: /*@
908:  MatNestSetSubMat - Set a single submatrix in the nest matrix.

910:  Logically collective on the submatrix communicator

912:  Input Parameters:
913: +   A  - nest matrix
914: .   idxm - index of the matrix within the nest matrix
915: .   jdxm - index of the matrix within the nest matrix
916: -   sub - matrix at index idxm,jdxm within the nest matrix

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

921:  This increments the reference count of the submatrix.

923:  Level: developer

925: .seealso: MatNestSetSubMats(), MatNestGetSubMats()
926: @*/
927: PetscErrorCode  MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
928: {

932:   PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
933:   return(0);
934: }

936: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
937: {
938:   Mat_Nest *bA = (Mat_Nest*)A->data;

941:   if (M)   *M   = bA->nr;
942:   if (N)   *N   = bA->nc;
943:   if (mat) *mat = bA->m;
944:   return(0);
945: }

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

950:  Not collective

952:  Input Parameters:
953: .   A  - nest matrix

955:  Output Parameter:
956: +   M - number of rows in the nest matrix
957: .   N - number of cols in the nest matrix
958: -   mat - 2d array of matrices

960:  Notes:

962:  The user should not free the array mat.

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

968:  Level: developer

970: .seealso: MatNestGetSize(), MatNestGetSubMat()
971: @*/
972: PetscErrorCode  MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
973: {

977:   PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
978:   return(0);
979: }

981: PetscErrorCode  MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
982: {
983:   Mat_Nest *bA = (Mat_Nest*)A->data;

986:   if (M) *M = bA->nr;
987:   if (N) *N = bA->nc;
988:   return(0);
989: }

991: /*@
992:  MatNestGetSize - Returns the size of the nest matrix.

994:  Not collective

996:  Input Parameters:
997: .   A  - nest matrix

999:  Output Parameter:
1000: +   M - number of rows in the nested mat
1001: -   N - number of cols in the nested mat

1003:  Notes:

1005:  Level: developer

1007: .seealso: MatNestGetSubMat(), MatNestGetSubMats()
1008: @*/
1009: PetscErrorCode  MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
1010: {

1014:   PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
1015:   return(0);
1016: }

1018: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
1019: {
1020:   Mat_Nest *vs = (Mat_Nest*)A->data;
1021:   PetscInt i;

1024:   if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
1025:   if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
1026:   return(0);
1027: }

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

1032:  Not collective

1034:  Input Parameters:
1035: .   A  - nest matrix

1037:  Output Parameter:
1038: +   rows - array of row index sets
1039: -   cols - array of column index sets

1041:  Level: advanced

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

1046: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs()
1047: @*/
1048: PetscErrorCode  MatNestGetISs(Mat A,IS rows[],IS cols[])
1049: {

1054:   PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1055:   return(0);
1056: }

1058: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
1059: {
1060:   Mat_Nest *vs = (Mat_Nest*)A->data;
1061:   PetscInt i;

1064:   if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
1065:   if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
1066:   return(0);
1067: }

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

1072:  Not collective

1074:  Input Parameters:
1075: .   A  - nest matrix

1077:  Output Parameter:
1078: +   rows - array of row index sets (or NULL to ignore)
1079: -   cols - array of column index sets (or NULL to ignore)

1081:  Level: advanced

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

1086: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs()
1087: @*/
1088: PetscErrorCode  MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
1089: {

1094:   PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1095:   return(0);
1096: }

1098: PetscErrorCode  MatNestSetVecType_Nest(Mat A,VecType vtype)
1099: {
1101:   PetscBool      flg;

1104:   PetscStrcmp(vtype,VECNEST,&flg);
1105:   /* In reality, this only distinguishes VECNEST and "other" */
1106:   if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1107:   else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
1108:   return(0);
1109: }

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

1114:  Not collective

1116:  Input Parameters:
1117: +  A  - nest matrix
1118: -  vtype - type to use for creating vectors

1120:  Notes:

1122:  Level: developer

1124: .seealso: MatCreateVecs()
1125: @*/
1126: PetscErrorCode  MatNestSetVecType(Mat A,VecType vtype)
1127: {

1131:   PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
1132:   return(0);
1133: }

1135: PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1136: {
1137:   Mat_Nest       *s = (Mat_Nest*)A->data;
1138:   PetscInt       i,j,m,n,M,N;
1140:   PetscBool      cong;

1143:   MatReset_Nest(A);

1145:   s->nr = nr;
1146:   s->nc = nc;

1148:   /* Create space for submatrices */
1149:   PetscMalloc1(nr,&s->m);
1150:   for (i=0; i<nr; i++) {
1151:     PetscMalloc1(nc,&s->m[i]);
1152:   }
1153:   for (i=0; i<nr; i++) {
1154:     for (j=0; j<nc; j++) {
1155:       s->m[i][j] = a[i*nc+j];
1156:       if (a[i*nc+j]) {
1157:         PetscObjectReference((PetscObject)a[i*nc+j]);
1158:       }
1159:     }
1160:   }

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

1164:   PetscMalloc1(nr,&s->row_len);
1165:   PetscMalloc1(nc,&s->col_len);
1166:   for (i=0; i<nr; i++) s->row_len[i]=-1;
1167:   for (j=0; j<nc; j++) s->col_len[j]=-1;

1169:   PetscCalloc1(nr*nc,&s->nnzstate);
1170:   for (i=0; i<nr; i++) {
1171:     for (j=0; j<nc; j++) {
1172:       if (s->m[i][j]) {
1173:         MatGetNonzeroState(s->m[i][j],&s->nnzstate[i*nc+j]);
1174:       }
1175:     }
1176:   }

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

1180:   PetscLayoutSetSize(A->rmap,M);
1181:   PetscLayoutSetLocalSize(A->rmap,m);
1182:   PetscLayoutSetSize(A->cmap,N);
1183:   PetscLayoutSetLocalSize(A->cmap,n);

1185:   PetscLayoutSetUp(A->rmap);
1186:   PetscLayoutSetUp(A->cmap);

1188:   /* disable operations that are not supported for non-square matrices,
1189:      or matrices for which is_row != is_col  */
1190:   MatHasCongruentLayouts(A,&cong);
1191:   if (cong && nr != nc) cong = PETSC_FALSE;
1192:   if (cong) {
1193:     for (i = 0; cong && i < nr; i++) {
1194:       ISEqualUnsorted(s->isglobal.row[i],s->isglobal.col[i],&cong);
1195:     }
1196:   }
1197:   if (!cong) {
1198:     A->ops->getdiagonal = NULL;
1199:     A->ops->shift       = NULL;
1200:     A->ops->diagonalset = NULL;
1201:   }

1203:   PetscCalloc2(nr,&s->left,nc,&s->right);
1204:   PetscObjectStateIncrease((PetscObject)A);
1205:   A->nonzerostate++;
1206:   return(0);
1207: }

1209: /*@
1210:    MatNestSetSubMats - Sets the nested submatrices

1212:    Collective on Mat

1214:    Input Parameter:
1215: +  A - nested matrix
1216: .  nr - number of nested row blocks
1217: .  is_row - index sets for each nested row block, or NULL to make contiguous
1218: .  nc - number of nested column blocks
1219: .  is_col - index sets for each nested column block, or NULL to make contiguous
1220: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

1222:    Notes: this always resets any submatrix information previously set

1224:    Level: advanced

1226: .seealso: MatCreateNest(), MATNEST
1227: @*/
1228: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1229: {
1231:   PetscInt       i;

1235:   if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1236:   if (nr && is_row) {
1239:   }
1240:   if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1241:   if (nc && is_col) {
1244:   }
1246:   PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1247:   return(0);
1248: }

1250: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1251: {
1253:   PetscBool      flg;
1254:   PetscInt       i,j,m,mi,*ix;

1257:   *ltog = NULL;
1258:   for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1259:     if (islocal[i]) {
1260:       ISGetLocalSize(islocal[i],&mi);
1261:       flg  = PETSC_TRUE;      /* We found a non-trivial entry */
1262:     } else {
1263:       ISGetLocalSize(isglobal[i],&mi);
1264:     }
1265:     m += mi;
1266:   }
1267:   if (!flg) return(0);

1269:   PetscMalloc1(m,&ix);
1270:   for (i=0,m=0; i<n; i++) {
1271:     ISLocalToGlobalMapping smap = NULL;
1272:     Mat                    sub = NULL;
1273:     PetscSF                sf;
1274:     PetscLayout            map;
1275:     const PetscInt         *ix2;

1277:     if (!colflg) {
1278:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1279:     } else {
1280:       MatNestFindNonzeroSubMatCol(A,i,&sub);
1281:     }
1282:     if (sub) {
1283:       if (!colflg) {
1284:         MatGetLocalToGlobalMapping(sub,&smap,NULL);
1285:       } else {
1286:         MatGetLocalToGlobalMapping(sub,NULL,&smap);
1287:       }
1288:     }
1289:     /*
1290:        Now we need to extract the monolithic global indices that correspond to the given split global indices.
1291:        In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1292:     */
1293:     ISGetIndices(isglobal[i],&ix2);
1294:     if (islocal[i]) {
1295:       PetscInt *ilocal,*iremote;
1296:       PetscInt mil,nleaves;

1298:       ISGetLocalSize(islocal[i],&mi);
1299:       if (!smap) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing local to global map");
1300:       for (j=0; j<mi; j++) ix[m+j] = j;
1301:       ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);

1303:       /* PetscSFSetGraphLayout does not like negative indices */
1304:       PetscMalloc2(mi,&ilocal,mi,&iremote);
1305:       for (j=0, nleaves = 0; j<mi; j++) {
1306:         if (ix[m+j] < 0) continue;
1307:         ilocal[nleaves]  = j;
1308:         iremote[nleaves] = ix[m+j];
1309:         nleaves++;
1310:       }
1311:       ISGetLocalSize(isglobal[i],&mil);
1312:       PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
1313:       PetscLayoutCreate(PetscObjectComm((PetscObject)A),&map);
1314:       PetscLayoutSetLocalSize(map,mil);
1315:       PetscLayoutSetUp(map);
1316:       PetscSFSetGraphLayout(sf,map,nleaves,ilocal,PETSC_USE_POINTER,iremote);
1317:       PetscLayoutDestroy(&map);
1318:       PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m);
1319:       PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m);
1320:       PetscSFDestroy(&sf);
1321:       PetscFree2(ilocal,iremote);
1322:     } else {
1323:       ISGetLocalSize(isglobal[i],&mi);
1324:       for (j=0; j<mi; j++) ix[m+j] = ix2[i];
1325:     }
1326:     ISRestoreIndices(isglobal[i],&ix2);
1327:     m   += mi;
1328:   }
1329:   ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);
1330:   return(0);
1331: }


1334: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1335: /*
1336:   nprocessors = NP
1337:   Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1338:        proc 0: => (g_0,h_0,)
1339:        proc 1: => (g_1,h_1,)
1340:        ...
1341:        proc nprocs-1: => (g_NP-1,h_NP-1,)

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

1346:             proc 0:
1347:     is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1348:             proc 1:
1349:     is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)

1351:             proc NP-1:
1352:     is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1353: */
1354: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1355: {
1356:   Mat_Nest       *vs = (Mat_Nest*)A->data;
1357:   PetscInt       i,j,offset,n,nsum,bs;
1359:   Mat            sub = NULL;

1362:   PetscMalloc1(nr,&vs->isglobal.row);
1363:   PetscMalloc1(nc,&vs->isglobal.col);
1364:   if (is_row) { /* valid IS is passed in */
1365:     /* refs on is[] are incremeneted */
1366:     for (i=0; i<vs->nr; i++) {
1367:       PetscObjectReference((PetscObject)is_row[i]);

1369:       vs->isglobal.row[i] = is_row[i];
1370:     }
1371:   } else {                      /* Create the ISs by inspecting sizes of a submatrix in each row */
1372:     nsum = 0;
1373:     for (i=0; i<vs->nr; i++) {  /* Add up the local sizes to compute the aggregate offset */
1374:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1375:       if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i);
1376:       MatGetLocalSize(sub,&n,NULL);
1377:       if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1378:       nsum += n;
1379:     }
1380:     MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1381:     offset -= nsum;
1382:     for (i=0; i<vs->nr; i++) {
1383:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1384:       MatGetLocalSize(sub,&n,NULL);
1385:       MatGetBlockSize(sub,&bs);
1386:       ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);
1387:       ISSetBlockSize(vs->isglobal.row[i],bs);
1388:       offset += n;
1389:     }
1390:   }

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

1397:       vs->isglobal.col[j] = is_col[j];
1398:     }
1399:   } else {                      /* Create the ISs by inspecting sizes of a submatrix in each column */
1400:     offset = A->cmap->rstart;
1401:     nsum   = 0;
1402:     for (j=0; j<vs->nc; j++) {
1403:       MatNestFindNonzeroSubMatCol(A,j,&sub);
1404:       if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i);
1405:       MatGetLocalSize(sub,NULL,&n);
1406:       if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1407:       nsum += n;
1408:     }
1409:     MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1410:     offset -= nsum;
1411:     for (j=0; j<vs->nc; j++) {
1412:       MatNestFindNonzeroSubMatCol(A,j,&sub);
1413:       MatGetLocalSize(sub,NULL,&n);
1414:       MatGetBlockSize(sub,&bs);
1415:       ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);
1416:       ISSetBlockSize(vs->isglobal.col[j],bs);
1417:       offset += n;
1418:     }
1419:   }

1421:   /* Set up the local ISs */
1422:   PetscMalloc1(vs->nr,&vs->islocal.row);
1423:   PetscMalloc1(vs->nc,&vs->islocal.col);
1424:   for (i=0,offset=0; i<vs->nr; i++) {
1425:     IS                     isloc;
1426:     ISLocalToGlobalMapping rmap = NULL;
1427:     PetscInt               nlocal,bs;
1428:     MatNestFindNonzeroSubMatRow(A,i,&sub);
1429:     if (sub) {MatGetLocalToGlobalMapping(sub,&rmap,NULL);}
1430:     if (rmap) {
1431:       MatGetBlockSize(sub,&bs);
1432:       ISLocalToGlobalMappingGetSize(rmap,&nlocal);
1433:       ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1434:       ISSetBlockSize(isloc,bs);
1435:     } else {
1436:       nlocal = 0;
1437:       isloc  = NULL;
1438:     }
1439:     vs->islocal.row[i] = isloc;
1440:     offset            += nlocal;
1441:   }
1442:   for (i=0,offset=0; i<vs->nc; i++) {
1443:     IS                     isloc;
1444:     ISLocalToGlobalMapping cmap = NULL;
1445:     PetscInt               nlocal,bs;
1446:     MatNestFindNonzeroSubMatCol(A,i,&sub);
1447:     if (sub) {MatGetLocalToGlobalMapping(sub,NULL,&cmap);}
1448:     if (cmap) {
1449:       MatGetBlockSize(sub,&bs);
1450:       ISLocalToGlobalMappingGetSize(cmap,&nlocal);
1451:       ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1452:       ISSetBlockSize(isloc,bs);
1453:     } else {
1454:       nlocal = 0;
1455:       isloc  = NULL;
1456:     }
1457:     vs->islocal.col[i] = isloc;
1458:     offset            += nlocal;
1459:   }

1461:   /* Set up the aggregate ISLocalToGlobalMapping */
1462:   {
1463:     ISLocalToGlobalMapping rmap,cmap;
1464:     MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);
1465:     MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);
1466:     if (rmap && cmap) {MatSetLocalToGlobalMapping(A,rmap,cmap);}
1467:     ISLocalToGlobalMappingDestroy(&rmap);
1468:     ISLocalToGlobalMappingDestroy(&cmap);
1469:   }

1471: #if defined(PETSC_USE_DEBUG)
1472:   for (i=0; i<vs->nr; i++) {
1473:     for (j=0; j<vs->nc; j++) {
1474:       PetscInt m,n,M,N,mi,ni,Mi,Ni;
1475:       Mat      B = vs->m[i][j];
1476:       if (!B) continue;
1477:       MatGetSize(B,&M,&N);
1478:       MatGetLocalSize(B,&m,&n);
1479:       ISGetSize(vs->isglobal.row[i],&Mi);
1480:       ISGetSize(vs->isglobal.col[j],&Ni);
1481:       ISGetLocalSize(vs->isglobal.row[i],&mi);
1482:       ISGetLocalSize(vs->isglobal.col[j],&ni);
1483:       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);
1484:       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);
1485:     }
1486:   }
1487: #endif

1489:   /* Set A->assembled if all non-null blocks are currently assembled */
1490:   for (i=0; i<vs->nr; i++) {
1491:     for (j=0; j<vs->nc; j++) {
1492:       if (vs->m[i][j] && !vs->m[i][j]->assembled) return(0);
1493:     }
1494:   }
1495:   A->assembled = PETSC_TRUE;
1496:   return(0);
1497: }

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

1502:    Collective on Mat

1504:    Input Parameter:
1505: +  comm - Communicator for the new Mat
1506: .  nr - number of nested row blocks
1507: .  is_row - index sets for each nested row block, or NULL to make contiguous
1508: .  nc - number of nested column blocks
1509: .  is_col - index sets for each nested column block, or NULL to make contiguous
1510: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

1512:    Output Parameter:
1513: .  B - new matrix

1515:    Level: advanced

1517: .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST
1518: @*/
1519: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1520: {
1521:   Mat            A;

1525:   *B   = 0;
1526:   MatCreate(comm,&A);
1527:   MatSetType(A,MATNEST);
1528:   A->preallocated = PETSC_TRUE;
1529:   MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1530:   *B   = A;
1531:   return(0);
1532: }

1534: static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1535: {
1536:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1537:   Mat            *trans;
1538:   PetscScalar    **avv;
1539:   PetscScalar    *vv;
1540:   PetscInt       **aii,**ajj;
1541:   PetscInt       *ii,*jj,*ci;
1542:   PetscInt       nr,nc,nnz,i,j;
1543:   PetscBool      done;

1547:   MatGetSize(A,&nr,&nc);
1548:   if (reuse == MAT_REUSE_MATRIX) {
1549:     PetscInt rnr;

1551:     MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);
1552:     if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ");
1553:     if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows");
1554:     MatSeqAIJGetArray(*newmat,&vv);
1555:   }
1556:   /* extract CSR for nested SeqAIJ matrices */
1557:   nnz  = 0;
1558:   PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);
1559:   for (i=0; i<nest->nr; ++i) {
1560:     for (j=0; j<nest->nc; ++j) {
1561:       Mat B = nest->m[i][j];
1562:       if (B) {
1563:         PetscScalar *naa;
1564:         PetscInt    *nii,*njj,nnr;
1565:         PetscBool   istrans;

1567:         PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1568:         if (istrans) {
1569:           Mat Bt;

1571:           MatTransposeGetMat(B,&Bt);
1572:           MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);
1573:           B    = trans[i*nest->nc+j];
1574:         }
1575:         MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);
1576:         if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ");
1577:         MatSeqAIJGetArray(B,&naa);
1578:         nnz += nii[nnr];

1580:         aii[i*nest->nc+j] = nii;
1581:         ajj[i*nest->nc+j] = njj;
1582:         avv[i*nest->nc+j] = naa;
1583:       }
1584:     }
1585:   }
1586:   if (reuse != MAT_REUSE_MATRIX) {
1587:     PetscMalloc1(nr+1,&ii);
1588:     PetscMalloc1(nnz,&jj);
1589:     PetscMalloc1(nnz,&vv);
1590:   } else {
1591:     if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros");
1592:   }

1594:   /* new row pointer */
1595:   PetscArrayzero(ii,nr+1);
1596:   for (i=0; i<nest->nr; ++i) {
1597:     PetscInt       ncr,rst;

1599:     ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1600:     ISGetLocalSize(nest->isglobal.row[i],&ncr);
1601:     for (j=0; j<nest->nc; ++j) {
1602:       if (aii[i*nest->nc+j]) {
1603:         PetscInt    *nii = aii[i*nest->nc+j];
1604:         PetscInt    ir;

1606:         for (ir=rst; ir<ncr+rst; ++ir) {
1607:           ii[ir+1] += nii[1]-nii[0];
1608:           nii++;
1609:         }
1610:       }
1611:     }
1612:   }
1613:   for (i=0; i<nr; i++) ii[i+1] += ii[i];

1615:   /* construct CSR for the new matrix */
1616:   PetscCalloc1(nr,&ci);
1617:   for (i=0; i<nest->nr; ++i) {
1618:     PetscInt       ncr,rst;

1620:     ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1621:     ISGetLocalSize(nest->isglobal.row[i],&ncr);
1622:     for (j=0; j<nest->nc; ++j) {
1623:       if (aii[i*nest->nc+j]) {
1624:         PetscScalar *nvv = avv[i*nest->nc+j];
1625:         PetscInt    *nii = aii[i*nest->nc+j];
1626:         PetscInt    *njj = ajj[i*nest->nc+j];
1627:         PetscInt    ir,cst;

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

1633:           for (ij=0;ij<rsize;ij++) {
1634:             jj[ist+ij] = *njj+cst;
1635:             vv[ist+ij] = *nvv;
1636:             njj++;
1637:             nvv++;
1638:           }
1639:           ci[ir] += rsize;
1640:           nii++;
1641:         }
1642:       }
1643:     }
1644:   }
1645:   PetscFree(ci);

1647:   /* restore info */
1648:   for (i=0; i<nest->nr; ++i) {
1649:     for (j=0; j<nest->nc; ++j) {
1650:       Mat B = nest->m[i][j];
1651:       if (B) {
1652:         PetscInt nnr = 0, k = i*nest->nc+j;

1654:         B    = (trans[k] ? trans[k] : B);
1655:         MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);
1656:         if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ");
1657:         MatSeqAIJRestoreArray(B,&avv[k]);
1658:         MatDestroy(&trans[k]);
1659:       }
1660:     }
1661:   }
1662:   PetscFree4(aii,ajj,avv,trans);

1664:   /* finalize newmat */
1665:   if (reuse == MAT_INITIAL_MATRIX) {
1666:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);
1667:   } else if (reuse == MAT_INPLACE_MATRIX) {
1668:     Mat B;

1670:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);
1671:     MatHeaderReplace(A,&B);
1672:   }
1673:   MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1674:   MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1675:   {
1676:     Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data);
1677:     a->free_a     = PETSC_TRUE;
1678:     a->free_ij    = PETSC_TRUE;
1679:   }
1680:   return(0);
1681: }

1683: PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1684: {
1686:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1687:   PetscInt       m,n,M,N,i,j,k,*dnnz,*onnz,rstart;
1688:   PetscInt       cstart,cend;
1689:   PetscMPIInt    size;
1690:   Mat            C;

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

1698:     PetscStrcmp(newtype,MATAIJ,&fast);
1699:     if (!fast) {
1700:       PetscStrcmp(newtype,MATSEQAIJ,&fast);
1701:     }
1702:     for (i=0; i<nest->nr && fast; ++i) {
1703:       for (j=0; j<nest->nc && fast; ++j) {
1704:         Mat B = nest->m[i][j];
1705:         if (B) {
1706:           PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);
1707:           if (!fast) {
1708:             PetscBool istrans;

1710:             PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1711:             if (istrans) {
1712:               Mat Bt;

1714:               MatTransposeGetMat(B,&Bt);
1715:               PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);
1716:             }
1717:           }
1718:         }
1719:       }
1720:     }
1721:     for (i=0, nf=0; i<nest->nr && fast; ++i) {
1722:       PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);
1723:       if (fast) {
1724:         PetscInt f,s;

1726:         ISStrideGetInfo(nest->isglobal.row[i],&f,&s);
1727:         if (f != nf || s != 1) { fast = PETSC_FALSE; }
1728:         else {
1729:           ISGetSize(nest->isglobal.row[i],&f);
1730:           nf  += f;
1731:         }
1732:       }
1733:     }
1734:     for (i=0, nf=0; i<nest->nc && fast; ++i) {
1735:       PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);
1736:       if (fast) {
1737:         PetscInt f,s;

1739:         ISStrideGetInfo(nest->isglobal.col[i],&f,&s);
1740:         if (f != nf || s != 1) { fast = PETSC_FALSE; }
1741:         else {
1742:           ISGetSize(nest->isglobal.col[i],&f);
1743:           nf  += f;
1744:         }
1745:       }
1746:     }
1747:     if (fast) {
1748:       MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);
1749:       return(0);
1750:     }
1751:   }
1752:   MatGetSize(A,&M,&N);
1753:   MatGetLocalSize(A,&m,&n);
1754:   MatGetOwnershipRangeColumn(A,&cstart,&cend);
1755:   switch (reuse) {
1756:   case MAT_INITIAL_MATRIX:
1757:     MatCreate(PetscObjectComm((PetscObject)A),&C);
1758:     MatSetType(C,newtype);
1759:     MatSetSizes(C,m,n,M,N);
1760:     *newmat = C;
1761:     break;
1762:   case MAT_REUSE_MATRIX:
1763:     C = *newmat;
1764:     break;
1765:   default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse");
1766:   }
1767:   PetscMalloc1(2*m,&dnnz);
1768:   onnz = dnnz + m;
1769:   for (k=0; k<m; k++) {
1770:     dnnz[k] = 0;
1771:     onnz[k] = 0;
1772:   }
1773:   for (j=0; j<nest->nc; ++j) {
1774:     IS             bNis;
1775:     PetscInt       bN;
1776:     const PetscInt *bNindices;
1777:     /* Using global column indices and ISAllGather() is not scalable. */
1778:     ISAllGather(nest->isglobal.col[j], &bNis);
1779:     ISGetSize(bNis, &bN);
1780:     ISGetIndices(bNis,&bNindices);
1781:     for (i=0; i<nest->nr; ++i) {
1782:       PetscSF        bmsf;
1783:       PetscSFNode    *iremote;
1784:       Mat            B;
1785:       PetscInt       bm, *sub_dnnz,*sub_onnz, br;
1786:       const PetscInt *bmindices;
1787:       B = nest->m[i][j];
1788:       if (!B) continue;
1789:       ISGetLocalSize(nest->isglobal.row[i],&bm);
1790:       ISGetIndices(nest->isglobal.row[i],&bmindices);
1791:       PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);
1792:       PetscMalloc1(bm,&iremote);
1793:       PetscMalloc1(bm,&sub_dnnz);
1794:       PetscMalloc1(bm,&sub_onnz);
1795:       for (k = 0; k < bm; ++k){
1796:             sub_dnnz[k] = 0;
1797:             sub_onnz[k] = 0;
1798:       }
1799:       /*
1800:        Locate the owners for all of the locally-owned global row indices for this row block.
1801:        These determine the roots of PetscSF used to communicate preallocation data to row owners.
1802:        The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
1803:        */
1804:       MatGetOwnershipRange(B,&rstart,NULL);
1805:       for (br = 0; br < bm; ++br) {
1806:         PetscInt       row = bmindices[br], rowowner = 0, brncols, col;
1807:         const PetscInt *brcols;
1808:         PetscInt       rowrel = 0; /* row's relative index on its owner rank */
1809:         PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);
1810:         /* how many roots  */
1811:         iremote[br].rank = rowowner; iremote[br].index = rowrel;           /* edge from bmdnnz to dnnz */
1812:         /* get nonzero pattern */
1813:         MatGetRow(B,br+rstart,&brncols,&brcols,NULL);
1814:         for (k=0; k<brncols; k++) {
1815:           col  = bNindices[brcols[k]];
1816:           if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) {
1817:             sub_dnnz[br]++;
1818:           } else {
1819:             sub_onnz[br]++;
1820:           }
1821:         }
1822:         MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);
1823:       }
1824:       ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1825:       /* bsf will have to take care of disposing of bedges. */
1826:       PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
1827:       PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
1828:       PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
1829:       PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
1830:       PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
1831:       PetscFree(sub_dnnz);
1832:       PetscFree(sub_onnz);
1833:       PetscSFDestroy(&bmsf);
1834:     }
1835:     ISRestoreIndices(bNis,&bNindices);
1836:     ISDestroy(&bNis);
1837:   }
1838:   /* Resize preallocation if overestimated */
1839:   for (i=0;i<m;i++) {
1840:     dnnz[i] = PetscMin(dnnz[i],A->cmap->n);
1841:     onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n);
1842:   }
1843:   MatSeqAIJSetPreallocation(C,0,dnnz);
1844:   MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);
1845:   PetscFree(dnnz);

1847:   /* Fill by row */
1848:   for (j=0; j<nest->nc; ++j) {
1849:     /* Using global column indices and ISAllGather() is not scalable. */
1850:     IS             bNis;
1851:     PetscInt       bN;
1852:     const PetscInt *bNindices;
1853:     ISAllGather(nest->isglobal.col[j], &bNis);
1854:     ISGetSize(bNis,&bN);
1855:     ISGetIndices(bNis,&bNindices);
1856:     for (i=0; i<nest->nr; ++i) {
1857:       Mat            B;
1858:       PetscInt       bm, br;
1859:       const PetscInt *bmindices;
1860:       B = nest->m[i][j];
1861:       if (!B) continue;
1862:       ISGetLocalSize(nest->isglobal.row[i],&bm);
1863:       ISGetIndices(nest->isglobal.row[i],&bmindices);
1864:       MatGetOwnershipRange(B,&rstart,NULL);
1865:       for (br = 0; br < bm; ++br) {
1866:         PetscInt          row = bmindices[br], brncols,  *cols;
1867:         const PetscInt    *brcols;
1868:         const PetscScalar *brcoldata;
1869:         MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1870:         PetscMalloc1(brncols,&cols);
1871:         for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
1872:         /*
1873:           Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
1874:           Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
1875:          */
1876:         MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);
1877:         MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1878:         PetscFree(cols);
1879:       }
1880:       ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1881:     }
1882:     ISRestoreIndices(bNis,&bNindices);
1883:     ISDestroy(&bNis);
1884:   }
1885:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
1886:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
1887:   return(0);
1888: }

1890: PetscErrorCode  MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool  *has)
1891: {
1893:   *has = PETSC_FALSE;
1894:   if (op == MATOP_MULT_TRANSPOSE) {
1895:     Mat_Nest       *bA = (Mat_Nest*)mat->data;
1896:     PetscInt       i,j,nr = bA->nr,nc = bA->nc;
1898:     PetscBool      flg;

1900:     for (j=0; j<nc; j++) {
1901:       for (i=0; i<nr; i++) {
1902:         if (!bA->m[i][j]) continue;
1903:         MatHasOperation(bA->m[i][j],MATOP_MULT_TRANSPOSE_ADD,&flg);
1904:         if (!flg) return(0);
1905:       }
1906:     }
1907:   }
1908:   if (((void**)mat->ops)[op]) *has =  PETSC_TRUE;
1909:   return(0);
1910: }

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

1915:   Level: intermediate

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

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

1926: .seealso: MatCreate(), MatType, MatCreateNest()
1927: M*/
1928: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
1929: {
1930:   Mat_Nest       *s;

1934:   PetscNewLog(A,&s);
1935:   A->data = (void*)s;

1937:   s->nr            = -1;
1938:   s->nc            = -1;
1939:   s->m             = NULL;
1940:   s->splitassembly = PETSC_FALSE;

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

1944:   A->ops->mult                  = MatMult_Nest;
1945:   A->ops->multadd               = MatMultAdd_Nest;
1946:   A->ops->multtranspose         = MatMultTranspose_Nest;
1947:   A->ops->multtransposeadd      = MatMultTransposeAdd_Nest;
1948:   A->ops->transpose             = MatTranspose_Nest;
1949:   A->ops->assemblybegin         = MatAssemblyBegin_Nest;
1950:   A->ops->assemblyend           = MatAssemblyEnd_Nest;
1951:   A->ops->zeroentries           = MatZeroEntries_Nest;
1952:   A->ops->copy                  = MatCopy_Nest;
1953:   A->ops->axpy                  = MatAXPY_Nest;
1954:   A->ops->duplicate             = MatDuplicate_Nest;
1955:   A->ops->createsubmatrix       = MatCreateSubMatrix_Nest;
1956:   A->ops->destroy               = MatDestroy_Nest;
1957:   A->ops->view                  = MatView_Nest;
1958:   A->ops->getvecs               = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
1959:   A->ops->getlocalsubmatrix     = MatGetLocalSubMatrix_Nest;
1960:   A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
1961:   A->ops->getdiagonal           = MatGetDiagonal_Nest;
1962:   A->ops->diagonalscale         = MatDiagonalScale_Nest;
1963:   A->ops->scale                 = MatScale_Nest;
1964:   A->ops->shift                 = MatShift_Nest;
1965:   A->ops->diagonalset           = MatDiagonalSet_Nest;
1966:   A->ops->setrandom             = MatSetRandom_Nest;
1967:   A->ops->hasoperation          = MatHasOperation_Nest;

1969:   A->spptr        = 0;
1970:   A->assembled    = PETSC_FALSE;

1972:   /* expose Nest api's */
1973:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",      MatNestGetSubMat_Nest);
1974:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",      MatNestSetSubMat_Nest);
1975:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",     MatNestGetSubMats_Nest);
1976:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",        MatNestGetSize_Nest);
1977:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",         MatNestGetISs_Nest);
1978:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",    MatNestGetLocalISs_Nest);
1979:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",     MatNestSetVecType_Nest);
1980:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",     MatNestSetSubMats_Nest);
1981:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",MatConvert_Nest_AIJ);
1982:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",MatConvert_Nest_AIJ);
1983:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",   MatConvert_Nest_AIJ);
1984:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",    MatConvert_Nest_IS);

1986:   PetscObjectChangeTypeName((PetscObject)A,MATNEST);
1987:   return(0);
1988: }