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_AB(Mat C)
251: {
252:   PetscFunctionBegin;
253:   C->ops->productsymbolic = MatProductSymbolic_Nest_Dense;
254:   PetscFunctionReturn(PETSC_SUCCESS);
255: }

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

261:   PetscFunctionBegin;
262:   if (product->type == MATPRODUCT_AB) PetscCall(MatProductSetFromOptions_Nest_Dense_AB(C));
263:   PetscFunctionReturn(PETSC_SUCCESS);
264: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1000:         PetscCheck(str == DIFFERENT_NONZERO_PATTERN || 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);
1001:         PetscCall(MatDuplicate(bX->m[i][j], MAT_COPY_VALUES, &M));
1002:         PetscCall(MatNestSetSubMat(Y, i, j, M));
1003:         PetscCall(MatDestroy(&M));
1004:       }
1005:       if (bY->m[i][j]) PetscCall(MatGetNonzeroState(bY->m[i][j], &subnnzstate));
1006:       nnzstate                 = (PetscBool)(nnzstate || bY->nnzstate[i * nc + j] != subnnzstate);
1007:       bY->nnzstate[i * nc + j] = subnnzstate;
1008:     }
1009:   }
1010:   if (nnzstate) Y->nonzerostate++;
1011:   PetscFunctionReturn(PETSC_SUCCESS);
1012: }

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

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

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

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

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

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

1056:   Not Collective

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

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

1066:   Level: developer

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

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

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

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

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

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

1110:   Logically Collective

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

1118:   Level: developer

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

1123:   This increments the reference count of the submatrix.

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

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

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

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

1149:   Not Collective

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

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

1159:   Level: developer

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

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

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

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

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

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

1193:   Not Collective

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

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

1202:   Level: developer

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

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

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

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

1230:   Not Collective

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

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

1239:   Level: advanced

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

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

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

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

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

1271:   Not Collective

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

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

1280:   Level: advanced

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

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

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

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

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

1311:   Not Collective

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

1317:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1420:   Collective

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

1430:   Level: advanced

1432:   Notes:
1433:   This always resets any submatrix information previously set

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

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

1444:   PetscFunctionBegin;
1446:   PetscCheck(nr >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Number of rows cannot be negative");
1447:   if (nr && is_row) {
1448:     PetscAssertPointer(is_row, 3);
1450:   }
1451:   PetscCheck(nc >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Number of columns cannot be negative");
1452:   if (nc && is_col) {
1453:     PetscAssertPointer(is_col, 5);
1455:   }
1456:   if (nr * nc > 0) PetscAssertPointer(a, 6);
1457:   PetscUseMethod(A, "MatNestSetSubMats_C", (Mat, PetscInt, const IS[], PetscInt, const IS[], const Mat[]), (A, nr, is_row, nc, is_col, a));
1458:   PetscFunctionReturn(PETSC_SUCCESS);
1459: }

1461: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A, PetscInt n, const IS islocal[], const IS isglobal[], PetscBool colflg, ISLocalToGlobalMapping *ltog)
1462: {
1463:   PetscBool flg;
1464:   PetscInt  i, j, m, mi, *ix;

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

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

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

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

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

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

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

1555:             proc 0:
1556:     is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1557:             proc 1:
1558:     is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)

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

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

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

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

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

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

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

1679:   if (PetscDefined(USE_DEBUG)) {
1680:     for (i = 0; i < vs->nr; i++) {
1681:       for (j = 0; j < vs->nc; j++) {
1682:         PetscInt m, n, M, N, mi, ni, Mi, Ni;
1683:         Mat      B = vs->m[i][j];
1684:         if (!B) continue;
1685:         PetscCall(MatGetSize(B, &M, &N));
1686:         PetscCall(MatGetLocalSize(B, &m, &n));
1687:         PetscCall(ISGetSize(vs->isglobal.row[i], &Mi));
1688:         PetscCall(ISGetSize(vs->isglobal.col[j], &Ni));
1689:         PetscCall(ISGetLocalSize(vs->isglobal.row[i], &mi));
1690:         PetscCall(ISGetLocalSize(vs->isglobal.col[j], &ni));
1691:         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);
1692:         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);
1693:       }
1694:     }
1695:   }

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

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

1710:   Collective

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

1720:   Output Parameter:
1721: . B - new matrix

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

1729:   Level: advanced

1731: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreate()`, `VecCreateNest()`, `DMCreateMatrix()`, `MatNestSetSubMat()`,
1732:           `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MatNestGetSize()`,
1733:           `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()`
1734: @*/
1735: PetscErrorCode MatCreateNest(MPI_Comm comm, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[], Mat *B)
1736: {
1737:   Mat A;

1739:   PetscFunctionBegin;
1740:   *B = NULL;
1741:   PetscCall(MatCreate(comm, &A));
1742:   PetscCall(MatSetType(A, MATNEST));
1743:   A->preallocated = PETSC_TRUE;
1744:   PetscCall(MatNestSetSubMats(A, nr, is_row, nc, is_col, a));
1745:   *B = A;
1746:   PetscFunctionReturn(PETSC_SUCCESS);
1747: }

1749: static PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1750: {
1751:   Mat_Nest     *nest = (Mat_Nest *)A->data;
1752:   Mat          *trans;
1753:   PetscScalar **avv;
1754:   PetscScalar  *vv;
1755:   PetscInt    **aii, **ajj;
1756:   PetscInt     *ii, *jj, *ci;
1757:   PetscInt      nr, nc, nnz, i, j;
1758:   PetscBool     done;

1760:   PetscFunctionBegin;
1761:   PetscCall(MatGetSize(A, &nr, &nc));
1762:   if (reuse == MAT_REUSE_MATRIX) {
1763:     PetscInt rnr;

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

1781:         PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &istrans));
1782:         if (istrans) {
1783:           Mat Bt;

1785:           PetscCall(MatTransposeGetMat(B, &Bt));
1786:           PetscCall(MatTranspose(Bt, MAT_INITIAL_MATRIX, &trans[i * nest->nc + j]));
1787:           B = trans[i * nest->nc + j];
1788:         } else {
1789:           PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &istrans));
1790:           if (istrans) {
1791:             Mat Bt;

1793:             PetscCall(MatHermitianTransposeGetMat(B, &Bt));
1794:             PetscCall(MatHermitianTranspose(Bt, MAT_INITIAL_MATRIX, &trans[i * nest->nc + j]));
1795:             B = trans[i * nest->nc + j];
1796:           }
1797:         }
1798:         PetscCall(MatGetRowIJ(B, 0, PETSC_FALSE, PETSC_FALSE, &nnr, (const PetscInt **)&nii, (const PetscInt **)&njj, &done));
1799:         PetscCheck(done, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "MatGetRowIJ");
1800:         PetscCall(MatSeqAIJGetArray(B, &naa));
1801:         nnz += nii[nnr];

1803:         aii[i * nest->nc + j] = nii;
1804:         ajj[i * nest->nc + j] = njj;
1805:         avv[i * nest->nc + j] = naa;
1806:       }
1807:     }
1808:   }
1809:   if (reuse != MAT_REUSE_MATRIX) {
1810:     PetscCall(PetscMalloc1(nr + 1, &ii));
1811:     PetscCall(PetscMalloc1(nnz, &jj));
1812:     PetscCall(PetscMalloc1(nnz, &vv));
1813:   } else {
1814:     PetscCheck(nnz == ii[nr], PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Cannot reuse matrix, wrong number of nonzeros");
1815:   }

1817:   /* new row pointer */
1818:   PetscCall(PetscArrayzero(ii, nr + 1));
1819:   for (i = 0; i < nest->nr; ++i) {
1820:     PetscInt ncr, rst;

1822:     PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &rst, NULL));
1823:     PetscCall(ISGetLocalSize(nest->isglobal.row[i], &ncr));
1824:     for (j = 0; j < nest->nc; ++j) {
1825:       if (aii[i * nest->nc + j]) {
1826:         PetscInt *nii = aii[i * nest->nc + j];
1827:         PetscInt  ir;

1829:         for (ir = rst; ir < ncr + rst; ++ir) {
1830:           ii[ir + 1] += nii[1] - nii[0];
1831:           nii++;
1832:         }
1833:       }
1834:     }
1835:   }
1836:   for (i = 0; i < nr; i++) ii[i + 1] += ii[i];

1838:   /* construct CSR for the new matrix */
1839:   PetscCall(PetscCalloc1(nr, &ci));
1840:   for (i = 0; i < nest->nr; ++i) {
1841:     PetscInt ncr, rst;

1843:     PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &rst, NULL));
1844:     PetscCall(ISGetLocalSize(nest->isglobal.row[i], &ncr));
1845:     for (j = 0; j < nest->nc; ++j) {
1846:       if (aii[i * nest->nc + j]) {
1847:         PetscScalar *nvv = avv[i * nest->nc + j];
1848:         PetscInt    *nii = aii[i * nest->nc + j];
1849:         PetscInt    *njj = ajj[i * nest->nc + j];
1850:         PetscInt     ir, cst;

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

1856:           for (ij = 0; ij < rsize; ij++) {
1857:             jj[ist + ij] = *njj + cst;
1858:             vv[ist + ij] = *nvv;
1859:             njj++;
1860:             nvv++;
1861:           }
1862:           ci[ir] += rsize;
1863:           nii++;
1864:         }
1865:       }
1866:     }
1867:   }
1868:   PetscCall(PetscFree(ci));

1870:   /* restore info */
1871:   for (i = 0; i < nest->nr; ++i) {
1872:     for (j = 0; j < nest->nc; ++j) {
1873:       Mat B = nest->m[i][j];
1874:       if (B) {
1875:         PetscInt nnr = 0, k = i * nest->nc + j;

1877:         B = (trans[k] ? trans[k] : B);
1878:         PetscCall(MatRestoreRowIJ(B, 0, PETSC_FALSE, PETSC_FALSE, &nnr, (const PetscInt **)&aii[k], (const PetscInt **)&ajj[k], &done));
1879:         PetscCheck(done, PetscObjectComm((PetscObject)B), PETSC_ERR_PLIB, "MatRestoreRowIJ");
1880:         PetscCall(MatSeqAIJRestoreArray(B, &avv[k]));
1881:         PetscCall(MatDestroy(&trans[k]));
1882:       }
1883:     }
1884:   }
1885:   PetscCall(PetscFree4(aii, ajj, avv, trans));

1887:   /* finalize newmat */
1888:   if (reuse == MAT_INITIAL_MATRIX) {
1889:     PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), nr, nc, ii, jj, vv, newmat));
1890:   } else if (reuse == MAT_INPLACE_MATRIX) {
1891:     Mat B;

1893:     PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), nr, nc, ii, jj, vv, &B));
1894:     PetscCall(MatHeaderReplace(A, &B));
1895:   }
1896:   PetscCall(MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY));
1897:   PetscCall(MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY));
1898:   {
1899:     Mat_SeqAIJ *a = (Mat_SeqAIJ *)((*newmat)->data);
1900:     a->free_a     = PETSC_TRUE;
1901:     a->free_ij    = PETSC_TRUE;
1902:   }
1903:   PetscFunctionReturn(PETSC_SUCCESS);
1904: }

1906: PETSC_INTERN PetscErrorCode MatAXPY_Dense_Nest(Mat Y, PetscScalar a, Mat X)
1907: {
1908:   Mat_Nest *nest = (Mat_Nest *)X->data;
1909:   PetscInt  i, j, k, rstart;
1910:   PetscBool flg;

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

1977: static PetscErrorCode MatConvert_Nest_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1978: {
1979:   Mat_Nest   *nest = (Mat_Nest *)A->data;
1980:   PetscInt    m, n, M, N, i, j, k, *dnnz, *onnz = NULL, rstart, cstart, cend;
1981:   PetscMPIInt size;
1982:   Mat         C;

1984:   PetscFunctionBegin;
1985:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1986:   if (size == 1) { /* look for a special case with SeqAIJ matrices and strided-1, contiguous, blocks */
1987:     PetscInt  nf;
1988:     PetscBool fast;

1990:     PetscCall(PetscStrcmp(newtype, MATAIJ, &fast));
1991:     if (!fast) PetscCall(PetscStrcmp(newtype, MATSEQAIJ, &fast));
1992:     for (i = 0; i < nest->nr && fast; ++i) {
1993:       for (j = 0; j < nest->nc && fast; ++j) {
1994:         Mat B = nest->m[i][j];
1995:         if (B) {
1996:           PetscCall(PetscObjectTypeCompare((PetscObject)B, MATSEQAIJ, &fast));
1997:           if (!fast) {
1998:             PetscBool istrans;

2000:             PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &istrans));
2001:             if (istrans) {
2002:               Mat Bt;

2004:               PetscCall(MatTransposeGetMat(B, &Bt));
2005:               PetscCall(PetscObjectTypeCompare((PetscObject)Bt, MATSEQAIJ, &fast));
2006:             } else {
2007:               PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &istrans));
2008:               if (istrans) {
2009:                 Mat Bt;

2011:                 PetscCall(MatHermitianTransposeGetMat(B, &Bt));
2012:                 PetscCall(PetscObjectTypeCompare((PetscObject)Bt, MATSEQAIJ, &fast));
2013:               }
2014:             }
2015:           }
2016:         }
2017:       }
2018:     }
2019:     for (i = 0, nf = 0; i < nest->nr && fast; ++i) {
2020:       PetscCall(PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i], ISSTRIDE, &fast));
2021:       if (fast) {
2022:         PetscInt f, s;

2024:         PetscCall(ISStrideGetInfo(nest->isglobal.row[i], &f, &s));
2025:         if (f != nf || s != 1) {
2026:           fast = PETSC_FALSE;
2027:         } else {
2028:           PetscCall(ISGetSize(nest->isglobal.row[i], &f));
2029:           nf += f;
2030:         }
2031:       }
2032:     }
2033:     for (i = 0, nf = 0; i < nest->nc && fast; ++i) {
2034:       PetscCall(PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i], ISSTRIDE, &fast));
2035:       if (fast) {
2036:         PetscInt f, s;

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

2165: static PetscErrorCode MatConvert_Nest_Dense(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
2166: {
2167:   Mat      B;
2168:   PetscInt m, n, M, N;

2170:   PetscFunctionBegin;
2171:   PetscCall(MatGetSize(A, &M, &N));
2172:   PetscCall(MatGetLocalSize(A, &m, &n));
2173:   if (reuse == MAT_REUSE_MATRIX) {
2174:     B = *newmat;
2175:     PetscCall(MatZeroEntries(B));
2176:   } else {
2177:     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A), m, PETSC_DECIDE, M, N, NULL, &B));
2178:   }
2179:   PetscCall(MatAXPY_Dense_Nest(B, 1.0, A));
2180:   if (reuse == MAT_INPLACE_MATRIX) {
2181:     PetscCall(MatHeaderReplace(A, &B));
2182:   } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
2183:   PetscFunctionReturn(PETSC_SUCCESS);
2184: }

2186: static PetscErrorCode MatHasOperation_Nest(Mat mat, MatOperation op, PetscBool *has)
2187: {
2188:   Mat_Nest    *bA = (Mat_Nest *)mat->data;
2189:   MatOperation opAdd;
2190:   PetscInt     i, j, nr = bA->nr, nc = bA->nc;
2191:   PetscBool    flg;
2192:   PetscFunctionBegin;

2194:   *has = PETSC_FALSE;
2195:   if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
2196:     opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD);
2197:     for (j = 0; j < nc; j++) {
2198:       for (i = 0; i < nr; i++) {
2199:         if (!bA->m[i][j]) continue;
2200:         PetscCall(MatHasOperation(bA->m[i][j], opAdd, &flg));
2201:         if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
2202:       }
2203:     }
2204:   }
2205:   if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
2206:   PetscFunctionReturn(PETSC_SUCCESS);
2207: }

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

2212:   Level: intermediate

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

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

2223: .seealso: [](ch_matrices), `Mat`, `MATNEST`, `MatCreate()`, `MatType`, `MatCreateNest()`, `MatNestSetSubMat()`, `MatNestGetSubMat()`,
2224:           `VecCreateNest()`, `DMCreateMatrix()`, `DMCOMPOSITE`, `MatNestSetVecType()`, `MatNestGetLocalISs()`,
2225:           `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()`
2226: M*/
2227: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2228: {
2229:   Mat_Nest *s;

2231:   PetscFunctionBegin;
2232:   PetscCall(PetscNew(&s));
2233:   A->data = (void *)s;

2235:   s->nr            = -1;
2236:   s->nc            = -1;
2237:   s->m             = NULL;
2238:   s->splitassembly = PETSC_FALSE;

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

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

2268:   A->spptr     = NULL;
2269:   A->assembled = PETSC_FALSE;

2271:   /* expose Nest api's */
2272:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMat_C", MatNestGetSubMat_Nest));
2273:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMat_C", MatNestSetSubMat_Nest));
2274:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMats_C", MatNestGetSubMats_Nest));
2275:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetSize_C", MatNestGetSize_Nest));
2276:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetISs_C", MatNestGetISs_Nest));
2277:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestGetLocalISs_C", MatNestGetLocalISs_Nest));
2278:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetVecType_C", MatNestSetVecType_Nest));
2279:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMats_C", MatNestSetSubMats_Nest));
2280:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ));
2281:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ));
2282:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_aij_C", MatConvert_Nest_AIJ));
2283:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_is_C", MatConvert_Nest_IS));
2284:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpidense_C", MatConvert_Nest_Dense));
2285:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqdense_C", MatConvert_Nest_Dense));
2286:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_seqdense_C", MatProductSetFromOptions_Nest_Dense));
2287:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_mpidense_C", MatProductSetFromOptions_Nest_Dense));
2288:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_dense_C", MatProductSetFromOptions_Nest_Dense));

2290:   PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATNEST));
2291:   PetscFunctionReturn(PETSC_SUCCESS);
2292: }