Actual source code: matnest.c

petsc-master 2020-07-09
Report Typos and Errors
  1:  #include <../src/mat/impls/nest/matnestimpl.h>
  2:  #include <../src/mat/impls/aij/seq/aij.h>
  3:  #include <petscsf.h>

  5: static PetscErrorCode MatSetUp_NestIS_Private(Mat,PetscInt,const IS[],PetscInt,const IS[]);
  6: static PetscErrorCode MatCreateVecs_Nest(Mat,Vec*,Vec*);
  7: static PetscErrorCode MatReset_Nest(Mat);

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

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

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

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

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

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

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

 89: typedef struct {
 90:   Mat          *workC;    /* array of Mat with specific containers depending on the underlying MatMatMult implementation */
 91:   PetscScalar  *tarray;   /* buffer for storing all temporary products A[i][j] B[j] */
 92:   PetscInt     *dm,*dn,k; /* displacements and number of submatrices */
 93: } Nest_Dense;

 95: PETSC_INTERN PetscErrorCode MatProductNumeric_Nest_Dense(Mat C)
 96: {
 97:   Mat_Nest          *bA;
 98:   Nest_Dense        *contents;
 99:   Mat               viewB,viewC,productB,workC;
100:   const PetscScalar *barray;
101:   PetscScalar       *carray;
102:   PetscInt          i,j,M,N,nr,nc,ldb,ldc;
103:   PetscErrorCode    ierr;
104:   Mat               A,B;

107:   MatCheckProduct(C,3);
108:   A    = C->product->A;
109:   B    = C->product->B;
110:   MatGetSize(B,NULL,&N);
111:   if (!N) {
112:     MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
113:     MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
114:     return(0);
115:   }
116:   contents = (Nest_Dense*)C->product->data;
117:   if (!contents) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data empty");
118:   bA   = (Mat_Nest*)A->data;
119:   nr   = bA->nr;
120:   nc   = bA->nc;
121:   MatDenseGetLDA(B,&ldb);
122:   MatDenseGetLDA(C,&ldc);
123:   MatZeroEntries(C);
124:   MatDenseGetArrayRead(B,&barray);
125:   MatDenseGetArrayWrite(C,&carray);
126:   for (i=0; i<nr; i++) {
127:     ISGetSize(bA->isglobal.row[i],&M);
128:     MatCreateDense(PetscObjectComm((PetscObject)A),contents->dm[i+1]-contents->dm[i],PETSC_DECIDE,M,N,carray+contents->dm[i],&viewC);
129:     MatDenseSetLDA(viewC,ldc);
130:     for (j=0; j<nc; j++) {
131:       if (!bA->m[i][j]) continue;
132:       ISGetSize(bA->isglobal.col[j],&M);
133:       MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);
134:       MatDenseSetLDA(viewB,ldb);

136:       /* MatMatMultNumeric(bA->m[i][j],viewB,contents->workC[i*nc + j]); */
137:       workC             = contents->workC[i*nc + j];
138:       productB          = workC->product->B;
139:       workC->product->B = viewB; /* use newly created dense matrix viewB */
140:       MatProductNumeric(workC);
141:       MatDestroy(&viewB);
142:       workC->product->B = productB; /* resume original B */

144:       /* C[i] <- workC + C[i] */
145:       MatAXPY(viewC,1.0,contents->workC[i*nc + j],SAME_NONZERO_PATTERN);
146:     }
147:     MatDestroy(&viewC);
148:   }
149:   MatDenseRestoreArrayWrite(C,&carray);
150:   MatDenseRestoreArrayRead(B,&barray);

152:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
153:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
154:   return(0);
155: }

157: PetscErrorCode MatNest_DenseDestroy(void *ctx)
158: {
159:   Nest_Dense     *contents = (Nest_Dense*)ctx;
160:   PetscInt       i;

164:   PetscFree(contents->tarray);
165:   for (i=0; i<contents->k; i++) {
166:     MatDestroy(contents->workC + i);
167:   }
168:   PetscFree3(contents->dm,contents->dn,contents->workC);
169:   PetscFree(contents);
170:   return(0);
171: }

173: PETSC_INTERN PetscErrorCode MatProductSymbolic_Nest_Dense(Mat C)
174: {
175:   Mat_Nest          *bA;
176:   Mat               viewB,workC;
177:   const PetscScalar *barray;
178:   PetscInt          i,j,M,N,m,n,nr,nc,maxm = 0,ldb;
179:   Nest_Dense        *contents=NULL;
180:   PetscBool         cisdense;
181:   PetscErrorCode    ierr;
182:   Mat               A,B;
183:   PetscReal         fill;

186:   MatCheckProduct(C,4);
187:   if (C->product->data) SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_PLIB,"Product data not empty");
188:   A    = C->product->A;
189:   B    = C->product->B;
190:   fill = C->product->fill;
191:   bA   = (Mat_Nest*)A->data;
192:   nr   = bA->nr;
193:   nc   = bA->nc;
194:   MatGetLocalSize(B,NULL,&n);
195:   MatGetSize(B,NULL,&N);
196:   MatGetLocalSize(A,&m,NULL);
197:   MatGetSize(A,&M,NULL);
198:   MatSetSizes(C,m,n,M,N);
199:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATMPIDENSE,MATSEQDENSECUDA,MATMPIDENSECUDA,"");
200:   if (!cisdense) {
201:     MatSetType(C,((PetscObject)B)->type_name);
202:   }
203:   MatSetUp(C);
204:   if (!N) {
205:     C->ops->productnumeric = MatProductNumeric_Nest_Dense;
206:     return(0);
207:   }

209:   PetscNew(&contents);
210:   C->product->data = contents;
211:   C->product->destroy = MatNest_DenseDestroy;
212:   PetscCalloc3(nr+1,&contents->dm,nc+1,&contents->dn,nr*nc,&contents->workC);
213:   contents->k = nr*nc;
214:   for (i=0; i<nr; i++) {
215:     ISGetLocalSize(bA->isglobal.row[i],contents->dm + i+1);
216:     maxm = PetscMax(maxm,contents->dm[i+1]);
217:     contents->dm[i+1] += contents->dm[i];
218:   }
219:   for (i=0; i<nc; i++) {
220:     ISGetLocalSize(bA->isglobal.col[i],contents->dn + i+1);
221:     contents->dn[i+1] += contents->dn[i];
222:   }
223:   PetscMalloc1(maxm*N,&contents->tarray);
224:   MatDenseGetLDA(B,&ldb);
225:   MatGetSize(B,NULL,&N);
226:   MatDenseGetArrayRead(B,&barray);
227:   /* loops are permuted compared to MatMatMultNumeric so that viewB is created only once per column of A */
228:   for (j=0; j<nc; j++) {
229:     ISGetSize(bA->isglobal.col[j],&M);
230:     MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);
231:     MatDenseSetLDA(viewB,ldb);
232:     for (i=0; i<nr; i++) {
233:       if (!bA->m[i][j]) continue;
234:       /* MatMatMultSymbolic may attach a specific container (depending on MatType of bA->m[i][j]) to workC[i][j] */

236:       MatProductCreate(bA->m[i][j],viewB,NULL,&contents->workC[i*nc + j]);
237:       workC = contents->workC[i*nc + j];
238:       MatProductSetType(workC,MATPRODUCT_AB);
239:       MatProductSetAlgorithm(workC,"default");
240:       MatProductSetFill(workC,fill);
241:       MatProductSetFromOptions(workC);
242:       MatProductSymbolic(workC);

244:       /* since tarray will be shared by all Mat */
245:       MatSeqDenseSetPreallocation(workC,contents->tarray);
246:       MatMPIDenseSetPreallocation(workC,contents->tarray);
247:     }
248:     MatDestroy(&viewB);
249:   }
250:   MatDenseRestoreArrayRead(B,&barray);

252:   C->ops->productnumeric = MatProductNumeric_Nest_Dense;
253:   return(0);
254: }

256: /* --------------------------------------------------------- */
257: static PetscErrorCode MatProductSetFromOptions_Nest_Dense_AB(Mat C)
258: {
260:   C->ops->productsymbolic = MatProductSymbolic_Nest_Dense;
261:   return(0);
262: }

264: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Nest_Dense(Mat C)
265: {
267:   Mat_Product    *product = C->product;

270:   if (product->type == MATPRODUCT_AB) {
271:     MatProductSetFromOptions_Nest_Dense_AB(C);
272:   }
273:   return(0);
274: }
275: /* --------------------------------------------------------- */

277: static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y)
278: {
279:   Mat_Nest       *bA = (Mat_Nest*)A->data;
280:   Vec            *bx = bA->left,*by = bA->right;
281:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

285:   for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
286:   for (i=0; i<nc; i++) {VecGetSubVector(y,bA->isglobal.col[i],&by[i]);}
287:   for (j=0; j<nc; j++) {
288:     VecZeroEntries(by[j]);
289:     for (i=0; i<nr; i++) {
290:       if (!bA->m[i][j]) continue;
291:       /* y[j] <- y[j] + (A[i][j])^T * x[i] */
292:       MatMultTransposeAdd(bA->m[i][j],bx[i],by[j],by[j]);
293:     }
294:   }
295:   for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
296:   for (i=0; i<nc; i++) {VecRestoreSubVector(y,bA->isglobal.col[i],&by[i]);}
297:   return(0);
298: }

300: static PetscErrorCode MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z)
301: {
302:   Mat_Nest       *bA = (Mat_Nest*)A->data;
303:   Vec            *bx = bA->left,*bz = bA->right;
304:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

308:   for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
309:   for (i=0; i<nc; i++) {VecGetSubVector(z,bA->isglobal.col[i],&bz[i]);}
310:   for (j=0; j<nc; j++) {
311:     if (y != z) {
312:       Vec by;
313:       VecGetSubVector(y,bA->isglobal.col[j],&by);
314:       VecCopy(by,bz[j]);
315:       VecRestoreSubVector(y,bA->isglobal.col[j],&by);
316:     }
317:     for (i=0; i<nr; i++) {
318:       if (!bA->m[i][j]) continue;
319:       /* z[j] <- y[j] + (A[i][j])^T * x[i] */
320:       MatMultTransposeAdd(bA->m[i][j],bx[i],bz[j],bz[j]);
321:     }
322:   }
323:   for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
324:   for (i=0; i<nc; i++) {VecRestoreSubVector(z,bA->isglobal.col[i],&bz[i]);}
325:   return(0);
326: }

328: static PetscErrorCode MatTranspose_Nest(Mat A,MatReuse reuse,Mat *B)
329: {
330:   Mat_Nest       *bA = (Mat_Nest*)A->data, *bC;
331:   Mat            C;
332:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

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

338:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
339:     Mat *subs;
340:     IS  *is_row,*is_col;

342:     PetscCalloc1(nr * nc,&subs);
343:     PetscMalloc2(nr,&is_row,nc,&is_col);
344:     MatNestGetISs(A,is_row,is_col);
345:     if (reuse == MAT_INPLACE_MATRIX) {
346:       for (i=0; i<nr; i++) {
347:         for (j=0; j<nc; j++) {
348:           subs[i + nr * j] = bA->m[i][j];
349:         }
350:       }
351:     }

353:     MatCreateNest(PetscObjectComm((PetscObject)A),nc,is_col,nr,is_row,subs,&C);
354:     PetscFree(subs);
355:     PetscFree2(is_row,is_col);
356:   } else {
357:     C = *B;
358:   }

360:   bC = (Mat_Nest*)C->data;
361:   for (i=0; i<nr; i++) {
362:     for (j=0; j<nc; j++) {
363:       if (bA->m[i][j]) {
364:         MatTranspose(bA->m[i][j], reuse, &(bC->m[j][i]));
365:       } else {
366:         bC->m[j][i] = NULL;
367:       }
368:     }
369:   }

371:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
372:     *B = C;
373:   } else {
374:     MatHeaderMerge(A, &C);
375:   }
376:   return(0);
377: }

379: static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list)
380: {
382:   IS             *lst = *list;
383:   PetscInt       i;

386:   if (!lst) return(0);
387:   for (i=0; i<n; i++) if (lst[i]) {ISDestroy(&lst[i]);}
388:   PetscFree(lst);
389:   *list = NULL;
390:   return(0);
391: }

393: static PetscErrorCode MatReset_Nest(Mat A)
394: {
395:   Mat_Nest       *vs = (Mat_Nest*)A->data;
396:   PetscInt       i,j;

400:   /* release the matrices and the place holders */
401:   MatNestDestroyISList(vs->nr,&vs->isglobal.row);
402:   MatNestDestroyISList(vs->nc,&vs->isglobal.col);
403:   MatNestDestroyISList(vs->nr,&vs->islocal.row);
404:   MatNestDestroyISList(vs->nc,&vs->islocal.col);

406:   PetscFree(vs->row_len);
407:   PetscFree(vs->col_len);
408:   PetscFree(vs->nnzstate);

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

412:   /* release the matrices and the place holders */
413:   if (vs->m) {
414:     for (i=0; i<vs->nr; i++) {
415:       for (j=0; j<vs->nc; j++) {
416:         MatDestroy(&vs->m[i][j]);
417:       }
418:       PetscFree(vs->m[i]);
419:     }
420:     PetscFree(vs->m);
421:   }

423:   /* restore defaults */
424:   vs->nr = 0;
425:   vs->nc = 0;
426:   vs->splitassembly = PETSC_FALSE;
427:   return(0);
428: }

430: static PetscErrorCode MatDestroy_Nest(Mat A)
431: {

434:   MatReset_Nest(A);
435:   PetscFree(A->data);

437:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",0);
438:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",0);
439:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",0);
440:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",0);
441:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",0);
442:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",0);
443:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",0);
444:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",0);
445:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",0);
446:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",0);
447:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",0);
448:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",0);
449:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",NULL);
450:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",NULL);
451:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",NULL);
452:   return(0);
453: }

455: static PetscErrorCode MatMissingDiagonal_Nest(Mat mat,PetscBool *missing,PetscInt *dd)
456: {
457:   Mat_Nest       *vs = (Mat_Nest*)mat->data;
458:   PetscInt       i;

462:   if (dd) *dd = 0;
463:   if (!vs->nr) {
464:     *missing = PETSC_TRUE;
465:     return(0);
466:   }
467:   *missing = PETSC_FALSE;
468:   for (i = 0; i < vs->nr && !(*missing); i++) {
469:     *missing = PETSC_TRUE;
470:     if (vs->m[i][i]) {
471:       MatMissingDiagonal(vs->m[i][i],missing,NULL);
472:       if (*missing && dd) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"First missing entry not yet implemented");
473:     }
474:   }
475:   return(0);
476: }

478: static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)
479: {
480:   Mat_Nest       *vs = (Mat_Nest*)A->data;
481:   PetscInt       i,j;
483:   PetscBool      nnzstate = PETSC_FALSE;

486:   for (i=0; i<vs->nr; i++) {
487:     for (j=0; j<vs->nc; j++) {
488:       PetscObjectState subnnzstate = 0;
489:       if (vs->m[i][j]) {
490:         MatAssemblyBegin(vs->m[i][j],type);
491:         if (!vs->splitassembly) {
492:           /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
493:            * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
494:            * already performing an assembly, but the result would by more complicated and appears to offer less
495:            * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
496:            * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
497:            */
498:           MatAssemblyEnd(vs->m[i][j],type);
499:           MatGetNonzeroState(vs->m[i][j],&subnnzstate);
500:         }
501:       }
502:       nnzstate = (PetscBool)(nnzstate || vs->nnzstate[i*vs->nc+j] != subnnzstate);
503:       vs->nnzstate[i*vs->nc+j] = subnnzstate;
504:     }
505:   }
506:   if (nnzstate) A->nonzerostate++;
507:   return(0);
508: }

510: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
511: {
512:   Mat_Nest       *vs = (Mat_Nest*)A->data;
513:   PetscInt       i,j;

517:   for (i=0; i<vs->nr; i++) {
518:     for (j=0; j<vs->nc; j++) {
519:       if (vs->m[i][j]) {
520:         if (vs->splitassembly) {
521:           MatAssemblyEnd(vs->m[i][j],type);
522:         }
523:       }
524:     }
525:   }
526:   return(0);
527: }

529: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B)
530: {
532:   Mat_Nest       *vs = (Mat_Nest*)A->data;
533:   PetscInt       j;
534:   Mat            sub;

537:   sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
538:   for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j];
539:   if (sub) {MatSetUp(sub);}       /* Ensure that the sizes are available */
540:   *B = sub;
541:   return(0);
542: }

544: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B)
545: {
547:   Mat_Nest       *vs = (Mat_Nest*)A->data;
548:   PetscInt       i;
549:   Mat            sub;

552:   sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
553:   for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col];
554:   if (sub) {MatSetUp(sub);}       /* Ensure that the sizes are available */
555:   *B = sub;
556:   return(0);
557: }

559: static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found)
560: {
562:   PetscInt       i;
563:   PetscBool      flg;

569:   *found = -1;
570:   for (i=0; i<n; i++) {
571:     if (!list[i]) continue;
572:     ISEqualUnsorted(list[i],is,&flg);
573:     if (flg) {
574:       *found = i;
575:       return(0);
576:     }
577:   }
578:   SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set");
579:   return(0);
580: }

582: /* Get a block row as a new MatNest */
583: static PetscErrorCode MatNestGetRow(Mat A,PetscInt row,Mat *B)
584: {
585:   Mat_Nest       *vs = (Mat_Nest*)A->data;
586:   char           keyname[256];

590:   *B   = NULL;
591:   PetscSNPrintf(keyname,sizeof(keyname),"NestRow_%D",row);
592:   PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);
593:   if (*B) return(0);

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

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

599:   PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);
600:   PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
601:   return(0);
602: }

604: static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
605: {
606:   Mat_Nest       *vs = (Mat_Nest*)A->data;
608:   PetscInt       row,col;
609:   PetscBool      same,isFullCol,isFullColGlobal;

612:   /* Check if full column space. This is a hack */
613:   isFullCol = PETSC_FALSE;
614:   PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&same);
615:   if (same) {
616:     PetscInt n,first,step,i,an,am,afirst,astep;
617:     ISStrideGetInfo(iscol,&first,&step);
618:     ISGetLocalSize(iscol,&n);
619:     isFullCol = PETSC_TRUE;
620:     for (i=0,an=A->cmap->rstart; i<vs->nc; i++) {
621:       PetscObjectTypeCompare((PetscObject)is->col[i],ISSTRIDE,&same);
622:       ISGetLocalSize(is->col[i],&am);
623:       if (same) {
624:         ISStrideGetInfo(is->col[i],&afirst,&astep);
625:         if (afirst != an || astep != step) isFullCol = PETSC_FALSE;
626:       } else isFullCol = PETSC_FALSE;
627:       an += am;
628:     }
629:     if (an != A->cmap->rstart+n) isFullCol = PETSC_FALSE;
630:   }
631:   MPIU_Allreduce(&isFullCol,&isFullColGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)iscol));

633:   if (isFullColGlobal && vs->nc > 1) {
634:     PetscInt row;
635:     MatNestFindIS(A,vs->nr,is->row,isrow,&row);
636:     MatNestGetRow(A,row,B);
637:   } else {
638:     MatNestFindIS(A,vs->nr,is->row,isrow,&row);
639:     MatNestFindIS(A,vs->nc,is->col,iscol,&col);
640:     if (!vs->m[row][col]) {
641:       PetscInt lr,lc;

643:       MatCreate(PetscObjectComm((PetscObject)A),&vs->m[row][col]);
644:       ISGetLocalSize(vs->isglobal.row[row],&lr);
645:       ISGetLocalSize(vs->isglobal.col[col],&lc);
646:       MatSetSizes(vs->m[row][col],lr,lc,PETSC_DECIDE,PETSC_DECIDE);
647:       MatSetType(vs->m[row][col],MATAIJ);
648:       MatSeqAIJSetPreallocation(vs->m[row][col],0,NULL);
649:       MatMPIAIJSetPreallocation(vs->m[row][col],0,NULL,0,NULL);
650:       MatSetUp(vs->m[row][col]);
651:       MatAssemblyBegin(vs->m[row][col],MAT_FINAL_ASSEMBLY);
652:       MatAssemblyEnd(vs->m[row][col],MAT_FINAL_ASSEMBLY);
653:     }
654:     *B = vs->m[row][col];
655:   }
656:   return(0);
657: }

659: /*
660:    TODO: This does not actually returns a submatrix we can modify
661: */
662: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
663: {
665:   Mat_Nest       *vs = (Mat_Nest*)A->data;
666:   Mat            sub;

669:   MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
670:   switch (reuse) {
671:   case MAT_INITIAL_MATRIX:
672:     if (sub) { PetscObjectReference((PetscObject)sub); }
673:     *B = sub;
674:     break;
675:   case MAT_REUSE_MATRIX:
676:     if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call");
677:     break;
678:   case MAT_IGNORE_MATRIX:       /* Nothing to do */
679:     break;
680:   case MAT_INPLACE_MATRIX:       /* Nothing to do */
681:     SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet");
682:     break;
683:   }
684:   return(0);
685: }

687: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
688: {
690:   Mat_Nest       *vs = (Mat_Nest*)A->data;
691:   Mat            sub;

694:   MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
695:   /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
696:   if (sub) {PetscObjectReference((PetscObject)sub);}
697:   *B = sub;
698:   return(0);
699: }

701: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
702: {
704:   Mat_Nest       *vs = (Mat_Nest*)A->data;
705:   Mat            sub;

708:   MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
709:   if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten");
710:   if (sub) {
711:     if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times");
712:     MatDestroy(B);
713:   }
714:   return(0);
715: }

717: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
718: {
719:   Mat_Nest       *bA = (Mat_Nest*)A->data;
720:   PetscInt       i;

724:   for (i=0; i<bA->nr; i++) {
725:     Vec bv;
726:     VecGetSubVector(v,bA->isglobal.row[i],&bv);
727:     if (bA->m[i][i]) {
728:       MatGetDiagonal(bA->m[i][i],bv);
729:     } else {
730:       VecSet(bv,0.0);
731:     }
732:     VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
733:   }
734:   return(0);
735: }

737: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
738: {
739:   Mat_Nest       *bA = (Mat_Nest*)A->data;
740:   Vec            bl,*br;
741:   PetscInt       i,j;

745:   PetscCalloc1(bA->nc,&br);
746:   if (r) {
747:     for (j=0; j<bA->nc; j++) {VecGetSubVector(r,bA->isglobal.col[j],&br[j]);}
748:   }
749:   bl = NULL;
750:   for (i=0; i<bA->nr; i++) {
751:     if (l) {
752:       VecGetSubVector(l,bA->isglobal.row[i],&bl);
753:     }
754:     for (j=0; j<bA->nc; j++) {
755:       if (bA->m[i][j]) {
756:         MatDiagonalScale(bA->m[i][j],bl,br[j]);
757:       }
758:     }
759:     if (l) {
760:       VecRestoreSubVector(l,bA->isglobal.row[i],&bl);
761:     }
762:   }
763:   if (r) {
764:     for (j=0; j<bA->nc; j++) {VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);}
765:   }
766:   PetscFree(br);
767:   return(0);
768: }

770: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
771: {
772:   Mat_Nest       *bA = (Mat_Nest*)A->data;
773:   PetscInt       i,j;

777:   for (i=0; i<bA->nr; i++) {
778:     for (j=0; j<bA->nc; j++) {
779:       if (bA->m[i][j]) {
780:         MatScale(bA->m[i][j],a);
781:       }
782:     }
783:   }
784:   return(0);
785: }

787: static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
788: {
789:   Mat_Nest       *bA = (Mat_Nest*)A->data;
790:   PetscInt       i;
792:   PetscBool      nnzstate = PETSC_FALSE;

795:   for (i=0; i<bA->nr; i++) {
796:     PetscObjectState subnnzstate = 0;
797:     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);
798:     MatShift(bA->m[i][i],a);
799:     MatGetNonzeroState(bA->m[i][i],&subnnzstate);
800:     nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
801:     bA->nnzstate[i*bA->nc+i] = subnnzstate;
802:   }
803:   if (nnzstate) A->nonzerostate++;
804:   return(0);
805: }

807: static PetscErrorCode MatDiagonalSet_Nest(Mat A,Vec D,InsertMode is)
808: {
809:   Mat_Nest       *bA = (Mat_Nest*)A->data;
810:   PetscInt       i;
812:   PetscBool      nnzstate = PETSC_FALSE;

815:   for (i=0; i<bA->nr; i++) {
816:     PetscObjectState subnnzstate = 0;
817:     Vec              bv;
818:     VecGetSubVector(D,bA->isglobal.row[i],&bv);
819:     if (bA->m[i][i]) {
820:       MatDiagonalSet(bA->m[i][i],bv,is);
821:       MatGetNonzeroState(bA->m[i][i],&subnnzstate);
822:     }
823:     VecRestoreSubVector(D,bA->isglobal.row[i],&bv);
824:     nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
825:     bA->nnzstate[i*bA->nc+i] = subnnzstate;
826:   }
827:   if (nnzstate) A->nonzerostate++;
828:   return(0);
829: }

831: static PetscErrorCode MatSetRandom_Nest(Mat A,PetscRandom rctx)
832: {
833:   Mat_Nest       *bA = (Mat_Nest*)A->data;
834:   PetscInt       i,j;

838:   for (i=0; i<bA->nr; i++) {
839:     for (j=0; j<bA->nc; j++) {
840:       if (bA->m[i][j]) {
841:         MatSetRandom(bA->m[i][j],rctx);
842:       }
843:     }
844:   }
845:   return(0);
846: }

848: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left)
849: {
850:   Mat_Nest       *bA = (Mat_Nest*)A->data;
851:   Vec            *L,*R;
852:   MPI_Comm       comm;
853:   PetscInt       i,j;

857:   PetscObjectGetComm((PetscObject)A,&comm);
858:   if (right) {
859:     /* allocate R */
860:     PetscMalloc1(bA->nc, &R);
861:     /* Create the right vectors */
862:     for (j=0; j<bA->nc; j++) {
863:       for (i=0; i<bA->nr; i++) {
864:         if (bA->m[i][j]) {
865:           MatCreateVecs(bA->m[i][j],&R[j],NULL);
866:           break;
867:         }
868:       }
869:       if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
870:     }
871:     VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
872:     /* hand back control to the nest vector */
873:     for (j=0; j<bA->nc; j++) {
874:       VecDestroy(&R[j]);
875:     }
876:     PetscFree(R);
877:   }

879:   if (left) {
880:     /* allocate L */
881:     PetscMalloc1(bA->nr, &L);
882:     /* Create the left vectors */
883:     for (i=0; i<bA->nr; i++) {
884:       for (j=0; j<bA->nc; j++) {
885:         if (bA->m[i][j]) {
886:           MatCreateVecs(bA->m[i][j],NULL,&L[i]);
887:           break;
888:         }
889:       }
890:       if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
891:     }

893:     VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
894:     for (i=0; i<bA->nr; i++) {
895:       VecDestroy(&L[i]);
896:     }

898:     PetscFree(L);
899:   }
900:   return(0);
901: }

903: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
904: {
905:   Mat_Nest       *bA = (Mat_Nest*)A->data;
906:   PetscBool      isascii,viewSub = PETSC_FALSE;
907:   PetscInt       i,j;

911:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
912:   if (isascii) {

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

919:     PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");
920:     for (i=0; i<bA->nr; i++) {
921:       for (j=0; j<bA->nc; j++) {
922:         MatType   type;
923:         char      name[256] = "",prefix[256] = "";
924:         PetscInt  NR,NC;
925:         PetscBool isNest = PETSC_FALSE;

927:         if (!bA->m[i][j]) {
928:           PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);
929:           continue;
930:         }
931:         MatGetSize(bA->m[i][j],&NR,&NC);
932:         MatGetType(bA->m[i][j], &type);
933:         if (((PetscObject)bA->m[i][j])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);}
934:         if (((PetscObject)bA->m[i][j])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);}
935:         PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);

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

939:         if (isNest || viewSub) {
940:           PetscViewerASCIIPushTab(viewer);  /* push1 */
941:           MatView(bA->m[i][j],viewer);
942:           PetscViewerASCIIPopTab(viewer);    /* pop1 */
943:         }
944:       }
945:     }
946:     PetscViewerASCIIPopTab(viewer);    /* pop0 */
947:   }
948:   return(0);
949: }

951: static PetscErrorCode MatZeroEntries_Nest(Mat A)
952: {
953:   Mat_Nest       *bA = (Mat_Nest*)A->data;
954:   PetscInt       i,j;

958:   for (i=0; i<bA->nr; i++) {
959:     for (j=0; j<bA->nc; j++) {
960:       if (!bA->m[i][j]) continue;
961:       MatZeroEntries(bA->m[i][j]);
962:     }
963:   }
964:   return(0);
965: }

967: static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str)
968: {
969:   Mat_Nest       *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data;
970:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;
972:   PetscBool      nnzstate = PETSC_FALSE;

975:   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);
976:   for (i=0; i<nr; i++) {
977:     for (j=0; j<nc; j++) {
978:       PetscObjectState subnnzstate = 0;
979:       if (bA->m[i][j] && bB->m[i][j]) {
980:         MatCopy(bA->m[i][j],bB->m[i][j],str);
981:       } 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);
982:       MatGetNonzeroState(bB->m[i][j],&subnnzstate);
983:       nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i*nc+j] != subnnzstate);
984:       bB->nnzstate[i*nc+j] = subnnzstate;
985:     }
986:   }
987:   if (nnzstate) B->nonzerostate++;
988:   return(0);
989: }

991: static PetscErrorCode MatAXPY_Nest(Mat Y,PetscScalar a,Mat X,MatStructure str)
992: {
993:   Mat_Nest       *bY = (Mat_Nest*)Y->data,*bX = (Mat_Nest*)X->data;
994:   PetscInt       i,j,nr = bY->nr,nc = bY->nc;
996:   PetscBool      nnzstate = PETSC_FALSE;

999:   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);
1000:   for (i=0; i<nr; i++) {
1001:     for (j=0; j<nc; j++) {
1002:       PetscObjectState subnnzstate = 0;
1003:       if (bY->m[i][j] && bX->m[i][j]) {
1004:         MatAXPY(bY->m[i][j],a,bX->m[i][j],str);
1005:       } else if (bX->m[i][j]) {
1006:         Mat M;

1008:         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);
1009:         MatDuplicate(bX->m[i][j],MAT_COPY_VALUES,&M);
1010:         MatNestSetSubMat(Y,i,j,M);
1011:         MatDestroy(&M);
1012:       }
1013:       if (bY->m[i][j]) { MatGetNonzeroState(bY->m[i][j],&subnnzstate); }
1014:       nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i*nc+j] != subnnzstate);
1015:       bY->nnzstate[i*nc+j] = subnnzstate;
1016:     }
1017:   }
1018:   if (nnzstate) Y->nonzerostate++;
1019:   return(0);
1020: }

1022: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
1023: {
1024:   Mat_Nest       *bA = (Mat_Nest*)A->data;
1025:   Mat            *b;
1026:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

1030:   PetscMalloc1(nr*nc,&b);
1031:   for (i=0; i<nr; i++) {
1032:     for (j=0; j<nc; j++) {
1033:       if (bA->m[i][j]) {
1034:         MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
1035:       } else {
1036:         b[i*nc+j] = NULL;
1037:       }
1038:     }
1039:   }
1040:   MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
1041:   /* Give the new MatNest exclusive ownership */
1042:   for (i=0; i<nr*nc; i++) {
1043:     MatDestroy(&b[i]);
1044:   }
1045:   PetscFree(b);

1047:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1048:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1049:   return(0);
1050: }

1052: /* nest api */
1053: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
1054: {
1055:   Mat_Nest *bA = (Mat_Nest*)A->data;

1058:   if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
1059:   if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
1060:   *mat = bA->m[idxm][jdxm];
1061:   return(0);
1062: }

1064: /*@
1065:  MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.

1067:  Not collective

1069:  Input Parameters:
1070: +   A  - nest matrix
1071: .   idxm - index of the matrix within the nest matrix
1072: -   jdxm - index of the matrix within the nest matrix

1074:  Output Parameter:
1075: .   sub - matrix at index idxm,jdxm within the nest matrix

1077:  Level: developer

1079: .seealso: MatNestGetSize(), MatNestGetSubMats(), MatCreateNest(), MATNEST, MatNestSetSubMat(),
1080:           MatNestGetLocalISs(), MatNestGetISs()
1081: @*/
1082: PetscErrorCode  MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
1083: {

1087:   PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
1088:   return(0);
1089: }

1091: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
1092: {
1093:   Mat_Nest       *bA = (Mat_Nest*)A->data;
1094:   PetscInt       m,n,M,N,mi,ni,Mi,Ni;

1098:   if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
1099:   if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
1100:   MatGetLocalSize(mat,&m,&n);
1101:   MatGetSize(mat,&M,&N);
1102:   ISGetLocalSize(bA->isglobal.row[idxm],&mi);
1103:   ISGetSize(bA->isglobal.row[idxm],&Mi);
1104:   ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
1105:   ISGetSize(bA->isglobal.col[jdxm],&Ni);
1106:   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);
1107:   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);

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

1112:   PetscObjectReference((PetscObject)mat);
1113:   MatDestroy(&bA->m[idxm][jdxm]);
1114:   bA->m[idxm][jdxm] = mat;
1115:   PetscObjectStateIncrease((PetscObject)A);
1116:   MatGetNonzeroState(mat,&bA->nnzstate[idxm*bA->nc+jdxm]);
1117:   A->nonzerostate++;
1118:   return(0);
1119: }

1121: /*@
1122:  MatNestSetSubMat - Set a single submatrix in the nest matrix.

1124:  Logically collective on the submatrix communicator

1126:  Input Parameters:
1127: +   A  - nest matrix
1128: .   idxm - index of the matrix within the nest matrix
1129: .   jdxm - index of the matrix within the nest matrix
1130: -   sub - matrix at index idxm,jdxm within the nest matrix

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

1135:  This increments the reference count of the submatrix.

1137:  Level: developer

1139: .seealso: MatNestSetSubMats(), MatNestGetSubMats(), MatNestGetLocalISs(), MATNEST, MatCreateNest(),
1140:           MatNestGetSubMat(), MatNestGetISs(), MatNestGetSize()
1141: @*/
1142: PetscErrorCode  MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
1143: {

1147:   PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
1148:   return(0);
1149: }

1151: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1152: {
1153:   Mat_Nest *bA = (Mat_Nest*)A->data;

1156:   if (M)   *M   = bA->nr;
1157:   if (N)   *N   = bA->nc;
1158:   if (mat) *mat = bA->m;
1159:   return(0);
1160: }

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

1165:  Not collective

1167:  Input Parameters:
1168: .   A  - nest matrix

1170:  Output Parameter:
1171: +   M - number of rows in the nest matrix
1172: .   N - number of cols in the nest matrix
1173: -   mat - 2d array of matrices

1175:  Notes:

1177:  The user should not free the array mat.

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

1183:  Level: developer

1185: .seealso: MatNestGetSize(), MatNestGetSubMat(), MatNestGetLocalISs(), MATNEST, MatCreateNest(),
1186:           MatNestSetSubMats(), MatNestGetISs(), MatNestSetSubMat()
1187: @*/
1188: PetscErrorCode  MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1189: {

1193:   PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
1194:   return(0);
1195: }

1197: PetscErrorCode  MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
1198: {
1199:   Mat_Nest *bA = (Mat_Nest*)A->data;

1202:   if (M) *M = bA->nr;
1203:   if (N) *N = bA->nc;
1204:   return(0);
1205: }

1207: /*@
1208:  MatNestGetSize - Returns the size of the nest matrix.

1210:  Not collective

1212:  Input Parameters:
1213: .   A  - nest matrix

1215:  Output Parameter:
1216: +   M - number of rows in the nested mat
1217: -   N - number of cols in the nested mat

1219:  Notes:

1221:  Level: developer

1223: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MATNEST, MatCreateNest(), MatNestGetLocalISs(),
1224:           MatNestGetISs()
1225: @*/
1226: PetscErrorCode  MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
1227: {

1231:   PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
1232:   return(0);
1233: }

1235: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
1236: {
1237:   Mat_Nest *vs = (Mat_Nest*)A->data;
1238:   PetscInt i;

1241:   if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
1242:   if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
1243:   return(0);
1244: }

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

1249:  Not collective

1251:  Input Parameters:
1252: .   A  - nest matrix

1254:  Output Parameter:
1255: +   rows - array of row index sets
1256: -   cols - array of column index sets

1258:  Level: advanced

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

1263: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs(), MATNEST,
1264:           MatCreateNest(), MatNestGetSubMats(), MatNestSetSubMats()
1265: @*/
1266: PetscErrorCode  MatNestGetISs(Mat A,IS rows[],IS cols[])
1267: {

1272:   PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1273:   return(0);
1274: }

1276: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
1277: {
1278:   Mat_Nest *vs = (Mat_Nest*)A->data;
1279:   PetscInt i;

1282:   if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
1283:   if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
1284:   return(0);
1285: }

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

1290:  Not collective

1292:  Input Parameters:
1293: .   A  - nest matrix

1295:  Output Parameter:
1296: +   rows - array of row index sets (or NULL to ignore)
1297: -   cols - array of column index sets (or NULL to ignore)

1299:  Level: advanced

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

1304: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs(), MatCreateNest(),
1305:           MATNEST, MatNestSetSubMats(), MatNestSetSubMat()
1306: @*/
1307: PetscErrorCode  MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
1308: {

1313:   PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1314:   return(0);
1315: }

1317: PetscErrorCode  MatNestSetVecType_Nest(Mat A,VecType vtype)
1318: {
1320:   PetscBool      flg;

1323:   PetscStrcmp(vtype,VECNEST,&flg);
1324:   /* In reality, this only distinguishes VECNEST and "other" */
1325:   if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1326:   else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
1327:   return(0);
1328: }

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

1333:  Not collective

1335:  Input Parameters:
1336: +  A  - nest matrix
1337: -  vtype - type to use for creating vectors

1339:  Notes:

1341:  Level: developer

1343: .seealso: MatCreateVecs(), MATNEST, MatCreateNest()
1344: @*/
1345: PetscErrorCode  MatNestSetVecType(Mat A,VecType vtype)
1346: {

1350:   PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
1351:   return(0);
1352: }

1354: PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1355: {
1356:   Mat_Nest       *s = (Mat_Nest*)A->data;
1357:   PetscInt       i,j,m,n,M,N;
1359:   PetscBool      cong;

1362:   MatReset_Nest(A);

1364:   s->nr = nr;
1365:   s->nc = nc;

1367:   /* Create space for submatrices */
1368:   PetscMalloc1(nr,&s->m);
1369:   for (i=0; i<nr; i++) {
1370:     PetscMalloc1(nc,&s->m[i]);
1371:   }
1372:   for (i=0; i<nr; i++) {
1373:     for (j=0; j<nc; j++) {
1374:       s->m[i][j] = a[i*nc+j];
1375:       if (a[i*nc+j]) {
1376:         PetscObjectReference((PetscObject)a[i*nc+j]);
1377:       }
1378:     }
1379:   }

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

1383:   PetscMalloc1(nr,&s->row_len);
1384:   PetscMalloc1(nc,&s->col_len);
1385:   for (i=0; i<nr; i++) s->row_len[i]=-1;
1386:   for (j=0; j<nc; j++) s->col_len[j]=-1;

1388:   PetscCalloc1(nr*nc,&s->nnzstate);
1389:   for (i=0; i<nr; i++) {
1390:     for (j=0; j<nc; j++) {
1391:       if (s->m[i][j]) {
1392:         MatGetNonzeroState(s->m[i][j],&s->nnzstate[i*nc+j]);
1393:       }
1394:     }
1395:   }

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

1399:   PetscLayoutSetSize(A->rmap,M);
1400:   PetscLayoutSetLocalSize(A->rmap,m);
1401:   PetscLayoutSetSize(A->cmap,N);
1402:   PetscLayoutSetLocalSize(A->cmap,n);

1404:   PetscLayoutSetUp(A->rmap);
1405:   PetscLayoutSetUp(A->cmap);

1407:   /* disable operations that are not supported for non-square matrices,
1408:      or matrices for which is_row != is_col  */
1409:   MatHasCongruentLayouts(A,&cong);
1410:   if (cong && nr != nc) cong = PETSC_FALSE;
1411:   if (cong) {
1412:     for (i = 0; cong && i < nr; i++) {
1413:       ISEqualUnsorted(s->isglobal.row[i],s->isglobal.col[i],&cong);
1414:     }
1415:   }
1416:   if (!cong) {
1417:     A->ops->missingdiagonal = NULL;
1418:     A->ops->getdiagonal     = NULL;
1419:     A->ops->shift           = NULL;
1420:     A->ops->diagonalset     = NULL;
1421:   }

1423:   PetscCalloc2(nr,&s->left,nc,&s->right);
1424:   PetscObjectStateIncrease((PetscObject)A);
1425:   A->nonzerostate++;
1426:   return(0);
1427: }

1429: /*@
1430:    MatNestSetSubMats - Sets the nested submatrices

1432:    Collective on Mat

1434:    Input Parameter:
1435: +  A - nested matrix
1436: .  nr - number of nested row blocks
1437: .  is_row - index sets for each nested row block, or NULL to make contiguous
1438: .  nc - number of nested column blocks
1439: .  is_col - index sets for each nested column block, or NULL to make contiguous
1440: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

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

1444:    Level: advanced

1446: .seealso: MatCreateNest(), MATNEST, MatNestSetSubMat(), MatNestGetSubMat(), MatNestGetSubMats()
1447: @*/
1448: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1449: {
1451:   PetscInt       i;

1455:   if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1456:   if (nr && is_row) {
1459:   }
1460:   if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1461:   if (nc && is_col) {
1464:   }
1466:   PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1467:   return(0);
1468: }

1470: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1471: {
1473:   PetscBool      flg;
1474:   PetscInt       i,j,m,mi,*ix;

1477:   *ltog = NULL;
1478:   for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1479:     if (islocal[i]) {
1480:       ISGetLocalSize(islocal[i],&mi);
1481:       flg  = PETSC_TRUE;      /* We found a non-trivial entry */
1482:     } else {
1483:       ISGetLocalSize(isglobal[i],&mi);
1484:     }
1485:     m += mi;
1486:   }
1487:   if (!flg) return(0);

1489:   PetscMalloc1(m,&ix);
1490:   for (i=0,m=0; i<n; i++) {
1491:     ISLocalToGlobalMapping smap = NULL;
1492:     Mat                    sub = NULL;
1493:     PetscSF                sf;
1494:     PetscLayout            map;
1495:     const PetscInt         *ix2;

1497:     if (!colflg) {
1498:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1499:     } else {
1500:       MatNestFindNonzeroSubMatCol(A,i,&sub);
1501:     }
1502:     if (sub) {
1503:       if (!colflg) {
1504:         MatGetLocalToGlobalMapping(sub,&smap,NULL);
1505:       } else {
1506:         MatGetLocalToGlobalMapping(sub,NULL,&smap);
1507:       }
1508:     }
1509:     /*
1510:        Now we need to extract the monolithic global indices that correspond to the given split global indices.
1511:        In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1512:     */
1513:     ISGetIndices(isglobal[i],&ix2);
1514:     if (islocal[i]) {
1515:       PetscInt *ilocal,*iremote;
1516:       PetscInt mil,nleaves;

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

1523:       /* PetscSFSetGraphLayout does not like negative indices */
1524:       PetscMalloc2(mi,&ilocal,mi,&iremote);
1525:       for (j=0, nleaves = 0; j<mi; j++) {
1526:         if (ix[m+j] < 0) continue;
1527:         ilocal[nleaves]  = j;
1528:         iremote[nleaves] = ix[m+j];
1529:         nleaves++;
1530:       }
1531:       ISGetLocalSize(isglobal[i],&mil);
1532:       PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
1533:       PetscLayoutCreate(PetscObjectComm((PetscObject)A),&map);
1534:       PetscLayoutSetLocalSize(map,mil);
1535:       PetscLayoutSetUp(map);
1536:       PetscSFSetGraphLayout(sf,map,nleaves,ilocal,PETSC_USE_POINTER,iremote);
1537:       PetscLayoutDestroy(&map);
1538:       PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m);
1539:       PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m);
1540:       PetscSFDestroy(&sf);
1541:       PetscFree2(ilocal,iremote);
1542:     } else {
1543:       ISGetLocalSize(isglobal[i],&mi);
1544:       for (j=0; j<mi; j++) ix[m+j] = ix2[i];
1545:     }
1546:     ISRestoreIndices(isglobal[i],&ix2);
1547:     m   += mi;
1548:   }
1549:   ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);
1550:   return(0);
1551: }


1554: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1555: /*
1556:   nprocessors = NP
1557:   Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1558:        proc 0: => (g_0,h_0,)
1559:        proc 1: => (g_1,h_1,)
1560:        ...
1561:        proc nprocs-1: => (g_NP-1,h_NP-1,)

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

1566:             proc 0:
1567:     is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1568:             proc 1:
1569:     is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)

1571:             proc NP-1:
1572:     is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1573: */
1574: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1575: {
1576:   Mat_Nest       *vs = (Mat_Nest*)A->data;
1577:   PetscInt       i,j,offset,n,nsum,bs;
1579:   Mat            sub = NULL;

1582:   PetscMalloc1(nr,&vs->isglobal.row);
1583:   PetscMalloc1(nc,&vs->isglobal.col);
1584:   if (is_row) { /* valid IS is passed in */
1585:     /* refs on is[] are incremeneted */
1586:     for (i=0; i<vs->nr; i++) {
1587:       PetscObjectReference((PetscObject)is_row[i]);

1589:       vs->isglobal.row[i] = is_row[i];
1590:     }
1591:   } else {                      /* Create the ISs by inspecting sizes of a submatrix in each row */
1592:     nsum = 0;
1593:     for (i=0; i<vs->nr; i++) {  /* Add up the local sizes to compute the aggregate offset */
1594:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1595:       if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i);
1596:       MatGetLocalSize(sub,&n,NULL);
1597:       if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1598:       nsum += n;
1599:     }
1600:     MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1601:     offset -= nsum;
1602:     for (i=0; i<vs->nr; i++) {
1603:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1604:       MatGetLocalSize(sub,&n,NULL);
1605:       MatGetBlockSizes(sub,&bs,NULL);
1606:       ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);
1607:       ISSetBlockSize(vs->isglobal.row[i],bs);
1608:       offset += n;
1609:     }
1610:   }

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

1617:       vs->isglobal.col[j] = is_col[j];
1618:     }
1619:   } else {                      /* Create the ISs by inspecting sizes of a submatrix in each column */
1620:     offset = A->cmap->rstart;
1621:     nsum   = 0;
1622:     for (j=0; j<vs->nc; j++) {
1623:       MatNestFindNonzeroSubMatCol(A,j,&sub);
1624:       if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i);
1625:       MatGetLocalSize(sub,NULL,&n);
1626:       if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1627:       nsum += n;
1628:     }
1629:     MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1630:     offset -= nsum;
1631:     for (j=0; j<vs->nc; j++) {
1632:       MatNestFindNonzeroSubMatCol(A,j,&sub);
1633:       MatGetLocalSize(sub,NULL,&n);
1634:       MatGetBlockSizes(sub,NULL,&bs);
1635:       ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);
1636:       ISSetBlockSize(vs->isglobal.col[j],bs);
1637:       offset += n;
1638:     }
1639:   }

1641:   /* Set up the local ISs */
1642:   PetscMalloc1(vs->nr,&vs->islocal.row);
1643:   PetscMalloc1(vs->nc,&vs->islocal.col);
1644:   for (i=0,offset=0; i<vs->nr; i++) {
1645:     IS                     isloc;
1646:     ISLocalToGlobalMapping rmap = NULL;
1647:     PetscInt               nlocal,bs;
1648:     MatNestFindNonzeroSubMatRow(A,i,&sub);
1649:     if (sub) {MatGetLocalToGlobalMapping(sub,&rmap,NULL);}
1650:     if (rmap) {
1651:       MatGetBlockSizes(sub,&bs,NULL);
1652:       ISLocalToGlobalMappingGetSize(rmap,&nlocal);
1653:       ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1654:       ISSetBlockSize(isloc,bs);
1655:     } else {
1656:       nlocal = 0;
1657:       isloc  = NULL;
1658:     }
1659:     vs->islocal.row[i] = isloc;
1660:     offset            += nlocal;
1661:   }
1662:   for (i=0,offset=0; i<vs->nc; i++) {
1663:     IS                     isloc;
1664:     ISLocalToGlobalMapping cmap = NULL;
1665:     PetscInt               nlocal,bs;
1666:     MatNestFindNonzeroSubMatCol(A,i,&sub);
1667:     if (sub) {MatGetLocalToGlobalMapping(sub,NULL,&cmap);}
1668:     if (cmap) {
1669:       MatGetBlockSizes(sub,NULL,&bs);
1670:       ISLocalToGlobalMappingGetSize(cmap,&nlocal);
1671:       ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1672:       ISSetBlockSize(isloc,bs);
1673:     } else {
1674:       nlocal = 0;
1675:       isloc  = NULL;
1676:     }
1677:     vs->islocal.col[i] = isloc;
1678:     offset            += nlocal;
1679:   }

1681:   /* Set up the aggregate ISLocalToGlobalMapping */
1682:   {
1683:     ISLocalToGlobalMapping rmap,cmap;
1684:     MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);
1685:     MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);
1686:     if (rmap && cmap) {MatSetLocalToGlobalMapping(A,rmap,cmap);}
1687:     ISLocalToGlobalMappingDestroy(&rmap);
1688:     ISLocalToGlobalMappingDestroy(&cmap);
1689:   }

1691:   if (PetscDefined(USE_DEBUG)) {
1692:     for (i=0; i<vs->nr; i++) {
1693:       for (j=0; j<vs->nc; j++) {
1694:         PetscInt m,n,M,N,mi,ni,Mi,Ni;
1695:         Mat      B = vs->m[i][j];
1696:         if (!B) continue;
1697:         MatGetSize(B,&M,&N);
1698:         MatGetLocalSize(B,&m,&n);
1699:         ISGetSize(vs->isglobal.row[i],&Mi);
1700:         ISGetSize(vs->isglobal.col[j],&Ni);
1701:         ISGetLocalSize(vs->isglobal.row[i],&mi);
1702:         ISGetLocalSize(vs->isglobal.col[j],&ni);
1703:         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);
1704:         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);
1705:       }
1706:     }
1707:   }

1709:   /* Set A->assembled if all non-null blocks are currently assembled */
1710:   for (i=0; i<vs->nr; i++) {
1711:     for (j=0; j<vs->nc; j++) {
1712:       if (vs->m[i][j] && !vs->m[i][j]->assembled) return(0);
1713:     }
1714:   }
1715:   A->assembled = PETSC_TRUE;
1716:   return(0);
1717: }

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

1722:    Collective on Mat

1724:    Input Parameter:
1725: +  comm - Communicator for the new Mat
1726: .  nr - number of nested row blocks
1727: .  is_row - index sets for each nested row block, or NULL to make contiguous
1728: .  nc - number of nested column blocks
1729: .  is_col - index sets for each nested column block, or NULL to make contiguous
1730: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

1732:    Output Parameter:
1733: .  B - new matrix

1735:    Level: advanced

1737: .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST, MatNestSetSubMat(),
1738:           MatNestGetSubMat(), MatNestGetLocalISs(), MatNestGetSize(),
1739:           MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
1740: @*/
1741: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1742: {
1743:   Mat            A;

1747:   *B   = 0;
1748:   MatCreate(comm,&A);
1749:   MatSetType(A,MATNEST);
1750:   A->preallocated = PETSC_TRUE;
1751:   MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1752:   *B   = A;
1753:   return(0);
1754: }

1756: static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1757: {
1758:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1759:   Mat            *trans;
1760:   PetscScalar    **avv;
1761:   PetscScalar    *vv;
1762:   PetscInt       **aii,**ajj;
1763:   PetscInt       *ii,*jj,*ci;
1764:   PetscInt       nr,nc,nnz,i,j;
1765:   PetscBool      done;

1769:   MatGetSize(A,&nr,&nc);
1770:   if (reuse == MAT_REUSE_MATRIX) {
1771:     PetscInt rnr;

1773:     MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);
1774:     if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ");
1775:     if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows");
1776:     MatSeqAIJGetArray(*newmat,&vv);
1777:   }
1778:   /* extract CSR for nested SeqAIJ matrices */
1779:   nnz  = 0;
1780:   PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);
1781:   for (i=0; i<nest->nr; ++i) {
1782:     for (j=0; j<nest->nc; ++j) {
1783:       Mat B = nest->m[i][j];
1784:       if (B) {
1785:         PetscScalar *naa;
1786:         PetscInt    *nii,*njj,nnr;
1787:         PetscBool   istrans;

1789:         PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1790:         if (istrans) {
1791:           Mat Bt;

1793:           MatTransposeGetMat(B,&Bt);
1794:           MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);
1795:           B    = trans[i*nest->nc+j];
1796:         }
1797:         MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);
1798:         if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ");
1799:         MatSeqAIJGetArray(B,&naa);
1800:         nnz += nii[nnr];

1802:         aii[i*nest->nc+j] = nii;
1803:         ajj[i*nest->nc+j] = njj;
1804:         avv[i*nest->nc+j] = naa;
1805:       }
1806:     }
1807:   }
1808:   if (reuse != MAT_REUSE_MATRIX) {
1809:     PetscMalloc1(nr+1,&ii);
1810:     PetscMalloc1(nnz,&jj);
1811:     PetscMalloc1(nnz,&vv);
1812:   } else {
1813:     if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros");
1814:   }

1816:   /* new row pointer */
1817:   PetscArrayzero(ii,nr+1);
1818:   for (i=0; i<nest->nr; ++i) {
1819:     PetscInt       ncr,rst;

1821:     ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1822:     ISGetLocalSize(nest->isglobal.row[i],&ncr);
1823:     for (j=0; j<nest->nc; ++j) {
1824:       if (aii[i*nest->nc+j]) {
1825:         PetscInt    *nii = aii[i*nest->nc+j];
1826:         PetscInt    ir;

1828:         for (ir=rst; ir<ncr+rst; ++ir) {
1829:           ii[ir+1] += nii[1]-nii[0];
1830:           nii++;
1831:         }
1832:       }
1833:     }
1834:   }
1835:   for (i=0; i<nr; i++) ii[i+1] += ii[i];

1837:   /* construct CSR for the new matrix */
1838:   PetscCalloc1(nr,&ci);
1839:   for (i=0; i<nest->nr; ++i) {
1840:     PetscInt       ncr,rst;

1842:     ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1843:     ISGetLocalSize(nest->isglobal.row[i],&ncr);
1844:     for (j=0; j<nest->nc; ++j) {
1845:       if (aii[i*nest->nc+j]) {
1846:         PetscScalar *nvv = avv[i*nest->nc+j];
1847:         PetscInt    *nii = aii[i*nest->nc+j];
1848:         PetscInt    *njj = ajj[i*nest->nc+j];
1849:         PetscInt    ir,cst;

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

1855:           for (ij=0;ij<rsize;ij++) {
1856:             jj[ist+ij] = *njj+cst;
1857:             vv[ist+ij] = *nvv;
1858:             njj++;
1859:             nvv++;
1860:           }
1861:           ci[ir] += rsize;
1862:           nii++;
1863:         }
1864:       }
1865:     }
1866:   }
1867:   PetscFree(ci);

1869:   /* restore info */
1870:   for (i=0; i<nest->nr; ++i) {
1871:     for (j=0; j<nest->nc; ++j) {
1872:       Mat B = nest->m[i][j];
1873:       if (B) {
1874:         PetscInt nnr = 0, k = i*nest->nc+j;

1876:         B    = (trans[k] ? trans[k] : B);
1877:         MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);
1878:         if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ");
1879:         MatSeqAIJRestoreArray(B,&avv[k]);
1880:         MatDestroy(&trans[k]);
1881:       }
1882:     }
1883:   }
1884:   PetscFree4(aii,ajj,avv,trans);

1886:   /* finalize newmat */
1887:   if (reuse == MAT_INITIAL_MATRIX) {
1888:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);
1889:   } else if (reuse == MAT_INPLACE_MATRIX) {
1890:     Mat B;

1892:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);
1893:     MatHeaderReplace(A,&B);
1894:   }
1895:   MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1896:   MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1897:   {
1898:     Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data);
1899:     a->free_a     = PETSC_TRUE;
1900:     a->free_ij    = PETSC_TRUE;
1901:   }
1902:   return(0);
1903: }

1905: PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1906: {
1908:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1909:   PetscInt       m,n,M,N,i,j,k,*dnnz,*onnz,rstart;
1910:   PetscInt       cstart,cend;
1911:   PetscMPIInt    size;
1912:   Mat            C;

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

1920:     PetscStrcmp(newtype,MATAIJ,&fast);
1921:     if (!fast) {
1922:       PetscStrcmp(newtype,MATSEQAIJ,&fast);
1923:     }
1924:     for (i=0; i<nest->nr && fast; ++i) {
1925:       for (j=0; j<nest->nc && fast; ++j) {
1926:         Mat B = nest->m[i][j];
1927:         if (B) {
1928:           PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);
1929:           if (!fast) {
1930:             PetscBool istrans;

1932:             PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1933:             if (istrans) {
1934:               Mat Bt;

1936:               MatTransposeGetMat(B,&Bt);
1937:               PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);
1938:             }
1939:           }
1940:         }
1941:       }
1942:     }
1943:     for (i=0, nf=0; i<nest->nr && fast; ++i) {
1944:       PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);
1945:       if (fast) {
1946:         PetscInt f,s;

1948:         ISStrideGetInfo(nest->isglobal.row[i],&f,&s);
1949:         if (f != nf || s != 1) { fast = PETSC_FALSE; }
1950:         else {
1951:           ISGetSize(nest->isglobal.row[i],&f);
1952:           nf  += f;
1953:         }
1954:       }
1955:     }
1956:     for (i=0, nf=0; i<nest->nc && fast; ++i) {
1957:       PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);
1958:       if (fast) {
1959:         PetscInt f,s;

1961:         ISStrideGetInfo(nest->isglobal.col[i],&f,&s);
1962:         if (f != nf || s != 1) { fast = PETSC_FALSE; }
1963:         else {
1964:           ISGetSize(nest->isglobal.col[i],&f);
1965:           nf  += f;
1966:         }
1967:       }
1968:     }
1969:     if (fast) {
1970:       MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);
1971:       return(0);
1972:     }
1973:   }
1974:   MatGetSize(A,&M,&N);
1975:   MatGetLocalSize(A,&m,&n);
1976:   MatGetOwnershipRangeColumn(A,&cstart,&cend);
1977:   switch (reuse) {
1978:   case MAT_INITIAL_MATRIX:
1979:     MatCreate(PetscObjectComm((PetscObject)A),&C);
1980:     MatSetType(C,newtype);
1981:     MatSetSizes(C,m,n,M,N);
1982:     *newmat = C;
1983:     break;
1984:   case MAT_REUSE_MATRIX:
1985:     C = *newmat;
1986:     break;
1987:   default: SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatReuse");
1988:   }
1989:   PetscMalloc1(2*m,&dnnz);
1990:   onnz = dnnz + m;
1991:   for (k=0; k<m; k++) {
1992:     dnnz[k] = 0;
1993:     onnz[k] = 0;
1994:   }
1995:   for (j=0; j<nest->nc; ++j) {
1996:     IS             bNis;
1997:     PetscInt       bN;
1998:     const PetscInt *bNindices;
1999:     /* Using global column indices and ISAllGather() is not scalable. */
2000:     ISAllGather(nest->isglobal.col[j], &bNis);
2001:     ISGetSize(bNis, &bN);
2002:     ISGetIndices(bNis,&bNindices);
2003:     for (i=0; i<nest->nr; ++i) {
2004:       PetscSF        bmsf;
2005:       PetscSFNode    *iremote;
2006:       Mat            B;
2007:       PetscInt       bm, *sub_dnnz,*sub_onnz, br;
2008:       const PetscInt *bmindices;
2009:       B = nest->m[i][j];
2010:       if (!B) continue;
2011:       ISGetLocalSize(nest->isglobal.row[i],&bm);
2012:       ISGetIndices(nest->isglobal.row[i],&bmindices);
2013:       PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);
2014:       PetscMalloc1(bm,&iremote);
2015:       PetscMalloc1(bm,&sub_dnnz);
2016:       PetscMalloc1(bm,&sub_onnz);
2017:       for (k = 0; k < bm; ++k){
2018:             sub_dnnz[k] = 0;
2019:             sub_onnz[k] = 0;
2020:       }
2021:       /*
2022:        Locate the owners for all of the locally-owned global row indices for this row block.
2023:        These determine the roots of PetscSF used to communicate preallocation data to row owners.
2024:        The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
2025:        */
2026:       MatGetOwnershipRange(B,&rstart,NULL);
2027:       for (br = 0; br < bm; ++br) {
2028:         PetscInt       row = bmindices[br], brncols, col;
2029:         const PetscInt *brcols;
2030:         PetscInt       rowrel = 0; /* row's relative index on its owner rank */
2031:         PetscMPIInt    rowowner = 0;
2032:         PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);
2033:         /* how many roots  */
2034:         iremote[br].rank = rowowner; iremote[br].index = rowrel;           /* edge from bmdnnz to dnnz */
2035:         /* get nonzero pattern */
2036:         MatGetRow(B,br+rstart,&brncols,&brcols,NULL);
2037:         for (k=0; k<brncols; k++) {
2038:           col  = bNindices[brcols[k]];
2039:           if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) {
2040:             sub_dnnz[br]++;
2041:           } else {
2042:             sub_onnz[br]++;
2043:           }
2044:         }
2045:         MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);
2046:       }
2047:       ISRestoreIndices(nest->isglobal.row[i],&bmindices);
2048:       /* bsf will have to take care of disposing of bedges. */
2049:       PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
2050:       PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
2051:       PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
2052:       PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
2053:       PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
2054:       PetscFree(sub_dnnz);
2055:       PetscFree(sub_onnz);
2056:       PetscSFDestroy(&bmsf);
2057:     }
2058:     ISRestoreIndices(bNis,&bNindices);
2059:     ISDestroy(&bNis);
2060:   }
2061:   /* Resize preallocation if overestimated */
2062:   for (i=0;i<m;i++) {
2063:     dnnz[i] = PetscMin(dnnz[i],A->cmap->n);
2064:     onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n);
2065:   }
2066:   MatSeqAIJSetPreallocation(C,0,dnnz);
2067:   MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);
2068:   PetscFree(dnnz);

2070:   /* Fill by row */
2071:   for (j=0; j<nest->nc; ++j) {
2072:     /* Using global column indices and ISAllGather() is not scalable. */
2073:     IS             bNis;
2074:     PetscInt       bN;
2075:     const PetscInt *bNindices;
2076:     ISAllGather(nest->isglobal.col[j], &bNis);
2077:     ISGetSize(bNis,&bN);
2078:     ISGetIndices(bNis,&bNindices);
2079:     for (i=0; i<nest->nr; ++i) {
2080:       Mat            B;
2081:       PetscInt       bm, br;
2082:       const PetscInt *bmindices;
2083:       B = nest->m[i][j];
2084:       if (!B) continue;
2085:       ISGetLocalSize(nest->isglobal.row[i],&bm);
2086:       ISGetIndices(nest->isglobal.row[i],&bmindices);
2087:       MatGetOwnershipRange(B,&rstart,NULL);
2088:       for (br = 0; br < bm; ++br) {
2089:         PetscInt          row = bmindices[br], brncols,  *cols;
2090:         const PetscInt    *brcols;
2091:         const PetscScalar *brcoldata;
2092:         MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);
2093:         PetscMalloc1(brncols,&cols);
2094:         for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
2095:         /*
2096:           Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
2097:           Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
2098:          */
2099:         MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);
2100:         MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);
2101:         PetscFree(cols);
2102:       }
2103:       ISRestoreIndices(nest->isglobal.row[i],&bmindices);
2104:     }
2105:     ISRestoreIndices(bNis,&bNindices);
2106:     ISDestroy(&bNis);
2107:   }
2108:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
2109:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
2110:   return(0);
2111: }

2113: PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has)
2114: {
2115:   Mat_Nest       *bA = (Mat_Nest*)mat->data;
2116:   MatOperation   opAdd;
2117:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;
2118:   PetscBool      flg;

2122:   *has = PETSC_FALSE;
2123:   if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
2124:     opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD);
2125:     for (j=0; j<nc; j++) {
2126:       for (i=0; i<nr; i++) {
2127:         if (!bA->m[i][j]) continue;
2128:         MatHasOperation(bA->m[i][j],opAdd,&flg);
2129:         if (!flg) return(0);
2130:       }
2131:     }
2132:   }
2133:   if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
2134:   return(0);
2135: }

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

2140:   Level: intermediate

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

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

2151: .seealso: MatCreate(), MatType, MatCreateNest(), MatNestSetSubMat(), MatNestGetSubMat(),
2152:           VecCreateNest(), DMCreateMatrix(), DMCOMPOSITE, MatNestSetVecType(), MatNestGetLocalISs(),
2153:           MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
2154: M*/
2155: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2156: {
2157:   Mat_Nest       *s;

2161:   PetscNewLog(A,&s);
2162:   A->data = (void*)s;

2164:   s->nr            = -1;
2165:   s->nc            = -1;
2166:   s->m             = NULL;
2167:   s->splitassembly = PETSC_FALSE;

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

2171:   A->ops->mult                  = MatMult_Nest;
2172:   A->ops->multadd               = MatMultAdd_Nest;
2173:   A->ops->multtranspose         = MatMultTranspose_Nest;
2174:   A->ops->multtransposeadd      = MatMultTransposeAdd_Nest;
2175:   A->ops->transpose             = MatTranspose_Nest;
2176:   A->ops->assemblybegin         = MatAssemblyBegin_Nest;
2177:   A->ops->assemblyend           = MatAssemblyEnd_Nest;
2178:   A->ops->zeroentries           = MatZeroEntries_Nest;
2179:   A->ops->copy                  = MatCopy_Nest;
2180:   A->ops->axpy                  = MatAXPY_Nest;
2181:   A->ops->duplicate             = MatDuplicate_Nest;
2182:   A->ops->createsubmatrix       = MatCreateSubMatrix_Nest;
2183:   A->ops->destroy               = MatDestroy_Nest;
2184:   A->ops->view                  = MatView_Nest;
2185:   A->ops->getvecs               = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2186:   A->ops->getlocalsubmatrix     = MatGetLocalSubMatrix_Nest;
2187:   A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2188:   A->ops->getdiagonal           = MatGetDiagonal_Nest;
2189:   A->ops->diagonalscale         = MatDiagonalScale_Nest;
2190:   A->ops->scale                 = MatScale_Nest;
2191:   A->ops->shift                 = MatShift_Nest;
2192:   A->ops->diagonalset           = MatDiagonalSet_Nest;
2193:   A->ops->setrandom             = MatSetRandom_Nest;
2194:   A->ops->hasoperation          = MatHasOperation_Nest;
2195:   A->ops->missingdiagonal       = MatMissingDiagonal_Nest;

2197:   A->spptr        = 0;
2198:   A->assembled    = PETSC_FALSE;

2200:   /* expose Nest api's */
2201:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",        MatNestGetSubMat_Nest);
2202:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",        MatNestSetSubMat_Nest);
2203:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",       MatNestGetSubMats_Nest);
2204:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",          MatNestGetSize_Nest);
2205:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",           MatNestGetISs_Nest);
2206:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",      MatNestGetLocalISs_Nest);
2207:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",       MatNestSetVecType_Nest);
2208:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",       MatNestSetSubMats_Nest);
2209:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",  MatConvert_Nest_AIJ);
2210:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",  MatConvert_Nest_AIJ);
2211:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",     MatConvert_Nest_AIJ);
2212:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",      MatConvert_Nest_IS);
2213:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",MatProductSetFromOptions_Nest_Dense);
2214:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",MatProductSetFromOptions_Nest_Dense);
2215:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",MatProductSetFromOptions_Nest_Dense);

2217:   PetscObjectChangeTypeName((PetscObject)A,MATNEST);
2218:   return(0);
2219: }