Actual source code: matnest.c

  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;

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

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

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

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

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

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

 89: PETSC_INTERN PetscErrorCode MatProductNumeric_Nest_Dense(Mat C)
 90: {
 91:   Mat_Nest          *bA;
 92:   Nest_Dense        *contents;
 93:   Mat               viewB,viewC,productB,workC;
 94:   const PetscScalar *barray;
 95:   PetscScalar       *carray;
 96:   PetscInt          i,j,M,N,nr,nc,ldb,ldc;
 97:   Mat               A,B;

 99:   MatCheckProduct(C,3);
100:   A    = C->product->A;
101:   B    = C->product->B;
102:   MatGetSize(B,NULL,&N);
103:   if (!N) {
104:     MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
105:     MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
106:     return 0;
107:   }
108:   contents = (Nest_Dense*)C->product->data;
110:   bA   = (Mat_Nest*)A->data;
111:   nr   = bA->nr;
112:   nc   = bA->nc;
113:   MatDenseGetLDA(B,&ldb);
114:   MatDenseGetLDA(C,&ldc);
115:   MatZeroEntries(C);
116:   MatDenseGetArrayRead(B,&barray);
117:   MatDenseGetArray(C,&carray);
118:   for (i=0; i<nr; i++) {
119:     ISGetSize(bA->isglobal.row[i],&M);
120:     MatCreateDense(PetscObjectComm((PetscObject)A),contents->dm[i+1]-contents->dm[i],PETSC_DECIDE,M,N,carray+contents->dm[i],&viewC);
121:     MatDenseSetLDA(viewC,ldc);
122:     for (j=0; j<nc; j++) {
123:       if (!bA->m[i][j]) continue;
124:       ISGetSize(bA->isglobal.col[j],&M);
125:       MatCreateDense(PetscObjectComm((PetscObject)A),contents->dn[j+1]-contents->dn[j],PETSC_DECIDE,M,N,(PetscScalar*)(barray+contents->dn[j]),&viewB);
126:       MatDenseSetLDA(viewB,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:       MatProductNumeric(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;

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

163: PETSC_INTERN PetscErrorCode MatProductSymbolic_Nest_Dense(Mat C)
164: {
165:   Mat_Nest          *bA;
166:   Mat               viewB,workC;
167:   const PetscScalar *barray;
168:   PetscInt          i,j,M,N,m,n,nr,nc,maxm = 0,ldb;
169:   Nest_Dense        *contents=NULL;
170:   PetscBool         cisdense;
171:   Mat               A,B;
172:   PetscReal         fill;

174:   MatCheckProduct(C,4);
176:   A    = C->product->A;
177:   B    = C->product->B;
178:   fill = C->product->fill;
179:   bA   = (Mat_Nest*)A->data;
180:   nr   = bA->nr;
181:   nc   = bA->nc;
182:   MatGetLocalSize(C,&m,&n);
183:   MatGetSize(C,&M,&N);
184:   if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
185:     MatGetLocalSize(B,NULL,&n);
186:     MatGetSize(B,NULL,&N);
187:     MatGetLocalSize(A,&m,NULL);
188:     MatGetSize(A,&M,NULL);
189:     MatSetSizes(C,m,n,M,N);
190:   }
191:   PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATMPIDENSE,MATSEQDENSECUDA,MATMPIDENSECUDA,"");
192:   if (!cisdense) {
193:     MatSetType(C,((PetscObject)B)->type_name);
194:   }
195:   MatSetUp(C);
196:   if (!N) {
197:     C->ops->productnumeric = MatProductNumeric_Nest_Dense;
198:     return 0;
199:   }

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

228:       MatProductCreate(bA->m[i][j],viewB,NULL,&contents->workC[i*nc + j]);
229:       workC = contents->workC[i*nc + j];
230:       MatProductSetType(workC,MATPRODUCT_AB);
231:       MatProductSetAlgorithm(workC,"default");
232:       MatProductSetFill(workC,fill);
233:       MatProductSetFromOptions(workC);
234:       MatProductSymbolic(workC);

236:       /* since tarray will be shared by all Mat */
237:       MatSeqDenseSetPreallocation(workC,contents->tarray);
238:       MatMPIDenseSetPreallocation(workC,contents->tarray);
239:     }
240:     MatDestroy(&viewB);
241:   }
242:   MatDenseRestoreArrayRead(B,&barray);

244:   C->ops->productnumeric = MatProductNumeric_Nest_Dense;
245:   return 0;
246: }

248: /* --------------------------------------------------------- */
249: static PetscErrorCode MatProductSetFromOptions_Nest_Dense_AB(Mat C)
250: {
251:   C->ops->productsymbolic = MatProductSymbolic_Nest_Dense;
252:   return 0;
253: }

255: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Nest_Dense(Mat C)
256: {
257:   Mat_Product    *product = C->product;

259:   if (product->type == MATPRODUCT_AB) {
260:     MatProductSetFromOptions_Nest_Dense_AB(C);
261:   }
262:   return 0;
263: }
264: /* --------------------------------------------------------- */

266: static PetscErrorCode MatMultTranspose_Nest(Mat A,Vec x,Vec y)
267: {
268:   Mat_Nest       *bA = (Mat_Nest*)A->data;
269:   Vec            *bx = bA->left,*by = bA->right;
270:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

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

287: static PetscErrorCode MatMultTransposeAdd_Nest(Mat A,Vec x,Vec y,Vec z)
288: {
289:   Mat_Nest       *bA = (Mat_Nest*)A->data;
290:   Vec            *bx = bA->left,*bz = bA->right;
291:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

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

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

319:   if (reuse == MAT_REUSE_MATRIX) MatTransposeCheckNonzeroState_Private(A,*B);

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

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

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

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

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

363: static PetscErrorCode MatNestDestroyISList(PetscInt n,IS **list)
364: {
365:   IS             *lst = *list;
366:   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;

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

386:   PetscFree(vs->row_len);
387:   PetscFree(vs->col_len);
388:   PetscFree(vs->nnzstate);

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

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

403:   /* restore defaults */
404:   vs->nr = 0;
405:   vs->nc = 0;
406:   vs->splitassembly = PETSC_FALSE;
407:   return 0;
408: }

410: static PetscErrorCode MatDestroy_Nest(Mat A)
411: {
412:   MatReset_Nest(A);
413:   PetscFree(A->data);
414:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",NULL);
415:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",NULL);
416:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",NULL);
417:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",NULL);
418:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",NULL);
419:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",NULL);
420:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",NULL);
421:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",NULL);
422:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",NULL);
423:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",NULL);
424:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",NULL);
425:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",NULL);
426:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpidense_C",NULL);
427:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqdense_C",NULL);
428:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",NULL);
429:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",NULL);
430:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",NULL);
431:   return 0;
432: }

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

439:   if (dd) *dd = 0;
440:   if (!vs->nr) {
441:     *missing = PETSC_TRUE;
442:     return 0;
443:   }
444:   *missing = PETSC_FALSE;
445:   for (i = 0; i < vs->nr && !(*missing); i++) {
446:     *missing = PETSC_TRUE;
447:     if (vs->m[i][i]) {
448:       MatMissingDiagonal(vs->m[i][i],missing,NULL);
450:     }
451:   }
452:   return 0;
453: }

455: static PetscErrorCode MatAssemblyBegin_Nest(Mat A,MatAssemblyType type)
456: {
457:   Mat_Nest       *vs = (Mat_Nest*)A->data;
458:   PetscInt       i,j;
459:   PetscBool      nnzstate = PETSC_FALSE;

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

485: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
486: {
487:   Mat_Nest       *vs = (Mat_Nest*)A->data;
488:   PetscInt       i,j;

490:   for (i=0; i<vs->nr; i++) {
491:     for (j=0; j<vs->nc; j++) {
492:       if (vs->m[i][j]) {
493:         if (vs->splitassembly) {
494:           MatAssemblyEnd(vs->m[i][j],type);
495:         }
496:       }
497:     }
498:   }
499:   return 0;
500: }

502: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A,PetscInt row,Mat *B)
503: {
504:   Mat_Nest       *vs = (Mat_Nest*)A->data;
505:   PetscInt       j;
506:   Mat            sub;

508:   sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
509:   for (j=0; !sub && j<vs->nc; j++) sub = vs->m[row][j];
510:   if (sub) MatSetUp(sub);       /* Ensure that the sizes are available */
511:   *B = sub;
512:   return 0;
513: }

515: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A,PetscInt col,Mat *B)
516: {
517:   Mat_Nest       *vs = (Mat_Nest*)A->data;
518:   PetscInt       i;
519:   Mat            sub;

521:   sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
522:   for (i=0; !sub && i<vs->nr; i++) sub = vs->m[i][col];
523:   if (sub) MatSetUp(sub);       /* Ensure that the sizes are available */
524:   *B = sub;
525:   return 0;
526: }

528: static PetscErrorCode MatNestFindISRange(Mat A,PetscInt n,const IS list[],IS is,PetscInt *begin,PetscInt *end)
529: {
530:   PetscInt       i,j,size,m;
531:   PetscBool      flg;
532:   IS             out,concatenate[2];

536:   if (begin) {
538:     *begin = -1;
539:   }
540:   if (end) {
542:     *end = -1;
543:   }
544:   for (i=0; i<n; i++) {
545:     if (!list[i]) continue;
546:     ISEqualUnsorted(list[i],is,&flg);
547:     if (flg) {
548:       if (begin) *begin = i;
549:       if (end) *end = i+1;
550:       return 0;
551:     }
552:   }
553:   ISGetSize(is,&size);
554:   for (i=0; i<n-1; i++) {
555:     if (!list[i]) continue;
556:     m = 0;
557:     ISConcatenate(PetscObjectComm((PetscObject)A),2,list+i,&out);
558:     ISGetSize(out,&m);
559:     for (j=i+2; j<n && m<size; j++) {
560:       if (list[j]) {
561:         concatenate[0] = out;
562:         concatenate[1] = list[j];
563:         ISConcatenate(PetscObjectComm((PetscObject)A),2,concatenate,&out);
564:         ISDestroy(concatenate);
565:         ISGetSize(out,&m);
566:       }
567:     }
568:     if (m == size) {
569:       ISEqualUnsorted(out,is,&flg);
570:       if (flg) {
571:         if (begin) *begin = i;
572:         if (end) *end = j;
573:         ISDestroy(&out);
574:         return 0;
575:       }
576:     }
577:     ISDestroy(&out);
578:   }
579:   return 0;
580: }

582: static PetscErrorCode MatNestFillEmptyMat_Private(Mat A,PetscInt i,PetscInt j,Mat *B)
583: {
584:   Mat_Nest       *vs = (Mat_Nest*)A->data;
585:   PetscInt       lr,lc;

587:   MatCreate(PetscObjectComm((PetscObject)A),B);
588:   ISGetLocalSize(vs->isglobal.row[i],&lr);
589:   ISGetLocalSize(vs->isglobal.col[j],&lc);
590:   MatSetSizes(*B,lr,lc,PETSC_DECIDE,PETSC_DECIDE);
591:   MatSetType(*B,MATAIJ);
592:   MatSeqAIJSetPreallocation(*B,0,NULL);
593:   MatMPIAIJSetPreallocation(*B,0,NULL,0,NULL);
594:   MatSetUp(*B);
595:   MatSetOption(*B,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);
596:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
597:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
598:   return 0;
599: }

601: static PetscErrorCode MatNestGetBlock_Private(Mat A,PetscInt rbegin,PetscInt rend,PetscInt cbegin,PetscInt cend,Mat *B)
602: {
603:   Mat_Nest       *vs = (Mat_Nest*)A->data;
604:   Mat            *a;
605:   PetscInt       i,j,k,l,nr=rend-rbegin,nc=cend-cbegin;
606:   char           keyname[256];
607:   PetscBool      *b;
608:   PetscBool      flg;

610:   *B   = NULL;
611:   PetscSNPrintf(keyname,sizeof(keyname),"NestBlock_%" PetscInt_FMT "-%" PetscInt_FMT "x%" PetscInt_FMT "-%" PetscInt_FMT,rbegin,rend,cbegin,cend);
612:   PetscObjectQuery((PetscObject)A,keyname,(PetscObject*)B);
613:   if (*B) return 0;

615:   PetscMalloc2(nr*nc,&a,nr*nc,&b);
616:   for (i=0; i<nr; i++) {
617:     for (j=0; j<nc; j++) {
618:       a[i*nc + j] = vs->m[rbegin+i][cbegin+j];
619:       b[i*nc + j] = PETSC_FALSE;
620:     }
621:   }
622:   if (nc!=vs->nc&&nr!=vs->nr) {
623:     for (i=0; i<nr; i++) {
624:       for (j=0; j<nc; j++) {
625:         flg = PETSC_FALSE;
626:         for (k=0; (k<nr&&!flg); k++) {
627:           if (a[j + k*nc]) flg = PETSC_TRUE;
628:         }
629:         if (flg) {
630:           flg = PETSC_FALSE;
631:           for (l=0; (l<nc&&!flg); l++) {
632:             if (a[i*nc + l]) flg = PETSC_TRUE;
633:           }
634:         }
635:         if (!flg) {
636:           b[i*nc + j] = PETSC_TRUE;
637:           MatNestFillEmptyMat_Private(A,rbegin+i,cbegin+j,a + i*nc + j);
638:         }
639:       }
640:     }
641:   }
642:   MatCreateNest(PetscObjectComm((PetscObject)A),nr,nr!=vs->nr?NULL:vs->isglobal.row,nc,nc!=vs->nc?NULL:vs->isglobal.col,a,B);
643:   for (i=0; i<nr; i++) {
644:     for (j=0; j<nc; j++) {
645:       if (b[i*nc + j]) {
646:         MatDestroy(a + i*nc + j);
647:       }
648:     }
649:   }
650:   PetscFree2(a,b);
651:   (*B)->assembled = A->assembled;
652:   PetscObjectCompose((PetscObject)A,keyname,(PetscObject)*B);
653:   PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
654:   return 0;
655: }

657: static PetscErrorCode MatNestFindSubMat(Mat A,struct MatNestISPair *is,IS isrow,IS iscol,Mat *B)
658: {
659:   Mat_Nest       *vs = (Mat_Nest*)A->data;
660:   PetscInt       rbegin,rend,cbegin,cend;

662:   MatNestFindISRange(A,vs->nr,is->row,isrow,&rbegin,&rend);
663:   MatNestFindISRange(A,vs->nc,is->col,iscol,&cbegin,&cend);
664:   if (rend == rbegin + 1 && cend == cbegin + 1) {
665:     if (!vs->m[rbegin][cbegin]) {
666:       MatNestFillEmptyMat_Private(A,rbegin,cbegin,vs->m[rbegin] + cbegin);
667:     }
668:     *B = vs->m[rbegin][cbegin];
669:   } else if (rbegin != -1 && cbegin != -1) {
670:     MatNestGetBlock_Private(A,rbegin,rend,cbegin,cend,B);
671:   } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_INCOMP,"Could not find index set");
672:   return 0;
673: }

675: /*
676:    TODO: This does not actually returns a submatrix we can modify
677: */
678: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A,IS isrow,IS iscol,MatReuse reuse,Mat *B)
679: {
680:   Mat_Nest       *vs = (Mat_Nest*)A->data;
681:   Mat            sub;

683:   MatNestFindSubMat(A,&vs->isglobal,isrow,iscol,&sub);
684:   switch (reuse) {
685:   case MAT_INITIAL_MATRIX:
686:     if (sub) PetscObjectReference((PetscObject)sub);
687:     *B = sub;
688:     break;
689:   case MAT_REUSE_MATRIX:
691:     break;
692:   case MAT_IGNORE_MATRIX:       /* Nothing to do */
693:     break;
694:   case MAT_INPLACE_MATRIX:       /* Nothing to do */
695:     SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_INPLACE_MATRIX is not supported yet");
696:   }
697:   return 0;
698: }

700: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
701: {
702:   Mat_Nest       *vs = (Mat_Nest*)A->data;
703:   Mat            sub;

705:   MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
706:   /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
707:   if (sub) PetscObjectReference((PetscObject)sub);
708:   *B = sub;
709:   return 0;
710: }

712: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A,IS isrow,IS iscol,Mat *B)
713: {
714:   Mat_Nest       *vs = (Mat_Nest*)A->data;
715:   Mat            sub;

717:   MatNestFindSubMat(A,&vs->islocal,isrow,iscol,&sub);
719:   if (sub) {
721:     MatDestroy(B);
722:   }
723:   return 0;
724: }

726: static PetscErrorCode MatGetDiagonal_Nest(Mat A,Vec v)
727: {
728:   Mat_Nest       *bA = (Mat_Nest*)A->data;
729:   PetscInt       i;

731:   for (i=0; i<bA->nr; i++) {
732:     Vec bv;
733:     VecGetSubVector(v,bA->isglobal.row[i],&bv);
734:     if (bA->m[i][i]) {
735:       MatGetDiagonal(bA->m[i][i],bv);
736:     } else {
737:       VecSet(bv,0.0);
738:     }
739:     VecRestoreSubVector(v,bA->isglobal.row[i],&bv);
740:   }
741:   return 0;
742: }

744: static PetscErrorCode MatDiagonalScale_Nest(Mat A,Vec l,Vec r)
745: {
746:   Mat_Nest       *bA = (Mat_Nest*)A->data;
747:   Vec            bl,*br;
748:   PetscInt       i,j;

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

775: static PetscErrorCode MatScale_Nest(Mat A,PetscScalar a)
776: {
777:   Mat_Nest       *bA = (Mat_Nest*)A->data;
778:   PetscInt       i,j;

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

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

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

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

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

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

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

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

852:   PetscObjectGetComm((PetscObject)A,&comm);
853:   if (right) {
854:     /* allocate R */
855:     PetscMalloc1(bA->nc, &R);
856:     /* Create the right vectors */
857:     for (j=0; j<bA->nc; j++) {
858:       for (i=0; i<bA->nr; i++) {
859:         if (bA->m[i][j]) {
860:           MatCreateVecs(bA->m[i][j],&R[j],NULL);
861:           break;
862:         }
863:       }
865:     }
866:     VecCreateNest(comm,bA->nc,bA->isglobal.col,R,right);
867:     /* hand back control to the nest vector */
868:     for (j=0; j<bA->nc; j++) {
869:       VecDestroy(&R[j]);
870:     }
871:     PetscFree(R);
872:   }

874:   if (left) {
875:     /* allocate L */
876:     PetscMalloc1(bA->nr, &L);
877:     /* Create the left vectors */
878:     for (i=0; i<bA->nr; i++) {
879:       for (j=0; j<bA->nc; j++) {
880:         if (bA->m[i][j]) {
881:           MatCreateVecs(bA->m[i][j],NULL,&L[i]);
882:           break;
883:         }
884:       }
886:     }

888:     VecCreateNest(comm,bA->nr,bA->isglobal.row,L,left);
889:     for (i=0; i<bA->nr; i++) {
890:       VecDestroy(&L[i]);
891:     }

893:     PetscFree(L);
894:   }
895:   return 0;
896: }

898: static PetscErrorCode MatView_Nest(Mat A,PetscViewer viewer)
899: {
900:   Mat_Nest       *bA = (Mat_Nest*)A->data;
901:   PetscBool      isascii,viewSub = PETSC_FALSE;
902:   PetscInt       i,j;

904:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
905:   if (isascii) {

907:     PetscOptionsGetBool(((PetscObject)A)->options,((PetscObject)A)->prefix,"-mat_view_nest_sub",&viewSub,NULL);
908:     PetscViewerASCIIPrintf(viewer,"Matrix object: \n");
909:     PetscViewerASCIIPushTab(viewer);
910:     PetscViewerASCIIPrintf(viewer, "type=nest, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT " \n",bA->nr,bA->nc);

912:     PetscViewerASCIIPrintf(viewer,"MatNest structure: \n");
913:     for (i=0; i<bA->nr; i++) {
914:       for (j=0; j<bA->nc; j++) {
915:         MatType   type;
916:         char      name[256] = "",prefix[256] = "";
917:         PetscInt  NR,NC;
918:         PetscBool isNest = PETSC_FALSE;

920:         if (!bA->m[i][j]) {
921:           PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ",%" PetscInt_FMT ") : NULL \n",i,j);
922:           continue;
923:         }
924:         MatGetSize(bA->m[i][j],&NR,&NC);
925:         MatGetType(bA->m[i][j], &type);
926:         if (((PetscObject)bA->m[i][j])->name) PetscSNPrintf(name,sizeof(name),"name=\"%s\", ",((PetscObject)bA->m[i][j])->name);
927:         if (((PetscObject)bA->m[i][j])->prefix) PetscSNPrintf(prefix,sizeof(prefix),"prefix=\"%s\", ",((PetscObject)bA->m[i][j])->prefix);
928:         PetscObjectTypeCompare((PetscObject)bA->m[i][j],MATNEST,&isNest);

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

932:         if (isNest || viewSub) {
933:           PetscViewerASCIIPushTab(viewer);  /* push1 */
934:           MatView(bA->m[i][j],viewer);
935:           PetscViewerASCIIPopTab(viewer);    /* pop1 */
936:         }
937:       }
938:     }
939:     PetscViewerASCIIPopTab(viewer);    /* pop0 */
940:   }
941:   return 0;
942: }

944: static PetscErrorCode MatZeroEntries_Nest(Mat A)
945: {
946:   Mat_Nest       *bA = (Mat_Nest*)A->data;
947:   PetscInt       i,j;

949:   for (i=0; i<bA->nr; i++) {
950:     for (j=0; j<bA->nc; j++) {
951:       if (!bA->m[i][j]) continue;
952:       MatZeroEntries(bA->m[i][j]);
953:     }
954:   }
955:   return 0;
956: }

958: static PetscErrorCode MatCopy_Nest(Mat A,Mat B,MatStructure str)
959: {
960:   Mat_Nest       *bA = (Mat_Nest*)A->data,*bB = (Mat_Nest*)B->data;
961:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;
962:   PetscBool      nnzstate = PETSC_FALSE;

965:   for (i=0; i<nr; i++) {
966:     for (j=0; j<nc; j++) {
967:       PetscObjectState subnnzstate = 0;
968:       if (bA->m[i][j] && bB->m[i][j]) {
969:         MatCopy(bA->m[i][j],bB->m[i][j],str);
971:       MatGetNonzeroState(bB->m[i][j],&subnnzstate);
972:       nnzstate = (PetscBool)(nnzstate || bB->nnzstate[i*nc+j] != subnnzstate);
973:       bB->nnzstate[i*nc+j] = subnnzstate;
974:     }
975:   }
976:   if (nnzstate) B->nonzerostate++;
977:   return 0;
978: }

980: static PetscErrorCode MatAXPY_Nest(Mat Y,PetscScalar a,Mat X,MatStructure str)
981: {
982:   Mat_Nest       *bY = (Mat_Nest*)Y->data,*bX = (Mat_Nest*)X->data;
983:   PetscInt       i,j,nr = bY->nr,nc = bY->nc;
984:   PetscBool      nnzstate = PETSC_FALSE;

987:   for (i=0; i<nr; i++) {
988:     for (j=0; j<nc; j++) {
989:       PetscObjectState subnnzstate = 0;
990:       if (bY->m[i][j] && bX->m[i][j]) {
991:         MatAXPY(bY->m[i][j],a,bX->m[i][j],str);
992:       } else if (bX->m[i][j]) {
993:         Mat M;

996:         MatDuplicate(bX->m[i][j],MAT_COPY_VALUES,&M);
997:         MatNestSetSubMat(Y,i,j,M);
998:         MatDestroy(&M);
999:       }
1000:       if (bY->m[i][j]) MatGetNonzeroState(bY->m[i][j],&subnnzstate);
1001:       nnzstate = (PetscBool)(nnzstate || bY->nnzstate[i*nc+j] != subnnzstate);
1002:       bY->nnzstate[i*nc+j] = subnnzstate;
1003:     }
1004:   }
1005:   if (nnzstate) Y->nonzerostate++;
1006:   return 0;
1007: }

1009: static PetscErrorCode MatDuplicate_Nest(Mat A,MatDuplicateOption op,Mat *B)
1010: {
1011:   Mat_Nest       *bA = (Mat_Nest*)A->data;
1012:   Mat            *b;
1013:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;

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

1032:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
1033:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
1034:   return 0;
1035: }

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

1044:   *mat = bA->m[idxm][jdxm];
1045:   return 0;
1046: }

1048: /*@
1049:  MatNestGetSubMat - Returns a single, sub-matrix from a nest matrix.

1051:  Not collective

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

1058:  Output Parameter:
1059: .   sub - matrix at index idxm,jdxm within the nest matrix

1061:  Level: developer

1063: .seealso: `MatNestGetSize()`, `MatNestGetSubMats()`, `MatCreateNest()`, `MATNEST`, `MatNestSetSubMat()`,
1064:           `MatNestGetLocalISs()`, `MatNestGetISs()`
1065: @*/
1066: PetscErrorCode  MatNestGetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat *sub)
1067: {
1068:   PetscUseMethod(A,"MatNestGetSubMat_C",(Mat,PetscInt,PetscInt,Mat*),(A,idxm,jdxm,sub));
1069:   return 0;
1070: }

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

1079:   MatGetLocalSize(mat,&m,&n);
1080:   MatGetSize(mat,&M,&N);
1081:   ISGetLocalSize(bA->isglobal.row[idxm],&mi);
1082:   ISGetSize(bA->isglobal.row[idxm],&Mi);
1083:   ISGetLocalSize(bA->isglobal.col[jdxm],&ni);
1084:   ISGetSize(bA->isglobal.col[jdxm],&Ni);

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

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

1100: /*@
1101:  MatNestSetSubMat - Set a single submatrix in the nest matrix.

1103:  Logically collective on the submatrix communicator

1105:  Input Parameters:
1106: +   A  - nest matrix
1107: .   idxm - index of the matrix within the nest matrix
1108: .   jdxm - index of the matrix within the nest matrix
1109: -   sub - matrix at index idxm,jdxm within the nest matrix

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

1114:  This increments the reference count of the submatrix.

1116:  Level: developer

1118: .seealso: `MatNestSetSubMats()`, `MatNestGetSubMats()`, `MatNestGetLocalISs()`, `MATNEST`, `MatCreateNest()`,
1119:           `MatNestGetSubMat()`, `MatNestGetISs()`, `MatNestGetSize()`
1120: @*/
1121: PetscErrorCode  MatNestSetSubMat(Mat A,PetscInt idxm,PetscInt jdxm,Mat sub)
1122: {
1123:   PetscUseMethod(A,"MatNestSetSubMat_C",(Mat,PetscInt,PetscInt,Mat),(A,idxm,jdxm,sub));
1124:   return 0;
1125: }

1127: PetscErrorCode MatNestGetSubMats_Nest(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1128: {
1129:   Mat_Nest *bA = (Mat_Nest*)A->data;

1131:   if (M)   *M   = bA->nr;
1132:   if (N)   *N   = bA->nc;
1133:   if (mat) *mat = bA->m;
1134:   return 0;
1135: }

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

1140:  Not collective

1142:  Input Parameter:
1143: .   A  - nest matrix

1145:  Output Parameters:
1146: +   M - number of rows in the nest matrix
1147: .   N - number of cols in the nest matrix
1148: -   mat - 2d array of matrices

1150:  Notes:

1152:  The user should not free the array mat.

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

1158:  Level: developer

1160: .seealso: `MatNestGetSize()`, `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MATNEST`, `MatCreateNest()`,
1161:           `MatNestSetSubMats()`, `MatNestGetISs()`, `MatNestSetSubMat()`
1162: @*/
1163: PetscErrorCode  MatNestGetSubMats(Mat A,PetscInt *M,PetscInt *N,Mat ***mat)
1164: {
1165:   PetscUseMethod(A,"MatNestGetSubMats_C",(Mat,PetscInt*,PetscInt*,Mat***),(A,M,N,mat));
1166:   return 0;
1167: }

1169: PetscErrorCode  MatNestGetSize_Nest(Mat A,PetscInt *M,PetscInt *N)
1170: {
1171:   Mat_Nest *bA = (Mat_Nest*)A->data;

1173:   if (M) *M = bA->nr;
1174:   if (N) *N = bA->nc;
1175:   return 0;
1176: }

1178: /*@
1179:  MatNestGetSize - Returns the size of the nest matrix.

1181:  Not collective

1183:  Input Parameter:
1184: .   A  - nest matrix

1186:  Output Parameters:
1187: +   M - number of rows in the nested mat
1188: -   N - number of cols in the nested mat

1190:  Notes:

1192:  Level: developer

1194: .seealso: `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MATNEST`, `MatCreateNest()`, `MatNestGetLocalISs()`,
1195:           `MatNestGetISs()`
1196: @*/
1197: PetscErrorCode  MatNestGetSize(Mat A,PetscInt *M,PetscInt *N)
1198: {
1199:   PetscUseMethod(A,"MatNestGetSize_C",(Mat,PetscInt*,PetscInt*),(A,M,N));
1200:   return 0;
1201: }

1203: static PetscErrorCode MatNestGetISs_Nest(Mat A,IS rows[],IS cols[])
1204: {
1205:   Mat_Nest *vs = (Mat_Nest*)A->data;
1206:   PetscInt i;

1208:   if (rows) for (i=0; i<vs->nr; i++) rows[i] = vs->isglobal.row[i];
1209:   if (cols) for (i=0; i<vs->nc; i++) cols[i] = vs->isglobal.col[i];
1210:   return 0;
1211: }

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

1216:  Not collective

1218:  Input Parameter:
1219: .   A  - nest matrix

1221:  Output Parameters:
1222: +   rows - array of row index sets
1223: -   cols - array of column index sets

1225:  Level: advanced

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

1230: .seealso: `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetLocalISs()`, `MATNEST`,
1231:           `MatCreateNest()`, `MatNestGetSubMats()`, `MatNestSetSubMats()`
1232: @*/
1233: PetscErrorCode  MatNestGetISs(Mat A,IS rows[],IS cols[])
1234: {
1236:   PetscUseMethod(A,"MatNestGetISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1237:   return 0;
1238: }

1240: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A,IS rows[],IS cols[])
1241: {
1242:   Mat_Nest *vs = (Mat_Nest*)A->data;
1243:   PetscInt i;

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

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

1253:  Not collective

1255:  Input Parameter:
1256: .   A  - nest matrix

1258:  Output Parameters:
1259: +   rows - array of row index sets (or NULL to ignore)
1260: -   cols - array of column index sets (or NULL to ignore)

1262:  Level: advanced

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

1267: .seealso: `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetISs()`, `MatCreateNest()`,
1268:           `MATNEST`, `MatNestSetSubMats()`, `MatNestSetSubMat()`
1269: @*/
1270: PetscErrorCode  MatNestGetLocalISs(Mat A,IS rows[],IS cols[])
1271: {
1273:   PetscUseMethod(A,"MatNestGetLocalISs_C",(Mat,IS[],IS[]),(A,rows,cols));
1274:   return 0;
1275: }

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

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

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

1291:  Not collective

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

1297:  Notes:

1299:  Level: developer

1301: .seealso: `MatCreateVecs()`, `MATNEST`, `MatCreateNest()`
1302: @*/
1303: PetscErrorCode  MatNestSetVecType(Mat A,VecType vtype)
1304: {
1305:   PetscTryMethod(A,"MatNestSetVecType_C",(Mat,VecType),(A,vtype));
1306:   return 0;
1307: }

1309: PetscErrorCode MatNestSetSubMats_Nest(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1310: {
1311:   Mat_Nest       *s = (Mat_Nest*)A->data;
1312:   PetscInt       i,j,m,n,M,N;
1313:   PetscBool      cong,isstd,sametype=PETSC_FALSE;
1314:   VecType        vtype,type;

1316:   MatReset_Nest(A);

1318:   s->nr = nr;
1319:   s->nc = nc;

1321:   /* Create space for submatrices */
1322:   PetscMalloc1(nr,&s->m);
1323:   for (i=0; i<nr; i++) {
1324:     PetscMalloc1(nc,&s->m[i]);
1325:   }
1326:   for (i=0; i<nr; i++) {
1327:     for (j=0; j<nc; j++) {
1328:       s->m[i][j] = a[i*nc+j];
1329:       if (a[i*nc+j]) {
1330:         PetscObjectReference((PetscObject)a[i*nc+j]);
1331:       }
1332:     }
1333:   }
1334:   MatGetVecType(A,&vtype);
1335:   PetscStrcmp(vtype,VECSTANDARD,&isstd);
1336:   if (isstd) {
1337:     /* check if all blocks have the same vectype */
1338:     vtype = NULL;
1339:     for (i=0; i<nr; i++) {
1340:       for (j=0; j<nc; j++) {
1341:         if (a[i*nc+j]) {
1342:           if (!vtype) {  /* first visited block */
1343:             MatGetVecType(a[i*nc+j],&vtype);
1344:             sametype = PETSC_TRUE;
1345:           } else if (sametype) {
1346:             MatGetVecType(a[i*nc+j],&type);
1347:             PetscStrcmp(vtype,type,&sametype);
1348:           }
1349:         }
1350:       }
1351:     }
1352:     if (sametype) {  /* propagate vectype */
1353:       MatSetVecType(A,vtype);
1354:     }
1355:   }

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

1359:   PetscMalloc1(nr,&s->row_len);
1360:   PetscMalloc1(nc,&s->col_len);
1361:   for (i=0; i<nr; i++) s->row_len[i]=-1;
1362:   for (j=0; j<nc; j++) s->col_len[j]=-1;

1364:   PetscCalloc1(nr*nc,&s->nnzstate);
1365:   for (i=0; i<nr; i++) {
1366:     for (j=0; j<nc; j++) {
1367:       if (s->m[i][j]) {
1368:         MatGetNonzeroState(s->m[i][j],&s->nnzstate[i*nc+j]);
1369:       }
1370:     }
1371:   }

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

1375:   PetscLayoutSetSize(A->rmap,M);
1376:   PetscLayoutSetLocalSize(A->rmap,m);
1377:   PetscLayoutSetSize(A->cmap,N);
1378:   PetscLayoutSetLocalSize(A->cmap,n);

1380:   PetscLayoutSetUp(A->rmap);
1381:   PetscLayoutSetUp(A->cmap);

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

1399:   PetscCalloc2(nr,&s->left,nc,&s->right);
1400:   PetscObjectStateIncrease((PetscObject)A);
1401:   A->nonzerostate++;
1402:   return 0;
1403: }

1405: /*@
1406:    MatNestSetSubMats - Sets the nested submatrices

1408:    Collective on Mat

1410:    Input Parameters:
1411: +  A - nested matrix
1412: .  nr - number of nested row blocks
1413: .  is_row - index sets for each nested row block, or NULL to make contiguous
1414: .  nc - number of nested column blocks
1415: .  is_col - index sets for each nested column block, or NULL to make contiguous
1416: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

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

1420:    Level: advanced

1422: .seealso: `MatCreateNest()`, `MATNEST`, `MatNestSetSubMat()`, `MatNestGetSubMat()`, `MatNestGetSubMats()`
1423: @*/
1424: PetscErrorCode MatNestSetSubMats(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[])
1425: {
1426:   PetscInt       i;

1430:   if (nr && is_row) {
1433:   }
1435:   if (nc && is_col) {
1438:   }
1440:   PetscUseMethod(A,"MatNestSetSubMats_C",(Mat,PetscInt,const IS[],PetscInt,const IS[],const Mat[]),(A,nr,is_row,nc,is_col,a));
1441:   return 0;
1442: }

1444: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A,PetscInt n,const IS islocal[],const IS isglobal[],PetscBool colflg,ISLocalToGlobalMapping *ltog)
1445: {
1446:   PetscBool      flg;
1447:   PetscInt       i,j,m,mi,*ix;

1449:   *ltog = NULL;
1450:   for (i=0,m=0,flg=PETSC_FALSE; i<n; i++) {
1451:     if (islocal[i]) {
1452:       ISGetLocalSize(islocal[i],&mi);
1453:       flg  = PETSC_TRUE;      /* We found a non-trivial entry */
1454:     } else {
1455:       ISGetLocalSize(isglobal[i],&mi);
1456:     }
1457:     m += mi;
1458:   }
1459:   if (!flg) return 0;

1461:   PetscMalloc1(m,&ix);
1462:   for (i=0,m=0; i<n; i++) {
1463:     ISLocalToGlobalMapping smap = NULL;
1464:     Mat                    sub = NULL;
1465:     PetscSF                sf;
1466:     PetscLayout            map;
1467:     const PetscInt         *ix2;

1469:     if (!colflg) {
1470:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1471:     } else {
1472:       MatNestFindNonzeroSubMatCol(A,i,&sub);
1473:     }
1474:     if (sub) {
1475:       if (!colflg) {
1476:         MatGetLocalToGlobalMapping(sub,&smap,NULL);
1477:       } else {
1478:         MatGetLocalToGlobalMapping(sub,NULL,&smap);
1479:       }
1480:     }
1481:     /*
1482:        Now we need to extract the monolithic global indices that correspond to the given split global indices.
1483:        In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1484:     */
1485:     ISGetIndices(isglobal[i],&ix2);
1486:     if (islocal[i]) {
1487:       PetscInt *ilocal,*iremote;
1488:       PetscInt mil,nleaves;

1490:       ISGetLocalSize(islocal[i],&mi);
1492:       for (j=0; j<mi; j++) ix[m+j] = j;
1493:       ISLocalToGlobalMappingApply(smap,mi,ix+m,ix+m);

1495:       /* PetscSFSetGraphLayout does not like negative indices */
1496:       PetscMalloc2(mi,&ilocal,mi,&iremote);
1497:       for (j=0, nleaves = 0; j<mi; j++) {
1498:         if (ix[m+j] < 0) continue;
1499:         ilocal[nleaves]  = j;
1500:         iremote[nleaves] = ix[m+j];
1501:         nleaves++;
1502:       }
1503:       ISGetLocalSize(isglobal[i],&mil);
1504:       PetscSFCreate(PetscObjectComm((PetscObject)A),&sf);
1505:       PetscLayoutCreate(PetscObjectComm((PetscObject)A),&map);
1506:       PetscLayoutSetLocalSize(map,mil);
1507:       PetscLayoutSetUp(map);
1508:       PetscSFSetGraphLayout(sf,map,nleaves,ilocal,PETSC_USE_POINTER,iremote);
1509:       PetscLayoutDestroy(&map);
1510:       PetscSFBcastBegin(sf,MPIU_INT,ix2,ix + m,MPI_REPLACE);
1511:       PetscSFBcastEnd(sf,MPIU_INT,ix2,ix + m,MPI_REPLACE);
1512:       PetscSFDestroy(&sf);
1513:       PetscFree2(ilocal,iremote);
1514:     } else {
1515:       ISGetLocalSize(isglobal[i],&mi);
1516:       for (j=0; j<mi; j++) ix[m+j] = ix2[i];
1517:     }
1518:     ISRestoreIndices(isglobal[i],&ix2);
1519:     m   += mi;
1520:   }
1521:   ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A),1,m,ix,PETSC_OWN_POINTER,ltog);
1522:   return 0;
1523: }

1525: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1526: /*
1527:   nprocessors = NP
1528:   Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1529:        proc 0: => (g_0,h_0,)
1530:        proc 1: => (g_1,h_1,)
1531:        ...
1532:        proc nprocs-1: => (g_NP-1,h_NP-1,)

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

1537:             proc 0:
1538:     is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1539:             proc 1:
1540:     is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)

1542:             proc NP-1:
1543:     is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1544: */
1545: static PetscErrorCode MatSetUp_NestIS_Private(Mat A,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[])
1546: {
1547:   Mat_Nest       *vs = (Mat_Nest*)A->data;
1548:   PetscInt       i,j,offset,n,nsum,bs;
1549:   Mat            sub = NULL;

1551:   PetscMalloc1(nr,&vs->isglobal.row);
1552:   PetscMalloc1(nc,&vs->isglobal.col);
1553:   if (is_row) { /* valid IS is passed in */
1554:     /* refs on is[] are incremented */
1555:     for (i=0; i<vs->nr; i++) {
1556:       PetscObjectReference((PetscObject)is_row[i]);

1558:       vs->isglobal.row[i] = is_row[i];
1559:     }
1560:   } else {                      /* Create the ISs by inspecting sizes of a submatrix in each row */
1561:     nsum = 0;
1562:     for (i=0; i<vs->nr; i++) {  /* Add up the local sizes to compute the aggregate offset */
1563:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1565:       MatGetLocalSize(sub,&n,NULL);
1567:       nsum += n;
1568:     }
1569:     MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1570:     offset -= nsum;
1571:     for (i=0; i<vs->nr; i++) {
1572:       MatNestFindNonzeroSubMatRow(A,i,&sub);
1573:       MatGetLocalSize(sub,&n,NULL);
1574:       MatGetBlockSizes(sub,&bs,NULL);
1575:       ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.row[i]);
1576:       ISSetBlockSize(vs->isglobal.row[i],bs);
1577:       offset += n;
1578:     }
1579:   }

1581:   if (is_col) { /* valid IS is passed in */
1582:     /* refs on is[] are incremented */
1583:     for (j=0; j<vs->nc; j++) {
1584:       PetscObjectReference((PetscObject)is_col[j]);

1586:       vs->isglobal.col[j] = is_col[j];
1587:     }
1588:   } else {                      /* Create the ISs by inspecting sizes of a submatrix in each column */
1589:     offset = A->cmap->rstart;
1590:     nsum   = 0;
1591:     for (j=0; j<vs->nc; j++) {
1592:       MatNestFindNonzeroSubMatCol(A,j,&sub);
1594:       MatGetLocalSize(sub,NULL,&n);
1596:       nsum += n;
1597:     }
1598:     MPI_Scan(&nsum,&offset,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)A));
1599:     offset -= nsum;
1600:     for (j=0; j<vs->nc; j++) {
1601:       MatNestFindNonzeroSubMatCol(A,j,&sub);
1602:       MatGetLocalSize(sub,NULL,&n);
1603:       MatGetBlockSizes(sub,NULL,&bs);
1604:       ISCreateStride(PetscObjectComm((PetscObject)sub),n,offset,1,&vs->isglobal.col[j]);
1605:       ISSetBlockSize(vs->isglobal.col[j],bs);
1606:       offset += n;
1607:     }
1608:   }

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

1650:   /* Set up the aggregate ISLocalToGlobalMapping */
1651:   {
1652:     ISLocalToGlobalMapping rmap,cmap;
1653:     MatNestCreateAggregateL2G_Private(A,vs->nr,vs->islocal.row,vs->isglobal.row,PETSC_FALSE,&rmap);
1654:     MatNestCreateAggregateL2G_Private(A,vs->nc,vs->islocal.col,vs->isglobal.col,PETSC_TRUE,&cmap);
1655:     if (rmap && cmap) MatSetLocalToGlobalMapping(A,rmap,cmap);
1656:     ISLocalToGlobalMappingDestroy(&rmap);
1657:     ISLocalToGlobalMappingDestroy(&cmap);
1658:   }

1660:   if (PetscDefined(USE_DEBUG)) {
1661:     for (i=0; i<vs->nr; i++) {
1662:       for (j=0; j<vs->nc; j++) {
1663:         PetscInt m,n,M,N,mi,ni,Mi,Ni;
1664:         Mat      B = vs->m[i][j];
1665:         if (!B) continue;
1666:         MatGetSize(B,&M,&N);
1667:         MatGetLocalSize(B,&m,&n);
1668:         ISGetSize(vs->isglobal.row[i],&Mi);
1669:         ISGetSize(vs->isglobal.col[j],&Ni);
1670:         ISGetLocalSize(vs->isglobal.row[i],&mi);
1671:         ISGetLocalSize(vs->isglobal.col[j],&ni);
1674:       }
1675:     }
1676:   }

1678:   /* Set A->assembled if all non-null blocks are currently assembled */
1679:   for (i=0; i<vs->nr; i++) {
1680:     for (j=0; j<vs->nc; j++) {
1681:       if (vs->m[i][j] && !vs->m[i][j]->assembled) return 0;
1682:     }
1683:   }
1684:   A->assembled = PETSC_TRUE;
1685:   return 0;
1686: }

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

1691:    Collective on Mat

1693:    Input Parameters:
1694: +  comm - Communicator for the new Mat
1695: .  nr - number of nested row blocks
1696: .  is_row - index sets for each nested row block, or NULL to make contiguous
1697: .  nc - number of nested column blocks
1698: .  is_col - index sets for each nested column block, or NULL to make contiguous
1699: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

1701:    Output Parameter:
1702: .  B - new matrix

1704:    Level: advanced

1706: .seealso: `MatCreate()`, `VecCreateNest()`, `DMCreateMatrix()`, `MATNEST`, `MatNestSetSubMat()`,
1707:           `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MatNestGetSize()`,
1708:           `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()`
1709: @*/
1710: PetscErrorCode MatCreateNest(MPI_Comm comm,PetscInt nr,const IS is_row[],PetscInt nc,const IS is_col[],const Mat a[],Mat *B)
1711: {
1712:   Mat            A;

1714:   *B   = NULL;
1715:   MatCreate(comm,&A);
1716:   MatSetType(A,MATNEST);
1717:   A->preallocated = PETSC_TRUE;
1718:   MatNestSetSubMats(A,nr,is_row,nc,is_col,a);
1719:   *B   = A;
1720:   return 0;
1721: }

1723: PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1724: {
1725:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1726:   Mat            *trans;
1727:   PetscScalar    **avv;
1728:   PetscScalar    *vv;
1729:   PetscInt       **aii,**ajj;
1730:   PetscInt       *ii,*jj,*ci;
1731:   PetscInt       nr,nc,nnz,i,j;
1732:   PetscBool      done;

1734:   MatGetSize(A,&nr,&nc);
1735:   if (reuse == MAT_REUSE_MATRIX) {
1736:     PetscInt rnr;

1738:     MatGetRowIJ(*newmat,0,PETSC_FALSE,PETSC_FALSE,&rnr,(const PetscInt**)&ii,(const PetscInt**)&jj,&done);
1741:     MatSeqAIJGetArray(*newmat,&vv);
1742:   }
1743:   /* extract CSR for nested SeqAIJ matrices */
1744:   nnz  = 0;
1745:   PetscCalloc4(nest->nr*nest->nc,&aii,nest->nr*nest->nc,&ajj,nest->nr*nest->nc,&avv,nest->nr*nest->nc,&trans);
1746:   for (i=0; i<nest->nr; ++i) {
1747:     for (j=0; j<nest->nc; ++j) {
1748:       Mat B = nest->m[i][j];
1749:       if (B) {
1750:         PetscScalar *naa;
1751:         PetscInt    *nii,*njj,nnr;
1752:         PetscBool   istrans;

1754:         PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1755:         if (istrans) {
1756:           Mat Bt;

1758:           MatTransposeGetMat(B,&Bt);
1759:           MatTranspose(Bt,MAT_INITIAL_MATRIX,&trans[i*nest->nc+j]);
1760:           B    = trans[i*nest->nc+j];
1761:         }
1762:         MatGetRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&nii,(const PetscInt**)&njj,&done);
1764:         MatSeqAIJGetArray(B,&naa);
1765:         nnz += nii[nnr];

1767:         aii[i*nest->nc+j] = nii;
1768:         ajj[i*nest->nc+j] = njj;
1769:         avv[i*nest->nc+j] = naa;
1770:       }
1771:     }
1772:   }
1773:   if (reuse != MAT_REUSE_MATRIX) {
1774:     PetscMalloc1(nr+1,&ii);
1775:     PetscMalloc1(nnz,&jj);
1776:     PetscMalloc1(nnz,&vv);
1777:   } else {
1779:   }

1781:   /* new row pointer */
1782:   PetscArrayzero(ii,nr+1);
1783:   for (i=0; i<nest->nr; ++i) {
1784:     PetscInt       ncr,rst;

1786:     ISStrideGetInfo(nest->isglobal.row[i],&rst,NULL);
1787:     ISGetLocalSize(nest->isglobal.row[i],&ncr);
1788:     for (j=0; j<nest->nc; ++j) {
1789:       if (aii[i*nest->nc+j]) {
1790:         PetscInt    *nii = aii[i*nest->nc+j];
1791:         PetscInt    ir;

1793:         for (ir=rst; ir<ncr+rst; ++ir) {
1794:           ii[ir+1] += nii[1]-nii[0];
1795:           nii++;
1796:         }
1797:       }
1798:     }
1799:   }
1800:   for (i=0; i<nr; i++) ii[i+1] += ii[i];

1802:   /* construct CSR for the new matrix */
1803:   PetscCalloc1(nr,&ci);
1804:   for (i=0; i<nest->nr; ++i) {
1805:     PetscInt       ncr,rst;

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

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

1820:           for (ij=0;ij<rsize;ij++) {
1821:             jj[ist+ij] = *njj+cst;
1822:             vv[ist+ij] = *nvv;
1823:             njj++;
1824:             nvv++;
1825:           }
1826:           ci[ir] += rsize;
1827:           nii++;
1828:         }
1829:       }
1830:     }
1831:   }
1832:   PetscFree(ci);

1834:   /* restore info */
1835:   for (i=0; i<nest->nr; ++i) {
1836:     for (j=0; j<nest->nc; ++j) {
1837:       Mat B = nest->m[i][j];
1838:       if (B) {
1839:         PetscInt nnr = 0, k = i*nest->nc+j;

1841:         B    = (trans[k] ? trans[k] : B);
1842:         MatRestoreRowIJ(B,0,PETSC_FALSE,PETSC_FALSE,&nnr,(const PetscInt**)&aii[k],(const PetscInt**)&ajj[k],&done);
1844:         MatSeqAIJRestoreArray(B,&avv[k]);
1845:         MatDestroy(&trans[k]);
1846:       }
1847:     }
1848:   }
1849:   PetscFree4(aii,ajj,avv,trans);

1851:   /* finalize newmat */
1852:   if (reuse == MAT_INITIAL_MATRIX) {
1853:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),nr,nc,ii,jj,vv,newmat);
1854:   } else if (reuse == MAT_INPLACE_MATRIX) {
1855:     Mat B;

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

1870: PETSC_INTERN PetscErrorCode MatAXPY_Dense_Nest(Mat Y,PetscScalar a,Mat X)
1871: {
1872:   Mat_Nest       *nest = (Mat_Nest*)X->data;
1873:   PetscInt       i,j,k,rstart;
1874:   PetscBool      flg;

1876:   /* Fill by row */
1877:   for (j=0; j<nest->nc; ++j) {
1878:     /* Using global column indices and ISAllGather() is not scalable. */
1879:     IS             bNis;
1880:     PetscInt       bN;
1881:     const PetscInt *bNindices;
1882:     ISAllGather(nest->isglobal.col[j], &bNis);
1883:     ISGetSize(bNis,&bN);
1884:     ISGetIndices(bNis,&bNindices);
1885:     for (i=0; i<nest->nr; ++i) {
1886:       Mat            B=nest->m[i][j],D=NULL;
1887:       PetscInt       bm,br;
1888:       const PetscInt *bmindices;
1889:       if (!B) continue;
1890:       PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&flg);
1891:       if (flg) {
1892:         PetscTryMethod(B,"MatTransposeGetMat_C",(Mat,Mat*),(B,&D));
1893:         PetscTryMethod(B,"MatHermitianTransposeGetMat_C",(Mat,Mat*),(B,&D));
1894:         MatConvert(B,((PetscObject)D)->type_name,MAT_INITIAL_MATRIX,&D);
1895:         B = D;
1896:       }
1897:       PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQSBAIJ,MATMPISBAIJ,"");
1898:       if (flg) {
1899:         if (D) MatConvert(D,MATBAIJ,MAT_INPLACE_MATRIX,&D);
1900:         else MatConvert(B,MATBAIJ,MAT_INITIAL_MATRIX,&D);
1901:         B = D;
1902:       }
1903:       ISGetLocalSize(nest->isglobal.row[i],&bm);
1904:       ISGetIndices(nest->isglobal.row[i],&bmindices);
1905:       MatGetOwnershipRange(B,&rstart,NULL);
1906:       for (br = 0; br < bm; ++br) {
1907:         PetscInt          row = bmindices[br], brncols, *cols;
1908:         const PetscInt    *brcols;
1909:         const PetscScalar *brcoldata;
1910:         PetscScalar       *vals = NULL;
1911:         MatGetRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1912:         PetscMalloc1(brncols,&cols);
1913:         for (k=0; k<brncols; k++) cols[k] = bNindices[brcols[k]];
1914:         /*
1915:           Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
1916:           Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
1917:          */
1918:         if (a != 1.0) {
1919:           PetscMalloc1(brncols,&vals);
1920:           for (k=0; k<brncols; k++) vals[k] = a * brcoldata[k];
1921:           MatSetValues(Y,1,&row,brncols,cols,vals,ADD_VALUES);
1922:           PetscFree(vals);
1923:         } else {
1924:           MatSetValues(Y,1,&row,brncols,cols,brcoldata,ADD_VALUES);
1925:         }
1926:         MatRestoreRow(B,br+rstart,&brncols,&brcols,&brcoldata);
1927:         PetscFree(cols);
1928:       }
1929:       if (D) MatDestroy(&D);
1930:       ISRestoreIndices(nest->isglobal.row[i],&bmindices);
1931:     }
1932:     ISRestoreIndices(bNis,&bNindices);
1933:     ISDestroy(&bNis);
1934:   }
1935:   MatAssemblyBegin(Y,MAT_FINAL_ASSEMBLY);
1936:   MatAssemblyEnd(Y,MAT_FINAL_ASSEMBLY);
1937:   return 0;
1938: }

1940: PetscErrorCode MatConvert_Nest_AIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
1941: {
1942:   Mat_Nest       *nest = (Mat_Nest*)A->data;
1943:   PetscInt       m,n,M,N,i,j,k,*dnnz,*onnz,rstart,cstart,cend;
1944:   PetscMPIInt    size;
1945:   Mat            C;

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

1952:     PetscStrcmp(newtype,MATAIJ,&fast);
1953:     if (!fast) {
1954:       PetscStrcmp(newtype,MATSEQAIJ,&fast);
1955:     }
1956:     for (i=0; i<nest->nr && fast; ++i) {
1957:       for (j=0; j<nest->nc && fast; ++j) {
1958:         Mat B = nest->m[i][j];
1959:         if (B) {
1960:           PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&fast);
1961:           if (!fast) {
1962:             PetscBool istrans;

1964:             PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&istrans);
1965:             if (istrans) {
1966:               Mat Bt;

1968:               MatTransposeGetMat(B,&Bt);
1969:               PetscObjectTypeCompare((PetscObject)Bt,MATSEQAIJ,&fast);
1970:             }
1971:           }
1972:         }
1973:       }
1974:     }
1975:     for (i=0, nf=0; i<nest->nr && fast; ++i) {
1976:       PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i],ISSTRIDE,&fast);
1977:       if (fast) {
1978:         PetscInt f,s;

1980:         ISStrideGetInfo(nest->isglobal.row[i],&f,&s);
1981:         if (f != nf || s != 1) { fast = PETSC_FALSE; }
1982:         else {
1983:           ISGetSize(nest->isglobal.row[i],&f);
1984:           nf  += f;
1985:         }
1986:       }
1987:     }
1988:     for (i=0, nf=0; i<nest->nc && fast; ++i) {
1989:       PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i],ISSTRIDE,&fast);
1990:       if (fast) {
1991:         PetscInt f,s;

1993:         ISStrideGetInfo(nest->isglobal.col[i],&f,&s);
1994:         if (f != nf || s != 1) { fast = PETSC_FALSE; }
1995:         else {
1996:           ISGetSize(nest->isglobal.col[i],&f);
1997:           nf  += f;
1998:         }
1999:       }
2000:     }
2001:     if (fast) {
2002:       MatConvert_Nest_SeqAIJ_fast(A,newtype,reuse,newmat);
2003:       return 0;
2004:     }
2005:   }
2006:   MatGetSize(A,&M,&N);
2007:   MatGetLocalSize(A,&m,&n);
2008:   MatGetOwnershipRangeColumn(A,&cstart,&cend);
2009:   if (reuse == MAT_REUSE_MATRIX) C = *newmat;
2010:   else {
2011:     MatCreate(PetscObjectComm((PetscObject)A),&C);
2012:     MatSetType(C,newtype);
2013:     MatSetSizes(C,m,n,M,N);
2014:   }
2015:   PetscMalloc1(2*m,&dnnz);
2016:   onnz = dnnz + m;
2017:   for (k=0; k<m; k++) {
2018:     dnnz[k] = 0;
2019:     onnz[k] = 0;
2020:   }
2021:   for (j=0; j<nest->nc; ++j) {
2022:     IS             bNis;
2023:     PetscInt       bN;
2024:     const PetscInt *bNindices;
2025:     PetscBool      flg;
2026:     /* Using global column indices and ISAllGather() is not scalable. */
2027:     ISAllGather(nest->isglobal.col[j], &bNis);
2028:     ISGetSize(bNis, &bN);
2029:     ISGetIndices(bNis,&bNindices);
2030:     for (i=0; i<nest->nr; ++i) {
2031:       PetscSF        bmsf;
2032:       PetscSFNode    *iremote;
2033:       Mat            B=nest->m[i][j],D=NULL;
2034:       PetscInt       bm,*sub_dnnz,*sub_onnz,br;
2035:       const PetscInt *bmindices;
2036:       if (!B) continue;
2037:       ISGetLocalSize(nest->isglobal.row[i],&bm);
2038:       ISGetIndices(nest->isglobal.row[i],&bmindices);
2039:       PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf);
2040:       PetscMalloc1(bm,&iremote);
2041:       PetscMalloc1(bm,&sub_dnnz);
2042:       PetscMalloc1(bm,&sub_onnz);
2043:       for (k = 0; k < bm; ++k) {
2044:         sub_dnnz[k] = 0;
2045:         sub_onnz[k] = 0;
2046:       }
2047:       PetscObjectTypeCompare((PetscObject)B,MATTRANSPOSEMAT,&flg);
2048:       if (flg) {
2049:         PetscTryMethod(B,"MatTransposeGetMat_C",(Mat,Mat*),(B,&D));
2050:         PetscTryMethod(B,"MatHermitianTransposeGetMat_C",(Mat,Mat*),(B,&D));
2051:         MatConvert(B,((PetscObject)D)->type_name,MAT_INITIAL_MATRIX,&D);
2052:         B = D;
2053:       }
2054:       PetscObjectTypeCompareAny((PetscObject)B,&flg,MATSEQSBAIJ,MATMPISBAIJ,"");
2055:       if (flg) {
2056:         if (D) MatConvert(D,MATBAIJ,MAT_INPLACE_MATRIX,&D);
2057:         else MatConvert(B,MATBAIJ,MAT_INITIAL_MATRIX,&D);
2058:         B = D;
2059:       }
2060:       /*
2061:        Locate the owners for all of the locally-owned global row indices for this row block.
2062:        These determine the roots of PetscSF used to communicate preallocation data to row owners.
2063:        The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
2064:        */
2065:       MatGetOwnershipRange(B,&rstart,NULL);
2066:       for (br = 0; br < bm; ++br) {
2067:         PetscInt       row = bmindices[br], brncols, col;
2068:         const PetscInt *brcols;
2069:         PetscInt       rowrel = 0; /* row's relative index on its owner rank */
2070:         PetscMPIInt    rowowner = 0;
2071:         PetscLayoutFindOwnerIndex(A->rmap,row,&rowowner,&rowrel);
2072:         /* how many roots  */
2073:         iremote[br].rank = rowowner; iremote[br].index = rowrel;           /* edge from bmdnnz to dnnz */
2074:         /* get nonzero pattern */
2075:         MatGetRow(B,br+rstart,&brncols,&brcols,NULL);
2076:         for (k=0; k<brncols; k++) {
2077:           col  = bNindices[brcols[k]];
2078:           if (col>=A->cmap->range[rowowner] && col<A->cmap->range[rowowner+1]) {
2079:             sub_dnnz[br]++;
2080:           } else {
2081:             sub_onnz[br]++;
2082:           }
2083:         }
2084:         MatRestoreRow(B,br+rstart,&brncols,&brcols,NULL);
2085:       }
2086:       if (D) MatDestroy(&D);
2087:       ISRestoreIndices(nest->isglobal.row[i],&bmindices);
2088:       /* bsf will have to take care of disposing of bedges. */
2089:       PetscSFSetGraph(bmsf,m,bm,NULL,PETSC_OWN_POINTER,iremote,PETSC_OWN_POINTER);
2090:       PetscSFReduceBegin(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
2091:       PetscSFReduceEnd(bmsf,MPIU_INT,sub_dnnz,dnnz,MPI_SUM);
2092:       PetscSFReduceBegin(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
2093:       PetscSFReduceEnd(bmsf,MPIU_INT,sub_onnz,onnz,MPI_SUM);
2094:       PetscFree(sub_dnnz);
2095:       PetscFree(sub_onnz);
2096:       PetscSFDestroy(&bmsf);
2097:     }
2098:     ISRestoreIndices(bNis,&bNindices);
2099:     ISDestroy(&bNis);
2100:   }
2101:   /* Resize preallocation if overestimated */
2102:   for (i=0;i<m;i++) {
2103:     dnnz[i] = PetscMin(dnnz[i],A->cmap->n);
2104:     onnz[i] = PetscMin(onnz[i],A->cmap->N - A->cmap->n);
2105:   }
2106:   MatSeqAIJSetPreallocation(C,0,dnnz);
2107:   MatMPIAIJSetPreallocation(C,0,dnnz,0,onnz);
2108:   PetscFree(dnnz);
2109:   MatAXPY_Dense_Nest(C,1.0,A);
2110:   if (reuse == MAT_INPLACE_MATRIX) {
2111:     MatHeaderReplace(A,&C);
2112:   } else *newmat = C;
2113:   return 0;
2114: }

2116: PetscErrorCode MatConvert_Nest_Dense(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
2117: {
2118:   Mat            B;
2119:   PetscInt       m,n,M,N;

2121:   MatGetSize(A,&M,&N);
2122:   MatGetLocalSize(A,&m,&n);
2123:   if (reuse == MAT_REUSE_MATRIX) {
2124:     B = *newmat;
2125:     MatZeroEntries(B);
2126:   } else {
2127:     MatCreateDense(PetscObjectComm((PetscObject)A),m,PETSC_DECIDE,M,N,NULL,&B);
2128:   }
2129:   MatAXPY_Dense_Nest(B,1.0,A);
2130:   if (reuse == MAT_INPLACE_MATRIX) {
2131:     MatHeaderReplace(A,&B);
2132:   } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
2133:   return 0;
2134: }

2136: PetscErrorCode MatHasOperation_Nest(Mat mat,MatOperation op,PetscBool *has)
2137: {
2138:   Mat_Nest       *bA = (Mat_Nest*)mat->data;
2139:   MatOperation   opAdd;
2140:   PetscInt       i,j,nr = bA->nr,nc = bA->nc;
2141:   PetscBool      flg;

2143:   *has = PETSC_FALSE;
2144:   if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
2145:     opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD);
2146:     for (j=0; j<nc; j++) {
2147:       for (i=0; i<nr; i++) {
2148:         if (!bA->m[i][j]) continue;
2149:         MatHasOperation(bA->m[i][j],opAdd,&flg);
2150:         if (!flg) return 0;
2151:       }
2152:     }
2153:   }
2154:   if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
2155:   return 0;
2156: }

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

2161:   Level: intermediate

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

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

2172: .seealso: `MatCreate()`, `MatType`, `MatCreateNest()`, `MatNestSetSubMat()`, `MatNestGetSubMat()`,
2173:           `VecCreateNest()`, `DMCreateMatrix()`, `DMCOMPOSITE`, `MatNestSetVecType()`, `MatNestGetLocalISs()`,
2174:           `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()`
2175: M*/
2176: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2177: {
2178:   Mat_Nest       *s;

2180:   PetscNewLog(A,&s);
2181:   A->data = (void*)s;

2183:   s->nr            = -1;
2184:   s->nc            = -1;
2185:   s->m             = NULL;
2186:   s->splitassembly = PETSC_FALSE;

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

2190:   A->ops->mult                  = MatMult_Nest;
2191:   A->ops->multadd               = MatMultAdd_Nest;
2192:   A->ops->multtranspose         = MatMultTranspose_Nest;
2193:   A->ops->multtransposeadd      = MatMultTransposeAdd_Nest;
2194:   A->ops->transpose             = MatTranspose_Nest;
2195:   A->ops->assemblybegin         = MatAssemblyBegin_Nest;
2196:   A->ops->assemblyend           = MatAssemblyEnd_Nest;
2197:   A->ops->zeroentries           = MatZeroEntries_Nest;
2198:   A->ops->copy                  = MatCopy_Nest;
2199:   A->ops->axpy                  = MatAXPY_Nest;
2200:   A->ops->duplicate             = MatDuplicate_Nest;
2201:   A->ops->createsubmatrix       = MatCreateSubMatrix_Nest;
2202:   A->ops->destroy               = MatDestroy_Nest;
2203:   A->ops->view                  = MatView_Nest;
2204:   A->ops->getvecs               = NULL; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2205:   A->ops->getlocalsubmatrix     = MatGetLocalSubMatrix_Nest;
2206:   A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2207:   A->ops->getdiagonal           = MatGetDiagonal_Nest;
2208:   A->ops->diagonalscale         = MatDiagonalScale_Nest;
2209:   A->ops->scale                 = MatScale_Nest;
2210:   A->ops->shift                 = MatShift_Nest;
2211:   A->ops->diagonalset           = MatDiagonalSet_Nest;
2212:   A->ops->setrandom             = MatSetRandom_Nest;
2213:   A->ops->hasoperation          = MatHasOperation_Nest;
2214:   A->ops->missingdiagonal       = MatMissingDiagonal_Nest;

2216:   A->spptr        = NULL;
2217:   A->assembled    = PETSC_FALSE;

2219:   /* expose Nest api's */
2220:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMat_C",        MatNestGetSubMat_Nest);
2221:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMat_C",        MatNestSetSubMat_Nest);
2222:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSubMats_C",       MatNestGetSubMats_Nest);
2223:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetSize_C",          MatNestGetSize_Nest);
2224:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetISs_C",           MatNestGetISs_Nest);
2225:   PetscObjectComposeFunction((PetscObject)A,"MatNestGetLocalISs_C",      MatNestGetLocalISs_Nest);
2226:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetVecType_C",       MatNestSetVecType_Nest);
2227:   PetscObjectComposeFunction((PetscObject)A,"MatNestSetSubMats_C",       MatNestSetSubMats_Nest);
2228:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpiaij_C",  MatConvert_Nest_AIJ);
2229:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqaij_C",  MatConvert_Nest_AIJ);
2230:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_aij_C",     MatConvert_Nest_AIJ);
2231:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_is_C",      MatConvert_Nest_IS);
2232:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_mpidense_C",MatConvert_Nest_Dense);
2233:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_nest_seqdense_C",MatConvert_Nest_Dense);
2234:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_seqdense_C",MatProductSetFromOptions_Nest_Dense);
2235:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_mpidense_C",MatProductSetFromOptions_Nest_Dense);
2236:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_nest_dense_C",MatProductSetFromOptions_Nest_Dense);

2238:   PetscObjectChangeTypeName((PetscObject)A,MATNEST);
2239:   return 0;
2240: }