Actual source code: matnest.c

petsc-3.13.4 2020-08-01
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 MatMatMultNumeric_Nest_Dense(Mat A,Mat B,Mat C)
 96: {
 97:   Mat_Nest          *bA = (Mat_Nest*)A->data;
 98:   PetscContainer    container;
 99:   Nest_Dense        *contents;
100:   Mat               viewB,viewC,seq,productB,workC;
101:   const PetscScalar *barray;
102:   PetscScalar       *carray;
103:   PetscInt          i,j,M,N,nr = bA->nr,nc = bA->nc,ldb,ldc;
104:   PetscErrorCode    ierr;

107:   PetscObjectQuery((PetscObject)C,"workC",(PetscObject*)&container);
108:   if (!container) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Container does not exist");
109:   PetscContainerGetPointer(container,(void**)&contents);
110:   MatDenseGetLDA(B,&ldb);
111:   MatDenseGetLDA(C,&ldc);
112:   MatGetSize(B,NULL,&N);
113:   MatZeroEntries(C);
114:   MatDenseGetArrayRead(B,&barray);
115:   MatDenseGetArray(C,&carray);
116:   for (i=0; i<nr; i++) {
117:     ISGetSize(bA->isglobal.row[i],&M);
118:     MatCreateDense(PetscObjectComm((PetscObject)A),contents->dm[i+1]-contents->dm[i],PETSC_DECIDE,M,N,carray+contents->dm[i],&viewC);
119:     MatDenseGetLocalMatrix(viewC,&seq);
120:     MatSeqDenseSetLDA(seq,ldc);
121:     for (j=0; j<nc; j++) {
122:       if (!bA->m[i][j]) continue;
123:       ISGetSize(bA->isglobal.col[j],&M);
124:       MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);
125:       MatDenseGetLocalMatrix(viewB,&seq);
126:       MatSeqDenseSetLDA(seq,ldb);

128:       /* MatMatMultNumeric(bA->m[i][j],viewB,contents->workC[i*nc + j]); */
129:       workC             = contents->workC[i*nc + j];
130:       productB          = workC->product->B;
131:       workC->product->B = viewB; /* use newly created dense matrix viewB */
132:       (*workC->ops->productnumeric)(workC);
133:       MatDestroy(&viewB);
134:       workC->product->B = productB; /* resume original B */

136:       /* C[i] <- workC + C[i] */
137:       MatAXPY(viewC,1.0,contents->workC[i*nc + j],SAME_NONZERO_PATTERN);
138:     }
139:     MatDestroy(&viewC);
140:   }
141:   MatDenseRestoreArray(C,&carray);
142:   MatDenseRestoreArrayRead(B,&barray);

144:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
145:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
146:   return(0);
147: }

149: PetscErrorCode MatNest_DenseDestroy(void *ctx)
150: {
151:   Nest_Dense     *contents = (Nest_Dense*)ctx;
152:   PetscInt       i;

156:   PetscFree(contents->tarray);
157:   for (i=0; i<contents->k; i++) {
158:     MatDestroy(contents->workC + i);
159:   }
160:   PetscFree3(contents->dm,contents->dn,contents->workC);
161:   PetscFree(contents);
162:   return(0);
163: }

165: PETSC_INTERN PetscErrorCode MatMatMultSymbolic_Nest_Dense(Mat A,Mat B,PetscReal fill,Mat C)
166: {
167:   Mat_Nest          *bA = (Mat_Nest*)A->data;
168:   Mat               viewB,viewSeq,workC;
169:   const PetscScalar *barray;
170:   PetscInt          i,j,M,N,m,nr = bA->nr,nc = bA->nc,maxm = 0,ldb;
171:   PetscContainer    container;
172:   Nest_Dense        *contents=NULL;
173:   PetscErrorCode    ierr;

176:   MatGetSize(B,NULL,&N);
177:   if (!C->assembled) {
178:     MatGetLocalSize(A,&m,NULL);
179:     MatGetSize(A,&M,NULL);

181:     MatSetSizes(C,m,PETSC_DECIDE,M,N);
182:     MatSetType(C,MATDENSE);
183:     MatSetUp(C);
184:   }

186:   PetscNew(&contents);
187:   PetscContainerCreate(PetscObjectComm((PetscObject)A),&container);
188:   PetscContainerSetPointer(container,contents);
189:   PetscContainerSetUserDestroy(container,MatNest_DenseDestroy);
190:   PetscObjectCompose((PetscObject)C,"workC",(PetscObject)container);
191:   PetscContainerDestroy(&container);
192:   PetscCalloc3(nr+1,&contents->dm,nc+1,&contents->dn,nr*nc,&contents->workC);
193:   contents->k = nr*nc;
194:   for (i=0; i<nr; i++) {
195:     ISGetLocalSize(bA->isglobal.row[i],contents->dm + i+1);
196:     maxm = PetscMax(maxm,contents->dm[i+1]);
197:     contents->dm[i+1] += contents->dm[i];
198:   }
199:   for (i=0; i<nc; i++) {
200:     ISGetLocalSize(bA->isglobal.col[i],contents->dn + i+1);
201:     contents->dn[i+1] += contents->dn[i];
202:   }
203:   PetscMalloc1(maxm*N,&contents->tarray);
204:   MatDenseGetLDA(B,&ldb);
205:   MatGetSize(B,NULL,&N);
206:   MatDenseGetArrayRead(B,&barray);
207:   /* loops are permuted compared to MatMatMultNumeric so that viewB is created only once per column of A */
208:   for (j=0; j<nc; j++) {
209:     ISGetSize(bA->isglobal.col[j],&M);
210:     MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);
211:     MatDenseGetLocalMatrix(viewB,&viewSeq);
212:     MatSeqDenseSetLDA(viewSeq,ldb);
213:     for (i=0; i<nr; i++) {
214:       if (!bA->m[i][j]) continue;
215:       /* MatMatMultSymbolic may attach a specific container (depending on MatType of bA->m[i][j]) to workC[i][j] */

217:       MatProductCreate(bA->m[i][j],viewB,NULL,&contents->workC[i*nc + j]);
218:       workC = contents->workC[i*nc + j];
219:       MatProductSetType(workC,MATPRODUCT_AB);
220:       MatProductSetAlgorithm(workC,"default");
221:       MatProductSetFill(workC,fill);
222:       MatProductSetFromOptions(workC);
223:       MatProductSymbolic(workC);

225:       MatDenseGetLocalMatrix(workC,&viewSeq);
226:       /* free the memory allocated in MatMatMultSymbolic, since tarray will be shared by all Mat */
227:       MatSeqDenseSetPreallocation(viewSeq,contents->tarray);
228:     }
229:     MatDestroy(&viewB);
230:   }
231:   MatDenseRestoreArrayRead(B,&barray);

233:   C->ops->matmultnumeric = MatMatMultNumeric_Nest_Dense;
234:   return(0);
235: }

237: /* --------------------------------------------------------- */
238: static PetscErrorCode MatProductSetFromOptions_Nest_Dense_AB(Mat C)
239: {
241:   C->ops->matmultsymbolic = MatMatMultSymbolic_Nest_Dense;
242:   C->ops->productsymbolic = MatProductSymbolic_AB;
243:   return(0);
244: }

246: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Nest_Dense(Mat C)
247: {
249:   Mat_Product    *product = C->product;

252:   if (product->type == MATPRODUCT_AB) {
253:     MatProductSetFromOptions_Nest_Dense_AB(C);
254:   } else SETERRQ1(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"MatProduct type %s is not supported for Nest and Dense matrices",MatProductTypes[product->type]);
255:   return(0);
256: }
257: /* --------------------------------------------------------- */

259: static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y)
260: {
261:   Mat_Nest       *bA = (Mat_Nest*)A->data;
262:   Vec            *bx = bA->left,*by = bA->right;
263:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

267:   for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
268:   for (i=0; i<nc; i++) {VecGetSubVector(y,bA->isglobal.col[i],&by[i]);}
269:   for (j=0; j<nc; j++) {
270:     VecZeroEntries(by[j]);
271:     for (i=0; i<nr; i++) {
272:       if (!bA->m[i][j]) continue;
273:       /* y[j] <- y[j] + (A[i][j])^T * x[i] */
274:       MatMultTransposeAdd(bA->m[i][j],bx[i],by[j],by[j]);
275:     }
276:   }
277:   for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
278:   for (i=0; i<nc; i++) {VecRestoreSubVector(y,bA->isglobal.col[i],&by[i]);}
279:   return(0);
280: }

282: static PetscErrorCode MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z)
283: {
284:   Mat_Nest       *bA = (Mat_Nest*)A->data;
285:   Vec            *bx = bA->left,*bz = bA->right;
286:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

290:   for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
291:   for (i=0; i<nc; i++) {VecGetSubVector(z,bA->isglobal.col[i],&bz[i]);}
292:   for (j=0; j<nc; j++) {
293:     if (y != z) {
294:       Vec by;
295:       VecGetSubVector(y,bA->isglobal.col[j],&by);
296:       VecCopy(by,bz[j]);
297:       VecRestoreSubVector(y,bA->isglobal.col[j],&by);
298:     }
299:     for (i=0; i<nr; i++) {
300:       if (!bA->m[i][j]) continue;
301:       /* z[j] <- y[j] + (A[i][j])^T * x[i] */
302:       MatMultTransposeAdd(bA->m[i][j],bx[i],bz[j],bz[j]);
303:     }
304:   }
305:   for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
306:   for (i=0; i<nc; i++) {VecRestoreSubVector(z,bA->isglobal.col[i],&bz[i]);}
307:   return(0);
308: }

310: static PetscErrorCode MatTranspose_Nest(Mat A,MatReuse reuse,Mat *B)
311: {
312:   Mat_Nest       *bA = (Mat_Nest*)A->data, *bC;
313:   Mat            C;
314:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

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

320:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
321:     Mat *subs;
322:     IS  *is_row,*is_col;

324:     PetscCalloc1(nr * nc,&subs);
325:     PetscMalloc2(nr,&is_row,nc,&is_col);
326:     MatNestGetISs(A,is_row,is_col);
327:     if (reuse == MAT_INPLACE_MATRIX) {
328:       for (i=0; i<nr; i++) {
329:         for (j=0; j<nc; j++) {
330:           subs[i + nr * j] = bA->m[i][j];
331:         }
332:       }
333:     }

335:     MatCreateNest(PetscObjectComm((PetscObject)A),nc,is_col,nr,is_row,subs,&C);
336:     PetscFree(subs);
337:     PetscFree2(is_row,is_col);
338:   } else {
339:     C = *B;
340:   }

342:   bC = (Mat_Nest*)C->data;
343:   for (i=0; i<nr; i++) {
344:     for (j=0; j<nc; j++) {
345:       if (bA->m[i][j]) {
346:         MatTranspose(bA->m[i][j], reuse, &(bC->m[j][i]));
347:       } else {
348:         bC->m[j][i] = NULL;
349:       }
350:     }
351:   }

353:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
354:     *B = C;
355:   } else {
356:     MatHeaderMerge(A, &C);
357:   }
358:   return(0);
359: }

361: static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list)
362: {
364:   IS             *lst = *list;
365:   PetscInt       i;

368:   if (!lst) return(0);
369:   for (i=0; i<n; i++) if (lst[i]) {ISDestroy(&lst[i]);}
370:   PetscFree(lst);
371:   *list = NULL;
372:   return(0);
373: }

375: static PetscErrorCode MatReset_Nest(Mat A)
376: {
377:   Mat_Nest       *vs = (Mat_Nest*)A->data;
378:   PetscInt       i,j;

382:   /* release the matrices and the place holders */
383:   MatNestDestroyISList(vs->nr,&vs->isglobal.row);
384:   MatNestDestroyISList(vs->nc,&vs->isglobal.col);
385:   MatNestDestroyISList(vs->nr,&vs->islocal.row);
386:   MatNestDestroyISList(vs->nc,&vs->islocal.col);

388:   PetscFree(vs->row_len);
389:   PetscFree(vs->col_len);
390:   PetscFree(vs->nnzstate);

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

394:   /* release the matrices and the place holders */
395:   if (vs->m) {
396:     for (i=0; i<vs->nr; i++) {
397:       for (j=0; j<vs->nc; j++) {
398:         MatDestroy(&vs->m[i][j]);
399:       }
400:       PetscFree(vs->m[i]);
401:     }
402:     PetscFree(vs->m);
403:   }

405:   /* restore defaults */
406:   vs->nr = 0;
407:   vs->nc = 0;
408:   vs->splitassembly = PETSC_FALSE;
409:   return(0);
410: }

412: static PetscErrorCode MatDestroy_Nest(Mat A)
413: {

416:   MatReset_Nest(A);
417:   PetscFree(A->data);

419:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",0);
420:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",0);
421:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",0);
422:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",0);
423:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",0);
424:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",0);
425:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",0);
426:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",0);
427:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",0);
428:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",0);
429:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",0);
430:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",0);
431:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",NULL);
432:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",NULL);
433:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",NULL);
434:   return(0);
435: }

437: static PetscErrorCode MatMissingDiagonal_Nest(Mat mat,PetscBool *missing,PetscInt *dd)
438: {
439:   Mat_Nest       *vs = (Mat_Nest*)mat->data;
440:   PetscInt       i;

444:   if (dd) *dd = 0;
445:   if (!vs->nr) {
446:     *missing = PETSC_TRUE;
447:     return(0);
448:   }
449:   *missing = PETSC_FALSE;
450:   for (i = 0; i < vs->nr && !(*missing); i++) {
451:     *missing = PETSC_TRUE;
452:     if (vs->m[i][i]) {
453:       MatMissingDiagonal(vs->m[i][i],missing,NULL);
454:       if (*missing && dd) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"First missing entry not yet implemented");
455:     }
456:   }
457:   return(0);
458: }

460: static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)
461: {
462:   Mat_Nest       *vs = (Mat_Nest*)A->data;
463:   PetscInt       i,j;
465:   PetscBool      nnzstate = PETSC_FALSE;

468:   for (i=0; i<vs->nr; i++) {
469:     for (j=0; j<vs->nc; j++) {
470:       PetscObjectState subnnzstate = 0;
471:       if (vs->m[i][j]) {
472:         MatAssemblyBegin(vs->m[i][j],type);
473:         if (!vs->splitassembly) {
474:           /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
475:            * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
476:            * already performing an assembly, but the result would by more complicated and appears to offer less
477:            * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
478:            * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
479:            */
480:           MatAssemblyEnd(vs->m[i][j],type);
481:           MatGetNonzeroState(vs->m[i][j],&subnnzstate);
482:         }
483:       }
484:       nnzstate = (PetscBool)(nnzstate || vs->nnzstate[i*vs->nc+j] != subnnzstate);
485:       vs->nnzstate[i*vs->nc+j] = subnnzstate;
486:     }
487:   }
488:   if (nnzstate) A->nonzerostate++;
489:   return(0);
490: }

492: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
493: {
494:   Mat_Nest       *vs = (Mat_Nest*)A->data;
495:   PetscInt       i,j;

499:   for (i=0; i<vs->nr; i++) {
500:     for (j=0; j<vs->nc; j++) {
501:       if (vs->m[i][j]) {
502:         if (vs->splitassembly) {
503:           MatAssemblyEnd(vs->m[i][j],type);
504:         }
505:       }
506:     }
507:   }
508:   return(0);
509: }

511: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B)
512: {
514:   Mat_Nest       *vs = (Mat_Nest*)A->data;
515:   PetscInt       j;
516:   Mat            sub;

519:   sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
520:   for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j];
521:   if (sub) {MatSetUp(sub);}       /* Ensure that the sizes are available */
522:   *B = sub;
523:   return(0);
524: }

526: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B)
527: {
529:   Mat_Nest       *vs = (Mat_Nest*)A->data;
530:   PetscInt       i;
531:   Mat            sub;

534:   sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
535:   for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col];
536:   if (sub) {MatSetUp(sub);}       /* Ensure that the sizes are available */
537:   *B = sub;
538:   return(0);
539: }

541: static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found)
542: {
544:   PetscInt       i;
545:   PetscBool      flg;

551:   *found = -1;
552:   for (i=0; i<n; i++) {
553:     if (!list[i]) continue;
554:     ISEqualUnsorted(list[i],is,&flg);
555:     if (flg) {
556:       *found = i;
557:       return(0);
558:     }
559:   }
560:   SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set");
561:   return(0);
562: }

564: /* Get a block row as a new MatNest */
565: static PetscErrorCode MatNestGetRow(Mat A,PetscInt row,Mat *B)
566: {
567:   Mat_Nest       *vs = (Mat_Nest*)A->data;
568:   char           keyname[256];

572:   *B   = NULL;
573:   PetscSNPrintf(keyname,sizeof(keyname),"NestRow_%D",row);
574:   PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);
575:   if (*B) return(0);

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

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

581:   PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);
582:   PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
583:   return(0);
584: }

586: static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
587: {
588:   Mat_Nest       *vs = (Mat_Nest*)A->data;
590:   PetscInt       row,col;
591:   PetscBool      same,isFullCol,isFullColGlobal;

594:   /* Check if full column space. This is a hack */
595:   isFullCol = PETSC_FALSE;
596:   PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&same);
597:   if (same) {
598:     PetscInt n,first,step,i,an,am,afirst,astep;
599:     ISStrideGetInfo(iscol,&first,&step);
600:     ISGetLocalSize(iscol,&n);
601:     isFullCol = PETSC_TRUE;
602:     for (i=0,an=A->cmap->rstart; i<vs->nc; i++) {
603:       PetscObjectTypeCompare((PetscObject)is->col[i],ISSTRIDE,&same);
604:       ISGetLocalSize(is->col[i],&am);
605:       if (same) {
606:         ISStrideGetInfo(is->col[i],&afirst,&astep);
607:         if (afirst != an || astep != step) isFullCol = PETSC_FALSE;
608:       } else isFullCol = PETSC_FALSE;
609:       an += am;
610:     }
611:     if (an != A->cmap->rstart+n) isFullCol = PETSC_FALSE;
612:   }
613:   MPIU_Allreduce(&isFullCol,&isFullColGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)iscol));

615:   if (isFullColGlobal && vs->nc > 1) {
616:     PetscInt row;
617:     MatNestFindIS(A,vs->nr,is->row,isrow,&row);
618:     MatNestGetRow(A,row,B);
619:   } else {
620:     MatNestFindIS(A,vs->nr,is->row,isrow,&row);
621:     MatNestFindIS(A,vs->nc,is->col,iscol,&col);
622:     if (!vs->m[row][col]) {
623:       PetscInt lr,lc;

625:       MatCreate(PetscObjectComm((PetscObject)A),&vs->m[row][col]);
626:       ISGetLocalSize(vs->isglobal.row[row],&lr);
627:       ISGetLocalSize(vs->isglobal.col[col],&lc);
628:       MatSetSizes(vs->m[row][col],lr,lc,PETSC_DECIDE,PETSC_DECIDE);
629:       MatSetType(vs->m[row][col],MATAIJ);
630:       MatSeqAIJSetPreallocation(vs->m[row][col],0,NULL);
631:       MatMPIAIJSetPreallocation(vs->m[row][col],0,NULL,0,NULL);
632:       MatSetUp(vs->m[row][col]);
633:       MatAssemblyBegin(vs->m[row][col],MAT_FINAL_ASSEMBLY);
634:       MatAssemblyEnd(vs->m[row][col],MAT_FINAL_ASSEMBLY);
635:     }
636:     *B = vs->m[row][col];
637:   }
638:   return(0);
639: }

641: /*
642:    TODO: This does not actually returns a submatrix we can modify
643: */
644: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
645: {
647:   Mat_Nest       *vs = (Mat_Nest*)A->data;
648:   Mat            sub;

651:   MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
652:   switch (reuse) {
653:   case MAT_INITIAL_MATRIX:
654:     if (sub) { PetscObjectReference((PetscObject)sub); }
655:     *B = sub;
656:     break;
657:   case MAT_REUSE_MATRIX:
658:     if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call");
659:     break;
660:   case MAT_IGNORE_MATRIX:       /* Nothing to do */
661:     break;
662:   case MAT_INPLACE_MATRIX:       /* Nothing to do */
663:     SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet");
664:     break;
665:   }
666:   return(0);
667: }

669: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
670: {
672:   Mat_Nest       *vs = (Mat_Nest*)A->data;
673:   Mat            sub;

676:   MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
677:   /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
678:   if (sub) {PetscObjectReference((PetscObject)sub);}
679:   *B = sub;
680:   return(0);
681: }

683: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
684: {
686:   Mat_Nest       *vs = (Mat_Nest*)A->data;
687:   Mat            sub;

690:   MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
691:   if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten");
692:   if (sub) {
693:     if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times");
694:     MatDestroy(B);
695:   }
696:   return(0);
697: }

699: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
700: {
701:   Mat_Nest       *bA = (Mat_Nest*)A->data;
702:   PetscInt       i;

706:   for (i=0; i<bA->nr; i++) {
707:     Vec bv;
708:     VecGetSubVector(v,bA->isglobal.row[i],&bv);
709:     if (bA->m[i][i]) {
710:       MatGetDiagonal(bA->m[i][i],bv);
711:     } else {
712:       VecSet(bv,0.0);
713:     }
714:     VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
715:   }
716:   return(0);
717: }

719: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
720: {
721:   Mat_Nest       *bA = (Mat_Nest*)A->data;
722:   Vec            bl,*br;
723:   PetscInt       i,j;

727:   PetscCalloc1(bA->nc,&br);
728:   if (r) {
729:     for (j=0; j<bA->nc; j++) {VecGetSubVector(r,bA->isglobal.col[j],&br[j]);}
730:   }
731:   bl = NULL;
732:   for (i=0; i<bA->nr; i++) {
733:     if (l) {
734:       VecGetSubVector(l,bA->isglobal.row[i],&bl);
735:     }
736:     for (j=0; j<bA->nc; j++) {
737:       if (bA->m[i][j]) {
738:         MatDiagonalScale(bA->m[i][j],bl,br[j]);
739:       }
740:     }
741:     if (l) {
742:       VecRestoreSubVector(l,bA->isglobal.row[i],&bl);
743:     }
744:   }
745:   if (r) {
746:     for (j=0; j<bA->nc; j++) {VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);}
747:   }
748:   PetscFree(br);
749:   return(0);
750: }

752: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
753: {
754:   Mat_Nest       *bA = (Mat_Nest*)A->data;
755:   PetscInt       i,j;

759:   for (i=0; i<bA->nr; i++) {
760:     for (j=0; j<bA->nc; j++) {
761:       if (bA->m[i][j]) {
762:         MatScale(bA->m[i][j],a);
763:       }
764:     }
765:   }
766:   return(0);
767: }

769: static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
770: {
771:   Mat_Nest       *bA = (Mat_Nest*)A->data;
772:   PetscInt       i;
774:   PetscBool      nnzstate = PETSC_FALSE;

777:   for (i=0; i<bA->nr; i++) {
778:     PetscObjectState subnnzstate = 0;
779:     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);
780:     MatShift(bA->m[i][i],a);
781:     MatGetNonzeroState(bA->m[i][i],&subnnzstate);
782:     nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
783:     bA->nnzstate[i*bA->nc+i] = subnnzstate;
784:   }
785:   if (nnzstate) A->nonzerostate++;
786:   return(0);
787: }

789: static PetscErrorCode MatDiagonalSet_Nest(Mat A,Vec D,InsertMode is)
790: {
791:   Mat_Nest       *bA = (Mat_Nest*)A->data;
792:   PetscInt       i;
794:   PetscBool      nnzstate = PETSC_FALSE;

797:   for (i=0; i<bA->nr; i++) {
798:     PetscObjectState subnnzstate = 0;
799:     Vec              bv;
800:     VecGetSubVector(D,bA->isglobal.row[i],&bv);
801:     if (bA->m[i][i]) {
802:       MatDiagonalSet(bA->m[i][i],bv,is);
803:       MatGetNonzeroState(bA->m[i][i],&subnnzstate);
804:     }
805:     VecRestoreSubVector(D,bA->isglobal.row[i],&bv);
806:     nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
807:     bA->nnzstate[i*bA->nc+i] = subnnzstate;
808:   }
809:   if (nnzstate) A->nonzerostate++;
810:   return(0);
811: }

813: static PetscErrorCode MatSetRandom_Nest(Mat A,PetscRandom rctx)
814: {
815:   Mat_Nest       *bA = (Mat_Nest*)A->data;
816:   PetscInt       i,j;

820:   for (i=0; i<bA->nr; i++) {
821:     for (j=0; j<bA->nc; j++) {
822:       if (bA->m[i][j]) {
823:         MatSetRandom(bA->m[i][j],rctx);
824:       }
825:     }
826:   }
827:   return(0);
828: }

830: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left)
831: {
832:   Mat_Nest       *bA = (Mat_Nest*)A->data;
833:   Vec            *L,*R;
834:   MPI_Comm       comm;
835:   PetscInt       i,j;

839:   PetscObjectGetComm((PetscObject)A,&comm);
840:   if (right) {
841:     /* allocate R */
842:     PetscMalloc1(bA->nc, &R);
843:     /* Create the right vectors */
844:     for (j=0; j<bA->nc; j++) {
845:       for (i=0; i<bA->nr; i++) {
846:         if (bA->m[i][j]) {
847:           MatCreateVecs(bA->m[i][j],&R[j],NULL);
848:           break;
849:         }
850:       }
851:       if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
852:     }
853:     VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
854:     /* hand back control to the nest vector */
855:     for (j=0; j<bA->nc; j++) {
856:       VecDestroy(&R[j]);
857:     }
858:     PetscFree(R);
859:   }

861:   if (left) {
862:     /* allocate L */
863:     PetscMalloc1(bA->nr, &L);
864:     /* Create the left vectors */
865:     for (i=0; i<bA->nr; i++) {
866:       for (j=0; j<bA->nc; j++) {
867:         if (bA->m[i][j]) {
868:           MatCreateVecs(bA->m[i][j],NULL,&L[i]);
869:           break;
870:         }
871:       }
872:       if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
873:     }

875:     VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
876:     for (i=0; i<bA->nr; i++) {
877:       VecDestroy(&L[i]);
878:     }

880:     PetscFree(L);
881:   }
882:   return(0);
883: }

885: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
886: {
887:   Mat_Nest       *bA = (Mat_Nest*)A->data;
888:   PetscBool      isascii,viewSub = PETSC_FALSE;
889:   PetscInt       i,j;

893:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
894:   if (isascii) {

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

901:     PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");
902:     for (i=0; i<bA->nr; i++) {
903:       for (j=0; j<bA->nc; j++) {
904:         MatType   type;
905:         char      name[256] = "",prefix[256] = "";
906:         PetscInt  NR,NC;
907:         PetscBool isNest = PETSC_FALSE;

909:         if (!bA->m[i][j]) {
910:           PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);
911:           continue;
912:         }
913:         MatGetSize(bA->m[i][j],&NR,&NC);
914:         MatGetType(bA->m[i][j], &type);
915:         if (((PetscObject)bA->m[i][j])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);}
916:         if (((PetscObject)bA->m[i][j])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);}
917:         PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);

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

921:         if (isNest || viewSub) {
922:           PetscViewerASCIIPushTab(viewer);  /* push1 */
923:           MatView(bA->m[i][j],viewer);
924:           PetscViewerASCIIPopTab(viewer);    /* pop1 */
925:         }
926:       }
927:     }
928:     PetscViewerASCIIPopTab(viewer);    /* pop0 */
929:   }
930:   return(0);
931: }

933: static PetscErrorCode MatZeroEntries_Nest(Mat A)
934: {
935:   Mat_Nest       *bA = (Mat_Nest*)A->data;
936:   PetscInt       i,j;

940:   for (i=0; i<bA->nr; i++) {
941:     for (j=0; j<bA->nc; j++) {
942:       if (!bA->m[i][j]) continue;
943:       MatZeroEntries(bA->m[i][j]);
944:     }
945:   }
946:   return(0);
947: }

949: static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str)
950: {
951:   Mat_Nest       *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data;
952:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;
954:   PetscBool      nnzstate = PETSC_FALSE;

957:   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);
958:   for (i=0; i<nr; i++) {
959:     for (j=0; j<nc; j++) {
960:       PetscObjectState subnnzstate = 0;
961:       if (bA->m[i][j] && bB->m[i][j]) {
962:         MatCopy(bA->m[i][j],bB->m[i][j],str);
963:       } 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);
964:       MatGetNonzeroState(bB->m[i][j],&subnnzstate);
965:       nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i*nc+j] != subnnzstate);
966:       bB->nnzstate[i*nc+j] = subnnzstate;
967:     }
968:   }
969:   if (nnzstate) B->nonzerostate++;
970:   return(0);
971: }

973: static PetscErrorCode MatAXPY_Nest(Mat Y,PetscScalar a,Mat X,MatStructure str)
974: {
975:   Mat_Nest       *bY = (Mat_Nest*)Y->data,*bX = (Mat_Nest*)X->data;
976:   PetscInt       i,j,nr = bY->nr,nc = bY->nc;
978:   PetscBool      nnzstate = PETSC_FALSE;

981:   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);
982:   for (i=0; i<nr; i++) {
983:     for (j=0; j<nc; j++) {
984:       PetscObjectState subnnzstate = 0;
985:       if (bY->m[i][j] && bX->m[i][j]) {
986:         MatAXPY(bY->m[i][j],a,bX->m[i][j],str);
987:       } else if (bX->m[i][j]) {
988:         Mat M;

990:         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);
991:         MatDuplicate(bX->m[i][j],MAT_COPY_VALUES,&M);
992:         MatNestSetSubMat(Y,i,j,M);
993:         MatDestroy(&M);
994:       }
995:       if (bY->m[i][j]) { MatGetNonzeroState(bY->m[i][j],&subnnzstate); }
996:       nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i*nc+j] != subnnzstate);
997:       bY->nnzstate[i*nc+j] = subnnzstate;
998:     }
999:   }
1000:   if (nnzstate) Y->nonzerostate++;
1001:   return(0);
1002: }

1004: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
1005: {
1006:   Mat_Nest       *bA = (Mat_Nest*)A->data;
1007:   Mat            *b;
1008:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

1012:   PetscMalloc1(nr*nc,&b);
1013:   for (i=0; i<nr; i++) {
1014:     for (j=0; j<nc; j++) {
1015:       if (bA->m[i][j]) {
1016:         MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
1017:       } else {
1018:         b[i*nc+j] = NULL;
1019:       }
1020:     }
1021:   }
1022:   MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
1023:   /* Give the new MatNest exclusive ownership */
1024:   for (i=0; i<nr*nc; i++) {
1025:     MatDestroy(&b[i]);
1026:   }
1027:   PetscFree(b);

1029:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1030:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1031:   return(0);
1032: }

1034: /* nest api */
1035: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
1036: {
1037:   Mat_Nest *bA = (Mat_Nest*)A->data;

1040:   if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
1041:   if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
1042:   *mat = bA->m[idxm][jdxm];
1043:   return(0);
1044: }

1046: /*@
1047:  MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.

1049:  Not collective

1051:  Input Parameters:
1052: +   A  - nest matrix
1053: .   idxm - index of the matrix within the nest matrix
1054: -   jdxm - index of the matrix within the nest matrix

1056:  Output Parameter:
1057: .   sub - matrix at index idxm,jdxm within the nest matrix

1059:  Level: developer

1061: .seealso: MatNestGetSize(), MatNestGetSubMats(), MatCreateNest(), MATNEST, MatNestSetSubMat(),
1062:           MatNestGetLocalISs(), MatNestGetISs()
1063: @*/
1064: PetscErrorCode  MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
1065: {

1069:   PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
1070:   return(0);
1071: }

1073: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
1074: {
1075:   Mat_Nest       *bA = (Mat_Nest*)A->data;
1076:   PetscInt       m,n,M,N,mi,ni,Mi,Ni;

1080:   if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
1081:   if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
1082:   MatGetLocalSize(mat,&m,&n);
1083:   MatGetSize(mat,&M,&N);
1084:   ISGetLocalSize(bA->isglobal.row[idxm],&mi);
1085:   ISGetSize(bA->isglobal.row[idxm],&Mi);
1086:   ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
1087:   ISGetSize(bA->isglobal.col[jdxm],&Ni);
1088:   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);
1089:   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);

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

1094:   PetscObjectReference((PetscObject)mat);
1095:   MatDestroy(&bA->m[idxm][jdxm]);
1096:   bA->m[idxm][jdxm] = mat;
1097:   PetscObjectStateIncrease((PetscObject)A);
1098:   MatGetNonzeroState(mat,&bA->nnzstate[idxm*bA->nc+jdxm]);
1099:   A->nonzerostate++;
1100:   return(0);
1101: }

1103: /*@
1104:  MatNestSetSubMat - Set a single submatrix in the nest matrix.

1106:  Logically collective on the submatrix communicator

1108:  Input Parameters:
1109: +   A  - nest matrix
1110: .   idxm - index of the matrix within the nest matrix
1111: .   jdxm - index of the matrix within the nest matrix
1112: -   sub - matrix at index idxm,jdxm within the nest matrix

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

1117:  This increments the reference count of the submatrix.

1119:  Level: developer

1121: .seealso: MatNestSetSubMats(), MatNestGetSubMats(), MatNestGetLocalISs(), MATNEST, MatCreateNest(),
1122:           MatNestGetSubMat(), MatNestGetISs(), MatNestGetSize()
1123: @*/
1124: PetscErrorCode  MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
1125: {

1129:   PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
1130:   return(0);
1131: }

1133: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1134: {
1135:   Mat_Nest *bA = (Mat_Nest*)A->data;

1138:   if (M)   *M   = bA->nr;
1139:   if (N)   *N   = bA->nc;
1140:   if (mat) *mat = bA->m;
1141:   return(0);
1142: }

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

1147:  Not collective

1149:  Input Parameters:
1150: .   A  - nest matrix

1152:  Output Parameter:
1153: +   M - number of rows in the nest matrix
1154: .   N - number of cols in the nest matrix
1155: -   mat - 2d array of matrices

1157:  Notes:

1159:  The user should not free the array mat.

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

1165:  Level: developer

1167: .seealso: MatNestGetSize(), MatNestGetSubMat(), MatNestGetLocalISs(), MATNEST, MatCreateNest(),
1168:           MatNestSetSubMats(), MatNestGetISs(), MatNestSetSubMat()
1169: @*/
1170: PetscErrorCode  MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1171: {

1175:   PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
1176:   return(0);
1177: }

1179: PetscErrorCode  MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
1180: {
1181:   Mat_Nest *bA = (Mat_Nest*)A->data;

1184:   if (M) *M = bA->nr;
1185:   if (N) *N = bA->nc;
1186:   return(0);
1187: }

1189: /*@
1190:  MatNestGetSize - Returns the size of the nest matrix.

1192:  Not collective

1194:  Input Parameters:
1195: .   A  - nest matrix

1197:  Output Parameter:
1198: +   M - number of rows in the nested mat
1199: -   N - number of cols in the nested mat

1201:  Notes:

1203:  Level: developer

1205: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MATNEST, MatCreateNest(), MatNestGetLocalISs(),
1206:           MatNestGetISs()
1207: @*/
1208: PetscErrorCode  MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
1209: {

1213:   PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
1214:   return(0);
1215: }

1217: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
1218: {
1219:   Mat_Nest *vs = (Mat_Nest*)A->data;
1220:   PetscInt i;

1223:   if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
1224:   if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
1225:   return(0);
1226: }

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

1231:  Not collective

1233:  Input Parameters:
1234: .   A  - nest matrix

1236:  Output Parameter:
1237: +   rows - array of row index sets
1238: -   cols - array of column index sets

1240:  Level: advanced

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

1245: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs(), MATNEST,
1246:           MatCreateNest(), MatNestGetSubMats(), MatNestSetSubMats()
1247: @*/
1248: PetscErrorCode  MatNestGetISs(Mat A,IS rows[],IS cols[])
1249: {

1254:   PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1255:   return(0);
1256: }

1258: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
1259: {
1260:   Mat_Nest *vs = (Mat_Nest*)A->data;
1261:   PetscInt i;

1264:   if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->islocal.row[i];
1265:   if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->islocal.col[i];
1266:   return(0);
1267: }

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

1272:  Not collective

1274:  Input Parameters:
1275: .   A  - nest matrix

1277:  Output Parameter:
1278: +   rows - array of row index sets (or NULL to ignore)
1279: -   cols - array of column index sets (or NULL to ignore)

1281:  Level: advanced

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

1286: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetISs(), MatCreateNest(),
1287:           MATNEST, MatNestSetSubMats(), MatNestSetSubMat()
1288: @*/
1289: PetscErrorCode  MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
1290: {

1295:   PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1296:   return(0);
1297: }

1299: PetscErrorCode  MatNestSetVecType_Nest(Mat A,VecType vtype)
1300: {
1302:   PetscBool      flg;

1305:   PetscStrcmp(vtype,VECNEST,&flg);
1306:   /* In reality, this only distinguishes VECNEST and "other" */
1307:   if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1308:   else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
1309:   return(0);
1310: }

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

1315:  Not collective

1317:  Input Parameters:
1318: +  A  - nest matrix
1319: -  vtype - type to use for creating vectors

1321:  Notes:

1323:  Level: developer

1325: .seealso: MatCreateVecs(), MATNEST, MatCreateNest()
1326: @*/
1327: PetscErrorCode  MatNestSetVecType(Mat A,VecType vtype)
1328: {

1332:   PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
1333:   return(0);
1334: }

1336: PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1337: {
1338:   Mat_Nest       *s = (Mat_Nest*)A->data;
1339:   PetscInt       i,j,m,n,M,N;
1341:   PetscBool      cong;

1344:   MatReset_Nest(A);

1346:   s->nr = nr;
1347:   s->nc = nc;

1349:   /* Create space for submatrices */
1350:   PetscMalloc1(nr,&s->m);
1351:   for (i=0; i<nr; i++) {
1352:     PetscMalloc1(nc,&s->m[i]);
1353:   }
1354:   for (i=0; i<nr; i++) {
1355:     for (j=0; j<nc; j++) {
1356:       s->m[i][j] = a[i*nc+j];
1357:       if (a[i*nc+j]) {
1358:         PetscObjectReference((PetscObject)a[i*nc+j]);
1359:       }
1360:     }
1361:   }

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

1365:   PetscMalloc1(nr,&s->row_len);
1366:   PetscMalloc1(nc,&s->col_len);
1367:   for (i=0; i<nr; i++) s->row_len[i]=-1;
1368:   for (j=0; j<nc; j++) s->col_len[j]=-1;

1370:   PetscCalloc1(nr*nc,&s->nnzstate);
1371:   for (i=0; i<nr; i++) {
1372:     for (j=0; j<nc; j++) {
1373:       if (s->m[i][j]) {
1374:         MatGetNonzeroState(s->m[i][j],&s->nnzstate[i*nc+j]);
1375:       }
1376:     }
1377:   }

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

1381:   PetscLayoutSetSize(A->rmap,M);
1382:   PetscLayoutSetLocalSize(A->rmap,m);
1383:   PetscLayoutSetSize(A->cmap,N);
1384:   PetscLayoutSetLocalSize(A->cmap,n);

1386:   PetscLayoutSetUp(A->rmap);
1387:   PetscLayoutSetUp(A->cmap);

1389:   /* disable operations that are not supported for non-square matrices,
1390:      or matrices for which is_row != is_col  */
1391:   MatHasCongruentLayouts(A,&cong);
1392:   if (cong && nr != nc) cong = PETSC_FALSE;
1393:   if (cong) {
1394:     for (i = 0; cong && i < nr; i++) {
1395:       ISEqualUnsorted(s->isglobal.row[i],s->isglobal.col[i],&cong);
1396:     }
1397:   }
1398:   if (!cong) {
1399:     A->ops->missingdiagonal = NULL;
1400:     A->ops->getdiagonal     = NULL;
1401:     A->ops->shift           = NULL;
1402:     A->ops->diagonalset     = NULL;
1403:   }

1405:   PetscCalloc2(nr,&s->left,nc,&s->right);
1406:   PetscObjectStateIncrease((PetscObject)A);
1407:   A->nonzerostate++;
1408:   return(0);
1409: }

1411: /*@
1412:    MatNestSetSubMats - Sets the nested submatrices

1414:    Collective on Mat

1416:    Input Parameter:
1417: +  A - nested matrix
1418: .  nr - number of nested row blocks
1419: .  is_row - index sets for each nested row block, or NULL to make contiguous
1420: .  nc - number of nested column blocks
1421: .  is_col - index sets for each nested column block, or NULL to make contiguous
1422: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

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

1426:    Level: advanced

1428: .seealso: MatCreateNest(), MATNEST, MatNestSetSubMat(), MatNestGetSubMat(), MatNestGetSubMats()
1429: @*/
1430: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1431: {
1433:   PetscInt       i;

1437:   if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1438:   if (nr && is_row) {
1441:   }
1442:   if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1443:   if (nc && is_col) {
1446:   }
1448:   PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1449:   return(0);
1450: }

1452: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1453: {
1455:   PetscBool      flg;
1456:   PetscInt       i,j,m,mi,*ix;

1459:   *ltog = NULL;
1460:   for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1461:     if (islocal[i]) {
1462:       ISGetLocalSize(islocal[i],&mi);
1463:       flg  = PETSC_TRUE;      /* We found a non-trivial entry */
1464:     } else {
1465:       ISGetLocalSize(isglobal[i],&mi);
1466:     }
1467:     m += mi;
1468:   }
1469:   if (!flg) return(0);

1471:   PetscMalloc1(m,&ix);
1472:   for (i=0,m=0; i<n; i++) {
1473:     ISLocalToGlobalMapping smap = NULL;
1474:     Mat                    sub = NULL;
1475:     PetscSF                sf;
1476:     PetscLayout            map;
1477:     const PetscInt         *ix2;

1479:     if (!colflg) {
1480:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1481:     } else {
1482:       MatNestFindNonzeroSubMatCol(A,i,&sub);
1483:     }
1484:     if (sub) {
1485:       if (!colflg) {
1486:         MatGetLocalToGlobalMapping(sub,&smap,NULL);
1487:       } else {
1488:         MatGetLocalToGlobalMapping(sub,NULL,&smap);
1489:       }
1490:     }
1491:     /*
1492:        Now we need to extract the monolithic global indices that correspond to the given split global indices.
1493:        In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1494:     */
1495:     ISGetIndices(isglobal[i],&ix2);
1496:     if (islocal[i]) {
1497:       PetscInt *ilocal,*iremote;
1498:       PetscInt mil,nleaves;

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

1505:       /* PetscSFSetGraphLayout does not like negative indices */
1506:       PetscMalloc2(mi,&ilocal,mi,&iremote);
1507:       for (j=0, nleaves = 0; j<mi; j++) {
1508:         if (ix[m+j] < 0) continue;
1509:         ilocal[nleaves]  = j;
1510:         iremote[nleaves] = ix[m+j];
1511:         nleaves++;
1512:       }
1513:       ISGetLocalSize(isglobal[i],&mil);
1514:       PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
1515:       PetscLayoutCreate(PetscObjectComm((PetscObject)A),&map);
1516:       PetscLayoutSetLocalSize(map,mil);
1517:       PetscLayoutSetUp(map);
1518:       PetscSFSetGraphLayout(sf,map,nleaves,ilocal,PETSC_USE_POINTER,iremote);
1519:       PetscLayoutDestroy(&map);
1520:       PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m);
1521:       PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m);
1522:       PetscSFDestroy(&sf);
1523:       PetscFree2(ilocal,iremote);
1524:     } else {
1525:       ISGetLocalSize(isglobal[i],&mi);
1526:       for (j=0; j<mi; j++) ix[m+j] = ix2[i];
1527:     }
1528:     ISRestoreIndices(isglobal[i],&ix2);
1529:     m   += mi;
1530:   }
1531:   ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);
1532:   return(0);
1533: }


1536: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1537: /*
1538:   nprocessors = NP
1539:   Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1540:        proc 0: => (g_0,h_0,)
1541:        proc 1: => (g_1,h_1,)
1542:        ...
1543:        proc nprocs-1: => (g_NP-1,h_NP-1,)

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

1548:             proc 0:
1549:     is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1550:             proc 1:
1551:     is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)

1553:             proc NP-1:
1554:     is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1555: */
1556: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1557: {
1558:   Mat_Nest       *vs = (Mat_Nest*)A->data;
1559:   PetscInt       i,j,offset,n,nsum,bs;
1561:   Mat            sub = NULL;

1564:   PetscMalloc1(nr,&vs->isglobal.row);
1565:   PetscMalloc1(nc,&vs->isglobal.col);
1566:   if (is_row) { /* valid IS is passed in */
1567:     /* refs on is[] are incremeneted */
1568:     for (i=0; i<vs->nr; i++) {
1569:       PetscObjectReference((PetscObject)is_row[i]);

1571:       vs->isglobal.row[i] = is_row[i];
1572:     }
1573:   } else {                      /* Create the ISs by inspecting sizes of a submatrix in each row */
1574:     nsum = 0;
1575:     for (i=0; i<vs->nr; i++) {  /* Add up the local sizes to compute the aggregate offset */
1576:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1577:       if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i);
1578:       MatGetLocalSize(sub,&n,NULL);
1579:       if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1580:       nsum += n;
1581:     }
1582:     MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1583:     offset -= nsum;
1584:     for (i=0; i<vs->nr; i++) {
1585:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1586:       MatGetLocalSize(sub,&n,NULL);
1587:       MatGetBlockSizes(sub,&bs,NULL);
1588:       ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);
1589:       ISSetBlockSize(vs->isglobal.row[i],bs);
1590:       offset += n;
1591:     }
1592:   }

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

1599:       vs->isglobal.col[j] = is_col[j];
1600:     }
1601:   } else {                      /* Create the ISs by inspecting sizes of a submatrix in each column */
1602:     offset = A->cmap->rstart;
1603:     nsum   = 0;
1604:     for (j=0; j<vs->nc; j++) {
1605:       MatNestFindNonzeroSubMatCol(A,j,&sub);
1606:       if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i);
1607:       MatGetLocalSize(sub,NULL,&n);
1608:       if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1609:       nsum += n;
1610:     }
1611:     MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1612:     offset -= nsum;
1613:     for (j=0; j<vs->nc; j++) {
1614:       MatNestFindNonzeroSubMatCol(A,j,&sub);
1615:       MatGetLocalSize(sub,NULL,&n);
1616:       MatGetBlockSizes(sub,NULL,&bs);
1617:       ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);
1618:       ISSetBlockSize(vs->isglobal.col[j],bs);
1619:       offset += n;
1620:     }
1621:   }

1623:   /* Set up the local ISs */
1624:   PetscMalloc1(vs->nr,&vs->islocal.row);
1625:   PetscMalloc1(vs->nc,&vs->islocal.col);
1626:   for (i=0,offset=0; i<vs->nr; i++) {
1627:     IS                     isloc;
1628:     ISLocalToGlobalMapping rmap = NULL;
1629:     PetscInt               nlocal,bs;
1630:     MatNestFindNonzeroSubMatRow(A,i,&sub);
1631:     if (sub) {MatGetLocalToGlobalMapping(sub,&rmap,NULL);}
1632:     if (rmap) {
1633:       MatGetBlockSizes(sub,&bs,NULL);
1634:       ISLocalToGlobalMappingGetSize(rmap,&nlocal);
1635:       ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1636:       ISSetBlockSize(isloc,bs);
1637:     } else {
1638:       nlocal = 0;
1639:       isloc  = NULL;
1640:     }
1641:     vs->islocal.row[i] = isloc;
1642:     offset            += nlocal;
1643:   }
1644:   for (i=0,offset=0; i<vs->nc; i++) {
1645:     IS                     isloc;
1646:     ISLocalToGlobalMapping cmap = NULL;
1647:     PetscInt               nlocal,bs;
1648:     MatNestFindNonzeroSubMatCol(A,i,&sub);
1649:     if (sub) {MatGetLocalToGlobalMapping(sub,NULL,&cmap);}
1650:     if (cmap) {
1651:       MatGetBlockSizes(sub,NULL,&bs);
1652:       ISLocalToGlobalMappingGetSize(cmap,&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.col[i] = isloc;
1660:     offset            += nlocal;
1661:   }

1663:   /* Set up the aggregate ISLocalToGlobalMapping */
1664:   {
1665:     ISLocalToGlobalMapping rmap,cmap;
1666:     MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);
1667:     MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);
1668:     if (rmap && cmap) {MatSetLocalToGlobalMapping(A,rmap,cmap);}
1669:     ISLocalToGlobalMappingDestroy(&rmap);
1670:     ISLocalToGlobalMappingDestroy(&cmap);
1671:   }

1673: #if defined(PETSC_USE_DEBUG)
1674:   for (i=0; i<vs->nr; i++) {
1675:     for (j=0; j<vs->nc; j++) {
1676:       PetscInt m,n,M,N,mi,ni,Mi,Ni;
1677:       Mat      B = vs->m[i][j];
1678:       if (!B) continue;
1679:       MatGetSize(B,&M,&N);
1680:       MatGetLocalSize(B,&m,&n);
1681:       ISGetSize(vs->isglobal.row[i],&Mi);
1682:       ISGetSize(vs->isglobal.col[j],&Ni);
1683:       ISGetLocalSize(vs->isglobal.row[i],&mi);
1684:       ISGetLocalSize(vs->isglobal.col[j],&ni);
1685:       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);
1686:       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);
1687:     }
1688:   }
1689: #endif

1691:   /* Set A->assembled if all non-null blocks are currently assembled */
1692:   for (i=0; i<vs->nr; i++) {
1693:     for (j=0; j<vs->nc; j++) {
1694:       if (vs->m[i][j] && !vs->m[i][j]->assembled) return(0);
1695:     }
1696:   }
1697:   A->assembled = PETSC_TRUE;
1698:   return(0);
1699: }

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

1704:    Collective on Mat

1706:    Input Parameter:
1707: +  comm - Communicator for the new Mat
1708: .  nr - number of nested row blocks
1709: .  is_row - index sets for each nested row block, or NULL to make contiguous
1710: .  nc - number of nested column blocks
1711: .  is_col - index sets for each nested column block, or NULL to make contiguous
1712: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

1714:    Output Parameter:
1715: .  B - new matrix

1717:    Level: advanced

1719: .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST, MatNestSetSubMat(),
1720:           MatNestGetSubMat(), MatNestGetLocalISs(), MatNestGetSize(),
1721:           MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
1722: @*/
1723: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1724: {
1725:   Mat            A;

1729:   *B   = 0;
1730:   MatCreate(comm,&A);
1731:   MatSetType(A,MATNEST);
1732:   A->preallocated = PETSC_TRUE;
1733:   MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1734:   *B   = A;
1735:   return(0);
1736: }

1738: static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1739: {
1740:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1741:   Mat            *trans;
1742:   PetscScalar    **avv;
1743:   PetscScalar    *vv;
1744:   PetscInt       **aii,**ajj;
1745:   PetscInt       *ii,*jj,*ci;
1746:   PetscInt       nr,nc,nnz,i,j;
1747:   PetscBool      done;

1751:   MatGetSize(A,&nr,&nc);
1752:   if (reuse == MAT_REUSE_MATRIX) {
1753:     PetscInt rnr;

1755:     MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);
1756:     if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ");
1757:     if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows");
1758:     MatSeqAIJGetArray(*newmat,&vv);
1759:   }
1760:   /* extract CSR for nested SeqAIJ matrices */
1761:   nnz  = 0;
1762:   PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);
1763:   for (i=0; i<nest->nr; ++i) {
1764:     for (j=0; j<nest->nc; ++j) {
1765:       Mat B = nest->m[i][j];
1766:       if (B) {
1767:         PetscScalar *naa;
1768:         PetscInt    *nii,*njj,nnr;
1769:         PetscBool   istrans;

1771:         PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1772:         if (istrans) {
1773:           Mat Bt;

1775:           MatTransposeGetMat(B,&Bt);
1776:           MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);
1777:           B    = trans[i*nest->nc+j];
1778:         }
1779:         MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);
1780:         if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ");
1781:         MatSeqAIJGetArray(B,&naa);
1782:         nnz += nii[nnr];

1784:         aii[i*nest->nc+j] = nii;
1785:         ajj[i*nest->nc+j] = njj;
1786:         avv[i*nest->nc+j] = naa;
1787:       }
1788:     }
1789:   }
1790:   if (reuse != MAT_REUSE_MATRIX) {
1791:     PetscMalloc1(nr+1,&ii);
1792:     PetscMalloc1(nnz,&jj);
1793:     PetscMalloc1(nnz,&vv);
1794:   } else {
1795:     if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros");
1796:   }

1798:   /* new row pointer */
1799:   PetscArrayzero(ii,nr+1);
1800:   for (i=0; i<nest->nr; ++i) {
1801:     PetscInt       ncr,rst;

1803:     ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1804:     ISGetLocalSize(nest->isglobal.row[i],&ncr);
1805:     for (j=0; j<nest->nc; ++j) {
1806:       if (aii[i*nest->nc+j]) {
1807:         PetscInt    *nii = aii[i*nest->nc+j];
1808:         PetscInt    ir;

1810:         for (ir=rst; ir<ncr+rst; ++ir) {
1811:           ii[ir+1] += nii[1]-nii[0];
1812:           nii++;
1813:         }
1814:       }
1815:     }
1816:   }
1817:   for (i=0; i<nr; i++) ii[i+1] += ii[i];

1819:   /* construct CSR for the new matrix */
1820:   PetscCalloc1(nr,&ci);
1821:   for (i=0; i<nest->nr; ++i) {
1822:     PetscInt       ncr,rst;

1824:     ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1825:     ISGetLocalSize(nest->isglobal.row[i],&ncr);
1826:     for (j=0; j<nest->nc; ++j) {
1827:       if (aii[i*nest->nc+j]) {
1828:         PetscScalar *nvv = avv[i*nest->nc+j];
1829:         PetscInt    *nii = aii[i*nest->nc+j];
1830:         PetscInt    *njj = ajj[i*nest->nc+j];
1831:         PetscInt    ir,cst;

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

1837:           for (ij=0;ij<rsize;ij++) {
1838:             jj[ist+ij] = *njj+cst;
1839:             vv[ist+ij] = *nvv;
1840:             njj++;
1841:             nvv++;
1842:           }
1843:           ci[ir] += rsize;
1844:           nii++;
1845:         }
1846:       }
1847:     }
1848:   }
1849:   PetscFree(ci);

1851:   /* restore info */
1852:   for (i=0; i<nest->nr; ++i) {
1853:     for (j=0; j<nest->nc; ++j) {
1854:       Mat B = nest->m[i][j];
1855:       if (B) {
1856:         PetscInt nnr = 0, k = i*nest->nc+j;

1858:         B    = (trans[k] ? trans[k] : B);
1859:         MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);
1860:         if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ");
1861:         MatSeqAIJRestoreArray(B,&avv[k]);
1862:         MatDestroy(&trans[k]);
1863:       }
1864:     }
1865:   }
1866:   PetscFree4(aii,ajj,avv,trans);

1868:   /* finalize newmat */
1869:   if (reuse == MAT_INITIAL_MATRIX) {
1870:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);
1871:   } else if (reuse == MAT_INPLACE_MATRIX) {
1872:     Mat B;

1874:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);
1875:     MatHeaderReplace(A,&B);
1876:   }
1877:   MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1878:   MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1879:   {
1880:     Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data);
1881:     a->free_a     = PETSC_TRUE;
1882:     a->free_ij    = PETSC_TRUE;
1883:   }
1884:   return(0);
1885: }

1887: PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1888: {
1890:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1891:   PetscInt       m,n,M,N,i,j,k,*dnnz,*onnz,rstart;
1892:   PetscInt       cstart,cend;
1893:   PetscMPIInt    size;
1894:   Mat            C;

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

1902:     PetscStrcmp(newtype,MATAIJ,&fast);
1903:     if (!fast) {
1904:       PetscStrcmp(newtype,MATSEQAIJ,&fast);
1905:     }
1906:     for (i=0; i<nest->nr && fast; ++i) {
1907:       for (j=0; j<nest->nc && fast; ++j) {
1908:         Mat B = nest->m[i][j];
1909:         if (B) {
1910:           PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);
1911:           if (!fast) {
1912:             PetscBool istrans;

1914:             PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1915:             if (istrans) {
1916:               Mat Bt;

1918:               MatTransposeGetMat(B,&Bt);
1919:               PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);
1920:             }
1921:           }
1922:         }
1923:       }
1924:     }
1925:     for (i=0, nf=0; i<nest->nr && fast; ++i) {
1926:       PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);
1927:       if (fast) {
1928:         PetscInt f,s;

1930:         ISStrideGetInfo(nest->isglobal.row[i],&f,&s);
1931:         if (f != nf || s != 1) { fast = PETSC_FALSE; }
1932:         else {
1933:           ISGetSize(nest->isglobal.row[i],&f);
1934:           nf  += f;
1935:         }
1936:       }
1937:     }
1938:     for (i=0, nf=0; i<nest->nc && fast; ++i) {
1939:       PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);
1940:       if (fast) {
1941:         PetscInt f,s;

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

2052:   /* Fill by row */
2053:   for (j=0; j<nest->nc; ++j) {
2054:     /* Using global column indices and ISAllGather() is not scalable. */
2055:     IS             bNis;
2056:     PetscInt       bN;
2057:     const PetscInt *bNindices;
2058:     ISAllGather(nest->isglobal.col[j], &bNis);
2059:     ISGetSize(bNis,&bN);
2060:     ISGetIndices(bNis,&bNindices);
2061:     for (i=0; i<nest->nr; ++i) {
2062:       Mat            B;
2063:       PetscInt       bm, br;
2064:       const PetscInt *bmindices;
2065:       B = nest->m[i][j];
2066:       if (!B) continue;
2067:       ISGetLocalSize(nest->isglobal.row[i],&bm);
2068:       ISGetIndices(nest->isglobal.row[i],&bmindices);
2069:       MatGetOwnershipRange(B,&rstart,NULL);
2070:       for (br = 0; br < bm; ++br) {
2071:         PetscInt          row = bmindices[br], brncols,  *cols;
2072:         const PetscInt    *brcols;
2073:         const PetscScalar *brcoldata;
2074:         MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);
2075:         PetscMalloc1(brncols,&cols);
2076:         for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
2077:         /*
2078:           Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
2079:           Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
2080:          */
2081:         MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);
2082:         MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);
2083:         PetscFree(cols);
2084:       }
2085:       ISRestoreIndices(nest->isglobal.row[i],&bmindices);
2086:     }
2087:     ISRestoreIndices(bNis,&bNindices);
2088:     ISDestroy(&bNis);
2089:   }
2090:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
2091:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
2092:   return(0);
2093: }

2095: PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has)
2096: {
2097:   Mat_Nest       *bA = (Mat_Nest*)mat->data;
2098:   MatOperation   opAdd;
2099:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;
2100:   PetscBool      flg;

2104:   *has = PETSC_FALSE;
2105:   if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
2106:     opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD);
2107:     for (j=0; j<nc; j++) {
2108:       for (i=0; i<nr; i++) {
2109:         if (!bA->m[i][j]) continue;
2110:         MatHasOperation(bA->m[i][j],opAdd,&flg);
2111:         if (!flg) return(0);
2112:       }
2113:     }
2114:   }
2115:   if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
2116:   return(0);
2117: }

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

2122:   Level: intermediate

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

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

2133: .seealso: MatCreate(), MatType, MatCreateNest(), MatNestSetSubMat(), MatNestGetSubMat(),
2134:           VecCreateNest(), DMCreateMatrix(), DMCOMPOSITE, MatNestSetVecType(), MatNestGetLocalISs(),
2135:           MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
2136: M*/
2137: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2138: {
2139:   Mat_Nest       *s;

2143:   PetscNewLog(A,&s);
2144:   A->data = (void*)s;

2146:   s->nr            = -1;
2147:   s->nc            = -1;
2148:   s->m             = NULL;
2149:   s->splitassembly = PETSC_FALSE;

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

2153:   A->ops->mult                  = MatMult_Nest;
2154:   A->ops->multadd               = MatMultAdd_Nest;
2155:   A->ops->multtranspose         = MatMultTranspose_Nest;
2156:   A->ops->multtransposeadd      = MatMultTransposeAdd_Nest;
2157:   A->ops->transpose             = MatTranspose_Nest;
2158:   A->ops->assemblybegin         = MatAssemblyBegin_Nest;
2159:   A->ops->assemblyend           = MatAssemblyEnd_Nest;
2160:   A->ops->zeroentries           = MatZeroEntries_Nest;
2161:   A->ops->copy                  = MatCopy_Nest;
2162:   A->ops->axpy                  = MatAXPY_Nest;
2163:   A->ops->duplicate             = MatDuplicate_Nest;
2164:   A->ops->createsubmatrix       = MatCreateSubMatrix_Nest;
2165:   A->ops->destroy               = MatDestroy_Nest;
2166:   A->ops->view                  = MatView_Nest;
2167:   A->ops->getvecs               = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2168:   A->ops->getlocalsubmatrix     = MatGetLocalSubMatrix_Nest;
2169:   A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2170:   A->ops->getdiagonal           = MatGetDiagonal_Nest;
2171:   A->ops->diagonalscale         = MatDiagonalScale_Nest;
2172:   A->ops->scale                 = MatScale_Nest;
2173:   A->ops->shift                 = MatShift_Nest;
2174:   A->ops->diagonalset           = MatDiagonalSet_Nest;
2175:   A->ops->setrandom             = MatSetRandom_Nest;
2176:   A->ops->hasoperation          = MatHasOperation_Nest;
2177:   A->ops->missingdiagonal       = MatMissingDiagonal_Nest;

2179:   A->spptr        = 0;
2180:   A->assembled    = PETSC_FALSE;

2182:   /* expose Nest api's */
2183:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",        MatNestGetSubMat_Nest);
2184:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",        MatNestSetSubMat_Nest);
2185:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",       MatNestGetSubMats_Nest);
2186:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",          MatNestGetSize_Nest);
2187:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",           MatNestGetISs_Nest);
2188:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",      MatNestGetLocalISs_Nest);
2189:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",       MatNestSetVecType_Nest);
2190:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",       MatNestSetSubMats_Nest);
2191:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",  MatConvert_Nest_AIJ);
2192:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",  MatConvert_Nest_AIJ);
2193:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",     MatConvert_Nest_AIJ);
2194:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",      MatConvert_Nest_IS);
2195:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",MatProductSetFromOptions_Nest_Dense);
2196:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",MatProductSetFromOptions_Nest_Dense);
2197:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",MatProductSetFromOptions_Nest_Dense);

2199:   PetscObjectChangeTypeName((PetscObject)A,MATNEST);
2200:   return(0);
2201: }