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

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

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

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

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

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

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

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

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

140:       /* C[i] <- workC + C[i] */
141:       PetscCall(MatAXPY(viewC, 1.0, contents->workC[i * nc + j], SAME_NONZERO_PATTERN));
142:     }
143:     PetscCall(MatDestroy(&viewC));
144:   }
145:   PetscCall(MatDenseRestoreArray(C, &carray));
146:   PetscCall(MatDenseRestoreArrayRead(B, &barray));

148:   PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
149:   PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
150:   PetscFunctionReturn(PETSC_SUCCESS);
151: }

153: PetscErrorCode MatNest_DenseDestroy(void *ctx)
154: {
155:   Nest_Dense *contents = (Nest_Dense *)ctx;
156:   PetscInt    i;

158:   PetscFunctionBegin;
159:   PetscCall(PetscFree(contents->tarray));
160:   for (i = 0; i < contents->k; i++) PetscCall(MatDestroy(contents->workC + i));
161:   PetscCall(PetscFree3(contents->dm, contents->dn, contents->workC));
162:   PetscCall(PetscFree(contents));
163:   PetscFunctionReturn(PETSC_SUCCESS);
164: }

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

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

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

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

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

246:   C->ops->productnumeric = MatProductNumeric_Nest_Dense;
247:   PetscFunctionReturn(PETSC_SUCCESS);
248: }

250: static PetscErrorCode MatProductSetFromOptions_Nest_Dense_AB(Mat C)
251: {
252:   PetscFunctionBegin;
253:   C->ops->productsymbolic = MatProductSymbolic_Nest_Dense;
254:   PetscFunctionReturn(PETSC_SUCCESS);
255: }

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

261:   PetscFunctionBegin;
262:   if (product->type == MATPRODUCT_AB) PetscCall(MatProductSetFromOptions_Nest_Dense_AB(C));
263:   PetscFunctionReturn(PETSC_SUCCESS);
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:   PetscFunctionBegin;
273:   for (i = 0; i < nr; i++) PetscCall(VecGetSubVector(x, bA->isglobal.row[i], &bx[i]));
274:   for (i = 0; i < nc; i++) PetscCall(VecGetSubVector(y, bA->isglobal.col[i], &by[i]));
275:   for (j = 0; j < nc; j++) {
276:     PetscCall(VecZeroEntries(by[j]));
277:     for (i = 0; i < nr; i++) {
278:       if (!bA->m[i][j]) continue;
279:       /* y[j] <- y[j] + (A[i][j])^T * x[i] */
280:       PetscCall(MatMultTransposeAdd(bA->m[i][j], bx[i], by[j], by[j]));
281:     }
282:   }
283:   for (i = 0; i < nr; i++) PetscCall(VecRestoreSubVector(x, bA->isglobal.row[i], &bx[i]));
284:   for (i = 0; i < nc; i++) PetscCall(VecRestoreSubVector(y, bA->isglobal.col[i], &by[i]));
285:   PetscFunctionReturn(PETSC_SUCCESS);
286: }

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

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

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

321:   PetscFunctionBegin;
322:   if (reuse == MAT_REUSE_MATRIX) PetscCall(MatTransposeCheckNonzeroState_Private(A, *B));
323:   PetscCheck(reuse != MAT_INPLACE_MATRIX || nr == nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_SIZ, "Square nested matrix only for in-place");

325:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) {
326:     Mat *subs;
327:     IS  *is_row, *is_col;

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

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

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

356:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
357:     *B = C;
358:   } else {
359:     PetscCall(MatHeaderMerge(A, &C));
360:   }
361:   PetscFunctionReturn(PETSC_SUCCESS);
362: }

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

369:   PetscFunctionBegin;
370:   if (!lst) PetscFunctionReturn(PETSC_SUCCESS);
371:   for (i = 0; i < n; i++)
372:     if (lst[i]) PetscCall(ISDestroy(&lst[i]));
373:   PetscCall(PetscFree(lst));
374:   *list = NULL;
375:   PetscFunctionReturn(PETSC_SUCCESS);
376: }

378: static PetscErrorCode MatReset_Nest(Mat A)
379: {
380:   Mat_Nest *vs = (Mat_Nest *)A->data;
381:   PetscInt  i, j;

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

390:   PetscCall(PetscFree(vs->row_len));
391:   PetscCall(PetscFree(vs->col_len));
392:   PetscCall(PetscFree(vs->nnzstate));

394:   PetscCall(PetscFree2(vs->left, vs->right));

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

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

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

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

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

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

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

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

495:   PetscFunctionBegin;
496:   for (i = 0; i < vs->nr; i++) {
497:     for (j = 0; j < vs->nc; j++) {
498:       if (vs->m[i][j]) {
499:         if (vs->splitassembly) PetscCall(MatAssemblyEnd(vs->m[i][j], type));
500:       }
501:     }
502:   }
503:   PetscFunctionReturn(PETSC_SUCCESS);
504: }

506: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A, PetscInt row, Mat *B)
507: {
508:   Mat_Nest *vs = (Mat_Nest *)A->data;
509:   PetscInt  j;
510:   Mat       sub;

512:   PetscFunctionBegin;
513:   sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
514:   for (j = 0; !sub && j < vs->nc; j++) sub = vs->m[row][j];
515:   if (sub) PetscCall(MatSetUp(sub)); /* Ensure that the sizes are available */
516:   *B = sub;
517:   PetscFunctionReturn(PETSC_SUCCESS);
518: }

520: static PetscErrorCode MatNestFindNonzeroSubMatCol(Mat A, PetscInt col, Mat *B)
521: {
522:   Mat_Nest *vs = (Mat_Nest *)A->data;
523:   PetscInt  i;
524:   Mat       sub;

526:   PetscFunctionBegin;
527:   sub = (col < vs->nr) ? vs->m[col][col] : (Mat)NULL; /* Prefer to find on the diagonal */
528:   for (i = 0; !sub && i < vs->nr; i++) sub = vs->m[i][col];
529:   if (sub) PetscCall(MatSetUp(sub)); /* Ensure that the sizes are available */
530:   *B = sub;
531:   PetscFunctionReturn(PETSC_SUCCESS);
532: }

534: static PetscErrorCode MatNestFindISRange(Mat A, PetscInt n, const IS list[], IS is, PetscInt *begin, PetscInt *end)
535: {
536:   PetscInt  i, j, size, m;
537:   PetscBool flg;
538:   IS        out, concatenate[2];

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

589: static PetscErrorCode MatNestFillEmptyMat_Private(Mat A, PetscInt i, PetscInt j, Mat *B)
590: {
591:   Mat_Nest *vs = (Mat_Nest *)A->data;
592:   PetscInt  lr, lc;

594:   PetscFunctionBegin;
595:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
596:   PetscCall(ISGetLocalSize(vs->isglobal.row[i], &lr));
597:   PetscCall(ISGetLocalSize(vs->isglobal.col[j], &lc));
598:   PetscCall(MatSetSizes(*B, lr, lc, PETSC_DECIDE, PETSC_DECIDE));
599:   PetscCall(MatSetType(*B, MATAIJ));
600:   PetscCall(MatSeqAIJSetPreallocation(*B, 0, NULL));
601:   PetscCall(MatMPIAIJSetPreallocation(*B, 0, NULL, 0, NULL));
602:   PetscCall(MatSetUp(*B));
603:   PetscCall(MatSetOption(*B, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
604:   PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
605:   PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
606:   PetscFunctionReturn(PETSC_SUCCESS);
607: }

609: static PetscErrorCode MatNestGetBlock_Private(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *B)
610: {
611:   Mat_Nest  *vs = (Mat_Nest *)A->data;
612:   Mat       *a;
613:   PetscInt   i, j, k, l, nr = rend - rbegin, nc = cend - cbegin;
614:   char       keyname[256];
615:   PetscBool *b;
616:   PetscBool  flg;

618:   PetscFunctionBegin;
619:   *B = NULL;
620:   PetscCall(PetscSNPrintf(keyname, sizeof(keyname), "NestBlock_%" PetscInt_FMT "-%" PetscInt_FMT "x%" PetscInt_FMT "-%" PetscInt_FMT, rbegin, rend, cbegin, cend));
621:   PetscCall(PetscObjectQuery((PetscObject)A, keyname, (PetscObject *)B));
622:   if (*B) PetscFunctionReturn(PETSC_SUCCESS);

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

664: static PetscErrorCode MatNestFindSubMat(Mat A, struct MatNestISPair *is, IS isrow, IS iscol, Mat *B)
665: {
666:   Mat_Nest *vs = (Mat_Nest *)A->data;
667:   PetscInt  rbegin, rend, cbegin, cend;

669:   PetscFunctionBegin;
670:   PetscCall(MatNestFindISRange(A, vs->nr, is->row, isrow, &rbegin, &rend));
671:   PetscCall(MatNestFindISRange(A, vs->nc, is->col, iscol, &cbegin, &cend));
672:   if (rend == rbegin + 1 && cend == cbegin + 1) {
673:     if (!vs->m[rbegin][cbegin]) PetscCall(MatNestFillEmptyMat_Private(A, rbegin, cbegin, vs->m[rbegin] + cbegin));
674:     *B = vs->m[rbegin][cbegin];
675:   } else if (rbegin != -1 && cbegin != -1) {
676:     PetscCall(MatNestGetBlock_Private(A, rbegin, rend, cbegin, cend, B));
677:   } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Could not find index set");
678:   PetscFunctionReturn(PETSC_SUCCESS);
679: }

681: /*
682:    TODO: This does not actually returns a submatrix we can modify
683: */
684: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A, IS isrow, IS iscol, MatReuse reuse, Mat *B)
685: {
686:   Mat_Nest *vs = (Mat_Nest *)A->data;
687:   Mat       sub;

689:   PetscFunctionBegin;
690:   PetscCall(MatNestFindSubMat(A, &vs->isglobal, isrow, iscol, &sub));
691:   switch (reuse) {
692:   case MAT_INITIAL_MATRIX:
693:     if (sub) PetscCall(PetscObjectReference((PetscObject)sub));
694:     *B = sub;
695:     break;
696:   case MAT_REUSE_MATRIX:
697:     PetscCheck(sub == *B, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Submatrix was not used before in this call");
698:     break;
699:   case MAT_IGNORE_MATRIX: /* Nothing to do */
700:     break;
701:   case MAT_INPLACE_MATRIX: /* Nothing to do */
702:     SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MAT_INPLACE_MATRIX is not supported yet");
703:   }
704:   PetscFunctionReturn(PETSC_SUCCESS);
705: }

707: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A, IS isrow, IS iscol, Mat *B)
708: {
709:   Mat_Nest *vs = (Mat_Nest *)A->data;
710:   Mat       sub;

712:   PetscFunctionBegin;
713:   PetscCall(MatNestFindSubMat(A, &vs->islocal, isrow, iscol, &sub));
714:   /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
715:   if (sub) PetscCall(PetscObjectReference((PetscObject)sub));
716:   *B = sub;
717:   PetscFunctionReturn(PETSC_SUCCESS);
718: }

720: static PetscErrorCode MatRestoreLocalSubMatrix_Nest(Mat A, IS isrow, IS iscol, Mat *B)
721: {
722:   Mat_Nest *vs = (Mat_Nest *)A->data;
723:   Mat       sub;

725:   PetscFunctionBegin;
726:   PetscCall(MatNestFindSubMat(A, &vs->islocal, isrow, iscol, &sub));
727:   PetscCheck(*B == sub, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Local submatrix has not been gotten");
728:   if (sub) {
729:     PetscCheck(((PetscObject)sub)->refct > 1, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Local submatrix has had reference count decremented too many times");
730:     PetscCall(MatDestroy(B));
731:   }
732:   PetscFunctionReturn(PETSC_SUCCESS);
733: }

735: static PetscErrorCode MatGetDiagonal_Nest(Mat A, Vec v)
736: {
737:   Mat_Nest *bA = (Mat_Nest *)A->data;
738:   PetscInt  i;

740:   PetscFunctionBegin;
741:   for (i = 0; i < bA->nr; i++) {
742:     Vec bv;
743:     PetscCall(VecGetSubVector(v, bA->isglobal.row[i], &bv));
744:     if (bA->m[i][i]) {
745:       PetscCall(MatGetDiagonal(bA->m[i][i], bv));
746:     } else {
747:       PetscCall(VecSet(bv, 0.0));
748:     }
749:     PetscCall(VecRestoreSubVector(v, bA->isglobal.row[i], &bv));
750:   }
751:   PetscFunctionReturn(PETSC_SUCCESS);
752: }

754: static PetscErrorCode MatDiagonalScale_Nest(Mat A, Vec l, Vec r)
755: {
756:   Mat_Nest *bA = (Mat_Nest *)A->data;
757:   Vec       bl, *br;
758:   PetscInt  i, j;

760:   PetscFunctionBegin;
761:   PetscCall(PetscCalloc1(bA->nc, &br));
762:   if (r) {
763:     for (j = 0; j < bA->nc; j++) PetscCall(VecGetSubVector(r, bA->isglobal.col[j], &br[j]));
764:   }
765:   bl = NULL;
766:   for (i = 0; i < bA->nr; i++) {
767:     if (l) PetscCall(VecGetSubVector(l, bA->isglobal.row[i], &bl));
768:     for (j = 0; j < bA->nc; j++) {
769:       if (bA->m[i][j]) PetscCall(MatDiagonalScale(bA->m[i][j], bl, br[j]));
770:     }
771:     if (l) PetscCall(VecRestoreSubVector(l, bA->isglobal.row[i], &bl));
772:   }
773:   if (r) {
774:     for (j = 0; j < bA->nc; j++) PetscCall(VecRestoreSubVector(r, bA->isglobal.col[j], &br[j]));
775:   }
776:   PetscCall(PetscFree(br));
777:   PetscFunctionReturn(PETSC_SUCCESS);
778: }

780: static PetscErrorCode MatScale_Nest(Mat A, PetscScalar a)
781: {
782:   Mat_Nest *bA = (Mat_Nest *)A->data;
783:   PetscInt  i, j;

785:   PetscFunctionBegin;
786:   for (i = 0; i < bA->nr; i++) {
787:     for (j = 0; j < bA->nc; j++) {
788:       if (bA->m[i][j]) PetscCall(MatScale(bA->m[i][j], a));
789:     }
790:   }
791:   PetscFunctionReturn(PETSC_SUCCESS);
792: }

794: static PetscErrorCode MatShift_Nest(Mat A, PetscScalar a)
795: {
796:   Mat_Nest *bA = (Mat_Nest *)A->data;
797:   PetscInt  i;
798:   PetscBool nnzstate = PETSC_FALSE;

800:   PetscFunctionBegin;
801:   for (i = 0; i < bA->nr; i++) {
802:     PetscObjectState subnnzstate = 0;
803:     PetscCheck(bA->m[i][i], PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for shifting an empty diagonal block, insert a matrix in block (%" PetscInt_FMT ",%" PetscInt_FMT ")", i, i);
804:     PetscCall(MatShift(bA->m[i][i], a));
805:     PetscCall(MatGetNonzeroState(bA->m[i][i], &subnnzstate));
806:     nnzstate                     = (PetscBool)(nnzstate || bA->nnzstate[i * bA->nc + i] != subnnzstate);
807:     bA->nnzstate[i * bA->nc + i] = subnnzstate;
808:   }
809:   if (nnzstate) A->nonzerostate++;
810:   PetscFunctionReturn(PETSC_SUCCESS);
811: }

813: static PetscErrorCode MatDiagonalSet_Nest(Mat A, Vec D, InsertMode is)
814: {
815:   Mat_Nest *bA = (Mat_Nest *)A->data;
816:   PetscInt  i;
817:   PetscBool nnzstate = PETSC_FALSE;

819:   PetscFunctionBegin;
820:   for (i = 0; i < bA->nr; i++) {
821:     PetscObjectState subnnzstate = 0;
822:     Vec              bv;
823:     PetscCall(VecGetSubVector(D, bA->isglobal.row[i], &bv));
824:     if (bA->m[i][i]) {
825:       PetscCall(MatDiagonalSet(bA->m[i][i], bv, is));
826:       PetscCall(MatGetNonzeroState(bA->m[i][i], &subnnzstate));
827:     }
828:     PetscCall(VecRestoreSubVector(D, bA->isglobal.row[i], &bv));
829:     nnzstate                     = (PetscBool)(nnzstate || bA->nnzstate[i * bA->nc + i] != subnnzstate);
830:     bA->nnzstate[i * bA->nc + i] = subnnzstate;
831:   }
832:   if (nnzstate) A->nonzerostate++;
833:   PetscFunctionReturn(PETSC_SUCCESS);
834: }

836: static PetscErrorCode MatSetRandom_Nest(Mat A, PetscRandom rctx)
837: {
838:   Mat_Nest *bA = (Mat_Nest *)A->data;
839:   PetscInt  i, j;

841:   PetscFunctionBegin;
842:   for (i = 0; i < bA->nr; i++) {
843:     for (j = 0; j < bA->nc; j++) {
844:       if (bA->m[i][j]) PetscCall(MatSetRandom(bA->m[i][j], rctx));
845:     }
846:   }
847:   PetscFunctionReturn(PETSC_SUCCESS);
848: }

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

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

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

892:     PetscCall(VecCreateNest(comm, bA->nr, bA->isglobal.row, L, left));
893:     for (i = 0; i < bA->nr; i++) PetscCall(VecDestroy(&L[i]));

895:     PetscCall(PetscFree(L));
896:   }
897:   PetscFunctionReturn(PETSC_SUCCESS);
898: }

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

906:   PetscFunctionBegin;
907:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
908:   if (isascii) {
909:     PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_view_nest_sub", &viewSub, NULL));
910:     PetscCall(PetscViewerASCIIPrintf(viewer, "Matrix object: \n"));
911:     PetscCall(PetscViewerASCIIPushTab(viewer));
912:     PetscCall(PetscViewerASCIIPrintf(viewer, "type=nest, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT " \n", bA->nr, bA->nc));

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

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

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

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

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

951:   PetscFunctionBegin;
952:   for (i = 0; i < bA->nr; i++) {
953:     for (j = 0; j < bA->nc; j++) {
954:       if (!bA->m[i][j]) continue;
955:       PetscCall(MatZeroEntries(bA->m[i][j]));
956:     }
957:   }
958:   PetscFunctionReturn(PETSC_SUCCESS);
959: }

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

967:   PetscFunctionBegin;
968:   PetscCheck(nr == bB->nr && nc == bB->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Cannot copy a Mat_Nest of block size (%" PetscInt_FMT ",%" PetscInt_FMT ") to a Mat_Nest of block size (%" PetscInt_FMT ",%" PetscInt_FMT ")", bB->nr, bB->nc, nr, nc);
969:   for (i = 0; i < nr; i++) {
970:     for (j = 0; j < nc; j++) {
971:       PetscObjectState subnnzstate = 0;
972:       if (bA->m[i][j] && bB->m[i][j]) {
973:         PetscCall(MatCopy(bA->m[i][j], bB->m[i][j], str));
974:       } else PetscCheck(!bA->m[i][j] && !bB->m[i][j], PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Matrix block does not exist at %" PetscInt_FMT ",%" PetscInt_FMT, i, j);
975:       PetscCall(MatGetNonzeroState(bB->m[i][j], &subnnzstate));
976:       nnzstate                 = (PetscBool)(nnzstate || bB->nnzstate[i * nc + j] != subnnzstate);
977:       bB->nnzstate[i * nc + j] = subnnzstate;
978:     }
979:   }
980:   if (nnzstate) B->nonzerostate++;
981:   PetscFunctionReturn(PETSC_SUCCESS);
982: }

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

990:   PetscFunctionBegin;
991:   PetscCheck(nr == bX->nr && nc == bX->nc, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_INCOMP, "Cannot AXPY a MatNest of block size (%" PetscInt_FMT ",%" PetscInt_FMT ") with a MatNest of block size (%" PetscInt_FMT ",%" PetscInt_FMT ")", bX->nr, bX->nc, nr, nc);
992:   for (i = 0; i < nr; i++) {
993:     for (j = 0; j < nc; j++) {
994:       PetscObjectState subnnzstate = 0;
995:       if (bY->m[i][j] && bX->m[i][j]) {
996:         PetscCall(MatAXPY(bY->m[i][j], a, bX->m[i][j], str));
997:       } else if (bX->m[i][j]) {
998:         Mat M;

1000:         PetscCheck(str == DIFFERENT_NONZERO_PATTERN, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_INCOMP, "Matrix block does not exist at %" PetscInt_FMT ",%" PetscInt_FMT ". Use DIFFERENT_NONZERO_PATTERN", i, j);
1001:         PetscCall(MatDuplicate(bX->m[i][j], MAT_COPY_VALUES, &M));
1002:         PetscCall(MatNestSetSubMat(Y, i, j, M));
1003:         PetscCall(MatDestroy(&M));
1004:       }
1005:       if (bY->m[i][j]) PetscCall(MatGetNonzeroState(bY->m[i][j], &subnnzstate));
1006:       nnzstate                 = (PetscBool)(nnzstate || bY->nnzstate[i * nc + j] != subnnzstate);
1007:       bY->nnzstate[i * nc + j] = subnnzstate;
1008:     }
1009:   }
1010:   if (nnzstate) Y->nonzerostate++;
1011:   PetscFunctionReturn(PETSC_SUCCESS);
1012: }

1014: static PetscErrorCode MatDuplicate_Nest(Mat A, MatDuplicateOption op, Mat *B)
1015: {
1016:   Mat_Nest *bA = (Mat_Nest *)A->data;
1017:   Mat      *b;
1018:   PetscInt  i, j, nr = bA->nr, nc = bA->nc;

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

1036:   PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
1037:   PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
1038:   PetscFunctionReturn(PETSC_SUCCESS);
1039: }

1041: /* nest api */
1042: PetscErrorCode MatNestGetSubMat_Nest(Mat A, PetscInt idxm, PetscInt jdxm, Mat *mat)
1043: {
1044:   Mat_Nest *bA = (Mat_Nest *)A->data;

1046:   PetscFunctionBegin;
1047:   PetscCheck(idxm < bA->nr, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, idxm, bA->nr - 1);
1048:   PetscCheck(jdxm < bA->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Col too large: row %" PetscInt_FMT " max %" PetscInt_FMT, jdxm, bA->nc - 1);
1049:   *mat = bA->m[idxm][jdxm];
1050:   PetscFunctionReturn(PETSC_SUCCESS);
1051: }

1053: /*@
1054:  MatNestGetSubMat - Returns a single, sub-matrix from a `MATNEST`

1056:  Not Collective

1058:  Input Parameters:
1059: +   A  - `MATNEST` matrix
1060: .   idxm - index of the matrix within the nest matrix
1061: -   jdxm - index of the matrix within the nest matrix

1063:  Output Parameter:
1064: .   sub - matrix at index `idxm`, `jdxm` within the nest matrix

1066:  Level: developer

1068: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSize()`, `MatNestGetSubMats()`, `MatCreateNest()`, `MATNEST`, `MatNestSetSubMat()`,
1069:           `MatNestGetLocalISs()`, `MatNestGetISs()`
1070: @*/
1071: PetscErrorCode MatNestGetSubMat(Mat A, PetscInt idxm, PetscInt jdxm, Mat *sub)
1072: {
1073:   PetscFunctionBegin;
1074:   PetscUseMethod(A, "MatNestGetSubMat_C", (Mat, PetscInt, PetscInt, Mat *), (A, idxm, jdxm, sub));
1075:   PetscFunctionReturn(PETSC_SUCCESS);
1076: }

1078: PetscErrorCode MatNestSetSubMat_Nest(Mat A, PetscInt idxm, PetscInt jdxm, Mat mat)
1079: {
1080:   Mat_Nest *bA = (Mat_Nest *)A->data;
1081:   PetscInt  m, n, M, N, mi, ni, Mi, Ni;

1083:   PetscFunctionBegin;
1084:   PetscCheck(idxm < bA->nr, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, idxm, bA->nr - 1);
1085:   PetscCheck(jdxm < bA->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Col too large: row %" PetscInt_FMT " max %" PetscInt_FMT, jdxm, bA->nc - 1);
1086:   PetscCall(MatGetLocalSize(mat, &m, &n));
1087:   PetscCall(MatGetSize(mat, &M, &N));
1088:   PetscCall(ISGetLocalSize(bA->isglobal.row[idxm], &mi));
1089:   PetscCall(ISGetSize(bA->isglobal.row[idxm], &Mi));
1090:   PetscCall(ISGetLocalSize(bA->isglobal.col[jdxm], &ni));
1091:   PetscCall(ISGetSize(bA->isglobal.col[jdxm], &Ni));
1092:   PetscCheck(M == Mi && N == Ni, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_INCOMP, "Submatrix dimension (%" PetscInt_FMT ",%" PetscInt_FMT ") incompatible with nest block (%" PetscInt_FMT ",%" PetscInt_FMT ")", M, N, Mi, Ni);
1093:   PetscCheck(m == mi && n == ni, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_INCOMP, "Submatrix local dimension (%" PetscInt_FMT ",%" PetscInt_FMT ") incompatible with nest block (%" PetscInt_FMT ",%" PetscInt_FMT ")", m, n, mi, ni);

1095:   /* do not increase object state */
1096:   if (mat == bA->m[idxm][jdxm]) PetscFunctionReturn(PETSC_SUCCESS);

1098:   PetscCall(PetscObjectReference((PetscObject)mat));
1099:   PetscCall(MatDestroy(&bA->m[idxm][jdxm]));
1100:   bA->m[idxm][jdxm] = mat;
1101:   PetscCall(PetscObjectStateIncrease((PetscObject)A));
1102:   PetscCall(MatGetNonzeroState(mat, &bA->nnzstate[idxm * bA->nc + jdxm]));
1103:   A->nonzerostate++;
1104:   PetscFunctionReturn(PETSC_SUCCESS);
1105: }

1107: /*@
1108:  MatNestSetSubMat - Set a single submatrix in the `MATNEST`

1110:  Logically Collective

1112:  Input Parameters:
1113: +   A  - `MATNEST` matrix
1114: .   idxm - index of the matrix within the nest matrix
1115: .   jdxm - index of the matrix within the nest matrix
1116: -   sub - matrix at index `idxm`, `jdxm` within the nest matrix

1118:  Level: developer

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

1123:  This increments the reference count of the submatrix.

1125: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestSetSubMats()`, `MatNestGetSubMats()`, `MatNestGetLocalISs()`, `MATNEST`, `MatCreateNest()`,
1126:           `MatNestGetSubMat()`, `MatNestGetISs()`, `MatNestGetSize()`
1127: @*/
1128: PetscErrorCode MatNestSetSubMat(Mat A, PetscInt idxm, PetscInt jdxm, Mat sub)
1129: {
1130:   PetscFunctionBegin;
1131:   PetscUseMethod(A, "MatNestSetSubMat_C", (Mat, PetscInt, PetscInt, Mat), (A, idxm, jdxm, sub));
1132:   PetscFunctionReturn(PETSC_SUCCESS);
1133: }

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

1139:   PetscFunctionBegin;
1140:   if (M) *M = bA->nr;
1141:   if (N) *N = bA->nc;
1142:   if (mat) *mat = bA->m;
1143:   PetscFunctionReturn(PETSC_SUCCESS);
1144: }

1146: /*@C
1147:  MatNestGetSubMats - Returns the entire two dimensional array of matrices defining a `MATNEST` matrix.

1149:  Not Collective

1151:  Input Parameter:
1152: .   A  - nest matrix

1154:  Output Parameters:
1155: +   M - number of rows in the nest matrix
1156: .   N - number of cols in the nest matrix
1157: -   mat - 2d array of matrices

1159:  Level: developer

1161:  Note:
1162:  The user should not free the array `mat`.

1164:  Fortran Note:
1165:  This routine has a calling sequence
1166: $   call MatNestGetSubMats(A, M, N, mat, ierr)
1167:  where the space allocated for the optional argument `mat` is assumed large enough (if provided).

1169: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSize()`, `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MATNEST`, `MatCreateNest()`,
1170:           `MatNestSetSubMats()`, `MatNestGetISs()`, `MatNestSetSubMat()`
1171: @*/
1172: PetscErrorCode MatNestGetSubMats(Mat A, PetscInt *M, PetscInt *N, Mat ***mat)
1173: {
1174:   PetscFunctionBegin;
1175:   PetscUseMethod(A, "MatNestGetSubMats_C", (Mat, PetscInt *, PetscInt *, Mat ***), (A, M, N, mat));
1176:   PetscFunctionReturn(PETSC_SUCCESS);
1177: }

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

1183:   PetscFunctionBegin;
1184:   if (M) *M = bA->nr;
1185:   if (N) *N = bA->nc;
1186:   PetscFunctionReturn(PETSC_SUCCESS);
1187: }

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

1192:  Not Collective

1194:  Input Parameter:
1195: .   A  - `MATNEST` matrix

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

1201:  Level: developer

1203: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MATNEST`, `MatCreateNest()`, `MatNestGetLocalISs()`,
1204:           `MatNestGetISs()`
1205: @*/
1206: PetscErrorCode MatNestGetSize(Mat A, PetscInt *M, PetscInt *N)
1207: {
1208:   PetscFunctionBegin;
1209:   PetscUseMethod(A, "MatNestGetSize_C", (Mat, PetscInt *, PetscInt *), (A, M, N));
1210:   PetscFunctionReturn(PETSC_SUCCESS);
1211: }

1213: static PetscErrorCode MatNestGetISs_Nest(Mat A, IS rows[], IS cols[])
1214: {
1215:   Mat_Nest *vs = (Mat_Nest *)A->data;
1216:   PetscInt  i;

1218:   PetscFunctionBegin;
1219:   if (rows)
1220:     for (i = 0; i < vs->nr; i++) rows[i] = vs->isglobal.row[i];
1221:   if (cols)
1222:     for (i = 0; i < vs->nc; i++) cols[i] = vs->isglobal.col[i];
1223:   PetscFunctionReturn(PETSC_SUCCESS);
1224: }

1226: /*@C
1227:  MatNestGetISs - Returns the index sets partitioning the row and column spaces of a `MATNEST`

1229:  Not Collective

1231:  Input Parameter:
1232: .   A  - `MATNEST` matrix

1234:  Output Parameters:
1235: +   rows - array of row index sets
1236: -   cols - array of column index sets

1238:  Level: advanced

1240:  Note:
1241:  The user must have allocated arrays of the correct size. The reference count is not increased on the returned `IS`s.

1243: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetLocalISs()`, `MATNEST`,
1244:           `MatCreateNest()`, `MatNestGetSubMats()`, `MatNestSetSubMats()`
1245: @*/
1246: PetscErrorCode MatNestGetISs(Mat A, IS rows[], IS cols[])
1247: {
1248:   PetscFunctionBegin;
1250:   PetscUseMethod(A, "MatNestGetISs_C", (Mat, IS[], IS[]), (A, rows, cols));
1251:   PetscFunctionReturn(PETSC_SUCCESS);
1252: }

1254: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A, IS rows[], IS cols[])
1255: {
1256:   Mat_Nest *vs = (Mat_Nest *)A->data;
1257:   PetscInt  i;

1259:   PetscFunctionBegin;
1260:   if (rows)
1261:     for (i = 0; i < vs->nr; i++) rows[i] = vs->islocal.row[i];
1262:   if (cols)
1263:     for (i = 0; i < vs->nc; i++) cols[i] = vs->islocal.col[i];
1264:   PetscFunctionReturn(PETSC_SUCCESS);
1265: }

1267: /*@C
1268:  MatNestGetLocalISs - Returns the index sets partitioning the row and column spaces of a `MATNEST`

1270:  Not Collective

1272:  Input Parameter:
1273: .   A  - `MATNEST` matrix

1275:  Output Parameters:
1276: +   rows - array of row index sets (or `NULL` to ignore)
1277: -   cols - array of column index sets (or `NULL` to ignore)

1279:  Level: advanced

1281:  Note:
1282:  The user must have allocated arrays of the correct size. The reference count is not increased on the returned `IS`s.

1284: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetISs()`, `MatCreateNest()`,
1285:           `MATNEST`, `MatNestSetSubMats()`, `MatNestSetSubMat()`
1286: @*/
1287: PetscErrorCode MatNestGetLocalISs(Mat A, IS rows[], IS cols[])
1288: {
1289:   PetscFunctionBegin;
1291:   PetscUseMethod(A, "MatNestGetLocalISs_C", (Mat, IS[], IS[]), (A, rows, cols));
1292:   PetscFunctionReturn(PETSC_SUCCESS);
1293: }

1295: PetscErrorCode MatNestSetVecType_Nest(Mat A, VecType vtype)
1296: {
1297:   PetscBool flg;

1299:   PetscFunctionBegin;
1300:   PetscCall(PetscStrcmp(vtype, VECNEST, &flg));
1301:   /* In reality, this only distinguishes VECNEST and "other" */
1302:   if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1303:   else A->ops->getvecs = (PetscErrorCode(*)(Mat, Vec *, Vec *))0;
1304:   PetscFunctionReturn(PETSC_SUCCESS);
1305: }

1307: /*@C
1308:  MatNestSetVecType - Sets the type of `Vec` returned by `MatCreateVecs()`

1310:  Not Collective

1312:  Input Parameters:
1313: +  A  - `MATNEST` matrix
1314: -  vtype - `VecType` to use for creating vectors

1316:  Level: developer

1318: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreateVecs()`, `MATNEST`, `MatCreateNest()`, `VecType`
1319: @*/
1320: PetscErrorCode MatNestSetVecType(Mat A, VecType vtype)
1321: {
1322:   PetscFunctionBegin;
1323:   PetscTryMethod(A, "MatNestSetVecType_C", (Mat, VecType), (A, vtype));
1324:   PetscFunctionReturn(PETSC_SUCCESS);
1325: }

1327: PetscErrorCode MatNestSetSubMats_Nest(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[])
1328: {
1329:   Mat_Nest *s = (Mat_Nest *)A->data;
1330:   PetscInt  i, j, m, n, M, N;
1331:   PetscBool cong, isstd, sametype = PETSC_FALSE;
1332:   VecType   vtype, type;

1334:   PetscFunctionBegin;
1335:   PetscCall(MatReset_Nest(A));

1337:   s->nr = nr;
1338:   s->nc = nc;

1340:   /* Create space for submatrices */
1341:   PetscCall(PetscMalloc1(nr, &s->m));
1342:   for (i = 0; i < nr; i++) PetscCall(PetscMalloc1(nc, &s->m[i]));
1343:   for (i = 0; i < nr; i++) {
1344:     for (j = 0; j < nc; j++) {
1345:       s->m[i][j] = a[i * nc + j];
1346:       if (a[i * nc + j]) PetscCall(PetscObjectReference((PetscObject)a[i * nc + j]));
1347:     }
1348:   }
1349:   PetscCall(MatGetVecType(A, &vtype));
1350:   PetscCall(PetscStrcmp(vtype, VECSTANDARD, &isstd));
1351:   if (isstd) {
1352:     /* check if all blocks have the same vectype */
1353:     vtype = NULL;
1354:     for (i = 0; i < nr; i++) {
1355:       for (j = 0; j < nc; j++) {
1356:         if (a[i * nc + j]) {
1357:           if (!vtype) { /* first visited block */
1358:             PetscCall(MatGetVecType(a[i * nc + j], &vtype));
1359:             sametype = PETSC_TRUE;
1360:           } else if (sametype) {
1361:             PetscCall(MatGetVecType(a[i * nc + j], &type));
1362:             PetscCall(PetscStrcmp(vtype, type, &sametype));
1363:           }
1364:         }
1365:       }
1366:     }
1367:     if (sametype) { /* propagate vectype */
1368:       PetscCall(MatSetVecType(A, vtype));
1369:     }
1370:   }

1372:   PetscCall(MatSetUp_NestIS_Private(A, nr, is_row, nc, is_col));

1374:   PetscCall(PetscMalloc1(nr, &s->row_len));
1375:   PetscCall(PetscMalloc1(nc, &s->col_len));
1376:   for (i = 0; i < nr; i++) s->row_len[i] = -1;
1377:   for (j = 0; j < nc; j++) s->col_len[j] = -1;

1379:   PetscCall(PetscCalloc1(nr * nc, &s->nnzstate));
1380:   for (i = 0; i < nr; i++) {
1381:     for (j = 0; j < nc; j++) {
1382:       if (s->m[i][j]) PetscCall(MatGetNonzeroState(s->m[i][j], &s->nnzstate[i * nc + j]));
1383:     }
1384:   }

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

1388:   PetscCall(PetscLayoutSetSize(A->rmap, M));
1389:   PetscCall(PetscLayoutSetLocalSize(A->rmap, m));
1390:   PetscCall(PetscLayoutSetSize(A->cmap, N));
1391:   PetscCall(PetscLayoutSetLocalSize(A->cmap, n));

1393:   PetscCall(PetscLayoutSetUp(A->rmap));
1394:   PetscCall(PetscLayoutSetUp(A->cmap));

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

1410:   PetscCall(PetscCalloc2(nr, &s->left, nc, &s->right));
1411:   PetscCall(PetscObjectStateIncrease((PetscObject)A));
1412:   A->nonzerostate++;
1413:   PetscFunctionReturn(PETSC_SUCCESS);
1414: }

1416: /*@
1417:    MatNestSetSubMats - Sets the nested submatrices in a `MATNEST`

1419:    Collective

1421:    Input Parameters:
1422: +  A - `MATNEST` matrix
1423: .  nr - number of nested row blocks
1424: .  is_row - index sets for each nested row block, or `NULL` to make contiguous
1425: .  nc - number of nested column blocks
1426: .  is_col - index sets for each nested column block, or `NULL` to make contiguous
1427: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using `NULL`

1429:    Level: advanced

1431:    Note:
1432:    This always resets any submatrix information previously set

1434: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreateNest()`, `MatNestSetSubMat()`, `MatNestGetSubMat()`, `MatNestGetSubMats()`
1435: @*/
1436: PetscErrorCode MatNestSetSubMats(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[])
1437: {
1438:   PetscInt i;

1440:   PetscFunctionBegin;
1442:   PetscCheck(nr >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Number of rows cannot be negative");
1443:   if (nr && is_row) {
1446:   }
1447:   PetscCheck(nc >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Number of columns cannot be negative");
1448:   if (nc && is_col) {
1451:   }
1453:   PetscUseMethod(A, "MatNestSetSubMats_C", (Mat, PetscInt, const IS[], PetscInt, const IS[], const Mat[]), (A, nr, is_row, nc, is_col, a));
1454:   PetscFunctionReturn(PETSC_SUCCESS);
1455: }

1457: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A, PetscInt n, const IS islocal[], const IS isglobal[], PetscBool colflg, ISLocalToGlobalMapping *ltog)
1458: {
1459:   PetscBool flg;
1460:   PetscInt  i, j, m, mi, *ix;

1462:   PetscFunctionBegin;
1463:   *ltog = NULL;
1464:   for (i = 0, m = 0, flg = PETSC_FALSE; i < n; i++) {
1465:     if (islocal[i]) {
1466:       PetscCall(ISGetLocalSize(islocal[i], &mi));
1467:       flg = PETSC_TRUE; /* We found a non-trivial entry */
1468:     } else {
1469:       PetscCall(ISGetLocalSize(isglobal[i], &mi));
1470:     }
1471:     m += mi;
1472:   }
1473:   if (!flg) PetscFunctionReturn(PETSC_SUCCESS);

1475:   PetscCall(PetscMalloc1(m, &ix));
1476:   for (i = 0, m = 0; i < n; i++) {
1477:     ISLocalToGlobalMapping smap = NULL;
1478:     Mat                    sub  = NULL;
1479:     PetscSF                sf;
1480:     PetscLayout            map;
1481:     const PetscInt        *ix2;

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

1504:       PetscCall(ISGetLocalSize(islocal[i], &mi));
1505:       PetscCheck(smap, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing local to global map");
1506:       for (j = 0; j < mi; j++) ix[m + j] = j;
1507:       PetscCall(ISLocalToGlobalMappingApply(smap, mi, ix + m, ix + m));

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

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

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

1551:             proc 0:
1552:     is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1553:             proc 1:
1554:     is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)

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

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

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

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

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

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

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

1675:   if (PetscDefined(USE_DEBUG)) {
1676:     for (i = 0; i < vs->nr; i++) {
1677:       for (j = 0; j < vs->nc; j++) {
1678:         PetscInt m, n, M, N, mi, ni, Mi, Ni;
1679:         Mat      B = vs->m[i][j];
1680:         if (!B) continue;
1681:         PetscCall(MatGetSize(B, &M, &N));
1682:         PetscCall(MatGetLocalSize(B, &m, &n));
1683:         PetscCall(ISGetSize(vs->isglobal.row[i], &Mi));
1684:         PetscCall(ISGetSize(vs->isglobal.col[j], &Ni));
1685:         PetscCall(ISGetLocalSize(vs->isglobal.row[i], &mi));
1686:         PetscCall(ISGetLocalSize(vs->isglobal.col[j], &ni));
1687:         PetscCheck(M == Mi && N == Ni, PetscObjectComm((PetscObject)sub), PETSC_ERR_ARG_INCOMP, "Global sizes (%" PetscInt_FMT ",%" PetscInt_FMT ") of nested submatrix (%" PetscInt_FMT ",%" PetscInt_FMT ") do not agree with space defined by index sets (%" PetscInt_FMT ",%" PetscInt_FMT ")", M, N, i, j, Mi, Ni);
1688:         PetscCheck(m == mi && n == ni, PetscObjectComm((PetscObject)sub), PETSC_ERR_ARG_INCOMP, "Local sizes (%" PetscInt_FMT ",%" PetscInt_FMT ") of nested submatrix (%" PetscInt_FMT ",%" PetscInt_FMT ") do not agree with space defined by index sets (%" PetscInt_FMT ",%" PetscInt_FMT ")", m, n, i, j, mi, ni);
1689:       }
1690:     }
1691:   }

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

1703: /*@C
1704:    MatCreateNest - Creates a new `MATNEST` matrix containing several nested submatrices, each stored separately

1706:    Collective

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

1716:    Output Parameter:
1717: .  B - new matrix

1719:    Level: advanced

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

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

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

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

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

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

1775:           PetscCall(MatTransposeGetMat(B, &Bt));
1776:           PetscCall(MatTranspose(Bt, MAT_INITIAL_MATRIX, &trans[i * nest->nc + j]));
1777:           B = trans[i * nest->nc + j];
1778:         } else {
1779:           PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &istrans));
1780:           if (istrans) {
1781:             Mat Bt;

1783:             PetscCall(MatHermitianTransposeGetMat(B, &Bt));
1784:             PetscCall(MatHermitianTranspose(Bt, MAT_INITIAL_MATRIX, &trans[i * nest->nc + j]));
1785:             B = trans[i * nest->nc + j];
1786:           }
1787:         }
1788:         PetscCall(MatGetRowIJ(B, 0, PETSC_FALSE, PETSC_FALSE, &nnr, (const PetscInt **)&nii, (const PetscInt **)&njj, &done));
1789:         PetscCheck(done, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "MatGetRowIJ");
1790:         PetscCall(MatSeqAIJGetArray(B, &naa));
1791:         nnz += nii[nnr];

1793:         aii[i * nest->nc + j] = nii;
1794:         ajj[i * nest->nc + j] = njj;
1795:         avv[i * nest->nc + j] = naa;
1796:       }
1797:     }
1798:   }
1799:   if (reuse != MAT_REUSE_MATRIX) {
1800:     PetscCall(PetscMalloc1(nr + 1, &ii));
1801:     PetscCall(PetscMalloc1(nnz, &jj));
1802:     PetscCall(PetscMalloc1(nnz, &vv));
1803:   } else {
1804:     PetscCheck(nnz == ii[nr], PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Cannot reuse matrix, wrong number of nonzeros");
1805:   }

1807:   /* new row pointer */
1808:   PetscCall(PetscArrayzero(ii, nr + 1));
1809:   for (i = 0; i < nest->nr; ++i) {
1810:     PetscInt ncr, rst;

1812:     PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &rst, NULL));
1813:     PetscCall(ISGetLocalSize(nest->isglobal.row[i], &ncr));
1814:     for (j = 0; j < nest->nc; ++j) {
1815:       if (aii[i * nest->nc + j]) {
1816:         PetscInt *nii = aii[i * nest->nc + j];
1817:         PetscInt  ir;

1819:         for (ir = rst; ir < ncr + rst; ++ir) {
1820:           ii[ir + 1] += nii[1] - nii[0];
1821:           nii++;
1822:         }
1823:       }
1824:     }
1825:   }
1826:   for (i = 0; i < nr; i++) ii[i + 1] += ii[i];

1828:   /* construct CSR for the new matrix */
1829:   PetscCall(PetscCalloc1(nr, &ci));
1830:   for (i = 0; i < nest->nr; ++i) {
1831:     PetscInt ncr, rst;

1833:     PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &rst, NULL));
1834:     PetscCall(ISGetLocalSize(nest->isglobal.row[i], &ncr));
1835:     for (j = 0; j < nest->nc; ++j) {
1836:       if (aii[i * nest->nc + j]) {
1837:         PetscScalar *nvv = avv[i * nest->nc + j];
1838:         PetscInt    *nii = aii[i * nest->nc + j];
1839:         PetscInt    *njj = ajj[i * nest->nc + j];
1840:         PetscInt     ir, cst;

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

1846:           for (ij = 0; ij < rsize; ij++) {
1847:             jj[ist + ij] = *njj + cst;
1848:             vv[ist + ij] = *nvv;
1849:             njj++;
1850:             nvv++;
1851:           }
1852:           ci[ir] += rsize;
1853:           nii++;
1854:         }
1855:       }
1856:     }
1857:   }
1858:   PetscCall(PetscFree(ci));

1860:   /* restore info */
1861:   for (i = 0; i < nest->nr; ++i) {
1862:     for (j = 0; j < nest->nc; ++j) {
1863:       Mat B = nest->m[i][j];
1864:       if (B) {
1865:         PetscInt nnr = 0, k = i * nest->nc + j;

1867:         B = (trans[k] ? trans[k] : B);
1868:         PetscCall(MatRestoreRowIJ(B, 0, PETSC_FALSE, PETSC_FALSE, &nnr, (const PetscInt **)&aii[k], (const PetscInt **)&ajj[k], &done));
1869:         PetscCheck(done, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "MatRestoreRowIJ");
1870:         PetscCall(MatSeqAIJRestoreArray(B, &avv[k]));
1871:         PetscCall(MatDestroy(&trans[k]));
1872:       }
1873:     }
1874:   }
1875:   PetscCall(PetscFree4(aii, ajj, avv, trans));

1877:   /* finalize newmat */
1878:   if (reuse == MAT_INITIAL_MATRIX) {
1879:     PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), nr, nc, ii, jj, vv, newmat));
1880:   } else if (reuse == MAT_INPLACE_MATRIX) {
1881:     Mat B;

1883:     PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), nr, nc, ii, jj, vv, &B));
1884:     PetscCall(MatHeaderReplace(A, &B));
1885:   }
1886:   PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
1887:   PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
1888:   {
1889:     Mat_SeqAIJ *a = (Mat_SeqAIJ *)((*newmat)->data);
1890:     a->free_a     = PETSC_TRUE;
1891:     a->free_ij    = PETSC_TRUE;
1892:   }
1893:   PetscFunctionReturn(PETSC_SUCCESS);
1894: }

1896: PETSC_INTERN PetscErrorCode MatAXPY_Dense_Nest(Mat Y, PetscScalar a, Mat X)
1897: {
1898:   Mat_Nest *nest = (Mat_Nest *)X->data;
1899:   PetscInt  i, j, k, rstart;
1900:   PetscBool flg;

1902:   PetscFunctionBegin;
1903:   /* Fill by row */
1904:   for (j = 0; j < nest->nc; ++j) {
1905:     /* Using global column indices and ISAllGather() is not scalable. */
1906:     IS              bNis;
1907:     PetscInt        bN;
1908:     const PetscInt *bNindices;
1909:     PetscCall(ISAllGather(nest->isglobal.col[j], &bNis));
1910:     PetscCall(ISGetSize(bNis, &bN));
1911:     PetscCall(ISGetIndices(bNis, &bNindices));
1912:     for (i = 0; i < nest->nr; ++i) {
1913:       Mat             B = nest->m[i][j], D = NULL;
1914:       PetscInt        bm, br;
1915:       const PetscInt *bmindices;
1916:       if (!B) continue;
1917:       PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, ""));
1918:       if (flg) {
1919:         PetscTryMethod(B, "MatTransposeGetMat_C", (Mat, Mat *), (B, &D));
1920:         PetscTryMethod(B, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (B, &D));
1921:         PetscCall(MatConvert(B, ((PetscObject)D)->type_name, MAT_INITIAL_MATRIX, &D));
1922:         B = D;
1923:       }
1924:       PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));
1925:       if (flg) {
1926:         if (D) PetscCall(MatConvert(D, MATBAIJ, MAT_INPLACE_MATRIX, &D));
1927:         else PetscCall(MatConvert(B, MATBAIJ, MAT_INITIAL_MATRIX, &D));
1928:         B = D;
1929:       }
1930:       PetscCall(ISGetLocalSize(nest->isglobal.row[i], &bm));
1931:       PetscCall(ISGetIndices(nest->isglobal.row[i], &bmindices));
1932:       PetscCall(MatGetOwnershipRange(B, &rstart, NULL));
1933:       for (br = 0; br < bm; ++br) {
1934:         PetscInt           row = bmindices[br], brncols, *cols;
1935:         const PetscInt    *brcols;
1936:         const PetscScalar *brcoldata;
1937:         PetscScalar       *vals = NULL;
1938:         PetscCall(MatGetRow(B, br + rstart, &brncols, &brcols, &brcoldata));
1939:         PetscCall(PetscMalloc1(brncols, &cols));
1940:         for (k = 0; k < brncols; k++) cols[k] = bNindices[brcols[k]];
1941:         /*
1942:           Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
1943:           Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
1944:          */
1945:         if (a != 1.0) {
1946:           PetscCall(PetscMalloc1(brncols, &vals));
1947:           for (k = 0; k < brncols; k++) vals[k] = a * brcoldata[k];
1948:           PetscCall(MatSetValues(Y, 1, &row, brncols, cols, vals, ADD_VALUES));
1949:           PetscCall(PetscFree(vals));
1950:         } else {
1951:           PetscCall(MatSetValues(Y, 1, &row, brncols, cols, brcoldata, ADD_VALUES));
1952:         }
1953:         PetscCall(MatRestoreRow(B, br + rstart, &brncols, &brcols, &brcoldata));
1954:         PetscCall(PetscFree(cols));
1955:       }
1956:       if (D) PetscCall(MatDestroy(&D));
1957:       PetscCall(ISRestoreIndices(nest->isglobal.row[i], &bmindices));
1958:     }
1959:     PetscCall(ISRestoreIndices(bNis, &bNindices));
1960:     PetscCall(ISDestroy(&bNis));
1961:   }
1962:   PetscCall(MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY));
1963:   PetscCall(MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY));
1964:   PetscFunctionReturn(PETSC_SUCCESS);
1965: }

1967: PetscErrorCode MatConvert_Nest_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1968: {
1969:   Mat_Nest   *nest = (Mat_Nest *)A->data;
1970:   PetscInt    m, n, M, N, i, j, k, *dnnz, *onnz, rstart, cstart, cend;
1971:   PetscMPIInt size;
1972:   Mat         C;

1974:   PetscFunctionBegin;
1975:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1976:   if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */
1977:     PetscInt  nf;
1978:     PetscBool fast;

1980:     PetscCall(PetscStrcmp(newtype, MATAIJ, &fast));
1981:     if (!fast) PetscCall(PetscStrcmp(newtype, MATSEQAIJ, &fast));
1982:     for (i = 0; i < nest->nr && fast; ++i) {
1983:       for (j = 0; j < nest->nc && fast; ++j) {
1984:         Mat B = nest->m[i][j];
1985:         if (B) {
1986:           PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQAIJ, &fast));
1987:           if (!fast) {
1988:             PetscBool istrans;

1990:             PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &istrans));
1991:             if (istrans) {
1992:               Mat Bt;

1994:               PetscCall(MatTransposeGetMat(B, &Bt));
1995:               PetscCall(PetscObjectTypeCompare((PetscObject)Bt, MATSEQAIJ, &fast));
1996:             } else {
1997:               PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &istrans));
1998:               if (istrans) {
1999:                 Mat Bt;

2001:                 PetscCall(MatHermitianTransposeGetMat(B, &Bt));
2002:                 PetscCall(PetscObjectTypeCompare((PetscObject)Bt, MATSEQAIJ, &fast));
2003:               }
2004:             }
2005:           }
2006:         }
2007:       }
2008:     }
2009:     for (i = 0, nf = 0; i < nest->nr && fast; ++i) {
2010:       PetscCall(PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i], ISSTRIDE, &fast));
2011:       if (fast) {
2012:         PetscInt f, s;

2014:         PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &f, &s));
2015:         if (f != nf || s != 1) {
2016:           fast = PETSC_FALSE;
2017:         } else {
2018:           PetscCall(ISGetSize(nest->isglobal.row[i], &f));
2019:           nf += f;
2020:         }
2021:       }
2022:     }
2023:     for (i = 0, nf = 0; i < nest->nc && fast; ++i) {
2024:       PetscCall(PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i], ISSTRIDE, &fast));
2025:       if (fast) {
2026:         PetscInt f, s;

2028:         PetscCall(ISStrideGetInfo(nest->isglobal.col[i], &f, &s));
2029:         if (f != nf || s != 1) {
2030:           fast = PETSC_FALSE;
2031:         } else {
2032:           PetscCall(ISGetSize(nest->isglobal.col[i], &f));
2033:           nf += f;
2034:         }
2035:       }
2036:     }
2037:     if (fast) {
2038:       PetscCall(MatConvert_Nest_SeqAIJ_fast(A, newtype, reuse, newmat));
2039:       PetscFunctionReturn(PETSC_SUCCESS);
2040:     }
2041:   }
2042:   PetscCall(MatGetSize(A, &M, &N));
2043:   PetscCall(MatGetLocalSize(A, &m, &n));
2044:   PetscCall(MatGetOwnershipRangeColumn(A, &cstart, &cend));
2045:   if (reuse == MAT_REUSE_MATRIX) C = *newmat;
2046:   else {
2047:     PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &C));
2048:     PetscCall(MatSetType(C, newtype));
2049:     PetscCall(MatSetSizes(C, m, n, M, N));
2050:   }
2051:   PetscCall(PetscMalloc1(2 * m, &dnnz));
2052:   onnz = dnnz + m;
2053:   for (k = 0; k < m; k++) {
2054:     dnnz[k] = 0;
2055:     onnz[k] = 0;
2056:   }
2057:   for (j = 0; j < nest->nc; ++j) {
2058:     IS              bNis;
2059:     PetscInt        bN;
2060:     const PetscInt *bNindices;
2061:     PetscBool       flg;
2062:     /* Using global column indices and ISAllGather() is not scalable. */
2063:     PetscCall(ISAllGather(nest->isglobal.col[j], &bNis));
2064:     PetscCall(ISGetSize(bNis, &bN));
2065:     PetscCall(ISGetIndices(bNis, &bNindices));
2066:     for (i = 0; i < nest->nr; ++i) {
2067:       PetscSF         bmsf;
2068:       PetscSFNode    *iremote;
2069:       Mat             B = nest->m[i][j], D = NULL;
2070:       PetscInt        bm, *sub_dnnz, *sub_onnz, br;
2071:       const PetscInt *bmindices;
2072:       if (!B) continue;
2073:       PetscCall(ISGetLocalSize(nest->isglobal.row[i], &bm));
2074:       PetscCall(ISGetIndices(nest->isglobal.row[i], &bmindices));
2075:       PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)A), &bmsf));
2076:       PetscCall(PetscMalloc1(bm, &iremote));
2077:       PetscCall(PetscMalloc1(bm, &sub_dnnz));
2078:       PetscCall(PetscMalloc1(bm, &sub_onnz));
2079:       for (k = 0; k < bm; ++k) {
2080:         sub_dnnz[k] = 0;
2081:         sub_onnz[k] = 0;
2082:       }
2083:       PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, ""));
2084:       if (flg) {
2085:         PetscTryMethod(B, "MatTransposeGetMat_C", (Mat, Mat *), (B, &D));
2086:         PetscTryMethod(B, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (B, &D));
2087:         PetscCall(MatConvert(B, ((PetscObject)D)->type_name, MAT_INITIAL_MATRIX, &D));
2088:         B = D;
2089:       }
2090:       PetscCall(PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));
2091:       if (flg) {
2092:         if (D) PetscCall(MatConvert(D, MATBAIJ, MAT_INPLACE_MATRIX, &D));
2093:         else PetscCall(MatConvert(B, MATBAIJ, MAT_INITIAL_MATRIX, &D));
2094:         B = D;
2095:       }
2096:       /*
2097:        Locate the owners for all of the locally-owned global row indices for this row block.
2098:        These determine the roots of PetscSF used to communicate preallocation data to row owners.
2099:        The roots correspond to the dnnz and onnz entries; thus, there are two roots per row.
2100:        */
2101:       PetscCall(MatGetOwnershipRange(B, &rstart, NULL));
2102:       for (br = 0; br < bm; ++br) {
2103:         PetscInt        row = bmindices[br], brncols, col;
2104:         const PetscInt *brcols;
2105:         PetscInt        rowrel   = 0; /* row's relative index on its owner rank */
2106:         PetscMPIInt     rowowner = 0;
2107:         PetscCall(PetscLayoutFindOwnerIndex(A->rmap, row, &rowowner, &rowrel));
2108:         /* how many roots  */
2109:         iremote[br].rank  = rowowner;
2110:         iremote[br].index = rowrel; /* edge from bmdnnz to dnnz */
2111:         /* get nonzero pattern */
2112:         PetscCall(MatGetRow(B, br + rstart, &brncols, &brcols, NULL));
2113:         for (k = 0; k < brncols; k++) {
2114:           col = bNindices[brcols[k]];
2115:           if (col >= A->cmap->range[rowowner] && col < A->cmap->range[rowowner + 1]) {
2116:             sub_dnnz[br]++;
2117:           } else {
2118:             sub_onnz[br]++;
2119:           }
2120:         }
2121:         PetscCall(MatRestoreRow(B, br + rstart, &brncols, &brcols, NULL));
2122:       }
2123:       if (D) PetscCall(MatDestroy(&D));
2124:       PetscCall(ISRestoreIndices(nest->isglobal.row[i], &bmindices));
2125:       /* bsf will have to take care of disposing of bedges. */
2126:       PetscCall(PetscSFSetGraph(bmsf, m, bm, NULL, PETSC_OWN_POINTER, iremote, PETSC_OWN_POINTER));
2127:       PetscCall(PetscSFReduceBegin(bmsf, MPIU_INT, sub_dnnz, dnnz, MPI_SUM));
2128:       PetscCall(PetscSFReduceEnd(bmsf, MPIU_INT, sub_dnnz, dnnz, MPI_SUM));
2129:       PetscCall(PetscSFReduceBegin(bmsf, MPIU_INT, sub_onnz, onnz, MPI_SUM));
2130:       PetscCall(PetscSFReduceEnd(bmsf, MPIU_INT, sub_onnz, onnz, MPI_SUM));
2131:       PetscCall(PetscFree(sub_dnnz));
2132:       PetscCall(PetscFree(sub_onnz));
2133:       PetscCall(PetscSFDestroy(&bmsf));
2134:     }
2135:     PetscCall(ISRestoreIndices(bNis, &bNindices));
2136:     PetscCall(ISDestroy(&bNis));
2137:   }
2138:   /* Resize preallocation if overestimated */
2139:   for (i = 0; i < m; i++) {
2140:     dnnz[i] = PetscMin(dnnz[i], A->cmap->n);
2141:     onnz[i] = PetscMin(onnz[i], A->cmap->N - A->cmap->n);
2142:   }
2143:   PetscCall(MatSeqAIJSetPreallocation(C, 0, dnnz));
2144:   PetscCall(MatMPIAIJSetPreallocation(C, 0, dnnz, 0, onnz));
2145:   PetscCall(PetscFree(dnnz));
2146:   PetscCall(MatAXPY_Dense_Nest(C, 1.0, A));
2147:   if (reuse == MAT_INPLACE_MATRIX) {
2148:     PetscCall(MatHeaderReplace(A, &C));
2149:   } else *newmat = C;
2150:   PetscFunctionReturn(PETSC_SUCCESS);
2151: }

2153: PetscErrorCode MatConvert_Nest_Dense(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
2154: {
2155:   Mat      B;
2156:   PetscInt m, n, M, N;

2158:   PetscFunctionBegin;
2159:   PetscCall(MatGetSize(A, &M, &N));
2160:   PetscCall(MatGetLocalSize(A, &m, &n));
2161:   if (reuse == MAT_REUSE_MATRIX) {
2162:     B = *newmat;
2163:     PetscCall(MatZeroEntries(B));
2164:   } else {
2165:     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), m, PETSC_DECIDE, M, N, NULL, &B));
2166:   }
2167:   PetscCall(MatAXPY_Dense_Nest(B, 1.0, A));
2168:   if (reuse == MAT_INPLACE_MATRIX) {
2169:     PetscCall(MatHeaderReplace(A, &B));
2170:   } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
2171:   PetscFunctionReturn(PETSC_SUCCESS);
2172: }

2174: PetscErrorCode MatHasOperation_Nest(Mat mat, MatOperation op, PetscBool *has)
2175: {
2176:   Mat_Nest    *bA = (Mat_Nest *)mat->data;
2177:   MatOperation opAdd;
2178:   PetscInt     i, j, nr = bA->nr, nc = bA->nc;
2179:   PetscBool    flg;
2180:   PetscFunctionBegin;

2182:   *has = PETSC_FALSE;
2183:   if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
2184:     opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD);
2185:     for (j = 0; j < nc; j++) {
2186:       for (i = 0; i < nr; i++) {
2187:         if (!bA->m[i][j]) continue;
2188:         PetscCall(MatHasOperation(bA->m[i][j], opAdd, &flg));
2189:         if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
2190:       }
2191:     }
2192:   }
2193:   if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
2194:   PetscFunctionReturn(PETSC_SUCCESS);
2195: }

2197: /*MC
2198:   MATNEST -  "nest" - Matrix type consisting of nested submatrices, each stored separately.

2200:   Level: intermediate

2202:   Notes:
2203:   This matrix type permits scalable use of `PCFIELDSPLIT` and avoids the large memory costs of extracting submatrices.
2204:   It allows the use of symmetric and block formats for parts of multi-physics simulations.
2205:   It is usually used with `DMCOMPOSITE` and `DMCreateMatrix()`

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

2211: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreate()`, `MatType`, `MatCreateNest()`, `MatNestSetSubMat()`, `MatNestGetSubMat()`,
2212:           `VecCreateNest()`, `DMCreateMatrix()`, `DMCOMPOSITE`, `MatNestSetVecType()`, `MatNestGetLocalISs()`,
2213:           `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()`
2214: M*/
2215: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2216: {
2217:   Mat_Nest *s;

2219:   PetscFunctionBegin;
2220:   PetscCall(PetscNew(&s));
2221:   A->data = (void *)s;

2223:   s->nr            = -1;
2224:   s->nc            = -1;
2225:   s->m             = NULL;
2226:   s->splitassembly = PETSC_FALSE;

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

2230:   A->ops->mult                  = MatMult_Nest;
2231:   A->ops->multadd               = MatMultAdd_Nest;
2232:   A->ops->multtranspose         = MatMultTranspose_Nest;
2233:   A->ops->multtransposeadd      = MatMultTransposeAdd_Nest;
2234:   A->ops->transpose             = MatTranspose_Nest;
2235:   A->ops->assemblybegin         = MatAssemblyBegin_Nest;
2236:   A->ops->assemblyend           = MatAssemblyEnd_Nest;
2237:   A->ops->zeroentries           = MatZeroEntries_Nest;
2238:   A->ops->copy                  = MatCopy_Nest;
2239:   A->ops->axpy                  = MatAXPY_Nest;
2240:   A->ops->duplicate             = MatDuplicate_Nest;
2241:   A->ops->createsubmatrix       = MatCreateSubMatrix_Nest;
2242:   A->ops->destroy               = MatDestroy_Nest;
2243:   A->ops->view                  = MatView_Nest;
2244:   A->ops->getvecs               = NULL; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2245:   A->ops->getlocalsubmatrix     = MatGetLocalSubMatrix_Nest;
2246:   A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2247:   A->ops->getdiagonal           = MatGetDiagonal_Nest;
2248:   A->ops->diagonalscale         = MatDiagonalScale_Nest;
2249:   A->ops->scale                 = MatScale_Nest;
2250:   A->ops->shift                 = MatShift_Nest;
2251:   A->ops->diagonalset           = MatDiagonalSet_Nest;
2252:   A->ops->setrandom             = MatSetRandom_Nest;
2253:   A->ops->hasoperation          = MatHasOperation_Nest;
2254:   A->ops->missingdiagonal       = MatMissingDiagonal_Nest;

2256:   A->spptr     = NULL;
2257:   A->assembled = PETSC_FALSE;

2259:   /* expose Nest api's */
2260:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMat_C", MatNestGetSubMat_Nest));
2261:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMat_C", MatNestSetSubMat_Nest));
2262:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMats_C", MatNestGetSubMats_Nest));
2263:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSize_C", MatNestGetSize_Nest));
2264:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetISs_C", MatNestGetISs_Nest));
2265:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetLocalISs_C", MatNestGetLocalISs_Nest));
2266:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetVecType_C", MatNestSetVecType_Nest));
2267:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMats_C", MatNestSetSubMats_Nest));
2268:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ));
2269:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ));
2270:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_aij_C", MatConvert_Nest_AIJ));
2271:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_is_C", MatConvert_Nest_IS));
2272:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpidense_C", MatConvert_Nest_Dense));
2273:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqdense_C", MatConvert_Nest_Dense));
2274:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_seqdense_C", MatProductSetFromOptions_Nest_Dense));
2275:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_mpidense_C", MatProductSetFromOptions_Nest_Dense));
2276:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_dense_C", MatProductSetFromOptions_Nest_Dense));

2278:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATNEST));
2279:   PetscFunctionReturn(PETSC_SUCCESS);
2280: }