Actual source code: matnest.c

petsc-master 2020-02-23
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;
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);
127:       /* workC <- A[i][j] * B[j] */
128:       MatMatMultNumeric(bA->m[i][j],viewB,contents->workC[i*nc + j]);
129:       MatDestroy(&viewB);
130:        /* C[i] <- workC + C[i] */
131:       MatAXPY(viewC,1.0,contents->workC[i*nc + j],SAME_NONZERO_PATTERN);
132:     }
133:     MatDestroy(&viewC);
134:   }
135:   MatDenseRestoreArray(C,&carray);
136:   MatDenseRestoreArrayRead(B,&barray);
137:   return(0);
138: }

140: PetscErrorCode MatNest_DenseDestroy(void *ctx)
141: {
142:   Nest_Dense     *contents = (Nest_Dense*)ctx;
143:   PetscInt       i;

147:   PetscFree(contents->tarray);
148:   for (i=0; i<contents->k; i++) {
149:     MatDestroy(contents->workC + i);
150:   }
151:   PetscFree3(contents->dm,contents->dn,contents->workC);
152:   PetscFree(contents);
153:   return(0);
154: }

156: PETSC_INTERN PetscErrorCode MatMatMultSymbolic_Nest_Dense(Mat A,Mat B,PetscReal fill,Mat *C)
157: {
158:   Mat_Nest          *bA = (Mat_Nest*)A->data;
159:   Mat               viewB,viewSeq;
160:   const PetscScalar *barray;
161:   PetscInt          i,j,M,N,m,nr = bA->nr,nc = bA->nc,maxm = 0,ldb;
162:   PetscContainer    container;
163:   Nest_Dense        *contents;
164:   PetscErrorCode    ierr;

167:   MatGetSize(B,NULL,&N);
168:   if (!(*C)) {
169:     MatGetLocalSize(A,&m,NULL);
170:     MatGetSize(A,&M,NULL);
171:     MatCreateDense(PetscObjectComm((PetscObject)A),m,PETSC_DECIDE,M,N,NULL,C);
172:   } else {
173:     if ((*C)->rmap->n != A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrix local row dimensions are incompatible, %D != %D",(*C)->rmap->n,A->rmap->n);
174:     if ((*C)->cmap->n != B->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Matrix local column dimensions are incompatible, %D != %D",(*C)->cmap->n,B->cmap->n);
175:   }

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

216:   (*C)->ops->matmultnumeric = MatMatMultNumeric_Nest_Dense;
217:   return(0);
218: }

220: PETSC_INTERN PetscErrorCode MatMatMult_Nest_Dense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
221: {

225:   if (scall == MAT_INITIAL_MATRIX) {
226:     *C = NULL;
227:     PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
228:     MatMatMultSymbolic_Nest_Dense(A,B,fill,C);
229:     PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
230:   }
231:   PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
232:   MatMatMultNumeric_Nest_Dense(A,B,*C);
233:   PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
234:   return(0);
235: }

237: static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y)
238: {
239:   Mat_Nest       *bA = (Mat_Nest*)A->data;
240:   Vec            *bx = bA->left,*by = bA->right;
241:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

245:   for (i=0; i<nr; i++) {VecGetSubVector(x,bA->isglobal.row[i],&bx[i]);}
246:   for (i=0; i<nc; i++) {VecGetSubVector(y,bA->isglobal.col[i],&by[i]);}
247:   for (j=0; j<nc; j++) {
248:     VecZeroEntries(by[j]);
249:     for (i=0; i<nr; i++) {
250:       if (!bA->m[i][j]) continue;
251:       /* y[j] <- y[j] + (A[i][j])^T * x[i] */
252:       MatMultTransposeAdd(bA->m[i][j],bx[i],by[j],by[j]);
253:     }
254:   }
255:   for (i=0; i<nr; i++) {VecRestoreSubVector(x,bA->isglobal.row[i],&bx[i]);}
256:   for (i=0; i<nc; i++) {VecRestoreSubVector(y,bA->isglobal.col[i],&by[i]);}
257:   return(0);
258: }

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

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

288: static PetscErrorCode MatTranspose_Nest(Mat A,MatReuse reuse,Mat *B)
289: {
290:   Mat_Nest       *bA = (Mat_Nest*)A->data, *bC;
291:   Mat            C;
292:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

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

298:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
299:     Mat *subs;
300:     IS  *is_row,*is_col;

302:     PetscCalloc1(nr * nc,&subs);
303:     PetscMalloc2(nr,&is_row,nc,&is_col);
304:     MatNestGetISs(A,is_row,is_col);
305:     if (reuse == MAT_INPLACE_MATRIX) {
306:       for (i=0; i<nr; i++) {
307:         for (j=0; j<nc; j++) {
308:           subs[i + nr * j] = bA->m[i][j];
309:         }
310:       }
311:     }

313:     MatCreateNest(PetscObjectComm((PetscObject)A),nc,is_col,nr,is_row,subs,&C);
314:     PetscFree(subs);
315:     PetscFree2(is_row,is_col);
316:   } else {
317:     C = *B;
318:   }

320:   bC = (Mat_Nest*)C->data;
321:   for (i=0; i<nr; i++) {
322:     for (j=0; j<nc; j++) {
323:       if (bA->m[i][j]) {
324:         MatTranspose(bA->m[i][j], reuse, &(bC->m[j][i]));
325:       } else {
326:         bC->m[j][i] = NULL;
327:       }
328:     }
329:   }

331:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
332:     *B = C;
333:   } else {
334:     MatHeaderMerge(A, &C);
335:   }
336:   return(0);
337: }

339: static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list)
340: {
342:   IS             *lst = *list;
343:   PetscInt       i;

346:   if (!lst) return(0);
347:   for (i=0; i<n; i++) if (lst[i]) {ISDestroy(&lst[i]);}
348:   PetscFree(lst);
349:   *list = NULL;
350:   return(0);
351: }

353: static PetscErrorCode MatReset_Nest(Mat A)
354: {
355:   Mat_Nest       *vs = (Mat_Nest*)A->data;
356:   PetscInt       i,j;

360:   /* release the matrices and the place holders */
361:   MatNestDestroyISList(vs->nr,&vs->isglobal.row);
362:   MatNestDestroyISList(vs->nc,&vs->isglobal.col);
363:   MatNestDestroyISList(vs->nr,&vs->islocal.row);
364:   MatNestDestroyISList(vs->nc,&vs->islocal.col);

366:   PetscFree(vs->row_len);
367:   PetscFree(vs->col_len);
368:   PetscFree(vs->nnzstate);

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

372:   /* release the matrices and the place holders */
373:   if (vs->m) {
374:     for (i=0; i<vs->nr; i++) {
375:       for (j=0; j<vs->nc; j++) {
376:         MatDestroy(&vs->m[i][j]);
377:       }
378:       PetscFree(vs->m[i]);
379:     }
380:     PetscFree(vs->m);
381:   }

383:   /* restore defaults */
384:   vs->nr = 0;
385:   vs->nc = 0;
386:   vs->splitassembly = PETSC_FALSE;
387:   return(0);
388: }

390: static PetscErrorCode MatDestroy_Nest(Mat A)
391: {

394:   MatReset_Nest(A);
395:   PetscFree(A->data);

397:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",0);
398:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",0);
399:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",0);
400:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",0);
401:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",0);
402:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",0);
403:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",0);
404:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",0);
405:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",0);
406:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",0);
407:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",0);
408:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",0);
409:   PetscObjectComposeFunction((PetscObject)A,"MatMatMult_nest_mpidense_C",0);
410:   PetscObjectComposeFunction((PetscObject)A,"MatMatMult_nest_seqdense_C",0);
411:   PetscObjectComposeFunction((PetscObject)A,"MatMatMult_nest_dense_C",0);
412:   return(0);
413: }

415: static PetscErrorCode MatMissingDiagonal_Nest(Mat mat,PetscBool *missing,PetscInt *dd)
416: {
417:   Mat_Nest       *vs = (Mat_Nest*)mat->data;
418:   PetscInt       i;

422:   if (dd) *dd = 0;
423:   if (!vs->nr) {
424:     *missing = PETSC_TRUE;
425:     return(0);
426:   }
427:   *missing = PETSC_FALSE;
428:   for (i = 0; i < vs->nr && !(*missing); i++) {
429:     *missing = PETSC_TRUE;
430:     if (vs->m[i][i]) {
431:       MatMissingDiagonal(vs->m[i][i],missing,NULL);
432:       if (*missing && dd) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"First missing entry not yet implemented");
433:     }
434:   }
435:   return(0);
436: }

438: static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)
439: {
440:   Mat_Nest       *vs = (Mat_Nest*)A->data;
441:   PetscInt       i,j;
443:   PetscBool      nnzstate = PETSC_FALSE;

446:   for (i=0; i<vs->nr; i++) {
447:     for (j=0; j<vs->nc; j++) {
448:       PetscObjectState subnnzstate = 0;
449:       if (vs->m[i][j]) {
450:         MatAssemblyBegin(vs->m[i][j],type);
451:         if (!vs->splitassembly) {
452:           /* Note: split assembly will fail if the same block appears more than once (even indirectly through a nested
453:            * sub-block). This could be fixed by adding a flag to Mat so that there was a way to check if a Mat was
454:            * already performing an assembly, but the result would by more complicated and appears to offer less
455:            * potential for diagnostics and correctness checking. Split assembly should be fixed once there is an
456:            * interface for libraries to make asynchronous progress in "user-defined non-blocking collectives".
457:            */
458:           MatAssemblyEnd(vs->m[i][j],type);
459:           MatGetNonzeroState(vs->m[i][j],&subnnzstate);
460:         }
461:       }
462:       nnzstate = (PetscBool)(nnzstate || vs->nnzstate[i*vs->nc+j] != subnnzstate);
463:       vs->nnzstate[i*vs->nc+j] = subnnzstate;
464:     }
465:   }
466:   if (nnzstate) A->nonzerostate++;
467:   return(0);
468: }

470: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
471: {
472:   Mat_Nest       *vs = (Mat_Nest*)A->data;
473:   PetscInt       i,j;

477:   for (i=0; i<vs->nr; i++) {
478:     for (j=0; j<vs->nc; j++) {
479:       if (vs->m[i][j]) {
480:         if (vs->splitassembly) {
481:           MatAssemblyEnd(vs->m[i][j],type);
482:         }
483:       }
484:     }
485:   }
486:   return(0);
487: }

489: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B)
490: {
492:   Mat_Nest       *vs = (Mat_Nest*)A->data;
493:   PetscInt       j;
494:   Mat            sub;

497:   sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
498:   for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j];
499:   if (sub) {MatSetUp(sub);}       /* Ensure that the sizes are available */
500:   *B = sub;
501:   return(0);
502: }

504: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B)
505: {
507:   Mat_Nest       *vs = (Mat_Nest*)A->data;
508:   PetscInt       i;
509:   Mat            sub;

512:   sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
513:   for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col];
514:   if (sub) {MatSetUp(sub);}       /* Ensure that the sizes are available */
515:   *B = sub;
516:   return(0);
517: }

519: static PetscErrorCode MatNestFindIS(Mat A,PetscInt n,const IS list[],IS is,PetscInt *found)
520: {
522:   PetscInt       i;
523:   PetscBool      flg;

529:   *found = -1;
530:   for (i=0; i<n; i++) {
531:     if (!list[i]) continue;
532:     ISEqualUnsorted(list[i],is,&flg);
533:     if (flg) {
534:       *found = i;
535:       return(0);
536:     }
537:   }
538:   SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set");
539:   return(0);
540: }

542: /* Get a block row as a new MatNest */
543: static PetscErrorCode MatNestGetRow(Mat A,PetscInt row,Mat *B)
544: {
545:   Mat_Nest       *vs = (Mat_Nest*)A->data;
546:   char           keyname[256];

550:   *B   = NULL;
551:   PetscSNPrintf(keyname,sizeof(keyname),"NestRow_%D",row);
552:   PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);
553:   if (*B) return(0);

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

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

559:   PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);
560:   PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
561:   return(0);
562: }

564: static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
565: {
566:   Mat_Nest       *vs = (Mat_Nest*)A->data;
568:   PetscInt       row,col;
569:   PetscBool      same,isFullCol,isFullColGlobal;

572:   /* Check if full column space. This is a hack */
573:   isFullCol = PETSC_FALSE;
574:   PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&same);
575:   if (same) {
576:     PetscInt n,first,step,i,an,am,afirst,astep;
577:     ISStrideGetInfo(iscol,&first,&step);
578:     ISGetLocalSize(iscol,&n);
579:     isFullCol = PETSC_TRUE;
580:     for (i=0,an=A->cmap->rstart; i<vs->nc; i++) {
581:       PetscObjectTypeCompare((PetscObject)is->col[i],ISSTRIDE,&same);
582:       ISGetLocalSize(is->col[i],&am);
583:       if (same) {
584:         ISStrideGetInfo(is->col[i],&afirst,&astep);
585:         if (afirst != an || astep != step) isFullCol = PETSC_FALSE;
586:       } else isFullCol = PETSC_FALSE;
587:       an += am;
588:     }
589:     if (an != A->cmap->rstart+n) isFullCol = PETSC_FALSE;
590:   }
591:   MPIU_Allreduce(&isFullCol,&isFullColGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)iscol));

593:   if (isFullColGlobal && vs->nc > 1) {
594:     PetscInt row;
595:     MatNestFindIS(A,vs->nr,is->row,isrow,&row);
596:     MatNestGetRow(A,row,B);
597:   } else {
598:     MatNestFindIS(A,vs->nr,is->row,isrow,&row);
599:     MatNestFindIS(A,vs->nc,is->col,iscol,&col);
600:     if (!vs->m[row][col]) {
601:       PetscInt lr,lc;

603:       MatCreate(PetscObjectComm((PetscObject)A),&vs->m[row][col]);
604:       ISGetLocalSize(vs->isglobal.row[row],&lr);
605:       ISGetLocalSize(vs->isglobal.col[col],&lc);
606:       MatSetSizes(vs->m[row][col],lr,lc,PETSC_DECIDE,PETSC_DECIDE);
607:       MatSetType(vs->m[row][col],MATAIJ);
608:       MatSeqAIJSetPreallocation(vs->m[row][col],0,NULL);
609:       MatMPIAIJSetPreallocation(vs->m[row][col],0,NULL,0,NULL);
610:       MatSetUp(vs->m[row][col]);
611:       MatAssemblyBegin(vs->m[row][col],MAT_FINAL_ASSEMBLY);
612:       MatAssemblyEnd(vs->m[row][col],MAT_FINAL_ASSEMBLY);
613:     }
614:     *B = vs->m[row][col];
615:   }
616:   return(0);
617: }

619: /*
620:    TODO: This does not actually returns a submatrix we can modify
621: */
622: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
623: {
625:   Mat_Nest       *vs = (Mat_Nest*)A->data;
626:   Mat            sub;

629:   MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
630:   switch (reuse) {
631:   case MAT_INITIAL_MATRIX:
632:     if (sub) { PetscObjectReference((PetscObject)sub); }
633:     *B = sub;
634:     break;
635:   case MAT_REUSE_MATRIX:
636:     if (sub != *B) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Submatrix was not used before in this call");
637:     break;
638:   case MAT_IGNORE_MATRIX:       /* Nothing to do */
639:     break;
640:   case MAT_INPLACE_MATRIX:       /* Nothing to do */
641:     SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet");
642:     break;
643:   }
644:   return(0);
645: }

647: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
648: {
650:   Mat_Nest       *vs = (Mat_Nest*)A->data;
651:   Mat            sub;

654:   MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
655:   /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
656:   if (sub) {PetscObjectReference((PetscObject)sub);}
657:   *B = sub;
658:   return(0);
659: }

661: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
662: {
664:   Mat_Nest       *vs = (Mat_Nest*)A->data;
665:   Mat            sub;

668:   MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
669:   if (*B != sub) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has not been gotten");
670:   if (sub) {
671:     if (((PetscObject)sub)->refct <= 1) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Local submatrix has had reference count decremented too many times");
672:     MatDestroy(B);
673:   }
674:   return(0);
675: }

677: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
678: {
679:   Mat_Nest       *bA = (Mat_Nest*)A->data;
680:   PetscInt       i;

684:   for (i=0; i<bA->nr; i++) {
685:     Vec bv;
686:     VecGetSubVector(v,bA->isglobal.row[i],&bv);
687:     if (bA->m[i][i]) {
688:       MatGetDiagonal(bA->m[i][i],bv);
689:     } else {
690:       VecSet(bv,0.0);
691:     }
692:     VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
693:   }
694:   return(0);
695: }

697: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
698: {
699:   Mat_Nest       *bA = (Mat_Nest*)A->data;
700:   Vec            bl,*br;
701:   PetscInt       i,j;

705:   PetscCalloc1(bA->nc,&br);
706:   if (r) {
707:     for (j=0; j<bA->nc; j++) {VecGetSubVector(r,bA->isglobal.col[j],&br[j]);}
708:   }
709:   bl = NULL;
710:   for (i=0; i<bA->nr; i++) {
711:     if (l) {
712:       VecGetSubVector(l,bA->isglobal.row[i],&bl);
713:     }
714:     for (j=0; j<bA->nc; j++) {
715:       if (bA->m[i][j]) {
716:         MatDiagonalScale(bA->m[i][j],bl,br[j]);
717:       }
718:     }
719:     if (l) {
720:       VecRestoreSubVector(l,bA->isglobal.row[i],&bl);
721:     }
722:   }
723:   if (r) {
724:     for (j=0; j<bA->nc; j++) {VecRestoreSubVector(r,bA->isglobal.col[j],&br[j]);}
725:   }
726:   PetscFree(br);
727:   return(0);
728: }

730: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
731: {
732:   Mat_Nest       *bA = (Mat_Nest*)A->data;
733:   PetscInt       i,j;

737:   for (i=0; i<bA->nr; i++) {
738:     for (j=0; j<bA->nc; j++) {
739:       if (bA->m[i][j]) {
740:         MatScale(bA->m[i][j],a);
741:       }
742:     }
743:   }
744:   return(0);
745: }

747: static PetscErrorCode MatShift_Nest(Mat A,PetscScalar a)
748: {
749:   Mat_Nest       *bA = (Mat_Nest*)A->data;
750:   PetscInt       i;
752:   PetscBool      nnzstate = PETSC_FALSE;

755:   for (i=0; i<bA->nr; i++) {
756:     PetscObjectState subnnzstate = 0;
757:     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);
758:     MatShift(bA->m[i][i],a);
759:     MatGetNonzeroState(bA->m[i][i],&subnnzstate);
760:     nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
761:     bA->nnzstate[i*bA->nc+i] = subnnzstate;
762:   }
763:   if (nnzstate) A->nonzerostate++;
764:   return(0);
765: }

767: static PetscErrorCode MatDiagonalSet_Nest(Mat A,Vec D,InsertMode is)
768: {
769:   Mat_Nest       *bA = (Mat_Nest*)A->data;
770:   PetscInt       i;
772:   PetscBool      nnzstate = PETSC_FALSE;

775:   for (i=0; i<bA->nr; i++) {
776:     PetscObjectState subnnzstate = 0;
777:     Vec              bv;
778:     VecGetSubVector(D,bA->isglobal.row[i],&bv);
779:     if (bA->m[i][i]) {
780:       MatDiagonalSet(bA->m[i][i],bv,is);
781:       MatGetNonzeroState(bA->m[i][i],&subnnzstate);
782:     }
783:     VecRestoreSubVector(D,bA->isglobal.row[i],&bv);
784:     nnzstate = (PetscBool)(nnzstate || bA->nnzstate[i*bA->nc+i] != subnnzstate);
785:     bA->nnzstate[i*bA->nc+i] = subnnzstate;
786:   }
787:   if (nnzstate) A->nonzerostate++;
788:   return(0);
789: }

791: static PetscErrorCode MatSetRandom_Nest(Mat A,PetscRandom rctx)
792: {
793:   Mat_Nest       *bA = (Mat_Nest*)A->data;
794:   PetscInt       i,j;

798:   for (i=0; i<bA->nr; i++) {
799:     for (j=0; j<bA->nc; j++) {
800:       if (bA->m[i][j]) {
801:         MatSetRandom(bA->m[i][j],rctx);
802:       }
803:     }
804:   }
805:   return(0);
806: }

808: static PetscErrorCode MatCreateVecs_Nest(Mat A,Vec *right,Vec *left)
809: {
810:   Mat_Nest       *bA = (Mat_Nest*)A->data;
811:   Vec            *L,*R;
812:   MPI_Comm       comm;
813:   PetscInt       i,j;

817:   PetscObjectGetComm((PetscObject)A,&comm);
818:   if (right) {
819:     /* allocate R */
820:     PetscMalloc1(bA->nc, &R);
821:     /* Create the right vectors */
822:     for (j=0; j<bA->nc; j++) {
823:       for (i=0; i<bA->nr; i++) {
824:         if (bA->m[i][j]) {
825:           MatCreateVecs(bA->m[i][j],&R[j],NULL);
826:           break;
827:         }
828:       }
829:       if (i==bA->nr) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null column.");
830:     }
831:     VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
832:     /* hand back control to the nest vector */
833:     for (j=0; j<bA->nc; j++) {
834:       VecDestroy(&R[j]);
835:     }
836:     PetscFree(R);
837:   }

839:   if (left) {
840:     /* allocate L */
841:     PetscMalloc1(bA->nr, &L);
842:     /* Create the left vectors */
843:     for (i=0; i<bA->nr; i++) {
844:       for (j=0; j<bA->nc; j++) {
845:         if (bA->m[i][j]) {
846:           MatCreateVecs(bA->m[i][j],NULL,&L[i]);
847:           break;
848:         }
849:       }
850:       if (j==bA->nc) SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
851:     }

853:     VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
854:     for (i=0; i<bA->nr; i++) {
855:       VecDestroy(&L[i]);
856:     }

858:     PetscFree(L);
859:   }
860:   return(0);
861: }

863: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
864: {
865:   Mat_Nest       *bA = (Mat_Nest*)A->data;
866:   PetscBool      isascii,viewSub = PETSC_FALSE;
867:   PetscInt       i,j;

871:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
872:   if (isascii) {

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

879:     PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");
880:     for (i=0; i<bA->nr; i++) {
881:       for (j=0; j<bA->nc; j++) {
882:         MatType   type;
883:         char      name[256] = "",prefix[256] = "";
884:         PetscInt  NR,NC;
885:         PetscBool isNest = PETSC_FALSE;

887:         if (!bA->m[i][j]) {
888:           PetscViewerASCIIPrintf(viewer, "(%D,%D) : NULL \n",i,j);
889:           continue;
890:         }
891:         MatGetSize(bA->m[i][j],&NR,&NC);
892:         MatGetType(bA->m[i][j], &type);
893:         if (((PetscObject)bA->m[i][j])->name) {PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);}
894:         if (((PetscObject)bA->m[i][j])->prefix) {PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);}
895:         PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);

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

899:         if (isNest || viewSub) {
900:           PetscViewerASCIIPushTab(viewer);  /* push1 */
901:           MatView(bA->m[i][j],viewer);
902:           PetscViewerASCIIPopTab(viewer);    /* pop1 */
903:         }
904:       }
905:     }
906:     PetscViewerASCIIPopTab(viewer);    /* pop0 */
907:   }
908:   return(0);
909: }

911: static PetscErrorCode MatZeroEntries_Nest(Mat A)
912: {
913:   Mat_Nest       *bA = (Mat_Nest*)A->data;
914:   PetscInt       i,j;

918:   for (i=0; i<bA->nr; i++) {
919:     for (j=0; j<bA->nc; j++) {
920:       if (!bA->m[i][j]) continue;
921:       MatZeroEntries(bA->m[i][j]);
922:     }
923:   }
924:   return(0);
925: }

927: static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str)
928: {
929:   Mat_Nest       *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data;
930:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;
932:   PetscBool      nnzstate = PETSC_FALSE;

935:   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);
936:   for (i=0; i<nr; i++) {
937:     for (j=0; j<nc; j++) {
938:       PetscObjectState subnnzstate = 0;
939:       if (bA->m[i][j] && bB->m[i][j]) {
940:         MatCopy(bA->m[i][j],bB->m[i][j],str);
941:       } 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);
942:       MatGetNonzeroState(bB->m[i][j],&subnnzstate);
943:       nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i*nc+j] != subnnzstate);
944:       bB->nnzstate[i*nc+j] = subnnzstate;
945:     }
946:   }
947:   if (nnzstate) B->nonzerostate++;
948:   return(0);
949: }

951: static PetscErrorCode MatAXPY_Nest(Mat Y,PetscScalar a,Mat X,MatStructure str)
952: {
953:   Mat_Nest       *bY = (Mat_Nest*)Y->data,*bX = (Mat_Nest*)X->data;
954:   PetscInt       i,j,nr = bY->nr,nc = bY->nc;
956:   PetscBool      nnzstate = PETSC_FALSE;

959:   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);
960:   for (i=0; i<nr; i++) {
961:     for (j=0; j<nc; j++) {
962:       PetscObjectState subnnzstate = 0;
963:       if (bY->m[i][j] && bX->m[i][j]) {
964:         MatAXPY(bY->m[i][j],a,bX->m[i][j],str);
965:       } else if (bX->m[i][j]) {
966:         Mat M;

968:         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);
969:         MatDuplicate(bX->m[i][j],MAT_COPY_VALUES,&M);
970:         MatNestSetSubMat(Y,i,j,M);
971:         MatDestroy(&M);
972:       }
973:       if (bY->m[i][j]) { MatGetNonzeroState(bY->m[i][j],&subnnzstate); }
974:       nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i*nc+j] != subnnzstate);
975:       bY->nnzstate[i*nc+j] = subnnzstate;
976:     }
977:   }
978:   if (nnzstate) Y->nonzerostate++;
979:   return(0);
980: }

982: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
983: {
984:   Mat_Nest       *bA = (Mat_Nest*)A->data;
985:   Mat            *b;
986:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

990:   PetscMalloc1(nr*nc,&b);
991:   for (i=0; i<nr; i++) {
992:     for (j=0; j<nc; j++) {
993:       if (bA->m[i][j]) {
994:         MatDuplicate(bA->m[i][j],op,&b[i*nc+j]);
995:       } else {
996:         b[i*nc+j] = NULL;
997:       }
998:     }
999:   }
1000:   MatCreateNest(PetscObjectComm((PetscObject)A),nr,bA->isglobal.row,nc,bA->isglobal.col,b,B);
1001:   /* Give the new MatNest exclusive ownership */
1002:   for (i=0; i<nr*nc; i++) {
1003:     MatDestroy(&b[i]);
1004:   }
1005:   PetscFree(b);

1007:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1008:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1009:   return(0);
1010: }

1012: /* nest api */
1013: PetscErrorCode MatNestGetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat *mat)
1014: {
1015:   Mat_Nest *bA = (Mat_Nest*)A->data;

1018:   if (idxm >= bA->nr) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",idxm,bA->nr-1);
1019:   if (jdxm >= bA->nc) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Col too large: row %D max %D",jdxm,bA->nc-1);
1020:   *mat = bA->m[idxm][jdxm];
1021:   return(0);
1022: }

1024: /*@
1025:  MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.

1027:  Not collective

1029:  Input Parameters:
1030: +   A  - nest matrix
1031: .   idxm - index of the matrix within the nest matrix
1032: -   jdxm - index of the matrix within the nest matrix

1034:  Output Parameter:
1035: .   sub - matrix at index idxm,jdxm within the nest matrix

1037:  Level: developer

1039: .seealso: MatNestGetSize(), MatNestGetSubMats(), MatNestCreate(), MATNEST, MatNestSetSubMat(),
1040:           MatNestGetLocalISs(), MatNestGetISs()
1041: @*/
1042: PetscErrorCode  MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
1043: {

1047:   PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
1048:   return(0);
1049: }

1051: PetscErrorCode MatNestSetSubMat_Nest(Mat A,PetscInt idxm,PetscInt jdxm,Mat mat)
1052: {
1053:   Mat_Nest       *bA = (Mat_Nest*)A->data;
1054:   PetscInt       m,n,M,N,mi,ni,Mi,Ni;

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:   MatGetLocalSize(mat,&m,&n);
1061:   MatGetSize(mat,&M,&N);
1062:   ISGetLocalSize(bA->isglobal.row[idxm],&mi);
1063:   ISGetSize(bA->isglobal.row[idxm],&Mi);
1064:   ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
1065:   ISGetSize(bA->isglobal.col[jdxm],&Ni);
1066:   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);
1067:   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);

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

1072:   PetscObjectReference((PetscObject)mat);
1073:   MatDestroy(&bA->m[idxm][jdxm]);
1074:   bA->m[idxm][jdxm] = mat;
1075:   PetscObjectStateIncrease((PetscObject)A);
1076:   MatGetNonzeroState(mat,&bA->nnzstate[idxm*bA->nc+jdxm]);
1077:   A->nonzerostate++;
1078:   return(0);
1079: }

1081: /*@
1082:  MatNestSetSubMat - Set a single submatrix in the nest matrix.

1084:  Logically collective on the submatrix communicator

1086:  Input Parameters:
1087: +   A  - nest matrix
1088: .   idxm - index of the matrix within the nest matrix
1089: .   jdxm - index of the matrix within the nest matrix
1090: -   sub - matrix at index idxm,jdxm within the nest matrix

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

1095:  This increments the reference count of the submatrix.

1097:  Level: developer

1099: .seealso: MatNestSetSubMats(), MatNestGetSubMats(), MatNestGetLocalISs(), MATNEST, MatNestCreate(),
1100:           MatNestGetSubMat(), MatNestGetISs(), MatNestGetSize()
1101: @*/
1102: PetscErrorCode  MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
1103: {

1107:   PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
1108:   return(0);
1109: }

1111: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1112: {
1113:   Mat_Nest *bA = (Mat_Nest*)A->data;

1116:   if (M)   *M   = bA->nr;
1117:   if (N)   *N   = bA->nc;
1118:   if (mat) *mat = bA->m;
1119:   return(0);
1120: }

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

1125:  Not collective

1127:  Input Parameters:
1128: .   A  - nest matrix

1130:  Output Parameter:
1131: +   M - number of rows in the nest matrix
1132: .   N - number of cols in the nest matrix
1133: -   mat - 2d array of matrices

1135:  Notes:

1137:  The user should not free the array mat.

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

1143:  Level: developer

1145: .seealso: MatNestGetSize(), MatNestGetSubMat(), MatNestGetLocalISs(), MATNEST, MatNestCreate(),
1146:           MatNestSetSubMats(), MatNestGetISs(), MatNestSetSubMat()
1147: @*/
1148: PetscErrorCode  MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1149: {

1153:   PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
1154:   return(0);
1155: }

1157: PetscErrorCode  MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
1158: {
1159:   Mat_Nest *bA = (Mat_Nest*)A->data;

1162:   if (M) *M = bA->nr;
1163:   if (N) *N = bA->nc;
1164:   return(0);
1165: }

1167: /*@
1168:  MatNestGetSize - Returns the size of the nest matrix.

1170:  Not collective

1172:  Input Parameters:
1173: .   A  - nest matrix

1175:  Output Parameter:
1176: +   M - number of rows in the nested mat
1177: -   N - number of cols in the nested mat

1179:  Notes:

1181:  Level: developer

1183: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MATNEST, MatNestCreate(), MatNestGetLocalISs(),
1184:           MatNestGetISs()
1185: @*/
1186: PetscErrorCode  MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
1187: {

1191:   PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
1192:   return(0);
1193: }

1195: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
1196: {
1197:   Mat_Nest *vs = (Mat_Nest*)A->data;
1198:   PetscInt i;

1201:   if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
1202:   if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
1203:   return(0);
1204: }

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

1209:  Not collective

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

1214:  Output Parameter:
1215: +   rows - array of row index sets
1216: -   cols - array of column index sets

1218:  Level: advanced

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

1223: .seealso: MatNestGetSubMat(), MatNestGetSubMats(), MatNestGetSize(), MatNestGetLocalISs(), MATNEST,
1224:           MatNestCreate(), MatNestGetSubMats(), MatNestSetSubMats()
1225: @*/
1226: PetscErrorCode  MatNestGetISs(Mat A,IS rows[],IS cols[])
1227: {

1232:   PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1233:   return(0);
1234: }

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

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

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

1250:  Not collective

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

1255:  Output Parameter:
1256: +   rows - array of row index sets (or NULL to ignore)
1257: -   cols - array of column index sets (or NULL to ignore)

1259:  Level: advanced

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

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

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

1277: PetscErrorCode  MatNestSetVecType_Nest(Mat A,VecType vtype)
1278: {
1280:   PetscBool      flg;

1283:   PetscStrcmp(vtype,VECNEST,&flg);
1284:   /* In reality, this only distinguishes VECNEST and "other" */
1285:   if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1286:   else A->ops->getvecs = (PetscErrorCode (*)(Mat,Vec*,Vec*)) 0;
1287:   return(0);
1288: }

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

1293:  Not collective

1295:  Input Parameters:
1296: +  A  - nest matrix
1297: -  vtype - type to use for creating vectors

1299:  Notes:

1301:  Level: developer

1303: .seealso: MatCreateVecs(), MATNEST, MatNestCreate()
1304: @*/
1305: PetscErrorCode  MatNestSetVecType(Mat A,VecType vtype)
1306: {

1310:   PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
1311:   return(0);
1312: }

1314: PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1315: {
1316:   Mat_Nest       *s = (Mat_Nest*)A->data;
1317:   PetscInt       i,j,m,n,M,N;
1319:   PetscBool      cong;

1322:   MatReset_Nest(A);

1324:   s->nr = nr;
1325:   s->nc = nc;

1327:   /* Create space for submatrices */
1328:   PetscMalloc1(nr,&s->m);
1329:   for (i=0; i<nr; i++) {
1330:     PetscMalloc1(nc,&s->m[i]);
1331:   }
1332:   for (i=0; i<nr; i++) {
1333:     for (j=0; j<nc; j++) {
1334:       s->m[i][j] = a[i*nc+j];
1335:       if (a[i*nc+j]) {
1336:         PetscObjectReference((PetscObject)a[i*nc+j]);
1337:       }
1338:     }
1339:   }

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

1343:   PetscMalloc1(nr,&s->row_len);
1344:   PetscMalloc1(nc,&s->col_len);
1345:   for (i=0; i<nr; i++) s->row_len[i]=-1;
1346:   for (j=0; j<nc; j++) s->col_len[j]=-1;

1348:   PetscCalloc1(nr*nc,&s->nnzstate);
1349:   for (i=0; i<nr; i++) {
1350:     for (j=0; j<nc; j++) {
1351:       if (s->m[i][j]) {
1352:         MatGetNonzeroState(s->m[i][j],&s->nnzstate[i*nc+j]);
1353:       }
1354:     }
1355:   }

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

1359:   PetscLayoutSetSize(A->rmap,M);
1360:   PetscLayoutSetLocalSize(A->rmap,m);
1361:   PetscLayoutSetSize(A->cmap,N);
1362:   PetscLayoutSetLocalSize(A->cmap,n);

1364:   PetscLayoutSetUp(A->rmap);
1365:   PetscLayoutSetUp(A->cmap);

1367:   /* disable operations that are not supported for non-square matrices,
1368:      or matrices for which is_row != is_col  */
1369:   MatHasCongruentLayouts(A,&cong);
1370:   if (cong && nr != nc) cong = PETSC_FALSE;
1371:   if (cong) {
1372:     for (i = 0; cong && i < nr; i++) {
1373:       ISEqualUnsorted(s->isglobal.row[i],s->isglobal.col[i],&cong);
1374:     }
1375:   }
1376:   if (!cong) {
1377:     A->ops->missingdiagonal = NULL;
1378:     A->ops->getdiagonal     = NULL;
1379:     A->ops->shift           = NULL;
1380:     A->ops->diagonalset     = NULL;
1381:   }

1383:   PetscCalloc2(nr,&s->left,nc,&s->right);
1384:   PetscObjectStateIncrease((PetscObject)A);
1385:   A->nonzerostate++;
1386:   return(0);
1387: }

1389: /*@
1390:    MatNestSetSubMats - Sets the nested submatrices

1392:    Collective on Mat

1394:    Input Parameter:
1395: +  A - nested matrix
1396: .  nr - number of nested row blocks
1397: .  is_row - index sets for each nested row block, or NULL to make contiguous
1398: .  nc - number of nested column blocks
1399: .  is_col - index sets for each nested column block, or NULL to make contiguous
1400: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

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

1404:    Level: advanced

1406: .seealso: MatCreateNest(), MATNEST, MatNestSetSubMat(), MatNestGetSubMat(), MatNestGetSubMats()
1407: @*/
1408: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1409: {
1411:   PetscInt       i;

1415:   if (nr < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of rows cannot be negative");
1416:   if (nr && is_row) {
1419:   }
1420:   if (nc < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Number of columns cannot be negative");
1421:   if (nc && is_col) {
1424:   }
1426:   PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1427:   return(0);
1428: }

1430: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1431: {
1433:   PetscBool      flg;
1434:   PetscInt       i,j,m,mi,*ix;

1437:   *ltog = NULL;
1438:   for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1439:     if (islocal[i]) {
1440:       ISGetLocalSize(islocal[i],&mi);
1441:       flg  = PETSC_TRUE;      /* We found a non-trivial entry */
1442:     } else {
1443:       ISGetLocalSize(isglobal[i],&mi);
1444:     }
1445:     m += mi;
1446:   }
1447:   if (!flg) return(0);

1449:   PetscMalloc1(m,&ix);
1450:   for (i=0,m=0; i<n; i++) {
1451:     ISLocalToGlobalMapping smap = NULL;
1452:     Mat                    sub = NULL;
1453:     PetscSF                sf;
1454:     PetscLayout            map;
1455:     const PetscInt         *ix2;

1457:     if (!colflg) {
1458:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1459:     } else {
1460:       MatNestFindNonzeroSubMatCol(A,i,&sub);
1461:     }
1462:     if (sub) {
1463:       if (!colflg) {
1464:         MatGetLocalToGlobalMapping(sub,&smap,NULL);
1465:       } else {
1466:         MatGetLocalToGlobalMapping(sub,NULL,&smap);
1467:       }
1468:     }
1469:     /*
1470:        Now we need to extract the monolithic global indices that correspond to the given split global indices.
1471:        In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1472:     */
1473:     ISGetIndices(isglobal[i],&ix2);
1474:     if (islocal[i]) {
1475:       PetscInt *ilocal,*iremote;
1476:       PetscInt mil,nleaves;

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

1483:       /* PetscSFSetGraphLayout does not like negative indices */
1484:       PetscMalloc2(mi,&ilocal,mi,&iremote);
1485:       for (j=0, nleaves = 0; j<mi; j++) {
1486:         if (ix[m+j] < 0) continue;
1487:         ilocal[nleaves]  = j;
1488:         iremote[nleaves] = ix[m+j];
1489:         nleaves++;
1490:       }
1491:       ISGetLocalSize(isglobal[i],&mil);
1492:       PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
1493:       PetscLayoutCreate(PetscObjectComm((PetscObject)A),&map);
1494:       PetscLayoutSetLocalSize(map,mil);
1495:       PetscLayoutSetUp(map);
1496:       PetscSFSetGraphLayout(sf,map,nleaves,ilocal,PETSC_USE_POINTER,iremote);
1497:       PetscLayoutDestroy(&map);
1498:       PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m);
1499:       PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m);
1500:       PetscSFDestroy(&sf);
1501:       PetscFree2(ilocal,iremote);
1502:     } else {
1503:       ISGetLocalSize(isglobal[i],&mi);
1504:       for (j=0; j<mi; j++) ix[m+j] = ix2[i];
1505:     }
1506:     ISRestoreIndices(isglobal[i],&ix2);
1507:     m   += mi;
1508:   }
1509:   ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);
1510:   return(0);
1511: }


1514: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1515: /*
1516:   nprocessors = NP
1517:   Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1518:        proc 0: => (g_0,h_0,)
1519:        proc 1: => (g_1,h_1,)
1520:        ...
1521:        proc nprocs-1: => (g_NP-1,h_NP-1,)

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

1526:             proc 0:
1527:     is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1528:             proc 1:
1529:     is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)

1531:             proc NP-1:
1532:     is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1533: */
1534: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1535: {
1536:   Mat_Nest       *vs = (Mat_Nest*)A->data;
1537:   PetscInt       i,j,offset,n,nsum,bs;
1539:   Mat            sub = NULL;

1542:   PetscMalloc1(nr,&vs->isglobal.row);
1543:   PetscMalloc1(nc,&vs->isglobal.col);
1544:   if (is_row) { /* valid IS is passed in */
1545:     /* refs on is[] are incremeneted */
1546:     for (i=0; i<vs->nr; i++) {
1547:       PetscObjectReference((PetscObject)is_row[i]);

1549:       vs->isglobal.row[i] = is_row[i];
1550:     }
1551:   } else {                      /* Create the ISs by inspecting sizes of a submatrix in each row */
1552:     nsum = 0;
1553:     for (i=0; i<vs->nr; i++) {  /* Add up the local sizes to compute the aggregate offset */
1554:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1555:       if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in row %D",i);
1556:       MatGetLocalSize(sub,&n,NULL);
1557:       if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1558:       nsum += n;
1559:     }
1560:     MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1561:     offset -= nsum;
1562:     for (i=0; i<vs->nr; i++) {
1563:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1564:       MatGetLocalSize(sub,&n,NULL);
1565:       MatGetBlockSize(sub,&bs);
1566:       ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);
1567:       ISSetBlockSize(vs->isglobal.row[i],bs);
1568:       offset += n;
1569:     }
1570:   }

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

1577:       vs->isglobal.col[j] = is_col[j];
1578:     }
1579:   } else {                      /* Create the ISs by inspecting sizes of a submatrix in each column */
1580:     offset = A->cmap->rstart;
1581:     nsum   = 0;
1582:     for (j=0; j<vs->nc; j++) {
1583:       MatNestFindNonzeroSubMatCol(A,j,&sub);
1584:       if (!sub) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"No nonzero submatrix in column %D",i);
1585:       MatGetLocalSize(sub,NULL,&n);
1586:       if (n < 0) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Sizes have not yet been set for submatrix");
1587:       nsum += n;
1588:     }
1589:     MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1590:     offset -= nsum;
1591:     for (j=0; j<vs->nc; j++) {
1592:       MatNestFindNonzeroSubMatCol(A,j,&sub);
1593:       MatGetLocalSize(sub,NULL,&n);
1594:       MatGetBlockSize(sub,&bs);
1595:       ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);
1596:       ISSetBlockSize(vs->isglobal.col[j],bs);
1597:       offset += n;
1598:     }
1599:   }

1601:   /* Set up the local ISs */
1602:   PetscMalloc1(vs->nr,&vs->islocal.row);
1603:   PetscMalloc1(vs->nc,&vs->islocal.col);
1604:   for (i=0,offset=0; i<vs->nr; i++) {
1605:     IS                     isloc;
1606:     ISLocalToGlobalMapping rmap = NULL;
1607:     PetscInt               nlocal,bs;
1608:     MatNestFindNonzeroSubMatRow(A,i,&sub);
1609:     if (sub) {MatGetLocalToGlobalMapping(sub,&rmap,NULL);}
1610:     if (rmap) {
1611:       MatGetBlockSize(sub,&bs);
1612:       ISLocalToGlobalMappingGetSize(rmap,&nlocal);
1613:       ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1614:       ISSetBlockSize(isloc,bs);
1615:     } else {
1616:       nlocal = 0;
1617:       isloc  = NULL;
1618:     }
1619:     vs->islocal.row[i] = isloc;
1620:     offset            += nlocal;
1621:   }
1622:   for (i=0,offset=0; i<vs->nc; i++) {
1623:     IS                     isloc;
1624:     ISLocalToGlobalMapping cmap = NULL;
1625:     PetscInt               nlocal,bs;
1626:     MatNestFindNonzeroSubMatCol(A,i,&sub);
1627:     if (sub) {MatGetLocalToGlobalMapping(sub,NULL,&cmap);}
1628:     if (cmap) {
1629:       MatGetBlockSize(sub,&bs);
1630:       ISLocalToGlobalMappingGetSize(cmap,&nlocal);
1631:       ISCreateStride(PETSC_COMM_SELF,nlocal,offset,1,&isloc);
1632:       ISSetBlockSize(isloc,bs);
1633:     } else {
1634:       nlocal = 0;
1635:       isloc  = NULL;
1636:     }
1637:     vs->islocal.col[i] = isloc;
1638:     offset            += nlocal;
1639:   }

1641:   /* Set up the aggregate ISLocalToGlobalMapping */
1642:   {
1643:     ISLocalToGlobalMapping rmap,cmap;
1644:     MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);
1645:     MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);
1646:     if (rmap && cmap) {MatSetLocalToGlobalMapping(A,rmap,cmap);}
1647:     ISLocalToGlobalMappingDestroy(&rmap);
1648:     ISLocalToGlobalMappingDestroy(&cmap);
1649:   }

1651: #if defined(PETSC_USE_DEBUG)
1652:   for (i=0; i<vs->nr; i++) {
1653:     for (j=0; j<vs->nc; j++) {
1654:       PetscInt m,n,M,N,mi,ni,Mi,Ni;
1655:       Mat      B = vs->m[i][j];
1656:       if (!B) continue;
1657:       MatGetSize(B,&M,&N);
1658:       MatGetLocalSize(B,&m,&n);
1659:       ISGetSize(vs->isglobal.row[i],&Mi);
1660:       ISGetSize(vs->isglobal.col[j],&Ni);
1661:       ISGetLocalSize(vs->isglobal.row[i],&mi);
1662:       ISGetLocalSize(vs->isglobal.col[j],&ni);
1663:       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);
1664:       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);
1665:     }
1666:   }
1667: #endif

1669:   /* Set A->assembled if all non-null blocks are currently assembled */
1670:   for (i=0; i<vs->nr; i++) {
1671:     for (j=0; j<vs->nc; j++) {
1672:       if (vs->m[i][j] && !vs->m[i][j]->assembled) return(0);
1673:     }
1674:   }
1675:   A->assembled = PETSC_TRUE;
1676:   return(0);
1677: }

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

1682:    Collective on Mat

1684:    Input Parameter:
1685: +  comm - Communicator for the new Mat
1686: .  nr - number of nested row blocks
1687: .  is_row - index sets for each nested row block, or NULL to make contiguous
1688: .  nc - number of nested column blocks
1689: .  is_col - index sets for each nested column block, or NULL to make contiguous
1690: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

1692:    Output Parameter:
1693: .  B - new matrix

1695:    Level: advanced

1697: .seealso: MatCreate(), VecCreateNest(), DMCreateMatrix(), MATNEST, MatNestSetSubMat(),
1698:           MatNestGetSubMat(), MatNestGetLocalISs(), MatNestGetSize(),
1699:           MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
1700: @*/
1701: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1702: {
1703:   Mat            A;

1707:   *B   = 0;
1708:   MatCreate(comm,&A);
1709:   MatSetType(A,MATNEST);
1710:   A->preallocated = PETSC_TRUE;
1711:   MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1712:   *B   = A;
1713:   return(0);
1714: }

1716: static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1717: {
1718:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1719:   Mat            *trans;
1720:   PetscScalar    **avv;
1721:   PetscScalar    *vv;
1722:   PetscInt       **aii,**ajj;
1723:   PetscInt       *ii,*jj,*ci;
1724:   PetscInt       nr,nc,nnz,i,j;
1725:   PetscBool      done;

1729:   MatGetSize(A,&nr,&nc);
1730:   if (reuse == MAT_REUSE_MATRIX) {
1731:     PetscInt rnr;

1733:     MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);
1734:     if (!done) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"MatGetRowIJ");
1735:     if (rnr != nr) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of rows");
1736:     MatSeqAIJGetArray(*newmat,&vv);
1737:   }
1738:   /* extract CSR for nested SeqAIJ matrices */
1739:   nnz  = 0;
1740:   PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);
1741:   for (i=0; i<nest->nr; ++i) {
1742:     for (j=0; j<nest->nc; ++j) {
1743:       Mat B = nest->m[i][j];
1744:       if (B) {
1745:         PetscScalar *naa;
1746:         PetscInt    *nii,*njj,nnr;
1747:         PetscBool   istrans;

1749:         PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1750:         if (istrans) {
1751:           Mat Bt;

1753:           MatTransposeGetMat(B,&Bt);
1754:           MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);
1755:           B    = trans[i*nest->nc+j];
1756:         }
1757:         MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);
1758:         if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatGetRowIJ");
1759:         MatSeqAIJGetArray(B,&naa);
1760:         nnz += nii[nnr];

1762:         aii[i*nest->nc+j] = nii;
1763:         ajj[i*nest->nc+j] = njj;
1764:         avv[i*nest->nc+j] = naa;
1765:       }
1766:     }
1767:   }
1768:   if (reuse != MAT_REUSE_MATRIX) {
1769:     PetscMalloc1(nr+1,&ii);
1770:     PetscMalloc1(nnz,&jj);
1771:     PetscMalloc1(nnz,&vv);
1772:   } else {
1773:     if (nnz != ii[nr]) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Cannot reuse matrix, wrong number of nonzeros");
1774:   }

1776:   /* new row pointer */
1777:   PetscArrayzero(ii,nr+1);
1778:   for (i=0; i<nest->nr; ++i) {
1779:     PetscInt       ncr,rst;

1781:     ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1782:     ISGetLocalSize(nest->isglobal.row[i],&ncr);
1783:     for (j=0; j<nest->nc; ++j) {
1784:       if (aii[i*nest->nc+j]) {
1785:         PetscInt    *nii = aii[i*nest->nc+j];
1786:         PetscInt    ir;

1788:         for (ir=rst; ir<ncr+rst; ++ir) {
1789:           ii[ir+1] += nii[1]-nii[0];
1790:           nii++;
1791:         }
1792:       }
1793:     }
1794:   }
1795:   for (i=0; i<nr; i++) ii[i+1] += ii[i];

1797:   /* construct CSR for the new matrix */
1798:   PetscCalloc1(nr,&ci);
1799:   for (i=0; i<nest->nr; ++i) {
1800:     PetscInt       ncr,rst;

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

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

1815:           for (ij=0;ij<rsize;ij++) {
1816:             jj[ist+ij] = *njj+cst;
1817:             vv[ist+ij] = *nvv;
1818:             njj++;
1819:             nvv++;
1820:           }
1821:           ci[ir] += rsize;
1822:           nii++;
1823:         }
1824:       }
1825:     }
1826:   }
1827:   PetscFree(ci);

1829:   /* restore info */
1830:   for (i=0; i<nest->nr; ++i) {
1831:     for (j=0; j<nest->nc; ++j) {
1832:       Mat B = nest->m[i][j];
1833:       if (B) {
1834:         PetscInt nnr = 0, k = i*nest->nc+j;

1836:         B    = (trans[k] ? trans[k] : B);
1837:         MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);
1838:         if (!done) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_PLIB,"MatRestoreRowIJ");
1839:         MatSeqAIJRestoreArray(B,&avv[k]);
1840:         MatDestroy(&trans[k]);
1841:       }
1842:     }
1843:   }
1844:   PetscFree4(aii,ajj,avv,trans);

1846:   /* finalize newmat */
1847:   if (reuse == MAT_INITIAL_MATRIX) {
1848:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);
1849:   } else if (reuse == MAT_INPLACE_MATRIX) {
1850:     Mat B;

1852:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,&B);
1853:     MatHeaderReplace(A,&B);
1854:   }
1855:   MatAssemblyBegin(*newmat,MAT_FINAL_ASSEMBLY);
1856:   MatAssemblyEnd(*newmat,MAT_FINAL_ASSEMBLY);
1857:   {
1858:     Mat_SeqAIJ *a = (Mat_SeqAIJ*)((*newmat)->data);
1859:     a->free_a     = PETSC_TRUE;
1860:     a->free_ij    = PETSC_TRUE;
1861:   }
1862:   return(0);
1863: }

1865: PETSC_INTERN PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1866: {
1868:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1869:   PetscInt       m,n,M,N,i,j,k,*dnnz,*onnz,rstart;
1870:   PetscInt       cstart,cend;
1871:   PetscMPIInt    size;
1872:   Mat            C;

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

1880:     PetscStrcmp(newtype,MATAIJ,&fast);
1881:     if (!fast) {
1882:       PetscStrcmp(newtype,MATSEQAIJ,&fast);
1883:     }
1884:     for (i=0; i<nest->nr && fast; ++i) {
1885:       for (j=0; j<nest->nc && fast; ++j) {
1886:         Mat B = nest->m[i][j];
1887:         if (B) {
1888:           PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);
1889:           if (!fast) {
1890:             PetscBool istrans;

1892:             PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1893:             if (istrans) {
1894:               Mat Bt;

1896:               MatTransposeGetMat(B,&Bt);
1897:               PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);
1898:             }
1899:           }
1900:         }
1901:       }
1902:     }
1903:     for (i=0, nf=0; i<nest->nr && fast; ++i) {
1904:       PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);
1905:       if (fast) {
1906:         PetscInt f,s;

1908:         ISStrideGetInfo(nest->isglobal.row[i],&f,&s);
1909:         if (f != nf || s != 1) { fast = PETSC_FALSE; }
1910:         else {
1911:           ISGetSize(nest->isglobal.row[i],&f);
1912:           nf  += f;
1913:         }
1914:       }
1915:     }
1916:     for (i=0, nf=0; i<nest->nc && fast; ++i) {
1917:       PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);
1918:       if (fast) {
1919:         PetscInt f,s;

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

2030:   /* Fill by row */
2031:   for (j=0; j<nest->nc; ++j) {
2032:     /* Using global column indices and ISAllGather() is not scalable. */
2033:     IS             bNis;
2034:     PetscInt       bN;
2035:     const PetscInt *bNindices;
2036:     ISAllGather(nest->isglobal.col[j], &bNis);
2037:     ISGetSize(bNis,&bN);
2038:     ISGetIndices(bNis,&bNindices);
2039:     for (i=0; i<nest->nr; ++i) {
2040:       Mat            B;
2041:       PetscInt       bm, br;
2042:       const PetscInt *bmindices;
2043:       B = nest->m[i][j];
2044:       if (!B) continue;
2045:       ISGetLocalSize(nest->isglobal.row[i],&bm);
2046:       ISGetIndices(nest->isglobal.row[i],&bmindices);
2047:       MatGetOwnershipRange(B,&rstart,NULL);
2048:       for (br = 0; br < bm; ++br) {
2049:         PetscInt          row = bmindices[br], brncols,  *cols;
2050:         const PetscInt    *brcols;
2051:         const PetscScalar *brcoldata;
2052:         MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);
2053:         PetscMalloc1(brncols,&cols);
2054:         for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
2055:         /*
2056:           Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
2057:           Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
2058:          */
2059:         MatSetValues(C,1,&row,brncols,cols,brcoldata,ADD_VALUES);
2060:         MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);
2061:         PetscFree(cols);
2062:       }
2063:       ISRestoreIndices(nest->isglobal.row[i],&bmindices);
2064:     }
2065:     ISRestoreIndices(bNis,&bNindices);
2066:     ISDestroy(&bNis);
2067:   }
2068:   MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
2069:   MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
2070:   return(0);
2071: }

2073: PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has)
2074: {
2075:   Mat_Nest       *bA = (Mat_Nest*)mat->data;
2076:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;
2077:   PetscBool      flg;

2081:   *has = PETSC_FALSE;
2082:   if (op == MATOP_MULT_TRANSPOSE || op == MATOP_MAT_MULT) {
2083:     for (j=0; j<nc; j++) {
2084:       for (i=0; i<nr; i++) {
2085:         if (!bA->m[i][j]) continue;
2086:         MatHasOperation(bA->m[i][j],op,&flg);
2087:         if (!flg) return(0);
2088:       }
2089:     }
2090:   }
2091:   if (((void**)mat->ops)[op] || (op == MATOP_MAT_MULT && flg)) *has = PETSC_TRUE;
2092:   return(0);
2093: }

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

2098:   Level: intermediate

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

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

2109: .seealso: MatCreate(), MatType, MatCreateNest(), MatNestSetSubMat(), MatNestGetSubMat(),
2110:           VecCreateNest(), DMCreateMatrix(), DMCOMPOSITE, MatNestSetVecType(), MatNestGetLocalISs(),
2111:           MatNestGetISs(), MatNestSetSubMats(), MatNestGetSubMats()
2112: M*/
2113: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2114: {
2115:   Mat_Nest       *s;

2119:   PetscNewLog(A,&s);
2120:   A->data = (void*)s;

2122:   s->nr            = -1;
2123:   s->nc            = -1;
2124:   s->m             = NULL;
2125:   s->splitassembly = PETSC_FALSE;

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

2129:   A->ops->mult                  = MatMult_Nest;
2130:   A->ops->multadd               = MatMultAdd_Nest;
2131:   A->ops->multtranspose         = MatMultTranspose_Nest;
2132:   A->ops->multtransposeadd      = MatMultTransposeAdd_Nest;
2133:   A->ops->transpose             = MatTranspose_Nest;
2134:   A->ops->assemblybegin         = MatAssemblyBegin_Nest;
2135:   A->ops->assemblyend           = MatAssemblyEnd_Nest;
2136:   A->ops->zeroentries           = MatZeroEntries_Nest;
2137:   A->ops->copy                  = MatCopy_Nest;
2138:   A->ops->axpy                  = MatAXPY_Nest;
2139:   A->ops->duplicate             = MatDuplicate_Nest;
2140:   A->ops->createsubmatrix       = MatCreateSubMatrix_Nest;
2141:   A->ops->destroy               = MatDestroy_Nest;
2142:   A->ops->view                  = MatView_Nest;
2143:   A->ops->getvecs               = 0; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2144:   A->ops->getlocalsubmatrix     = MatGetLocalSubMatrix_Nest;
2145:   A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2146:   A->ops->getdiagonal           = MatGetDiagonal_Nest;
2147:   A->ops->diagonalscale         = MatDiagonalScale_Nest;
2148:   A->ops->scale                 = MatScale_Nest;
2149:   A->ops->shift                 = MatShift_Nest;
2150:   A->ops->diagonalset           = MatDiagonalSet_Nest;
2151:   A->ops->setrandom             = MatSetRandom_Nest;
2152:   A->ops->hasoperation          = MatHasOperation_Nest;
2153:   A->ops->missingdiagonal       = MatMissingDiagonal_Nest;

2155:   A->spptr        = 0;
2156:   A->assembled    = PETSC_FALSE;

2158:   /* expose Nest api's */
2159:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",        MatNestGetSubMat_Nest);
2160:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",        MatNestSetSubMat_Nest);
2161:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",       MatNestGetSubMats_Nest);
2162:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",          MatNestGetSize_Nest);
2163:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",           MatNestGetISs_Nest);
2164:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",      MatNestGetLocalISs_Nest);
2165:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",       MatNestSetVecType_Nest);
2166:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",       MatNestSetSubMats_Nest);
2167:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",  MatConvert_Nest_AIJ);
2168:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",  MatConvert_Nest_AIJ);
2169:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",     MatConvert_Nest_AIJ);
2170:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",      MatConvert_Nest_IS);
2171:   PetscObjectComposeFunction((PetscObject)A,"MatMatMult_nest_mpidense_C",MatMatMult_Nest_Dense);
2172:   PetscObjectComposeFunction((PetscObject)A,"MatMatMult_nest_seqdense_C",MatMatMult_Nest_Dense);
2173:   PetscObjectComposeFunction((PetscObject)A,"MatMatMult_nest_dense_C",   MatMatMult_Nest_Dense);

2175:   PetscObjectChangeTypeName((PetscObject)A,MATNEST);
2176:   return(0);
2177: }