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: static 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, 1);
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: static 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: static 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, 1);
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(Mat C)
251: {
252:   Mat_Product *product = C->product;

254:   PetscFunctionBegin;
255:   if (product->type == MATPRODUCT_AB) C->ops->productsymbolic = MatProductSymbolic_Nest_Dense;
256:   PetscFunctionReturn(PETSC_SUCCESS);
257: }

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

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

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

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

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

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

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

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

331:     PetscCall(MatCreateNest(PetscObjectComm((PetscObject)A), nc, is_col, nr, is_row, subs, &C));
332:     PetscCall(PetscFree(subs));
333:     PetscCall(PetscFree2(is_row, is_col));
334:   } else {
335:     C = *B;
336:   }

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

349:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
350:     *B = C;
351:   } else {
352:     PetscCall(MatHeaderMerge(A, &C));
353:   }
354:   PetscFunctionReturn(PETSC_SUCCESS);
355: }

357: static PetscErrorCode MatNestDestroyISList(PetscInt n, IS **list)
358: {
359:   IS      *lst = *list;
360:   PetscInt i;

362:   PetscFunctionBegin;
363:   if (!lst) PetscFunctionReturn(PETSC_SUCCESS);
364:   for (i = 0; i < n; i++)
365:     if (lst[i]) PetscCall(ISDestroy(&lst[i]));
366:   PetscCall(PetscFree(lst));
367:   *list = NULL;
368:   PetscFunctionReturn(PETSC_SUCCESS);
369: }

371: static PetscErrorCode MatReset_Nest(Mat A)
372: {
373:   Mat_Nest *vs = (Mat_Nest *)A->data;
374:   PetscInt  i, j;

376:   PetscFunctionBegin;
377:   /* release the matrices and the place holders */
378:   PetscCall(MatNestDestroyISList(vs->nr, &vs->isglobal.row));
379:   PetscCall(MatNestDestroyISList(vs->nc, &vs->isglobal.col));
380:   PetscCall(MatNestDestroyISList(vs->nr, &vs->islocal.row));
381:   PetscCall(MatNestDestroyISList(vs->nc, &vs->islocal.col));

383:   PetscCall(PetscFree(vs->row_len));
384:   PetscCall(PetscFree(vs->col_len));
385:   PetscCall(PetscFree(vs->nnzstate));

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

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

398:   /* restore defaults */
399:   vs->nr            = 0;
400:   vs->nc            = 0;
401:   vs->splitassembly = PETSC_FALSE;
402:   PetscFunctionReturn(PETSC_SUCCESS);
403: }

405: static PetscErrorCode MatDestroy_Nest(Mat A)
406: {
407:   PetscFunctionBegin;
408:   PetscCall(MatReset_Nest(A));
409:   PetscCall(PetscFree(A->data));
410:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMat_C", NULL));
411:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMat_C", NULL));
412:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMats_C", NULL));
413:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSize_C", NULL));
414:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetISs_C", NULL));
415:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetLocalISs_C", NULL));
416:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetVecType_C", NULL));
417:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMats_C", NULL));
418:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpiaij_C", NULL));
419:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqaij_C", NULL));
420:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_aij_C", NULL));
421:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_is_C", NULL));
422:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpidense_C", NULL));
423:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqdense_C", NULL));
424:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_seqdense_C", NULL));
425:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_mpidense_C", NULL));
426:   PetscFunctionReturn(PETSC_SUCCESS);
427: }

429: static PetscErrorCode MatMissingDiagonal_Nest(Mat mat, PetscBool *missing, PetscInt *dd)
430: {
431:   Mat_Nest *vs = (Mat_Nest *)mat->data;
432:   PetscInt  i;

434:   PetscFunctionBegin;
435:   if (dd) *dd = 0;
436:   if (!vs->nr) {
437:     *missing = PETSC_TRUE;
438:     PetscFunctionReturn(PETSC_SUCCESS);
439:   }
440:   *missing = PETSC_FALSE;
441:   for (i = 0; i < vs->nr && !(*missing); i++) {
442:     *missing = PETSC_TRUE;
443:     if (vs->m[i][i]) {
444:       PetscCall(MatMissingDiagonal(vs->m[i][i], missing, NULL));
445:       PetscCheck(!*missing || !dd, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "First missing entry not yet implemented");
446:     }
447:   }
448:   PetscFunctionReturn(PETSC_SUCCESS);
449: }

451: static PetscErrorCode MatAssemblyBegin_Nest(Mat A, MatAssemblyType type)
452: {
453:   Mat_Nest *vs = (Mat_Nest *)A->data;
454:   PetscInt  i, j;
455:   PetscBool nnzstate = PETSC_FALSE;

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

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

487:   PetscFunctionBegin;
488:   for (i = 0; i < vs->nr; i++) {
489:     for (j = 0; j < vs->nc; j++) {
490:       if (vs->m[i][j]) {
491:         if (vs->splitassembly) PetscCall(MatAssemblyEnd(vs->m[i][j], type));
492:       }
493:     }
494:   }
495:   PetscFunctionReturn(PETSC_SUCCESS);
496: }

498: static PetscErrorCode MatNestFindNonzeroSubMatRow(Mat A, PetscInt row, Mat *B)
499: {
500:   Mat_Nest *vs = (Mat_Nest *)A->data;
501:   PetscInt  j;
502:   Mat       sub;

504:   PetscFunctionBegin;
505:   sub = (row < vs->nc) ? vs->m[row][row] : (Mat)NULL; /* Prefer to find on the diagonal */
506:   for (j = 0; !sub && j < vs->nc; j++) sub = vs->m[row][j];
507:   if (sub) PetscCall(MatSetUp(sub)); /* Ensure that the sizes are available */
508:   *B = sub;
509:   PetscFunctionReturn(PETSC_SUCCESS);
510: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

717:   PetscFunctionBegin;
718:   PetscCall(MatNestFindSubMat(A, &vs->islocal, isrow, iscol, &sub));
719:   PetscCheck(*B == sub, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Local submatrix has not been gotten");
720:   if (sub) {
721:     PetscCheck(((PetscObject)sub)->refct > 1, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Local submatrix has had reference count decremented too many times");
722:     PetscCall(MatDestroy(B));
723:   }
724:   PetscFunctionReturn(PETSC_SUCCESS);
725: }

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

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

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

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

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

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

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

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

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

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

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

833:   PetscFunctionBegin;
834:   for (i = 0; i < bA->nr; i++) {
835:     for (j = 0; j < bA->nc; j++) {
836:       if (bA->m[i][j]) PetscCall(MatSetRandom(bA->m[i][j], rctx));
837:     }
838:   }
839:   PetscFunctionReturn(PETSC_SUCCESS);
840: }

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

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

870:   if (left) {
871:     /* allocate L */
872:     PetscCall(PetscMalloc1(bA->nr, &L));
873:     /* Create the left vectors */
874:     for (i = 0; i < bA->nr; i++) {
875:       for (j = 0; j < bA->nc; j++) {
876:         if (bA->m[i][j]) {
877:           PetscCall(MatCreateVecs(bA->m[i][j], NULL, &L[i]));
878:           break;
879:         }
880:       }
881:       PetscCheck(j != bA->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Mat(Nest) contains a null row.");
882:     }

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

887:     PetscCall(PetscFree(L));
888:   }
889:   PetscFunctionReturn(PETSC_SUCCESS);
890: }

892: static PetscErrorCode MatView_Nest(Mat A, PetscViewer viewer)
893: {
894:   Mat_Nest *bA = (Mat_Nest *)A->data;
895:   PetscBool isascii, viewSub = PETSC_FALSE;
896:   PetscInt  i, j;

898:   PetscFunctionBegin;
899:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
900:   if (isascii) {
901:     PetscCall(PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_view_nest_sub", &viewSub, NULL));
902:     PetscCall(PetscViewerASCIIPrintf(viewer, "Matrix object:\n"));
903:     PetscCall(PetscViewerASCIIPushTab(viewer));
904:     PetscCall(PetscViewerASCIIPrintf(viewer, "type=nest, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT "\n", bA->nr, bA->nc));

906:     PetscCall(PetscViewerASCIIPrintf(viewer, "MatNest structure:\n"));
907:     for (i = 0; i < bA->nr; i++) {
908:       for (j = 0; j < bA->nc; j++) {
909:         MatType   type;
910:         char      name[256] = "", prefix[256] = "";
911:         PetscInt  NR, NC;
912:         PetscBool isNest = PETSC_FALSE;

914:         if (!bA->m[i][j]) {
915:           PetscCall(PetscViewerASCIIPrintf(viewer, "(%" PetscInt_FMT ",%" PetscInt_FMT ") : NULL\n", i, j));
916:           continue;
917:         }
918:         PetscCall(MatGetSize(bA->m[i][j], &NR, &NC));
919:         PetscCall(MatGetType(bA->m[i][j], &type));
920:         if (((PetscObject)bA->m[i][j])->name) PetscCall(PetscSNPrintf(name, sizeof(name), "name=\"%s\", ", ((PetscObject)bA->m[i][j])->name));
921:         if (((PetscObject)bA->m[i][j])->prefix) PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "prefix=\"%s\", ", ((PetscObject)bA->m[i][j])->prefix));
922:         PetscCall(PetscObjectTypeCompare((PetscObject)bA->m[i][j], MATNEST, &isNest));

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

926:         if (isNest || viewSub) {
927:           PetscCall(PetscViewerASCIIPushTab(viewer)); /* push1 */
928:           PetscCall(MatView(bA->m[i][j], viewer));
929:           PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop1 */
930:         }
931:       }
932:     }
933:     PetscCall(PetscViewerASCIIPopTab(viewer)); /* pop0 */
934:   }
935:   PetscFunctionReturn(PETSC_SUCCESS);
936: }

938: static PetscErrorCode MatZeroEntries_Nest(Mat A)
939: {
940:   Mat_Nest *bA = (Mat_Nest *)A->data;
941:   PetscInt  i, j;

943:   PetscFunctionBegin;
944:   for (i = 0; i < bA->nr; i++) {
945:     for (j = 0; j < bA->nc; j++) {
946:       if (!bA->m[i][j]) continue;
947:       PetscCall(MatZeroEntries(bA->m[i][j]));
948:     }
949:   }
950:   PetscFunctionReturn(PETSC_SUCCESS);
951: }

953: static PetscErrorCode MatCopy_Nest(Mat A, Mat B, MatStructure str)
954: {
955:   Mat_Nest *bA = (Mat_Nest *)A->data, *bB = (Mat_Nest *)B->data;
956:   PetscInt  i, j, nr = bA->nr, nc = bA->nc;
957:   PetscBool nnzstate = PETSC_FALSE;

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

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

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

992:         PetscCheck(str == DIFFERENT_NONZERO_PATTERN || str == UNKNOWN_NONZERO_PATTERN, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_INCOMP, "Matrix block does not exist at %" PetscInt_FMT ",%" PetscInt_FMT ". Use DIFFERENT_NONZERO_PATTERN or UNKNOWN_NONZERO_PATTERN", i, j);
993:         PetscCall(MatDuplicate(bX->m[i][j], MAT_COPY_VALUES, &M));
994:         PetscCall(MatNestSetSubMat(Y, i, j, M));
995:         PetscCall(MatDestroy(&M));
996:       }
997:       if (bY->m[i][j]) PetscCall(MatGetNonzeroState(bY->m[i][j], &subnnzstate));
998:       nnzstate                 = (PetscBool)(nnzstate || bY->nnzstate[i * nc + j] != subnnzstate);
999:       bY->nnzstate[i * nc + j] = subnnzstate;
1000:     }
1001:   }
1002:   if (nnzstate) Y->nonzerostate++;
1003:   PetscFunctionReturn(PETSC_SUCCESS);
1004: }

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

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

1028:   PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
1029:   PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
1030:   PetscFunctionReturn(PETSC_SUCCESS);
1031: }

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

1038:   PetscFunctionBegin;
1039:   PetscCheck(idxm < bA->nr, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, idxm, bA->nr - 1);
1040:   PetscCheck(jdxm < bA->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Col too large: row %" PetscInt_FMT " max %" PetscInt_FMT, jdxm, bA->nc - 1);
1041:   *mat = bA->m[idxm][jdxm];
1042:   PetscFunctionReturn(PETSC_SUCCESS);
1043: }

1045: /*@
1046:   MatNestGetSubMat - Returns a single, sub-matrix from a `MATNEST`

1048:   Not Collective

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

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

1058:   Level: developer

1060: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSize()`, `MatNestGetSubMats()`, `MatCreateNest()`, `MatNestSetSubMat()`,
1061:           `MatNestGetLocalISs()`, `MatNestGetISs()`
1062: @*/
1063: PetscErrorCode MatNestGetSubMat(Mat A, PetscInt idxm, PetscInt jdxm, Mat *sub)
1064: {
1065:   PetscFunctionBegin;
1066:   PetscUseMethod(A, "MatNestGetSubMat_C", (Mat, PetscInt, PetscInt, Mat *), (A, idxm, jdxm, sub));
1067:   PetscFunctionReturn(PETSC_SUCCESS);
1068: }

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

1075:   PetscFunctionBegin;
1076:   PetscCheck(idxm < bA->nr, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, idxm, bA->nr - 1);
1077:   PetscCheck(jdxm < bA->nc, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Col too large: row %" PetscInt_FMT " max %" PetscInt_FMT, jdxm, bA->nc - 1);
1078:   PetscCall(MatGetLocalSize(mat, &m, &n));
1079:   PetscCall(MatGetSize(mat, &M, &N));
1080:   PetscCall(ISGetLocalSize(bA->isglobal.row[idxm], &mi));
1081:   PetscCall(ISGetSize(bA->isglobal.row[idxm], &Mi));
1082:   PetscCall(ISGetLocalSize(bA->isglobal.col[jdxm], &ni));
1083:   PetscCall(ISGetSize(bA->isglobal.col[jdxm], &Ni));
1084:   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);
1085:   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);

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

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

1099: /*@
1100:   MatNestSetSubMat - Set a single submatrix in the `MATNEST`

1102:   Logically Collective

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

1110:   Level: developer

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

1115:   This increments the reference count of the submatrix.

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

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

1131:   PetscFunctionBegin;
1132:   if (M) *M = bA->nr;
1133:   if (N) *N = bA->nc;
1134:   if (mat) *mat = bA->m;
1135:   PetscFunctionReturn(PETSC_SUCCESS);
1136: }

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

1141:   Not Collective

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

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

1151:   Level: developer

1153:   Note:
1154:   The user should not free the array `mat`.

1156:   Fortran Notes:
1157:   This routine has a calling sequence
1158: $   call MatNestGetSubMats(A, M, N, mat, ierr)
1159:   where the space allocated for the optional argument `mat` is assumed large enough (if provided).
1160:   Matrices in `mat` are returned in row-major order, see `MatCreateNest()` for an example.

1162: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSize()`, `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MatCreateNest()`,
1163:           `MatNestSetSubMats()`, `MatNestGetISs()`, `MatNestSetSubMat()`
1164: @*/
1165: PetscErrorCode MatNestGetSubMats(Mat A, PetscInt *M, PetscInt *N, Mat ***mat)
1166: {
1167:   PetscFunctionBegin;
1168:   PetscUseMethod(A, "MatNestGetSubMats_C", (Mat, PetscInt *, PetscInt *, Mat ***), (A, M, N, mat));
1169:   PetscFunctionReturn(PETSC_SUCCESS);
1170: }

1172: static PetscErrorCode MatNestGetSize_Nest(Mat A, PetscInt *M, PetscInt *N)
1173: {
1174:   Mat_Nest *bA = (Mat_Nest *)A->data;

1176:   PetscFunctionBegin;
1177:   if (M) *M = bA->nr;
1178:   if (N) *N = bA->nc;
1179:   PetscFunctionReturn(PETSC_SUCCESS);
1180: }

1182: /*@
1183:   MatNestGetSize - Returns the size of the `MATNEST` matrix.

1185:   Not Collective

1187:   Input Parameter:
1188: . A - `MATNEST` matrix

1190:   Output Parameters:
1191: + M - number of rows in the nested mat
1192: - N - number of cols in the nested mat

1194:   Level: developer

1196: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatCreateNest()`, `MatNestGetLocalISs()`,
1197:           `MatNestGetISs()`
1198: @*/
1199: PetscErrorCode MatNestGetSize(Mat A, PetscInt *M, PetscInt *N)
1200: {
1201:   PetscFunctionBegin;
1202:   PetscUseMethod(A, "MatNestGetSize_C", (Mat, PetscInt *, PetscInt *), (A, M, N));
1203:   PetscFunctionReturn(PETSC_SUCCESS);
1204: }

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

1211:   PetscFunctionBegin;
1212:   if (rows)
1213:     for (i = 0; i < vs->nr; i++) rows[i] = vs->isglobal.row[i];
1214:   if (cols)
1215:     for (i = 0; i < vs->nc; i++) cols[i] = vs->isglobal.col[i];
1216:   PetscFunctionReturn(PETSC_SUCCESS);
1217: }

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

1222:   Not Collective

1224:   Input Parameter:
1225: . A - `MATNEST` matrix

1227:   Output Parameters:
1228: + rows - array of row index sets
1229: - cols - array of column index sets

1231:   Level: advanced

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

1236: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetLocalISs()`,
1237:           `MatCreateNest()`, `MatNestSetSubMats()`
1238: @*/
1239: PetscErrorCode MatNestGetISs(Mat A, IS rows[], IS cols[])
1240: {
1241:   PetscFunctionBegin;
1243:   PetscUseMethod(A, "MatNestGetISs_C", (Mat, IS[], IS[]), (A, rows, cols));
1244:   PetscFunctionReturn(PETSC_SUCCESS);
1245: }

1247: static PetscErrorCode MatNestGetLocalISs_Nest(Mat A, IS rows[], IS cols[])
1248: {
1249:   Mat_Nest *vs = (Mat_Nest *)A->data;
1250:   PetscInt  i;

1252:   PetscFunctionBegin;
1253:   if (rows)
1254:     for (i = 0; i < vs->nr; i++) rows[i] = vs->islocal.row[i];
1255:   if (cols)
1256:     for (i = 0; i < vs->nc; i++) cols[i] = vs->islocal.col[i];
1257:   PetscFunctionReturn(PETSC_SUCCESS);
1258: }

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

1263:   Not Collective

1265:   Input Parameter:
1266: . A - `MATNEST` matrix

1268:   Output Parameters:
1269: + rows - array of row index sets (or `NULL` to ignore)
1270: - cols - array of column index sets (or `NULL` to ignore)

1272:   Level: advanced

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

1277: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetISs()`, `MatCreateNest()`,
1278:           `MatNestSetSubMats()`, `MatNestSetSubMat()`
1279: @*/
1280: PetscErrorCode MatNestGetLocalISs(Mat A, IS rows[], IS cols[])
1281: {
1282:   PetscFunctionBegin;
1284:   PetscUseMethod(A, "MatNestGetLocalISs_C", (Mat, IS[], IS[]), (A, rows, cols));
1285:   PetscFunctionReturn(PETSC_SUCCESS);
1286: }

1288: static PetscErrorCode MatNestSetVecType_Nest(Mat A, VecType vtype)
1289: {
1290:   PetscBool flg;

1292:   PetscFunctionBegin;
1293:   PetscCall(PetscStrcmp(vtype, VECNEST, &flg));
1294:   /* In reality, this only distinguishes VECNEST and "other" */
1295:   if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1296:   else A->ops->getvecs = (PetscErrorCode(*)(Mat, Vec *, Vec *))0;
1297:   PetscFunctionReturn(PETSC_SUCCESS);
1298: }

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

1303:   Not Collective

1305:   Input Parameters:
1306: + A     - `MATNEST` matrix
1307: - vtype - `VecType` to use for creating vectors

1309:   Level: developer

1311: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreateVecs()`, `MatCreateNest()`, `VecType`
1312: @*/
1313: PetscErrorCode MatNestSetVecType(Mat A, VecType vtype)
1314: {
1315:   PetscFunctionBegin;
1316:   PetscTryMethod(A, "MatNestSetVecType_C", (Mat, VecType), (A, vtype));
1317:   PetscFunctionReturn(PETSC_SUCCESS);
1318: }

1320: static PetscErrorCode MatNestSetSubMats_Nest(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[])
1321: {
1322:   Mat_Nest *s = (Mat_Nest *)A->data;
1323:   PetscInt  i, j, m, n, M, N;
1324:   PetscBool cong, isstd, sametype = PETSC_FALSE;
1325:   VecType   vtype, type;

1327:   PetscFunctionBegin;
1328:   PetscCall(MatReset_Nest(A));

1330:   s->nr = nr;
1331:   s->nc = nc;

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

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

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

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

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

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

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

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

1403:   PetscCall(PetscCalloc2(nr, &s->left, nc, &s->right));
1404:   PetscCall(PetscObjectStateIncrease((PetscObject)A));
1405:   A->nonzerostate++;
1406:   PetscFunctionReturn(PETSC_SUCCESS);
1407: }

1409: /*@
1410:   MatNestSetSubMats - Sets the nested submatrices in a `MATNEST`

1412:   Collective

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

1422:   Level: advanced

1424:   Notes:
1425:   This always resets any submatrix information previously set

1427:   In both C and Fortran, `a` must be a row-major order array containing the matrices. See
1428:   `MatCreateNest()` for an example.

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

1436:   PetscFunctionBegin;
1438:   PetscCheck(nr >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Number of rows cannot be negative");
1439:   if (nr && is_row) {
1440:     PetscAssertPointer(is_row, 3);
1442:   }
1443:   PetscCheck(nc >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Number of columns cannot be negative");
1444:   if (nc && is_col) {
1445:     PetscAssertPointer(is_col, 5);
1447:   }
1448:   if (nr * nc > 0) PetscAssertPointer(a, 6);
1449:   PetscUseMethod(A, "MatNestSetSubMats_C", (Mat, PetscInt, const IS[], PetscInt, const IS[], const Mat[]), (A, nr, is_row, nc, is_col, a));
1450:   PetscFunctionReturn(PETSC_SUCCESS);
1451: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1702:   Collective

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

1712:   Output Parameter:
1713: . B - new matrix

1715:   Note:
1716:   In both C and Fortran, `a` must be a row-major order array holding references to the matrices.
1717:   For instance, to represent the matrix
1718:   $\begin{bmatrix} A_{11} & A_{12} \\ A_{21} & A_{22}\end{bmatrix}$
1719:   one should use `Mat a[4]={A11,A12,A21,A22}`.

1721:   Level: advanced

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2157: static PetscErrorCode MatConvert_Nest_Dense(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
2158: {
2159:   Mat      B;
2160:   PetscInt m, n, M, N;

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

2178: static PetscErrorCode MatHasOperation_Nest(Mat mat, MatOperation op, PetscBool *has)
2179: {
2180:   Mat_Nest    *bA = (Mat_Nest *)mat->data;
2181:   MatOperation opAdd;
2182:   PetscInt     i, j, nr = bA->nr, nc = bA->nc;
2183:   PetscBool    flg;
2184:   PetscFunctionBegin;

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

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

2204:   Level: intermediate

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

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

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

2223:   PetscFunctionBegin;
2224:   PetscCall(PetscNew(&s));
2225:   A->data = (void *)s;

2227:   s->nr            = -1;
2228:   s->nc            = -1;
2229:   s->m             = NULL;
2230:   s->splitassembly = PETSC_FALSE;

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

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

2260:   A->spptr     = NULL;
2261:   A->assembled = PETSC_FALSE;

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

2281:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATNEST));
2282:   PetscFunctionReturn(PETSC_SUCCESS);
2283: }