Actual source code: matnest.c

  1: #include <../src/mat/impls/nest/matnestimpl.h>
  2: #include <../src/mat/impls/aij/seq/aij.h>
  3: #include <petscsf.h>

  5: static PetscErrorCode MatSetUp_NestIS_Private(Mat, PetscInt, const IS[], PetscInt, const IS[]);
  6: static PetscErrorCode MatCreateVecs_Nest(Mat, Vec *, Vec *);
  7: static PetscErrorCode MatReset_Nest(Mat);

  9: PETSC_INTERN PetscErrorCode MatConvert_Nest_IS(Mat, MatType, MatReuse, Mat *);

 11: /* private functions */
 12: static PetscErrorCode MatNestGetSizes_Private(Mat A, PetscInt *m, PetscInt *n, PetscInt *M, PetscInt *N)
 13: {
 14:   Mat_Nest *bA = (Mat_Nest *)A->data;
 15:   PetscInt  i, j;

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

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

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

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

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

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

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

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

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

136:       /* C[i] <- workC + C[i] */
137:       MatAXPY(viewC, 1.0, contents->workC[i * nc + j], SAME_NONZERO_PATTERN);
138:     }
139:     MatDestroy(&viewC);
140:   }
141:   MatDenseRestoreArray(C, &carray);
142:   MatDenseRestoreArrayRead(B, &barray);

144:   MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
145:   MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
146:   return 0;
147: }

149: PetscErrorCode MatNest_DenseDestroy(void *ctx)
150: {
151:   Nest_Dense *contents = (Nest_Dense *)ctx;
152:   PetscInt    i;

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

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

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

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

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

232:       /* since tarray will be shared by all Mat */
233:       MatSeqDenseSetPreallocation(workC, contents->tarray);
234:       MatMPIDenseSetPreallocation(workC, contents->tarray);
235:     }
236:     MatDestroy(&viewB);
237:   }
238:   MatDenseRestoreArrayRead(B, &barray);

240:   C->ops->productnumeric = MatProductNumeric_Nest_Dense;
241:   return 0;
242: }

244: /* --------------------------------------------------------- */
245: static PetscErrorCode MatProductSetFromOptions_Nest_Dense_AB(Mat C)
246: {
247:   C->ops->productsymbolic = MatProductSymbolic_Nest_Dense;
248:   return 0;
249: }

251: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Nest_Dense(Mat C)
252: {
253:   Mat_Product *product = C->product;

255:   if (product->type == MATPRODUCT_AB) MatProductSetFromOptions_Nest_Dense_AB(C);
256:   return 0;
257: }
258: /* --------------------------------------------------------- */

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

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

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

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

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

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

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

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

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

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

347:   if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) {
348:     *B = C;
349:   } else {
350:     MatHeaderMerge(A, &C);
351:   }
352:   return 0;
353: }

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

360:   if (!lst) return 0;
361:   for (i = 0; i < n; i++)
362:     if (lst[i]) ISDestroy(&lst[i]);
363:   PetscFree(lst);
364:   *list = NULL;
365:   return 0;
366: }

368: static PetscErrorCode MatReset_Nest(Mat A)
369: {
370:   Mat_Nest *vs = (Mat_Nest *)A->data;
371:   PetscInt  i, j;

373:   /* release the matrices and the place holders */
374:   MatNestDestroyISList(vs->nr, &vs->isglobal.row);
375:   MatNestDestroyISList(vs->nc, &vs->isglobal.col);
376:   MatNestDestroyISList(vs->nr, &vs->islocal.row);
377:   MatNestDestroyISList(vs->nc, &vs->islocal.col);

379:   PetscFree(vs->row_len);
380:   PetscFree(vs->col_len);
381:   PetscFree(vs->nnzstate);

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

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

394:   /* restore defaults */
395:   vs->nr            = 0;
396:   vs->nc            = 0;
397:   vs->splitassembly = PETSC_FALSE;
398:   return 0;
399: }

401: static PetscErrorCode MatDestroy_Nest(Mat A)
402: {
403:   MatReset_Nest(A);
404:   PetscFree(A->data);
405:   PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMat_C", NULL);
406:   PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMat_C", NULL);
407:   PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMats_C", NULL);
408:   PetscObjectComposeFunction((PetscObject)A, "MatNestGetSize_C", NULL);
409:   PetscObjectComposeFunction((PetscObject)A, "MatNestGetISs_C", NULL);
410:   PetscObjectComposeFunction((PetscObject)A, "MatNestGetLocalISs_C", NULL);
411:   PetscObjectComposeFunction((PetscObject)A, "MatNestSetVecType_C", NULL);
412:   PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMats_C", NULL);
413:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpiaij_C", NULL);
414:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqaij_C", NULL);
415:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_aij_C", NULL);
416:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_is_C", NULL);
417:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpidense_C", NULL);
418:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqdense_C", NULL);
419:   PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_seqdense_C", NULL);
420:   PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_mpidense_C", NULL);
421:   PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_dense_C", NULL);
422:   return 0;
423: }

425: static PetscErrorCode MatMissingDiagonal_Nest(Mat mat, PetscBool *missing, PetscInt *dd)
426: {
427:   Mat_Nest *vs = (Mat_Nest *)mat->data;
428:   PetscInt  i;

430:   if (dd) *dd = 0;
431:   if (!vs->nr) {
432:     *missing = PETSC_TRUE;
433:     return 0;
434:   }
435:   *missing = PETSC_FALSE;
436:   for (i = 0; i < vs->nr && !(*missing); i++) {
437:     *missing = PETSC_TRUE;
438:     if (vs->m[i][i]) {
439:       MatMissingDiagonal(vs->m[i][i], missing, NULL);
441:     }
442:   }
443:   return 0;
444: }

446: static PetscErrorCode MatAssemblyBegin_Nest(Mat A, MatAssemblyType type)
447: {
448:   Mat_Nest *vs = (Mat_Nest *)A->data;
449:   PetscInt  i, j;
450:   PetscBool nnzstate = PETSC_FALSE;

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

476: static PetscErrorCode MatAssemblyEnd_Nest(Mat A, MatAssemblyType type)
477: {
478:   Mat_Nest *vs = (Mat_Nest *)A->data;
479:   PetscInt  i, j;

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

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

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

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

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

517: static PetscErrorCode MatNestFindISRange(Mat A, PetscInt n, const IS list[], IS is, PetscInt *begin, PetscInt *end)
518: {
519:   PetscInt  i, j, size, m;
520:   PetscBool flg;
521:   IS        out, concatenate[2];

525:   if (begin) {
527:     *begin = -1;
528:   }
529:   if (end) {
531:     *end = -1;
532:   }
533:   for (i = 0; i < n; i++) {
534:     if (!list[i]) continue;
535:     ISEqualUnsorted(list[i], is, &flg);
536:     if (flg) {
537:       if (begin) *begin = i;
538:       if (end) *end = i + 1;
539:       return 0;
540:     }
541:   }
542:   ISGetSize(is, &size);
543:   for (i = 0; i < n - 1; i++) {
544:     if (!list[i]) continue;
545:     m = 0;
546:     ISConcatenate(PetscObjectComm((PetscObject)A), 2, list + i, &out);
547:     ISGetSize(out, &m);
548:     for (j = i + 2; j < n && m < size; j++) {
549:       if (list[j]) {
550:         concatenate[0] = out;
551:         concatenate[1] = list[j];
552:         ISConcatenate(PetscObjectComm((PetscObject)A), 2, concatenate, &out);
553:         ISDestroy(concatenate);
554:         ISGetSize(out, &m);
555:       }
556:     }
557:     if (m == size) {
558:       ISEqualUnsorted(out, is, &flg);
559:       if (flg) {
560:         if (begin) *begin = i;
561:         if (end) *end = j;
562:         ISDestroy(&out);
563:         return 0;
564:       }
565:     }
566:     ISDestroy(&out);
567:   }
568:   return 0;
569: }

571: static PetscErrorCode MatNestFillEmptyMat_Private(Mat A, PetscInt i, PetscInt j, Mat *B)
572: {
573:   Mat_Nest *vs = (Mat_Nest *)A->data;
574:   PetscInt  lr, lc;

576:   MatCreate(PetscObjectComm((PetscObject)A), B);
577:   ISGetLocalSize(vs->isglobal.row[i], &lr);
578:   ISGetLocalSize(vs->isglobal.col[j], &lc);
579:   MatSetSizes(*B, lr, lc, PETSC_DECIDE, PETSC_DECIDE);
580:   MatSetType(*B, MATAIJ);
581:   MatSeqAIJSetPreallocation(*B, 0, NULL);
582:   MatMPIAIJSetPreallocation(*B, 0, NULL, 0, NULL);
583:   MatSetUp(*B);
584:   MatSetOption(*B, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
585:   MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
586:   MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY);
587:   return 0;
588: }

590: static PetscErrorCode MatNestGetBlock_Private(Mat A, PetscInt rbegin, PetscInt rend, PetscInt cbegin, PetscInt cend, Mat *B)
591: {
592:   Mat_Nest  *vs = (Mat_Nest *)A->data;
593:   Mat       *a;
594:   PetscInt   i, j, k, l, nr = rend - rbegin, nc = cend - cbegin;
595:   char       keyname[256];
596:   PetscBool *b;
597:   PetscBool  flg;

599:   *B = NULL;
600:   PetscSNPrintf(keyname, sizeof(keyname), "NestBlock_%" PetscInt_FMT "-%" PetscInt_FMT "x%" PetscInt_FMT "-%" PetscInt_FMT, rbegin, rend, cbegin, cend);
601:   PetscObjectQuery((PetscObject)A, keyname, (PetscObject *)B);
602:   if (*B) return 0;

604:   PetscMalloc2(nr * nc, &a, nr * nc, &b);
605:   for (i = 0; i < nr; i++) {
606:     for (j = 0; j < nc; j++) {
607:       a[i * nc + j] = vs->m[rbegin + i][cbegin + j];
608:       b[i * nc + j] = PETSC_FALSE;
609:     }
610:   }
611:   if (nc != vs->nc && nr != vs->nr) {
612:     for (i = 0; i < nr; i++) {
613:       for (j = 0; j < nc; j++) {
614:         flg = PETSC_FALSE;
615:         for (k = 0; (k < nr && !flg); k++) {
616:           if (a[j + k * nc]) flg = PETSC_TRUE;
617:         }
618:         if (flg) {
619:           flg = PETSC_FALSE;
620:           for (l = 0; (l < nc && !flg); l++) {
621:             if (a[i * nc + l]) flg = PETSC_TRUE;
622:           }
623:         }
624:         if (!flg) {
625:           b[i * nc + j] = PETSC_TRUE;
626:           MatNestFillEmptyMat_Private(A, rbegin + i, cbegin + j, a + i * nc + j);
627:         }
628:       }
629:     }
630:   }
631:   MatCreateNest(PetscObjectComm((PetscObject)A), nr, nr != vs->nr ? NULL : vs->isglobal.row, nc, nc != vs->nc ? NULL : vs->isglobal.col, a, B);
632:   for (i = 0; i < nr; i++) {
633:     for (j = 0; j < nc; j++) {
634:       if (b[i * nc + j]) MatDestroy(a + i * nc + j);
635:     }
636:   }
637:   PetscFree2(a, b);
638:   (*B)->assembled = A->assembled;
639:   PetscObjectCompose((PetscObject)A, keyname, (PetscObject)*B);
640:   PetscObjectDereference((PetscObject)*B); /* Leave the only remaining reference in the composition */
641:   return 0;
642: }

644: static PetscErrorCode MatNestFindSubMat(Mat A, struct MatNestISPair *is, IS isrow, IS iscol, Mat *B)
645: {
646:   Mat_Nest *vs = (Mat_Nest *)A->data;
647:   PetscInt  rbegin, rend, cbegin, cend;

649:   MatNestFindISRange(A, vs->nr, is->row, isrow, &rbegin, &rend);
650:   MatNestFindISRange(A, vs->nc, is->col, iscol, &cbegin, &cend);
651:   if (rend == rbegin + 1 && cend == cbegin + 1) {
652:     if (!vs->m[rbegin][cbegin]) MatNestFillEmptyMat_Private(A, rbegin, cbegin, vs->m[rbegin] + cbegin);
653:     *B = vs->m[rbegin][cbegin];
654:   } else if (rbegin != -1 && cbegin != -1) {
655:     MatNestGetBlock_Private(A, rbegin, rend, cbegin, cend, B);
656:   } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Could not find index set");
657:   return 0;
658: }

660: /*
661:    TODO: This does not actually returns a submatrix we can modify
662: */
663: static PetscErrorCode MatCreateSubMatrix_Nest(Mat A, IS isrow, IS iscol, MatReuse reuse, Mat *B)
664: {
665:   Mat_Nest *vs = (Mat_Nest *)A->data;
666:   Mat       sub;

668:   MatNestFindSubMat(A, &vs->isglobal, isrow, iscol, &sub);
669:   switch (reuse) {
670:   case MAT_INITIAL_MATRIX:
671:     if (sub) PetscObjectReference((PetscObject)sub);
672:     *B = sub;
673:     break;
674:   case MAT_REUSE_MATRIX:
676:     break;
677:   case MAT_IGNORE_MATRIX: /* Nothing to do */
678:     break;
679:   case MAT_INPLACE_MATRIX: /* Nothing to do */
680:     SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MAT_INPLACE_MATRIX is not supported yet");
681:   }
682:   return 0;
683: }

685: PetscErrorCode MatGetLocalSubMatrix_Nest(Mat A, IS isrow, IS iscol, Mat *B)
686: {
687:   Mat_Nest *vs = (Mat_Nest *)A->data;
688:   Mat       sub;

690:   MatNestFindSubMat(A, &vs->islocal, isrow, iscol, &sub);
691:   /* We allow the submatrix to be NULL, perhaps it would be better for the user to return an empty matrix instead */
692:   if (sub) PetscObjectReference((PetscObject)sub);
693:   *B = sub;
694:   return 0;
695: }

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

702:   MatNestFindSubMat(A, &vs->islocal, isrow, iscol, &sub);
704:   if (sub) {
706:     MatDestroy(B);
707:   }
708:   return 0;
709: }

711: static PetscErrorCode MatGetDiagonal_Nest(Mat A, Vec v)
712: {
713:   Mat_Nest *bA = (Mat_Nest *)A->data;
714:   PetscInt  i;

716:   for (i = 0; i < bA->nr; i++) {
717:     Vec bv;
718:     VecGetSubVector(v, bA->isglobal.row[i], &bv);
719:     if (bA->m[i][i]) {
720:       MatGetDiagonal(bA->m[i][i], bv);
721:     } else {
722:       VecSet(bv, 0.0);
723:     }
724:     VecRestoreSubVector(v, bA->isglobal.row[i], &bv);
725:   }
726:   return 0;
727: }

729: static PetscErrorCode MatDiagonalScale_Nest(Mat A, Vec l, Vec r)
730: {
731:   Mat_Nest *bA = (Mat_Nest *)A->data;
732:   Vec       bl, *br;
733:   PetscInt  i, j;

735:   PetscCalloc1(bA->nc, &br);
736:   if (r) {
737:     for (j = 0; j < bA->nc; j++) VecGetSubVector(r, bA->isglobal.col[j], &br[j]);
738:   }
739:   bl = NULL;
740:   for (i = 0; i < bA->nr; i++) {
741:     if (l) VecGetSubVector(l, bA->isglobal.row[i], &bl);
742:     for (j = 0; j < bA->nc; j++) {
743:       if (bA->m[i][j]) MatDiagonalScale(bA->m[i][j], bl, br[j]);
744:     }
745:     if (l) VecRestoreSubVector(l, bA->isglobal.row[i], &bl);
746:   }
747:   if (r) {
748:     for (j = 0; j < bA->nc; j++) VecRestoreSubVector(r, bA->isglobal.col[j], &br[j]);
749:   }
750:   PetscFree(br);
751:   return 0;
752: }

754: static PetscErrorCode MatScale_Nest(Mat A, PetscScalar a)
755: {
756:   Mat_Nest *bA = (Mat_Nest *)A->data;
757:   PetscInt  i, j;

759:   for (i = 0; i < bA->nr; i++) {
760:     for (j = 0; j < bA->nc; j++) {
761:       if (bA->m[i][j]) MatScale(bA->m[i][j], a);
762:     }
763:   }
764:   return 0;
765: }

767: static PetscErrorCode MatShift_Nest(Mat A, PetscScalar a)
768: {
769:   Mat_Nest *bA = (Mat_Nest *)A->data;
770:   PetscInt  i;
771:   PetscBool nnzstate = PETSC_FALSE;

773:   for (i = 0; i < bA->nr; i++) {
774:     PetscObjectState subnnzstate = 0;
776:     MatShift(bA->m[i][i], a);
777:     MatGetNonzeroState(bA->m[i][i], &subnnzstate);
778:     nnzstate                     = (PetscBool)(nnzstate || bA->nnzstate[i * bA->nc + i] != subnnzstate);
779:     bA->nnzstate[i * bA->nc + i] = subnnzstate;
780:   }
781:   if (nnzstate) A->nonzerostate++;
782:   return 0;
783: }

785: static PetscErrorCode MatDiagonalSet_Nest(Mat A, Vec D, InsertMode is)
786: {
787:   Mat_Nest *bA = (Mat_Nest *)A->data;
788:   PetscInt  i;
789:   PetscBool nnzstate = PETSC_FALSE;

791:   for (i = 0; i < bA->nr; i++) {
792:     PetscObjectState subnnzstate = 0;
793:     Vec              bv;
794:     VecGetSubVector(D, bA->isglobal.row[i], &bv);
795:     if (bA->m[i][i]) {
796:       MatDiagonalSet(bA->m[i][i], bv, is);
797:       MatGetNonzeroState(bA->m[i][i], &subnnzstate);
798:     }
799:     VecRestoreSubVector(D, bA->isglobal.row[i], &bv);
800:     nnzstate                     = (PetscBool)(nnzstate || bA->nnzstate[i * bA->nc + i] != subnnzstate);
801:     bA->nnzstate[i * bA->nc + i] = subnnzstate;
802:   }
803:   if (nnzstate) A->nonzerostate++;
804:   return 0;
805: }

807: static PetscErrorCode MatSetRandom_Nest(Mat A, PetscRandom rctx)
808: {
809:   Mat_Nest *bA = (Mat_Nest *)A->data;
810:   PetscInt  i, j;

812:   for (i = 0; i < bA->nr; i++) {
813:     for (j = 0; j < bA->nc; j++) {
814:       if (bA->m[i][j]) MatSetRandom(bA->m[i][j], rctx);
815:     }
816:   }
817:   return 0;
818: }

820: static PetscErrorCode MatCreateVecs_Nest(Mat A, Vec *right, Vec *left)
821: {
822:   Mat_Nest *bA = (Mat_Nest *)A->data;
823:   Vec      *L, *R;
824:   MPI_Comm  comm;
825:   PetscInt  i, j;

827:   PetscObjectGetComm((PetscObject)A, &comm);
828:   if (right) {
829:     /* allocate R */
830:     PetscMalloc1(bA->nc, &R);
831:     /* Create the right vectors */
832:     for (j = 0; j < bA->nc; j++) {
833:       for (i = 0; i < bA->nr; i++) {
834:         if (bA->m[i][j]) {
835:           MatCreateVecs(bA->m[i][j], &R[j], NULL);
836:           break;
837:         }
838:       }
840:     }
841:     VecCreateNest(comm, bA->nc, bA->isglobal.col, R, right);
842:     /* hand back control to the nest vector */
843:     for (j = 0; j < bA->nc; j++) VecDestroy(&R[j]);
844:     PetscFree(R);
845:   }

847:   if (left) {
848:     /* allocate L */
849:     PetscMalloc1(bA->nr, &L);
850:     /* Create the left vectors */
851:     for (i = 0; i < bA->nr; i++) {
852:       for (j = 0; j < bA->nc; j++) {
853:         if (bA->m[i][j]) {
854:           MatCreateVecs(bA->m[i][j], NULL, &L[i]);
855:           break;
856:         }
857:       }
859:     }

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

864:     PetscFree(L);
865:   }
866:   return 0;
867: }

869: static PetscErrorCode MatView_Nest(Mat A, PetscViewer viewer)
870: {
871:   Mat_Nest *bA = (Mat_Nest *)A->data;
872:   PetscBool isascii, viewSub = PETSC_FALSE;
873:   PetscInt  i, j;

875:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
876:   if (isascii) {
877:     PetscOptionsGetBool(((PetscObject)A)->options, ((PetscObject)A)->prefix, "-mat_view_nest_sub", &viewSub, NULL);
878:     PetscViewerASCIIPrintf(viewer, "Matrix object: \n");
879:     PetscViewerASCIIPushTab(viewer);
880:     PetscViewerASCIIPrintf(viewer, "type=nest, rows=%" PetscInt_FMT ", cols=%" PetscInt_FMT " \n", bA->nr, bA->nc);

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

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

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

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

914: static PetscErrorCode MatZeroEntries_Nest(Mat A)
915: {
916:   Mat_Nest *bA = (Mat_Nest *)A->data;
917:   PetscInt  i, j;

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

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

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

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

957:   for (i = 0; i < nr; i++) {
958:     for (j = 0; j < nc; j++) {
959:       PetscObjectState subnnzstate = 0;
960:       if (bY->m[i][j] && bX->m[i][j]) {
961:         MatAXPY(bY->m[i][j], a, bX->m[i][j], str);
962:       } else if (bX->m[i][j]) {
963:         Mat M;

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

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

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

1000:   MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
1001:   MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY);
1002:   return 0;
1003: }

1005: /* nest api */
1006: PetscErrorCode MatNestGetSubMat_Nest(Mat A, PetscInt idxm, PetscInt jdxm, Mat *mat)
1007: {
1008:   Mat_Nest *bA = (Mat_Nest *)A->data;

1012:   *mat = bA->m[idxm][jdxm];
1013:   return 0;
1014: }

1016: /*@
1017:  MatNestGetSubMat - Returns a single, sub-matrix from a `MATNEST`

1019:  Not collective

1021:  Input Parameters:
1022: +   A  - `MATNEST` matrix
1023: .   idxm - index of the matrix within the nest matrix
1024: -   jdxm - index of the matrix within the nest matrix

1026:  Output Parameter:
1027: .   sub - matrix at index idxm,jdxm within the nest matrix

1029:  Level: developer

1031: .seealso: `MATNEST`, `MatNestGetSize()`, `MatNestGetSubMats()`, `MatCreateNest()`, `MATNEST`, `MatNestSetSubMat()`,
1032:           `MatNestGetLocalISs()`, `MatNestGetISs()`
1033: @*/
1034: PetscErrorCode MatNestGetSubMat(Mat A, PetscInt idxm, PetscInt jdxm, Mat *sub)
1035: {
1036:   PetscUseMethod(A, "MatNestGetSubMat_C", (Mat, PetscInt, PetscInt, Mat *), (A, idxm, jdxm, sub));
1037:   return 0;
1038: }

1040: PetscErrorCode MatNestSetSubMat_Nest(Mat A, PetscInt idxm, PetscInt jdxm, Mat mat)
1041: {
1042:   Mat_Nest *bA = (Mat_Nest *)A->data;
1043:   PetscInt  m, n, M, N, mi, ni, Mi, Ni;

1047:   MatGetLocalSize(mat, &m, &n);
1048:   MatGetSize(mat, &M, &N);
1049:   ISGetLocalSize(bA->isglobal.row[idxm], &mi);
1050:   ISGetSize(bA->isglobal.row[idxm], &Mi);
1051:   ISGetLocalSize(bA->isglobal.col[jdxm], &ni);
1052:   ISGetSize(bA->isglobal.col[jdxm], &Ni);

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

1059:   PetscObjectReference((PetscObject)mat);
1060:   MatDestroy(&bA->m[idxm][jdxm]);
1061:   bA->m[idxm][jdxm] = mat;
1062:   PetscObjectStateIncrease((PetscObject)A);
1063:   MatGetNonzeroState(mat, &bA->nnzstate[idxm * bA->nc + jdxm]);
1064:   A->nonzerostate++;
1065:   return 0;
1066: }

1068: /*@
1069:  MatNestSetSubMat - Set a single submatrix in the `MATNEST`

1071:  Logically collective on the submatrix communicator

1073:  Input Parameters:
1074: +   A  - `MATNEST` matrix
1075: .   idxm - index of the matrix within the nest matrix
1076: .   jdxm - index of the matrix within the nest matrix
1077: -   sub - matrix at index idxm,jdxm within the nest matrix

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

1082:  This increments the reference count of the submatrix.

1084:  Level: developer

1086: .seealso: `MATNEST`, `MatNestSetSubMats()`, `MatNestGetSubMats()`, `MatNestGetLocalISs()`, `MATNEST`, `MatCreateNest()`,
1087:           `MatNestGetSubMat()`, `MatNestGetISs()`, `MatNestGetSize()`
1088: @*/
1089: PetscErrorCode MatNestSetSubMat(Mat A, PetscInt idxm, PetscInt jdxm, Mat sub)
1090: {
1091:   PetscUseMethod(A, "MatNestSetSubMat_C", (Mat, PetscInt, PetscInt, Mat), (A, idxm, jdxm, sub));
1092:   return 0;
1093: }

1095: PetscErrorCode MatNestGetSubMats_Nest(Mat A, PetscInt *M, PetscInt *N, Mat ***mat)
1096: {
1097:   Mat_Nest *bA = (Mat_Nest *)A->data;

1099:   if (M) *M = bA->nr;
1100:   if (N) *N = bA->nc;
1101:   if (mat) *mat = bA->m;
1102:   return 0;
1103: }

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

1108:  Not collective

1110:  Input Parameter:
1111: .   A  - nest matrix

1113:  Output Parameters:
1114: +   M - number of rows in the nest matrix
1115: .   N - number of cols in the nest matrix
1116: -   mat - 2d array of matrices

1118:  Note:
1119:  The user should not free the array mat.

1121:  Fortran Note:
1122:  This routine has a calling sequence
1123: $   call MatNestGetSubMats(A, M, N, mat, ierr)
1124:  where the space allocated for the optional argument mat is assumed large enough (if provided).

1126:  Level: developer

1128: .seealso: `MATNEST`, `MatNestGetSize()`, `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MATNEST`, `MatCreateNest()`,
1129:           `MatNestSetSubMats()`, `MatNestGetISs()`, `MatNestSetSubMat()`
1130: @*/
1131: PetscErrorCode MatNestGetSubMats(Mat A, PetscInt *M, PetscInt *N, Mat ***mat)
1132: {
1133:   PetscUseMethod(A, "MatNestGetSubMats_C", (Mat, PetscInt *, PetscInt *, Mat ***), (A, M, N, mat));
1134:   return 0;
1135: }

1137: PetscErrorCode MatNestGetSize_Nest(Mat A, PetscInt *M, PetscInt *N)
1138: {
1139:   Mat_Nest *bA = (Mat_Nest *)A->data;

1141:   if (M) *M = bA->nr;
1142:   if (N) *N = bA->nc;
1143:   return 0;
1144: }

1146: /*@
1147:  MatNestGetSize - Returns the size of the `MATNEST` matrix.

1149:  Not collective

1151:  Input Parameter:
1152: .   A  - `MATNEST` matrix

1154:  Output Parameters:
1155: +   M - number of rows in the nested mat
1156: -   N - number of cols in the nested mat

1158:  Level: developer

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

1169: static PetscErrorCode MatNestGetISs_Nest(Mat A, IS rows[], IS cols[])
1170: {
1171:   Mat_Nest *vs = (Mat_Nest *)A->data;
1172:   PetscInt  i;

1174:   if (rows)
1175:     for (i = 0; i < vs->nr; i++) rows[i] = vs->isglobal.row[i];
1176:   if (cols)
1177:     for (i = 0; i < vs->nc; i++) cols[i] = vs->isglobal.col[i];
1178:   return 0;
1179: }

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

1184:  Not collective

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

1189:  Output Parameters:
1190: +   rows - array of row index sets
1191: -   cols - array of column index sets

1193:  Level: advanced

1195:  Note:
1196:  The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.

1198: .seealso: `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetLocalISs()`, `MATNEST`,
1199:           `MatCreateNest()`, `MatNestGetSubMats()`, `MatNestSetSubMats()`
1200: @*/
1201: PetscErrorCode MatNestGetISs(Mat A, IS rows[], IS cols[])
1202: {
1204:   PetscUseMethod(A, "MatNestGetISs_C", (Mat, IS[], IS[]), (A, rows, cols));
1205:   return 0;
1206: }

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

1213:   if (rows)
1214:     for (i = 0; i < vs->nr; i++) rows[i] = vs->islocal.row[i];
1215:   if (cols)
1216:     for (i = 0; i < vs->nc; i++) cols[i] = vs->islocal.col[i];
1217:   return 0;
1218: }

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

1223:  Not collective

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

1228:  Output Parameters:
1229: +   rows - array of row index sets (or NULL to ignore)
1230: -   cols - array of column index sets (or NULL to ignore)

1232:  Level: advanced

1234:  Note:
1235:  The user must have allocated arrays of the correct size. The reference count is not increased on the returned ISs.

1237: .seealso: `MATNEST`, `MatNestGetSubMat()`, `MatNestGetSubMats()`, `MatNestGetSize()`, `MatNestGetISs()`, `MatCreateNest()`,
1238:           `MATNEST`, `MatNestSetSubMats()`, `MatNestSetSubMat()`
1239: @*/
1240: PetscErrorCode MatNestGetLocalISs(Mat A, IS rows[], IS cols[])
1241: {
1243:   PetscUseMethod(A, "MatNestGetLocalISs_C", (Mat, IS[], IS[]), (A, rows, cols));
1244:   return 0;
1245: }

1247: PetscErrorCode MatNestSetVecType_Nest(Mat A, VecType vtype)
1248: {
1249:   PetscBool flg;

1251:   PetscStrcmp(vtype, VECNEST, &flg);
1252:   /* In reality, this only distinguishes VECNEST and "other" */
1253:   if (flg) A->ops->getvecs = MatCreateVecs_Nest;
1254:   else A->ops->getvecs = (PetscErrorCode(*)(Mat, Vec *, Vec *))0;
1255:   return 0;
1256: }

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

1261:  Not collective

1263:  Input Parameters:
1264: +  A  - `MATNEST` matrix
1265: -  vtype - `VecType` to use for creating vectors

1267:  Level: developer

1269: .seealso: `MATNEST`, `MatCreateVecs()`, `MATNEST`, `MatCreateNest()`, `VecType`
1270: @*/
1271: PetscErrorCode MatNestSetVecType(Mat A, VecType vtype)
1272: {
1273:   PetscTryMethod(A, "MatNestSetVecType_C", (Mat, VecType), (A, vtype));
1274:   return 0;
1275: }

1277: PetscErrorCode MatNestSetSubMats_Nest(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[])
1278: {
1279:   Mat_Nest *s = (Mat_Nest *)A->data;
1280:   PetscInt  i, j, m, n, M, N;
1281:   PetscBool cong, isstd, sametype = PETSC_FALSE;
1282:   VecType   vtype, type;

1284:   MatReset_Nest(A);

1286:   s->nr = nr;
1287:   s->nc = nc;

1289:   /* Create space for submatrices */
1290:   PetscMalloc1(nr, &s->m);
1291:   for (i = 0; i < nr; i++) PetscMalloc1(nc, &s->m[i]);
1292:   for (i = 0; i < nr; i++) {
1293:     for (j = 0; j < nc; j++) {
1294:       s->m[i][j] = a[i * nc + j];
1295:       if (a[i * nc + j]) PetscObjectReference((PetscObject)a[i * nc + j]);
1296:     }
1297:   }
1298:   MatGetVecType(A, &vtype);
1299:   PetscStrcmp(vtype, VECSTANDARD, &isstd);
1300:   if (isstd) {
1301:     /* check if all blocks have the same vectype */
1302:     vtype = NULL;
1303:     for (i = 0; i < nr; i++) {
1304:       for (j = 0; j < nc; j++) {
1305:         if (a[i * nc + j]) {
1306:           if (!vtype) { /* first visited block */
1307:             MatGetVecType(a[i * nc + j], &vtype);
1308:             sametype = PETSC_TRUE;
1309:           } else if (sametype) {
1310:             MatGetVecType(a[i * nc + j], &type);
1311:             PetscStrcmp(vtype, type, &sametype);
1312:           }
1313:         }
1314:       }
1315:     }
1316:     if (sametype) { /* propagate vectype */
1317:       MatSetVecType(A, vtype);
1318:     }
1319:   }

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

1323:   PetscMalloc1(nr, &s->row_len);
1324:   PetscMalloc1(nc, &s->col_len);
1325:   for (i = 0; i < nr; i++) s->row_len[i] = -1;
1326:   for (j = 0; j < nc; j++) s->col_len[j] = -1;

1328:   PetscCalloc1(nr * nc, &s->nnzstate);
1329:   for (i = 0; i < nr; i++) {
1330:     for (j = 0; j < nc; j++) {
1331:       if (s->m[i][j]) MatGetNonzeroState(s->m[i][j], &s->nnzstate[i * nc + j]);
1332:     }
1333:   }

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

1337:   PetscLayoutSetSize(A->rmap, M);
1338:   PetscLayoutSetLocalSize(A->rmap, m);
1339:   PetscLayoutSetSize(A->cmap, N);
1340:   PetscLayoutSetLocalSize(A->cmap, n);

1342:   PetscLayoutSetUp(A->rmap);
1343:   PetscLayoutSetUp(A->cmap);

1345:   /* disable operations that are not supported for non-square matrices,
1346:      or matrices for which is_row != is_col  */
1347:   MatHasCongruentLayouts(A, &cong);
1348:   if (cong && nr != nc) cong = PETSC_FALSE;
1349:   if (cong) {
1350:     for (i = 0; cong && i < nr; i++) ISEqualUnsorted(s->isglobal.row[i], s->isglobal.col[i], &cong);
1351:   }
1352:   if (!cong) {
1353:     A->ops->missingdiagonal = NULL;
1354:     A->ops->getdiagonal     = NULL;
1355:     A->ops->shift           = NULL;
1356:     A->ops->diagonalset     = NULL;
1357:   }

1359:   PetscCalloc2(nr, &s->left, nc, &s->right);
1360:   PetscObjectStateIncrease((PetscObject)A);
1361:   A->nonzerostate++;
1362:   return 0;
1363: }

1365: /*@
1366:    MatNestSetSubMats - Sets the nested submatrices in a `MATNEST`

1368:    Collective

1370:    Input Parameters:
1371: +  A - `MATNEST` matrix
1372: .  nr - number of nested row blocks
1373: .  is_row - index sets for each nested row block, or NULL to make contiguous
1374: .  nc - number of nested column blocks
1375: .  is_col - index sets for each nested column block, or NULL to make contiguous
1376: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

1378:    Note:
1379:    This always resets any submatrix information previously set

1381:    Level: advanced

1383: .seealso: `MATNEST`, `MatCreateNest()`, `MATNEST`, `MatNestSetSubMat()`, `MatNestGetSubMat()`, `MatNestGetSubMats()`
1384: @*/
1385: PetscErrorCode MatNestSetSubMats(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[])
1386: {
1387:   PetscInt i;

1391:   if (nr && is_row) {
1394:   }
1396:   if (nc && is_col) {
1399:   }
1401:   PetscUseMethod(A, "MatNestSetSubMats_C", (Mat, PetscInt, const IS[], PetscInt, const IS[], const Mat[]), (A, nr, is_row, nc, is_col, a));
1402:   return 0;
1403: }

1405: static PetscErrorCode MatNestCreateAggregateL2G_Private(Mat A, PetscInt n, const IS islocal[], const IS isglobal[], PetscBool colflg, ISLocalToGlobalMapping *ltog)
1406: {
1407:   PetscBool flg;
1408:   PetscInt  i, j, m, mi, *ix;

1410:   *ltog = NULL;
1411:   for (i = 0, m = 0, flg = PETSC_FALSE; i < n; i++) {
1412:     if (islocal[i]) {
1413:       ISGetLocalSize(islocal[i], &mi);
1414:       flg = PETSC_TRUE; /* We found a non-trivial entry */
1415:     } else {
1416:       ISGetLocalSize(isglobal[i], &mi);
1417:     }
1418:     m += mi;
1419:   }
1420:   if (!flg) return 0;

1422:   PetscMalloc1(m, &ix);
1423:   for (i = 0, m = 0; i < n; i++) {
1424:     ISLocalToGlobalMapping smap = NULL;
1425:     Mat                    sub  = NULL;
1426:     PetscSF                sf;
1427:     PetscLayout            map;
1428:     const PetscInt        *ix2;

1430:     if (!colflg) {
1431:       MatNestFindNonzeroSubMatRow(A, i, &sub);
1432:     } else {
1433:       MatNestFindNonzeroSubMatCol(A, i, &sub);
1434:     }
1435:     if (sub) {
1436:       if (!colflg) {
1437:         MatGetLocalToGlobalMapping(sub, &smap, NULL);
1438:       } else {
1439:         MatGetLocalToGlobalMapping(sub, NULL, &smap);
1440:       }
1441:     }
1442:     /*
1443:        Now we need to extract the monolithic global indices that correspond to the given split global indices.
1444:        In many/most cases, we only want MatGetLocalSubMatrix() to work, in which case we only need to know the size of the local spaces.
1445:     */
1446:     ISGetIndices(isglobal[i], &ix2);
1447:     if (islocal[i]) {
1448:       PetscInt *ilocal, *iremote;
1449:       PetscInt  mil, nleaves;

1451:       ISGetLocalSize(islocal[i], &mi);
1453:       for (j = 0; j < mi; j++) ix[m + j] = j;
1454:       ISLocalToGlobalMappingApply(smap, mi, ix + m, ix + m);

1456:       /* PetscSFSetGraphLayout does not like negative indices */
1457:       PetscMalloc2(mi, &ilocal, mi, &iremote);
1458:       for (j = 0, nleaves = 0; j < mi; j++) {
1459:         if (ix[m + j] < 0) continue;
1460:         ilocal[nleaves]  = j;
1461:         iremote[nleaves] = ix[m + j];
1462:         nleaves++;
1463:       }
1464:       ISGetLocalSize(isglobal[i], &mil);
1465:       PetscSFCreate(PetscObjectComm((PetscObject)A), &sf);
1466:       PetscLayoutCreate(PetscObjectComm((PetscObject)A), &map);
1467:       PetscLayoutSetLocalSize(map, mil);
1468:       PetscLayoutSetUp(map);
1469:       PetscSFSetGraphLayout(sf, map, nleaves, ilocal, PETSC_USE_POINTER, iremote);
1470:       PetscLayoutDestroy(&map);
1471:       PetscSFBcastBegin(sf, MPIU_INT, ix2, ix + m, MPI_REPLACE);
1472:       PetscSFBcastEnd(sf, MPIU_INT, ix2, ix + m, MPI_REPLACE);
1473:       PetscSFDestroy(&sf);
1474:       PetscFree2(ilocal, iremote);
1475:     } else {
1476:       ISGetLocalSize(isglobal[i], &mi);
1477:       for (j = 0; j < mi; j++) ix[m + j] = ix2[i];
1478:     }
1479:     ISRestoreIndices(isglobal[i], &ix2);
1480:     m += mi;
1481:   }
1482:   ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)A), 1, m, ix, PETSC_OWN_POINTER, ltog);
1483:   return 0;
1484: }

1486: /* If an IS was provided, there is nothing Nest needs to do, otherwise Nest will build a strided IS */
1487: /*
1488:   nprocessors = NP
1489:   Nest x^T = ((g_0,g_1,...g_nprocs-1), (h_0,h_1,...h_NP-1))
1490:        proc 0: => (g_0,h_0,)
1491:        proc 1: => (g_1,h_1,)
1492:        ...
1493:        proc nprocs-1: => (g_NP-1,h_NP-1,)

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

1498:             proc 0:
1499:     is[1] = (nlocal(g_0),nlocal(g_0)+1,...,nlocal(g_0)+nlocal(h_0)-1)
1500:             proc 1:
1501:     is[1] = (nlocal(g_1),nlocal(g_1)+1,...,nlocal(g_1)+nlocal(h_1)-1)

1503:             proc NP-1:
1504:     is[1] = (nlocal(g_NP-1),nlocal(g_NP-1)+1,...,nlocal(g_NP-1)+nlocal(h_NP-1)-1)
1505: */
1506: static PetscErrorCode MatSetUp_NestIS_Private(Mat A, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[])
1507: {
1508:   Mat_Nest *vs = (Mat_Nest *)A->data;
1509:   PetscInt  i, j, offset, n, nsum, bs;
1510:   Mat       sub = NULL;

1512:   PetscMalloc1(nr, &vs->isglobal.row);
1513:   PetscMalloc1(nc, &vs->isglobal.col);
1514:   if (is_row) { /* valid IS is passed in */
1515:     /* refs on is[] are incremented */
1516:     for (i = 0; i < vs->nr; i++) {
1517:       PetscObjectReference((PetscObject)is_row[i]);

1519:       vs->isglobal.row[i] = is_row[i];
1520:     }
1521:   } else { /* Create the ISs by inspecting sizes of a submatrix in each row */
1522:     nsum = 0;
1523:     for (i = 0; i < vs->nr; i++) { /* Add up the local sizes to compute the aggregate offset */
1524:       MatNestFindNonzeroSubMatRow(A, i, &sub);
1526:       MatGetLocalSize(sub, &n, NULL);
1528:       nsum += n;
1529:     }
1530:     MPI_Scan(&nsum, &offset, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A));
1531:     offset -= nsum;
1532:     for (i = 0; i < vs->nr; i++) {
1533:       MatNestFindNonzeroSubMatRow(A, i, &sub);
1534:       MatGetLocalSize(sub, &n, NULL);
1535:       MatGetBlockSizes(sub, &bs, NULL);
1536:       ISCreateStride(PetscObjectComm((PetscObject)sub), n, offset, 1, &vs->isglobal.row[i]);
1537:       ISSetBlockSize(vs->isglobal.row[i], bs);
1538:       offset += n;
1539:     }
1540:   }

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

1547:       vs->isglobal.col[j] = is_col[j];
1548:     }
1549:   } else { /* Create the ISs by inspecting sizes of a submatrix in each column */
1550:     offset = A->cmap->rstart;
1551:     nsum   = 0;
1552:     for (j = 0; j < vs->nc; j++) {
1553:       MatNestFindNonzeroSubMatCol(A, j, &sub);
1555:       MatGetLocalSize(sub, NULL, &n);
1557:       nsum += n;
1558:     }
1559:     MPI_Scan(&nsum, &offset, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)A));
1560:     offset -= nsum;
1561:     for (j = 0; j < vs->nc; j++) {
1562:       MatNestFindNonzeroSubMatCol(A, j, &sub);
1563:       MatGetLocalSize(sub, NULL, &n);
1564:       MatGetBlockSizes(sub, NULL, &bs);
1565:       ISCreateStride(PetscObjectComm((PetscObject)sub), n, offset, 1, &vs->isglobal.col[j]);
1566:       ISSetBlockSize(vs->isglobal.col[j], bs);
1567:       offset += n;
1568:     }
1569:   }

1571:   /* Set up the local ISs */
1572:   PetscMalloc1(vs->nr, &vs->islocal.row);
1573:   PetscMalloc1(vs->nc, &vs->islocal.col);
1574:   for (i = 0, offset = 0; i < vs->nr; i++) {
1575:     IS                     isloc;
1576:     ISLocalToGlobalMapping rmap = NULL;
1577:     PetscInt               nlocal, bs;
1578:     MatNestFindNonzeroSubMatRow(A, i, &sub);
1579:     if (sub) MatGetLocalToGlobalMapping(sub, &rmap, NULL);
1580:     if (rmap) {
1581:       MatGetBlockSizes(sub, &bs, NULL);
1582:       ISLocalToGlobalMappingGetSize(rmap, &nlocal);
1583:       ISCreateStride(PETSC_COMM_SELF, nlocal, offset, 1, &isloc);
1584:       ISSetBlockSize(isloc, bs);
1585:     } else {
1586:       nlocal = 0;
1587:       isloc  = NULL;
1588:     }
1589:     vs->islocal.row[i] = isloc;
1590:     offset += nlocal;
1591:   }
1592:   for (i = 0, offset = 0; i < vs->nc; i++) {
1593:     IS                     isloc;
1594:     ISLocalToGlobalMapping cmap = NULL;
1595:     PetscInt               nlocal, bs;
1596:     MatNestFindNonzeroSubMatCol(A, i, &sub);
1597:     if (sub) MatGetLocalToGlobalMapping(sub, NULL, &cmap);
1598:     if (cmap) {
1599:       MatGetBlockSizes(sub, NULL, &bs);
1600:       ISLocalToGlobalMappingGetSize(cmap, &nlocal);
1601:       ISCreateStride(PETSC_COMM_SELF, nlocal, offset, 1, &isloc);
1602:       ISSetBlockSize(isloc, bs);
1603:     } else {
1604:       nlocal = 0;
1605:       isloc  = NULL;
1606:     }
1607:     vs->islocal.col[i] = isloc;
1608:     offset += nlocal;
1609:   }

1611:   /* Set up the aggregate ISLocalToGlobalMapping */
1612:   {
1613:     ISLocalToGlobalMapping rmap, cmap;
1614:     MatNestCreateAggregateL2G_Private(A, vs->nr, vs->islocal.row, vs->isglobal.row, PETSC_FALSE, &rmap);
1615:     MatNestCreateAggregateL2G_Private(A, vs->nc, vs->islocal.col, vs->isglobal.col, PETSC_TRUE, &cmap);
1616:     if (rmap && cmap) MatSetLocalToGlobalMapping(A, rmap, cmap);
1617:     ISLocalToGlobalMappingDestroy(&rmap);
1618:     ISLocalToGlobalMappingDestroy(&cmap);
1619:   }

1621:   if (PetscDefined(USE_DEBUG)) {
1622:     for (i = 0; i < vs->nr; i++) {
1623:       for (j = 0; j < vs->nc; j++) {
1624:         PetscInt m, n, M, N, mi, ni, Mi, Ni;
1625:         Mat      B = vs->m[i][j];
1626:         if (!B) continue;
1627:         MatGetSize(B, &M, &N);
1628:         MatGetLocalSize(B, &m, &n);
1629:         ISGetSize(vs->isglobal.row[i], &Mi);
1630:         ISGetSize(vs->isglobal.col[j], &Ni);
1631:         ISGetLocalSize(vs->isglobal.row[i], &mi);
1632:         ISGetLocalSize(vs->isglobal.col[j], &ni);
1635:       }
1636:     }
1637:   }

1639:   /* Set A->assembled if all non-null blocks are currently assembled */
1640:   for (i = 0; i < vs->nr; i++) {
1641:     for (j = 0; j < vs->nc; j++) {
1642:       if (vs->m[i][j] && !vs->m[i][j]->assembled) return 0;
1643:     }
1644:   }
1645:   A->assembled = PETSC_TRUE;
1646:   return 0;
1647: }

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

1652:    Collective

1654:    Input Parameters:
1655: +  comm - Communicator for the new `MATNEST`
1656: .  nr - number of nested row blocks
1657: .  is_row - index sets for each nested row block, or NULL to make contiguous
1658: .  nc - number of nested column blocks
1659: .  is_col - index sets for each nested column block, or NULL to make contiguous
1660: -  a - row-aligned array of nr*nc submatrices, empty submatrices can be passed using NULL

1662:    Output Parameter:
1663: .  B - new matrix

1665:    Level: advanced

1667: .seealso: `MATNEST`, `MatCreate()`, `VecCreateNest()`, `DMCreateMatrix()`, `MATNEST`, `MatNestSetSubMat()`,
1668:           `MatNestGetSubMat()`, `MatNestGetLocalISs()`, `MatNestGetSize()`,
1669:           `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()`
1670: @*/
1671: PetscErrorCode MatCreateNest(MPI_Comm comm, PetscInt nr, const IS is_row[], PetscInt nc, const IS is_col[], const Mat a[], Mat *B)
1672: {
1673:   Mat A;

1675:   *B = NULL;
1676:   MatCreate(comm, &A);
1677:   MatSetType(A, MATNEST);
1678:   A->preallocated = PETSC_TRUE;
1679:   MatNestSetSubMats(A, nr, is_row, nc, is_col, a);
1680:   *B = A;
1681:   return 0;
1682: }

1684: PetscErrorCode MatConvert_Nest_SeqAIJ_fast(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1685: {
1686:   Mat_Nest     *nest = (Mat_Nest *)A->data;
1687:   Mat          *trans;
1688:   PetscScalar **avv;
1689:   PetscScalar  *vv;
1690:   PetscInt    **aii, **ajj;
1691:   PetscInt     *ii, *jj, *ci;
1692:   PetscInt      nr, nc, nnz, i, j;
1693:   PetscBool     done;

1695:   MatGetSize(A, &nr, &nc);
1696:   if (reuse == MAT_REUSE_MATRIX) {
1697:     PetscInt rnr;

1699:     MatGetRowIJ(*newmat, 0, PETSC_FALSE, PETSC_FALSE, &rnr, (const PetscInt **)&ii, (const PetscInt **)&jj, &done);
1702:     MatSeqAIJGetArray(*newmat, &vv);
1703:   }
1704:   /* extract CSR for nested SeqAIJ matrices */
1705:   nnz = 0;
1706:   PetscCalloc4(nest->nr * nest->nc, &aii, nest->nr * nest->nc, &ajj, nest->nr * nest->nc, &avv, nest->nr * nest->nc, &trans);
1707:   for (i = 0; i < nest->nr; ++i) {
1708:     for (j = 0; j < nest->nc; ++j) {
1709:       Mat B = nest->m[i][j];
1710:       if (B) {
1711:         PetscScalar *naa;
1712:         PetscInt    *nii, *njj, nnr;
1713:         PetscBool    istrans;

1715:         PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &istrans);
1716:         if (istrans) {
1717:           Mat Bt;

1719:           MatTransposeGetMat(B, &Bt);
1720:           MatTranspose(Bt, MAT_INITIAL_MATRIX, &trans[i * nest->nc + j]);
1721:           B = trans[i * nest->nc + j];
1722:         } else {
1723:           PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &istrans);
1724:           if (istrans) {
1725:             Mat Bt;

1727:             MatHermitianTransposeGetMat(B, &Bt);
1728:             MatHermitianTranspose(Bt, MAT_INITIAL_MATRIX, &trans[i * nest->nc + j]);
1729:             B = trans[i * nest->nc + j];
1730:           }
1731:         }
1732:         MatGetRowIJ(B, 0, PETSC_FALSE, PETSC_FALSE, &nnr, (const PetscInt **)&nii, (const PetscInt **)&njj, &done);
1734:         MatSeqAIJGetArray(B, &naa);
1735:         nnz += nii[nnr];

1737:         aii[i * nest->nc + j] = nii;
1738:         ajj[i * nest->nc + j] = njj;
1739:         avv[i * nest->nc + j] = naa;
1740:       }
1741:     }
1742:   }
1743:   if (reuse != MAT_REUSE_MATRIX) {
1744:     PetscMalloc1(nr + 1, &ii);
1745:     PetscMalloc1(nnz, &jj);
1746:     PetscMalloc1(nnz, &vv);
1747:   } else {
1749:   }

1751:   /* new row pointer */
1752:   PetscArrayzero(ii, nr + 1);
1753:   for (i = 0; i < nest->nr; ++i) {
1754:     PetscInt ncr, rst;

1756:     ISStrideGetInfo(nest->isglobal.row[i], &rst, NULL);
1757:     ISGetLocalSize(nest->isglobal.row[i], &ncr);
1758:     for (j = 0; j < nest->nc; ++j) {
1759:       if (aii[i * nest->nc + j]) {
1760:         PetscInt *nii = aii[i * nest->nc + j];
1761:         PetscInt  ir;

1763:         for (ir = rst; ir < ncr + rst; ++ir) {
1764:           ii[ir + 1] += nii[1] - nii[0];
1765:           nii++;
1766:         }
1767:       }
1768:     }
1769:   }
1770:   for (i = 0; i < nr; i++) ii[i + 1] += ii[i];

1772:   /* construct CSR for the new matrix */
1773:   PetscCalloc1(nr, &ci);
1774:   for (i = 0; i < nest->nr; ++i) {
1775:     PetscInt ncr, rst;

1777:     ISStrideGetInfo(nest->isglobal.row[i], &rst, NULL);
1778:     ISGetLocalSize(nest->isglobal.row[i], &ncr);
1779:     for (j = 0; j < nest->nc; ++j) {
1780:       if (aii[i * nest->nc + j]) {
1781:         PetscScalar *nvv = avv[i * nest->nc + j];
1782:         PetscInt    *nii = aii[i * nest->nc + j];
1783:         PetscInt    *njj = ajj[i * nest->nc + j];
1784:         PetscInt     ir, cst;

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

1790:           for (ij = 0; ij < rsize; ij++) {
1791:             jj[ist + ij] = *njj + cst;
1792:             vv[ist + ij] = *nvv;
1793:             njj++;
1794:             nvv++;
1795:           }
1796:           ci[ir] += rsize;
1797:           nii++;
1798:         }
1799:       }
1800:     }
1801:   }
1802:   PetscFree(ci);

1804:   /* restore info */
1805:   for (i = 0; i < nest->nr; ++i) {
1806:     for (j = 0; j < nest->nc; ++j) {
1807:       Mat B = nest->m[i][j];
1808:       if (B) {
1809:         PetscInt nnr = 0, k = i * nest->nc + j;

1811:         B = (trans[k] ? trans[k] : B);
1812:         MatRestoreRowIJ(B, 0, PETSC_FALSE, PETSC_FALSE, &nnr, (const PetscInt **)&aii[k], (const PetscInt **)&ajj[k], &done);
1814:         MatSeqAIJRestoreArray(B, &avv[k]);
1815:         MatDestroy(&trans[k]);
1816:       }
1817:     }
1818:   }
1819:   PetscFree4(aii, ajj, avv, trans);

1821:   /* finalize newmat */
1822:   if (reuse == MAT_INITIAL_MATRIX) {
1823:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), nr, nc, ii, jj, vv, newmat);
1824:   } else if (reuse == MAT_INPLACE_MATRIX) {
1825:     Mat B;

1827:     MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A), nr, nc, ii, jj, vv, &B);
1828:     MatHeaderReplace(A, &B);
1829:   }
1830:   MatAssemblyBegin(*newmat, MAT_FINAL_ASSEMBLY);
1831:   MatAssemblyEnd(*newmat, MAT_FINAL_ASSEMBLY);
1832:   {
1833:     Mat_SeqAIJ *a = (Mat_SeqAIJ *)((*newmat)->data);
1834:     a->free_a     = PETSC_TRUE;
1835:     a->free_ij    = PETSC_TRUE;
1836:   }
1837:   return 0;
1838: }

1840: PETSC_INTERN PetscErrorCode MatAXPY_Dense_Nest(Mat Y, PetscScalar a, Mat X)
1841: {
1842:   Mat_Nest *nest = (Mat_Nest *)X->data;
1843:   PetscInt  i, j, k, rstart;
1844:   PetscBool flg;

1846:   /* Fill by row */
1847:   for (j = 0; j < nest->nc; ++j) {
1848:     /* Using global column indices and ISAllGather() is not scalable. */
1849:     IS              bNis;
1850:     PetscInt        bN;
1851:     const PetscInt *bNindices;
1852:     ISAllGather(nest->isglobal.col[j], &bNis);
1853:     ISGetSize(bNis, &bN);
1854:     ISGetIndices(bNis, &bNindices);
1855:     for (i = 0; i < nest->nr; ++i) {
1856:       Mat             B = nest->m[i][j], D = NULL;
1857:       PetscInt        bm, br;
1858:       const PetscInt *bmindices;
1859:       if (!B) continue;
1860:       PetscObjectTypeCompareAny((PetscObject)B, &flg, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, "");
1861:       if (flg) {
1862:         PetscTryMethod(B, "MatTransposeGetMat_C", (Mat, Mat *), (B, &D));
1863:         PetscTryMethod(B, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (B, &D));
1864:         MatConvert(B, ((PetscObject)D)->type_name, MAT_INITIAL_MATRIX, &D);
1865:         B = D;
1866:       }
1867:       PetscObjectTypeCompareAny((PetscObject)B, &flg, MATSEQSBAIJ, MATMPISBAIJ, "");
1868:       if (flg) {
1869:         if (D) MatConvert(D, MATBAIJ, MAT_INPLACE_MATRIX, &D);
1870:         else MatConvert(B, MATBAIJ, MAT_INITIAL_MATRIX, &D);
1871:         B = D;
1872:       }
1873:       ISGetLocalSize(nest->isglobal.row[i], &bm);
1874:       ISGetIndices(nest->isglobal.row[i], &bmindices);
1875:       MatGetOwnershipRange(B, &rstart, NULL);
1876:       for (br = 0; br < bm; ++br) {
1877:         PetscInt           row = bmindices[br], brncols, *cols;
1878:         const PetscInt    *brcols;
1879:         const PetscScalar *brcoldata;
1880:         PetscScalar       *vals = NULL;
1881:         MatGetRow(B, br + rstart, &brncols, &brcols, &brcoldata);
1882:         PetscMalloc1(brncols, &cols);
1883:         for (k = 0; k < brncols; k++) cols[k] = bNindices[brcols[k]];
1884:         /*
1885:           Nest blocks are required to be nonoverlapping -- otherwise nest and monolithic index layouts wouldn't match.
1886:           Thus, we could use INSERT_VALUES, but I prefer ADD_VALUES.
1887:          */
1888:         if (a != 1.0) {
1889:           PetscMalloc1(brncols, &vals);
1890:           for (k = 0; k < brncols; k++) vals[k] = a * brcoldata[k];
1891:           MatSetValues(Y, 1, &row, brncols, cols, vals, ADD_VALUES);
1892:           PetscFree(vals);
1893:         } else {
1894:           MatSetValues(Y, 1, &row, brncols, cols, brcoldata, ADD_VALUES);
1895:         }
1896:         MatRestoreRow(B, br + rstart, &brncols, &brcols, &brcoldata);
1897:         PetscFree(cols);
1898:       }
1899:       if (D) MatDestroy(&D);
1900:       ISRestoreIndices(nest->isglobal.row[i], &bmindices);
1901:     }
1902:     ISRestoreIndices(bNis, &bNindices);
1903:     ISDestroy(&bNis);
1904:   }
1905:   MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY);
1906:   MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY);
1907:   return 0;
1908: }

1910: PetscErrorCode MatConvert_Nest_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
1911: {
1912:   Mat_Nest   *nest = (Mat_Nest *)A->data;
1913:   PetscInt    m, n, M, N, i, j, k, *dnnz, *onnz, rstart, cstart, cend;
1914:   PetscMPIInt size;
1915:   Mat         C;

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

1922:     PetscStrcmp(newtype, MATAIJ, &fast);
1923:     if (!fast) PetscStrcmp(newtype, MATSEQAIJ, &fast);
1924:     for (i = 0; i < nest->nr && fast; ++i) {
1925:       for (j = 0; j < nest->nc && fast; ++j) {
1926:         Mat B = nest->m[i][j];
1927:         if (B) {
1928:           PetscObjectTypeCompare((PetscObject)B, MATSEQAIJ, &fast);
1929:           if (!fast) {
1930:             PetscBool istrans;

1932:             PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &istrans);
1933:             if (istrans) {
1934:               Mat Bt;

1936:               MatTransposeGetMat(B, &Bt);
1937:               PetscObjectTypeCompare((PetscObject)Bt, MATSEQAIJ, &fast);
1938:             } else {
1939:               PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &istrans);
1940:               if (istrans) {
1941:                 Mat Bt;

1943:                 MatHermitianTransposeGetMat(B, &Bt);
1944:                 PetscObjectTypeCompare((PetscObject)Bt, MATSEQAIJ, &fast);
1945:               }
1946:             }
1947:           }
1948:         }
1949:       }
1950:     }
1951:     for (i = 0, nf = 0; i < nest->nr && fast; ++i) {
1952:       PetscObjectTypeCompare((PetscObject)nest->isglobal.row[i], ISSTRIDE, &fast);
1953:       if (fast) {
1954:         PetscInt f, s;

1956:         ISStrideGetInfo(nest->isglobal.row[i], &f, &s);
1957:         if (f != nf || s != 1) {
1958:           fast = PETSC_FALSE;
1959:         } else {
1960:           ISGetSize(nest->isglobal.row[i], &f);
1961:           nf += f;
1962:         }
1963:       }
1964:     }
1965:     for (i = 0, nf = 0; i < nest->nc && fast; ++i) {
1966:       PetscObjectTypeCompare((PetscObject)nest->isglobal.col[i], ISSTRIDE, &fast);
1967:       if (fast) {
1968:         PetscInt f, s;

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

2095: PetscErrorCode MatConvert_Nest_Dense(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
2096: {
2097:   Mat      B;
2098:   PetscInt m, n, M, N;

2100:   MatGetSize(A, &M, &N);
2101:   MatGetLocalSize(A, &m, &n);
2102:   if (reuse == MAT_REUSE_MATRIX) {
2103:     B = *newmat;
2104:     MatZeroEntries(B);
2105:   } else {
2106:     MatCreateDense(PetscObjectComm((PetscObject)A), m, PETSC_DECIDE, M, N, NULL, &B);
2107:   }
2108:   MatAXPY_Dense_Nest(B, 1.0, A);
2109:   if (reuse == MAT_INPLACE_MATRIX) {
2110:     MatHeaderReplace(A, &B);
2111:   } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
2112:   return 0;
2113: }

2115: PetscErrorCode MatHasOperation_Nest(Mat mat, MatOperation op, PetscBool *has)
2116: {
2117:   Mat_Nest    *bA = (Mat_Nest *)mat->data;
2118:   MatOperation opAdd;
2119:   PetscInt     i, j, nr = bA->nr, nc = bA->nc;
2120:   PetscBool    flg;

2122:   *has = PETSC_FALSE;
2123:   if (op == MATOP_MULT || op == MATOP_MULT_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
2124:     opAdd = (op == MATOP_MULT || op == MATOP_MULT_ADD ? MATOP_MULT_ADD : MATOP_MULT_TRANSPOSE_ADD);
2125:     for (j = 0; j < nc; j++) {
2126:       for (i = 0; i < nr; i++) {
2127:         if (!bA->m[i][j]) continue;
2128:         MatHasOperation(bA->m[i][j], opAdd, &flg);
2129:         if (!flg) return 0;
2130:       }
2131:     }
2132:   }
2133:   if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
2134:   return 0;
2135: }

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

2140:   Level: intermediate

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

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

2151: .seealso: `MATNEST`, `MatCreate()`, `MatType`, `MatCreateNest()`, `MatNestSetSubMat()`, `MatNestGetSubMat()`,
2152:           `VecCreateNest()`, `DMCreateMatrix()`, `DMCOMPOSITE`, `MatNestSetVecType()`, `MatNestGetLocalISs()`,
2153:           `MatNestGetISs()`, `MatNestSetSubMats()`, `MatNestGetSubMats()`
2154: M*/
2155: PETSC_EXTERN PetscErrorCode MatCreate_Nest(Mat A)
2156: {
2157:   Mat_Nest *s;

2159:   PetscNew(&s);
2160:   A->data = (void *)s;

2162:   s->nr            = -1;
2163:   s->nc            = -1;
2164:   s->m             = NULL;
2165:   s->splitassembly = PETSC_FALSE;

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

2169:   A->ops->mult                  = MatMult_Nest;
2170:   A->ops->multadd               = MatMultAdd_Nest;
2171:   A->ops->multtranspose         = MatMultTranspose_Nest;
2172:   A->ops->multtransposeadd      = MatMultTransposeAdd_Nest;
2173:   A->ops->transpose             = MatTranspose_Nest;
2174:   A->ops->assemblybegin         = MatAssemblyBegin_Nest;
2175:   A->ops->assemblyend           = MatAssemblyEnd_Nest;
2176:   A->ops->zeroentries           = MatZeroEntries_Nest;
2177:   A->ops->copy                  = MatCopy_Nest;
2178:   A->ops->axpy                  = MatAXPY_Nest;
2179:   A->ops->duplicate             = MatDuplicate_Nest;
2180:   A->ops->createsubmatrix       = MatCreateSubMatrix_Nest;
2181:   A->ops->destroy               = MatDestroy_Nest;
2182:   A->ops->view                  = MatView_Nest;
2183:   A->ops->getvecs               = NULL; /* Use VECNEST by calling MatNestSetVecType(A,VECNEST) */
2184:   A->ops->getlocalsubmatrix     = MatGetLocalSubMatrix_Nest;
2185:   A->ops->restorelocalsubmatrix = MatRestoreLocalSubMatrix_Nest;
2186:   A->ops->getdiagonal           = MatGetDiagonal_Nest;
2187:   A->ops->diagonalscale         = MatDiagonalScale_Nest;
2188:   A->ops->scale                 = MatScale_Nest;
2189:   A->ops->shift                 = MatShift_Nest;
2190:   A->ops->diagonalset           = MatDiagonalSet_Nest;
2191:   A->ops->setrandom             = MatSetRandom_Nest;
2192:   A->ops->hasoperation          = MatHasOperation_Nest;
2193:   A->ops->missingdiagonal       = MatMissingDiagonal_Nest;

2195:   A->spptr     = NULL;
2196:   A->assembled = PETSC_FALSE;

2198:   /* expose Nest api's */
2199:   PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMat_C", MatNestGetSubMat_Nest);
2200:   PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMat_C", MatNestSetSubMat_Nest);
2201:   PetscObjectComposeFunction((PetscObject)A, "MatNestGetSubMats_C", MatNestGetSubMats_Nest);
2202:   PetscObjectComposeFunction((PetscObject)A, "MatNestGetSize_C", MatNestGetSize_Nest);
2203:   PetscObjectComposeFunction((PetscObject)A, "MatNestGetISs_C", MatNestGetISs_Nest);
2204:   PetscObjectComposeFunction((PetscObject)A, "MatNestGetLocalISs_C", MatNestGetLocalISs_Nest);
2205:   PetscObjectComposeFunction((PetscObject)A, "MatNestSetVecType_C", MatNestSetVecType_Nest);
2206:   PetscObjectComposeFunction((PetscObject)A, "MatNestSetSubMats_C", MatNestSetSubMats_Nest);
2207:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpiaij_C", MatConvert_Nest_AIJ);
2208:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqaij_C", MatConvert_Nest_AIJ);
2209:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_aij_C", MatConvert_Nest_AIJ);
2210:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_is_C", MatConvert_Nest_IS);
2211:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_mpidense_C", MatConvert_Nest_Dense);
2212:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_nest_seqdense_C", MatConvert_Nest_Dense);
2213:   PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_seqdense_C", MatProductSetFromOptions_Nest_Dense);
2214:   PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_mpidense_C", MatProductSetFromOptions_Nest_Dense);
2215:   PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_nest_dense_C", MatProductSetFromOptions_Nest_Dense);

2217:   PetscObjectChangeTypeName((PetscObject)A, MATNEST);
2218:   return 0;
2219: }