Actual source code: bjacobi.c

  1: /*
  2:    Defines a block Jacobi preconditioner.
  3: */

  5: #include <../src/ksp/pc/impls/bjacobi/bjacobi.h>

  7: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC, Mat, Mat);
  8: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC, Mat, Mat);
  9: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC);

 11: static PetscErrorCode PCSetUp_BJacobi(PC pc)
 12: {
 13:   PC_BJacobi *jac = (PC_BJacobi *)pc->data;
 14:   Mat         mat = pc->mat, pmat = pc->pmat;
 15:   PetscBool   hasop;
 16:   PetscInt    N, M, start, i, sum, end;
 17:   PetscInt    bs, i_start = -1, i_end = -1;
 18:   PetscMPIInt rank, size;

 20:   PetscFunctionBegin;
 21:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
 22:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
 23:   PetscCall(MatGetLocalSize(pc->pmat, &M, &N));
 24:   PetscCall(MatGetBlockSize(pc->pmat, &bs));

 26:   if (jac->n > 0 && jac->n < size) {
 27:     PetscCall(PCSetUp_BJacobi_Multiproc(pc));
 28:     PetscFunctionReturn(PETSC_SUCCESS);
 29:   }

 31:   /*    Determines the number of blocks assigned to each processor */
 32:   /*   local block count  given */
 33:   if (jac->n_local > 0 && jac->n < 0) {
 34:     PetscCallMPI(MPIU_Allreduce(&jac->n_local, &jac->n, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)pc)));
 35:     if (jac->l_lens) { /* check that user set these correctly */
 36:       sum = 0;
 37:       for (i = 0; i < jac->n_local; i++) {
 38:         PetscCheck(jac->l_lens[i] / bs * bs == jac->l_lens[i], PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Mat blocksize doesn't match block Jacobi layout");
 39:         sum += jac->l_lens[i];
 40:       }
 41:       PetscCheck(sum == M, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local lens set incorrectly");
 42:     } else {
 43:       PetscCall(PetscMalloc1(jac->n_local, &jac->l_lens));
 44:       for (i = 0; i < jac->n_local; i++) jac->l_lens[i] = bs * ((M / bs) / jac->n_local + (((M / bs) % jac->n_local) > i));
 45:     }
 46:   } else if (jac->n > 0 && jac->n_local < 0) { /* global block count given */
 47:     /* global blocks given: determine which ones are local */
 48:     if (jac->g_lens) {
 49:       /* check if the g_lens is has valid entries */
 50:       for (i = 0; i < jac->n; i++) {
 51:         PetscCheck(jac->g_lens[i], PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Zero block not allowed");
 52:         PetscCheck(jac->g_lens[i] / bs * bs == jac->g_lens[i], PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Mat blocksize doesn't match block Jacobi layout");
 53:       }
 54:       if (size == 1) {
 55:         jac->n_local = jac->n;
 56:         PetscCall(PetscMalloc1(jac->n_local, &jac->l_lens));
 57:         PetscCall(PetscArraycpy(jac->l_lens, jac->g_lens, jac->n_local));
 58:         /* check that user set these correctly */
 59:         sum = 0;
 60:         for (i = 0; i < jac->n_local; i++) sum += jac->l_lens[i];
 61:         PetscCheck(sum == M, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Global lens set incorrectly");
 62:       } else {
 63:         PetscCall(MatGetOwnershipRange(pc->pmat, &start, &end));
 64:         /* loop over blocks determining first one owned by me */
 65:         sum = 0;
 66:         for (i = 0; i < jac->n + 1; i++) {
 67:           if (sum == start) {
 68:             i_start = i;
 69:             goto start_1;
 70:           }
 71:           if (i < jac->n) sum += jac->g_lens[i];
 72:         }
 73:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Block sizes used in PCBJacobiSetTotalBlocks()\nare not compatible with parallel matrix layout");
 74:       start_1:
 75:         for (i = i_start; i < jac->n + 1; i++) {
 76:           if (sum == end) {
 77:             i_end = i;
 78:             goto end_1;
 79:           }
 80:           if (i < jac->n) sum += jac->g_lens[i];
 81:         }
 82:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Block sizes used in PCBJacobiSetTotalBlocks()\nare not compatible with parallel matrix layout");
 83:       end_1:
 84:         jac->n_local = i_end - i_start;
 85:         PetscCall(PetscMalloc1(jac->n_local, &jac->l_lens));
 86:         PetscCall(PetscArraycpy(jac->l_lens, jac->g_lens + i_start, jac->n_local));
 87:       }
 88:     } else { /* no global blocks given, determine then using default layout */
 89:       jac->n_local = jac->n / size + ((jac->n % size) > rank);
 90:       PetscCall(PetscMalloc1(jac->n_local, &jac->l_lens));
 91:       for (i = 0; i < jac->n_local; i++) {
 92:         jac->l_lens[i] = ((M / bs) / jac->n_local + (((M / bs) % jac->n_local) > i)) * bs;
 93:         PetscCheck(jac->l_lens[i], PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Too many blocks given");
 94:       }
 95:     }
 96:   } else if (jac->n < 0 && jac->n_local < 0) { /* no blocks given */
 97:     jac->n       = size;
 98:     jac->n_local = 1;
 99:     PetscCall(PetscMalloc1(1, &jac->l_lens));
100:     jac->l_lens[0] = M;
101:   } else { /* jac->n > 0 && jac->n_local > 0 */
102:     if (!jac->l_lens) {
103:       PetscCall(PetscMalloc1(jac->n_local, &jac->l_lens));
104:       for (i = 0; i < jac->n_local; i++) jac->l_lens[i] = bs * ((M / bs) / jac->n_local + (((M / bs) % jac->n_local) > i));
105:     }
106:   }
107:   PetscCheck(jac->n_local >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of blocks is less than number of processors");

109:   /*    Determines mat and pmat */
110:   PetscCall(MatHasOperation(pc->mat, MATOP_GET_DIAGONAL_BLOCK, &hasop));
111:   if (!hasop && size == 1) {
112:     mat  = pc->mat;
113:     pmat = pc->pmat;
114:   } else {
115:     if (pc->useAmat) {
116:       /* use block from Amat matrix, not Pmat for local MatMult() */
117:       PetscCall(MatGetDiagonalBlock(pc->mat, &mat));
118:     }
119:     if (pc->pmat != pc->mat || !pc->useAmat) PetscCall(MatGetDiagonalBlock(pc->pmat, &pmat));
120:     else pmat = mat;
121:   }

123:   /*
124:      Setup code depends on the number of blocks
125:   */
126:   if (jac->n_local == 1) {
127:     PetscCall(PCSetUp_BJacobi_Singleblock(pc, mat, pmat));
128:   } else {
129:     PetscCall(PCSetUp_BJacobi_Multiblock(pc, mat, pmat));
130:   }
131:   PetscFunctionReturn(PETSC_SUCCESS);
132: }

134: /* Default destroy, if it has never been setup */
135: static PetscErrorCode PCDestroy_BJacobi(PC pc)
136: {
137:   PC_BJacobi *jac = (PC_BJacobi *)pc->data;

139:   PetscFunctionBegin;
140:   PetscCall(PetscFree(jac->g_lens));
141:   PetscCall(PetscFree(jac->l_lens));
142:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetSubKSP_C", NULL));
143:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetTotalBlocks_C", NULL));
144:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetTotalBlocks_C", NULL));
145:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetLocalBlocks_C", NULL));
146:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetLocalBlocks_C", NULL));
147:   PetscCall(PetscFree(pc->data));
148:   PetscFunctionReturn(PETSC_SUCCESS);
149: }

151: static PetscErrorCode PCSetFromOptions_BJacobi(PC pc, PetscOptionItems PetscOptionsObject)
152: {
153:   PC_BJacobi *jac = (PC_BJacobi *)pc->data;
154:   PetscInt    blocks, i;
155:   PetscBool   flg;

157:   PetscFunctionBegin;
158:   PetscOptionsHeadBegin(PetscOptionsObject, "Block Jacobi options");
159:   PetscCall(PetscOptionsInt("-pc_bjacobi_blocks", "Total number of blocks", "PCBJacobiSetTotalBlocks", jac->n, &blocks, &flg));
160:   if (flg) PetscCall(PCBJacobiSetTotalBlocks(pc, blocks, NULL));
161:   PetscCall(PetscOptionsInt("-pc_bjacobi_local_blocks", "Local number of blocks", "PCBJacobiSetLocalBlocks", jac->n_local, &blocks, &flg));
162:   if (flg) PetscCall(PCBJacobiSetLocalBlocks(pc, blocks, NULL));
163:   if (jac->ksp) {
164:     /* The sub-KSP has already been set up (e.g., PCSetUp_BJacobi_Singleblock), but KSPSetFromOptions was not called
165:      * unless we had already been called. */
166:     for (i = 0; i < jac->n_local; i++) PetscCall(KSPSetFromOptions(jac->ksp[i]));
167:   }
168:   PetscOptionsHeadEnd();
169:   PetscFunctionReturn(PETSC_SUCCESS);
170: }

172: #include <petscdraw.h>
173: static PetscErrorCode PCView_BJacobi(PC pc, PetscViewer viewer)
174: {
175:   PC_BJacobi           *jac   = (PC_BJacobi *)pc->data;
176:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
177:   PetscMPIInt           rank;
178:   PetscInt              i;
179:   PetscBool             isascii, isstring, isdraw;
180:   PetscViewer           sviewer;
181:   PetscViewerFormat     format;
182:   const char           *prefix;

184:   PetscFunctionBegin;
185:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
186:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
187:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
188:   if (isascii) {
189:     if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, "  using Amat local matrix, number of blocks = %" PetscInt_FMT "\n", jac->n));
190:     PetscCall(PetscViewerASCIIPrintf(viewer, "  number of blocks = %" PetscInt_FMT "\n", jac->n));
191:     PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
192:     PetscCall(PetscViewerGetFormat(viewer, &format));
193:     if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
194:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Local solver information for first block is in the following KSP and PC objects on rank 0:\n"));
195:       PetscCall(PCGetOptionsPrefix(pc, &prefix));
196:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Use -%sksp_view ::ascii_info_detail to display information for all blocks\n", prefix ? prefix : ""));
197:       if (jac->ksp && !jac->psubcomm) {
198:         PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
199:         if (rank == 0) {
200:           PetscCall(PetscViewerASCIIPushTab(sviewer));
201:           PetscCall(KSPView(jac->ksp[0], sviewer));
202:           PetscCall(PetscViewerASCIIPopTab(sviewer));
203:         }
204:         PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
205:         /*  extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
206:         PetscCall(PetscViewerASCIIPopSynchronized(viewer));
207:       } else if (mpjac && jac->ksp && mpjac->psubcomm) {
208:         PetscCall(PetscViewerGetSubViewer(viewer, mpjac->psubcomm->child, &sviewer));
209:         if (!mpjac->psubcomm->color) {
210:           PetscCall(PetscViewerASCIIPushTab(sviewer));
211:           PetscCall(KSPView(*jac->ksp, sviewer));
212:           PetscCall(PetscViewerASCIIPopTab(sviewer));
213:         }
214:         PetscCall(PetscViewerRestoreSubViewer(viewer, mpjac->psubcomm->child, &sviewer));
215:         /*  extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
216:         PetscCall(PetscViewerASCIIPopSynchronized(viewer));
217:       }
218:     } else {
219:       PetscInt n_global;
220:       PetscCallMPI(MPIU_Allreduce(&jac->n_local, &n_global, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc)));
221:       PetscCall(PetscViewerASCIIPushSynchronized(viewer));
222:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Local solver information for each block is in the following KSP and PC objects:\n"));
223:       PetscCall(PetscViewerASCIIPushTab(viewer));
224:       PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
225:       PetscCall(PetscViewerASCIIPrintf(sviewer, "[%d] number of local blocks = %" PetscInt_FMT ", first local block number = %" PetscInt_FMT "\n", rank, jac->n_local, jac->first_local));
226:       for (i = 0; i < jac->n_local; i++) {
227:         PetscCall(PetscViewerASCIIPrintf(sviewer, "[%d] local block number %" PetscInt_FMT "\n", rank, i));
228:         PetscCall(KSPView(jac->ksp[i], sviewer));
229:         PetscCall(PetscViewerASCIIPrintf(sviewer, "- - - - - - - - - - - - - - - - - -\n"));
230:       }
231:       PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
232:       PetscCall(PetscViewerASCIIPopTab(viewer));
233:       PetscCall(PetscViewerASCIIPopSynchronized(viewer));
234:     }
235:   } else if (isstring) {
236:     PetscCall(PetscViewerStringSPrintf(viewer, " blks=%" PetscInt_FMT, jac->n));
237:     PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
238:     if (jac->ksp) PetscCall(KSPView(jac->ksp[0], sviewer));
239:     PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
240:   } else if (isdraw) {
241:     PetscDraw draw;
242:     char      str[25];
243:     PetscReal x, y, bottom, h;

245:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
246:     PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
247:     PetscCall(PetscSNPrintf(str, 25, "Number blocks %" PetscInt_FMT, jac->n));
248:     PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
249:     bottom = y - h;
250:     PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
251:     /* warning the communicator on viewer is different then on ksp in parallel */
252:     if (jac->ksp) PetscCall(KSPView(jac->ksp[0], viewer));
253:     PetscCall(PetscDrawPopCurrentPoint(draw));
254:   }
255:   PetscFunctionReturn(PETSC_SUCCESS);
256: }

258: static PetscErrorCode PCBJacobiGetSubKSP_BJacobi(PC pc, PetscInt *n_local, PetscInt *first_local, KSP **ksp)
259: {
260:   PC_BJacobi *jac = (PC_BJacobi *)pc->data;

262:   PetscFunctionBegin;
263:   PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() first");

265:   if (n_local) *n_local = jac->n_local;
266:   if (first_local) *first_local = jac->first_local;
267:   if (ksp) *ksp = jac->ksp;
268:   PetscFunctionReturn(PETSC_SUCCESS);
269: }

271: static PetscErrorCode PCBJacobiSetTotalBlocks_BJacobi(PC pc, PetscInt blocks, const PetscInt *lens)
272: {
273:   PC_BJacobi *jac = (PC_BJacobi *)pc->data;

275:   PetscFunctionBegin;
276:   PetscCheck(!pc->setupcalled || jac->n == blocks, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Cannot alter number of blocks after PCSetUp()/KSPSetUp() has been called");
277:   jac->n = blocks;
278:   if (!lens) jac->g_lens = NULL;
279:   else {
280:     PetscCall(PetscMalloc1(blocks, &jac->g_lens));
281:     PetscCall(PetscArraycpy(jac->g_lens, lens, blocks));
282:   }
283:   PetscFunctionReturn(PETSC_SUCCESS);
284: }

286: static PetscErrorCode PCBJacobiGetTotalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
287: {
288:   PC_BJacobi *jac = (PC_BJacobi *)pc->data;

290:   PetscFunctionBegin;
291:   *blocks = jac->n;
292:   if (lens) *lens = jac->g_lens;
293:   PetscFunctionReturn(PETSC_SUCCESS);
294: }

296: static PetscErrorCode PCBJacobiSetLocalBlocks_BJacobi(PC pc, PetscInt blocks, const PetscInt lens[])
297: {
298:   PC_BJacobi *jac;

300:   PetscFunctionBegin;
301:   jac = (PC_BJacobi *)pc->data;

303:   jac->n_local = blocks;
304:   if (!lens) jac->l_lens = NULL;
305:   else {
306:     PetscCall(PetscMalloc1(blocks, &jac->l_lens));
307:     PetscCall(PetscArraycpy(jac->l_lens, lens, blocks));
308:   }
309:   PetscFunctionReturn(PETSC_SUCCESS);
310: }

312: static PetscErrorCode PCBJacobiGetLocalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
313: {
314:   PC_BJacobi *jac = (PC_BJacobi *)pc->data;

316:   PetscFunctionBegin;
317:   *blocks = jac->n_local;
318:   if (lens) *lens = jac->l_lens;
319:   PetscFunctionReturn(PETSC_SUCCESS);
320: }

322: /*@C
323:   PCBJacobiGetSubKSP - Gets the local `KSP` contexts for all blocks on
324:   this processor.

326:   Not Collective

328:   Input Parameter:
329: . pc - the preconditioner context

331:   Output Parameters:
332: + n_local     - the number of blocks on this processor, or NULL
333: . first_local - the global number of the first block on this processor, or NULL
334: - ksp         - the array of KSP contexts

336:   Level: advanced

338:   Notes:
339:   After `PCBJacobiGetSubKSP()` the array of `KSP` contexts is not to be freed.

341:   Currently for some matrix implementations only 1 block per processor
342:   is supported.

344:   You must call `KSPSetUp()` or `PCSetUp()` before calling `PCBJacobiGetSubKSP()`.

346:   Fortran Note:
347:   Call `PCBJacobiRestoreSubKSP()` when you no longer need access to the array of `KSP`

349: .seealso: [](ch_ksp), `PCBJACOBI`, `PCASM`, `PCASMGetSubKSP()`
350: @*/
351: PetscErrorCode PCBJacobiGetSubKSP(PC pc, PetscInt *n_local, PetscInt *first_local, KSP *ksp[])
352: {
353:   PetscFunctionBegin;
355:   PetscUseMethod(pc, "PCBJacobiGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (pc, n_local, first_local, ksp));
356:   PetscFunctionReturn(PETSC_SUCCESS);
357: }

359: /*@
360:   PCBJacobiSetTotalBlocks - Sets the global number of blocks for the block
361:   Jacobi preconditioner.

363:   Collective

365:   Input Parameters:
366: + pc     - the preconditioner context
367: . blocks - the number of blocks
368: - lens   - [optional] integer array containing the size of each block

370:   Options Database Key:
371: . -pc_bjacobi_blocks <blocks> - Sets the number of global blocks

373:   Level: intermediate

375:   Note:
376:   Currently only a limited number of blocking configurations are supported.
377:   All processors sharing the `PC` must call this routine with the same data.

379: .seealso: [](ch_ksp), `PCBJACOBI`, `PCSetUseAmat()`, `PCBJacobiSetLocalBlocks()`
380: @*/
381: PetscErrorCode PCBJacobiSetTotalBlocks(PC pc, PetscInt blocks, const PetscInt lens[])
382: {
383:   PetscFunctionBegin;
385:   PetscCheck(blocks > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Must have positive blocks");
386:   PetscTryMethod(pc, "PCBJacobiSetTotalBlocks_C", (PC, PetscInt, const PetscInt[]), (pc, blocks, lens));
387:   PetscFunctionReturn(PETSC_SUCCESS);
388: }

390: /*@C
391:   PCBJacobiGetTotalBlocks - Gets the global number of blocks for the block
392:   Jacobi, `PCBJACOBI`, preconditioner.

394:   Not Collective

396:   Input Parameter:
397: . pc - the preconditioner context

399:   Output Parameters:
400: + blocks - the number of blocks
401: - lens   - integer array containing the size of each block

403:   Level: intermediate

405: .seealso: [](ch_ksp), `PCBJACOBI`, `PCSetUseAmat()`, `PCBJacobiGetLocalBlocks()`
406: @*/
407: PetscErrorCode PCBJacobiGetTotalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
408: {
409:   PetscFunctionBegin;
411:   PetscAssertPointer(blocks, 2);
412:   PetscUseMethod(pc, "PCBJacobiGetTotalBlocks_C", (PC, PetscInt *, const PetscInt *[]), (pc, blocks, lens));
413:   PetscFunctionReturn(PETSC_SUCCESS);
414: }

416: /*@
417:   PCBJacobiSetLocalBlocks - Sets the local number of blocks for the block
418:   Jacobi, `PCBJACOBI`,  preconditioner.

420:   Not Collective

422:   Input Parameters:
423: + pc     - the preconditioner context
424: . blocks - the number of blocks
425: - lens   - [optional] integer array containing size of each block

427:   Options Database Key:
428: . -pc_bjacobi_local_blocks <blocks> - Sets the number of local blocks

430:   Level: intermediate

432:   Note:
433:   Currently only a limited number of blocking configurations are supported.

435: .seealso: [](ch_ksp), `PCBJACOBI`, `PCSetUseAmat()`, `PCBJacobiSetTotalBlocks()`
436: @*/
437: PetscErrorCode PCBJacobiSetLocalBlocks(PC pc, PetscInt blocks, const PetscInt lens[])
438: {
439:   PetscFunctionBegin;
441:   PetscCheck(blocks >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have nonegative blocks");
442:   PetscTryMethod(pc, "PCBJacobiSetLocalBlocks_C", (PC, PetscInt, const PetscInt[]), (pc, blocks, lens));
443:   PetscFunctionReturn(PETSC_SUCCESS);
444: }

446: /*@C
447:   PCBJacobiGetLocalBlocks - Gets the local number of blocks for the block
448:   Jacobi, `PCBJACOBI`, preconditioner.

450:   Not Collective

452:   Input Parameters:
453: + pc     - the preconditioner context
454: . blocks - the number of blocks
455: - lens   - [optional] integer array containing size of each block

457:   Level: intermediate

459:   Note:
460:   Currently only a limited number of blocking configurations are supported.

462: .seealso: [](ch_ksp), `PCBJACOBI`, `PCSetUseAmat()`, `PCBJacobiGetTotalBlocks()`
463: @*/
464: PetscErrorCode PCBJacobiGetLocalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
465: {
466:   PetscFunctionBegin;
468:   PetscAssertPointer(blocks, 2);
469:   PetscUseMethod(pc, "PCBJacobiGetLocalBlocks_C", (PC, PetscInt *, const PetscInt *[]), (pc, blocks, lens));
470:   PetscFunctionReturn(PETSC_SUCCESS);
471: }

473: /*MC
474:    PCBJACOBI - Use block Jacobi preconditioning, each block is (approximately) solved with
475:            its own `KSP` object.

477:    Options Database Keys:
478: +  -pc_use_amat - use Amat to apply block of operator in inner Krylov method
479: -  -pc_bjacobi_blocks <n> - use n total blocks

481:    Level: beginner

483:    Notes:
484:     See `PCJACOBI` for diagonal Jacobi, `PCVPBJACOBI` for variable point block, and `PCPBJACOBI` for fixed size point block

486:     Each processor can have one or more blocks, or a single block can be shared by several processes. Defaults to one block per processor.

488:      To set options on the solvers for each block append -sub_ to all the `KSP` and `PC`
489:         options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly

491:      To set the options on the solvers separate for each block call `PCBJacobiGetSubKSP()`
492:          and set the options directly on the resulting `KSP` object (you can access its `PC`
493:          `KSPGetPC()`)

495:      For GPU-based vectors (`VECCUDA`, `VECViennaCL`) it is recommended to use exactly one block per MPI process for best
496:          performance.  Different block partitioning may lead to additional data transfers
497:          between host and GPU that lead to degraded performance.

499:      When multiple processes share a single block, each block encompasses exactly all the unknowns owned its set of processes.

501: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCType`,
502:           `PCASM`, `PCSetUseAmat()`, `PCGetUseAmat()`, `PCBJacobiGetSubKSP()`, `PCBJacobiSetTotalBlocks()`,
503:           `PCBJacobiSetLocalBlocks()`, `PCSetModifySubMatrices()`, `PCJACOBI`, `PCVPBJACOBI`, `PCPBJACOBI`
504: M*/

506: PETSC_EXTERN PetscErrorCode PCCreate_BJacobi(PC pc)
507: {
508:   PetscMPIInt rank;
509:   PC_BJacobi *jac;

511:   PetscFunctionBegin;
512:   PetscCall(PetscNew(&jac));
513:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));

515:   pc->ops->apply             = NULL;
516:   pc->ops->matapply          = NULL;
517:   pc->ops->applytranspose    = NULL;
518:   pc->ops->matapplytranspose = NULL;
519:   pc->ops->setup             = PCSetUp_BJacobi;
520:   pc->ops->destroy           = PCDestroy_BJacobi;
521:   pc->ops->setfromoptions    = PCSetFromOptions_BJacobi;
522:   pc->ops->view              = PCView_BJacobi;
523:   pc->ops->applyrichardson   = NULL;

525:   pc->data         = (void *)jac;
526:   jac->n           = -1;
527:   jac->n_local     = -1;
528:   jac->first_local = rank;
529:   jac->ksp         = NULL;
530:   jac->g_lens      = NULL;
531:   jac->l_lens      = NULL;
532:   jac->psubcomm    = NULL;

534:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetSubKSP_C", PCBJacobiGetSubKSP_BJacobi));
535:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetTotalBlocks_C", PCBJacobiSetTotalBlocks_BJacobi));
536:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetTotalBlocks_C", PCBJacobiGetTotalBlocks_BJacobi));
537:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetLocalBlocks_C", PCBJacobiSetLocalBlocks_BJacobi));
538:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetLocalBlocks_C", PCBJacobiGetLocalBlocks_BJacobi));
539:   PetscFunctionReturn(PETSC_SUCCESS);
540: }

542: /*
543:         These are for a single block per processor; works for AIJ, BAIJ; Seq and MPI
544: */
545: static PetscErrorCode PCReset_BJacobi_Singleblock(PC pc)
546: {
547:   PC_BJacobi             *jac  = (PC_BJacobi *)pc->data;
548:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;

550:   PetscFunctionBegin;
551:   PetscCall(KSPReset(jac->ksp[0]));
552:   PetscCall(VecDestroy(&bjac->x));
553:   PetscCall(VecDestroy(&bjac->y));
554:   PetscFunctionReturn(PETSC_SUCCESS);
555: }

557: static PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
558: {
559:   PC_BJacobi             *jac  = (PC_BJacobi *)pc->data;
560:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;

562:   PetscFunctionBegin;
563:   PetscCall(PCReset_BJacobi_Singleblock(pc));
564:   PetscCall(KSPDestroy(&jac->ksp[0]));
565:   PetscCall(PetscFree(jac->ksp));
566:   PetscCall(PetscFree(bjac));
567:   PetscCall(PCDestroy_BJacobi(pc));
568:   PetscFunctionReturn(PETSC_SUCCESS);
569: }

571: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
572: {
573:   PC_BJacobi        *jac    = (PC_BJacobi *)pc->data;
574:   KSP                subksp = jac->ksp[0];
575:   KSPConvergedReason reason;

577:   PetscFunctionBegin;
578:   PetscCall(KSPSetUp(subksp));
579:   PetscCall(KSPGetConvergedReason(subksp, &reason));
580:   if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
581:   PetscFunctionReturn(PETSC_SUCCESS);
582: }

584: static PetscErrorCode PCApply_BJacobi_Singleblock(PC pc, Vec x, Vec y)
585: {
586:   PC_BJacobi             *jac  = (PC_BJacobi *)pc->data;
587:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;

589:   PetscFunctionBegin;
590:   PetscCall(VecGetLocalVectorRead(x, bjac->x));
591:   PetscCall(VecGetLocalVector(y, bjac->y));
592:   /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
593:      matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
594:      of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
595:   PetscCall(KSPSetReusePreconditioner(jac->ksp[0], pc->reusepreconditioner));
596:   PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
597:   PetscCall(KSPSolve(jac->ksp[0], bjac->x, bjac->y));
598:   PetscCall(KSPCheckSolve(jac->ksp[0], pc, bjac->y));
599:   PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
600:   PetscCall(VecRestoreLocalVectorRead(x, bjac->x));
601:   PetscCall(VecRestoreLocalVector(y, bjac->y));
602:   PetscFunctionReturn(PETSC_SUCCESS);
603: }

605: static PetscErrorCode PCMatApply_BJacobi_Singleblock_Private(PC pc, Mat X, Mat Y, PetscBool transpose)
606: {
607:   PC_BJacobi *jac = (PC_BJacobi *)pc->data;
608:   Mat         sX, sY;

610:   PetscFunctionBegin;
611:   /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
612:      matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
613:      of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
614:   PetscCall(KSPSetReusePreconditioner(jac->ksp[0], pc->reusepreconditioner));
615:   PetscCall(MatDenseGetLocalMatrix(X, &sX));
616:   PetscCall(MatDenseGetLocalMatrix(Y, &sY));
617:   if (!transpose) {
618:     PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], sX, sY, 0));
619:     PetscCall(KSPMatSolve(jac->ksp[0], sX, sY));
620:     PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], sX, sY, 0));
621:   } else {
622:     PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, jac->ksp[0], sX, sY, 0));
623:     PetscCall(KSPMatSolveTranspose(jac->ksp[0], sX, sY));
624:     PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, jac->ksp[0], sX, sY, 0));
625:   }
626:   PetscFunctionReturn(PETSC_SUCCESS);
627: }

629: static PetscErrorCode PCMatApply_BJacobi_Singleblock(PC pc, Mat X, Mat Y)
630: {
631:   PetscFunctionBegin;
632:   PetscCall(PCMatApply_BJacobi_Singleblock_Private(pc, X, Y, PETSC_FALSE));
633:   PetscFunctionReturn(PETSC_SUCCESS);
634: }

636: static PetscErrorCode PCMatApplyTranspose_BJacobi_Singleblock(PC pc, Mat X, Mat Y)
637: {
638:   PetscFunctionBegin;
639:   PetscCall(PCMatApply_BJacobi_Singleblock_Private(pc, X, Y, PETSC_TRUE));
640:   PetscFunctionReturn(PETSC_SUCCESS);
641: }

643: static PetscErrorCode PCApplySymmetricLeft_BJacobi_Singleblock(PC pc, Vec x, Vec y)
644: {
645:   PC_BJacobi             *jac  = (PC_BJacobi *)pc->data;
646:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
647:   PetscScalar            *y_array;
648:   const PetscScalar      *x_array;
649:   PC                      subpc;

651:   PetscFunctionBegin;
652:   /*
653:       The VecPlaceArray() is to avoid having to copy the
654:     y vector into the bjac->x vector. The reason for
655:     the bjac->x vector is that we need a sequential vector
656:     for the sequential solve.
657:   */
658:   PetscCall(VecGetArrayRead(x, &x_array));
659:   PetscCall(VecGetArray(y, &y_array));
660:   PetscCall(VecPlaceArray(bjac->x, x_array));
661:   PetscCall(VecPlaceArray(bjac->y, y_array));
662:   /* apply the symmetric left portion of the inner PC operator */
663:   /* note this bypasses the inner KSP and its options completely */
664:   PetscCall(KSPGetPC(jac->ksp[0], &subpc));
665:   PetscCall(PCApplySymmetricLeft(subpc, bjac->x, bjac->y));
666:   PetscCall(VecResetArray(bjac->x));
667:   PetscCall(VecResetArray(bjac->y));
668:   PetscCall(VecRestoreArrayRead(x, &x_array));
669:   PetscCall(VecRestoreArray(y, &y_array));
670:   PetscFunctionReturn(PETSC_SUCCESS);
671: }

673: static PetscErrorCode PCApplySymmetricRight_BJacobi_Singleblock(PC pc, Vec x, Vec y)
674: {
675:   PC_BJacobi             *jac  = (PC_BJacobi *)pc->data;
676:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
677:   PetscScalar            *y_array;
678:   const PetscScalar      *x_array;
679:   PC                      subpc;

681:   PetscFunctionBegin;
682:   /*
683:       The VecPlaceArray() is to avoid having to copy the
684:     y vector into the bjac->x vector. The reason for
685:     the bjac->x vector is that we need a sequential vector
686:     for the sequential solve.
687:   */
688:   PetscCall(VecGetArrayRead(x, &x_array));
689:   PetscCall(VecGetArray(y, &y_array));
690:   PetscCall(VecPlaceArray(bjac->x, x_array));
691:   PetscCall(VecPlaceArray(bjac->y, y_array));

693:   /* apply the symmetric right portion of the inner PC operator */
694:   /* note this bypasses the inner KSP and its options completely */

696:   PetscCall(KSPGetPC(jac->ksp[0], &subpc));
697:   PetscCall(PCApplySymmetricRight(subpc, bjac->x, bjac->y));

699:   PetscCall(VecResetArray(bjac->x));
700:   PetscCall(VecResetArray(bjac->y));
701:   PetscCall(VecRestoreArrayRead(x, &x_array));
702:   PetscCall(VecRestoreArray(y, &y_array));
703:   PetscFunctionReturn(PETSC_SUCCESS);
704: }

706: static PetscErrorCode PCApplyTranspose_BJacobi_Singleblock(PC pc, Vec x, Vec y)
707: {
708:   PC_BJacobi             *jac  = (PC_BJacobi *)pc->data;
709:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
710:   PetscScalar            *y_array;
711:   const PetscScalar      *x_array;

713:   PetscFunctionBegin;
714:   /*
715:       The VecPlaceArray() is to avoid having to copy the
716:     y vector into the bjac->x vector. The reason for
717:     the bjac->x vector is that we need a sequential vector
718:     for the sequential solve.
719:   */
720:   PetscCall(VecGetArrayRead(x, &x_array));
721:   PetscCall(VecGetArray(y, &y_array));
722:   PetscCall(VecPlaceArray(bjac->x, x_array));
723:   PetscCall(VecPlaceArray(bjac->y, y_array));
724:   PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
725:   PetscCall(KSPSolveTranspose(jac->ksp[0], bjac->x, bjac->y));
726:   PetscCall(KSPCheckSolve(jac->ksp[0], pc, bjac->y));
727:   PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
728:   PetscCall(VecResetArray(bjac->x));
729:   PetscCall(VecResetArray(bjac->y));
730:   PetscCall(VecRestoreArrayRead(x, &x_array));
731:   PetscCall(VecRestoreArray(y, &y_array));
732:   PetscFunctionReturn(PETSC_SUCCESS);
733: }

735: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC pc, Mat mat, Mat pmat)
736: {
737:   PC_BJacobi             *jac = (PC_BJacobi *)pc->data;
738:   PetscInt                m;
739:   KSP                     ksp;
740:   PC_BJacobi_Singleblock *bjac;
741:   PetscBool               wasSetup = PETSC_TRUE;
742:   VecType                 vectype;
743:   const char             *prefix;

745:   PetscFunctionBegin;
746:   if (!pc->setupcalled) {
747:     if (!jac->ksp) {
748:       PetscInt nestlevel;

750:       wasSetup = PETSC_FALSE;

752:       PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
753:       PetscCall(PCGetKSPNestLevel(pc, &nestlevel));
754:       PetscCall(KSPSetNestLevel(ksp, nestlevel + 1));
755:       PetscCall(KSPSetErrorIfNotConverged(ksp, pc->erroriffailure));
756:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1));
757:       PetscCall(KSPSetType(ksp, KSPPREONLY));
758:       PetscCall(PCGetOptionsPrefix(pc, &prefix));
759:       PetscCall(KSPSetOptionsPrefix(ksp, prefix));
760:       PetscCall(KSPAppendOptionsPrefix(ksp, "sub_"));

762:       pc->ops->reset               = PCReset_BJacobi_Singleblock;
763:       pc->ops->destroy             = PCDestroy_BJacobi_Singleblock;
764:       pc->ops->apply               = PCApply_BJacobi_Singleblock;
765:       pc->ops->matapply            = PCMatApply_BJacobi_Singleblock;
766:       pc->ops->matapplytranspose   = PCMatApplyTranspose_BJacobi_Singleblock;
767:       pc->ops->applysymmetricleft  = PCApplySymmetricLeft_BJacobi_Singleblock;
768:       pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Singleblock;
769:       pc->ops->applytranspose      = PCApplyTranspose_BJacobi_Singleblock;
770:       pc->ops->setuponblocks       = PCSetUpOnBlocks_BJacobi_Singleblock;

772:       PetscCall(PetscMalloc1(1, &jac->ksp));
773:       jac->ksp[0] = ksp;

775:       PetscCall(PetscNew(&bjac));
776:       jac->data = (void *)bjac;
777:     } else {
778:       ksp  = jac->ksp[0];
779:       bjac = (PC_BJacobi_Singleblock *)jac->data;
780:     }

782:     /*
783:       The reason we need to generate these vectors is to serve
784:       as the right-hand side and solution vector for the solve on the
785:       block. We do not need to allocate space for the vectors since
786:       that is provided via VecPlaceArray() just before the call to
787:       KSPSolve() on the block.
788:     */
789:     PetscCall(MatGetSize(pmat, &m, &m));
790:     PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, m, NULL, &bjac->x));
791:     PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, m, NULL, &bjac->y));
792:     PetscCall(MatGetVecType(pmat, &vectype));
793:     PetscCall(VecSetType(bjac->x, vectype));
794:     PetscCall(VecSetType(bjac->y, vectype));
795:   } else {
796:     ksp  = jac->ksp[0];
797:     bjac = (PC_BJacobi_Singleblock *)jac->data;
798:   }
799:   PetscCall(KSPGetOptionsPrefix(ksp, &prefix));
800:   if (pc->useAmat) {
801:     PetscCall(KSPSetOperators(ksp, mat, pmat));
802:     PetscCall(MatSetOptionsPrefix(mat, prefix));
803:   } else {
804:     PetscCall(KSPSetOperators(ksp, pmat, pmat));
805:   }
806:   PetscCall(MatSetOptionsPrefix(pmat, prefix));
807:   if (!wasSetup && pc->setfromoptionscalled) {
808:     /* If PCSetFromOptions_BJacobi is called later, KSPSetFromOptions will be called at that time. */
809:     PetscCall(KSPSetFromOptions(ksp));
810:   }
811:   PetscFunctionReturn(PETSC_SUCCESS);
812: }

814: static PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
815: {
816:   PC_BJacobi            *jac  = (PC_BJacobi *)pc->data;
817:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
818:   PetscInt               i;

820:   PetscFunctionBegin;
821:   if (bjac && bjac->pmat) {
822:     PetscCall(MatDestroyMatrices(jac->n_local, &bjac->pmat));
823:     if (pc->useAmat) PetscCall(MatDestroyMatrices(jac->n_local, &bjac->mat));
824:   }

826:   for (i = 0; i < jac->n_local; i++) {
827:     PetscCall(KSPReset(jac->ksp[i]));
828:     if (bjac && bjac->x) {
829:       PetscCall(VecDestroy(&bjac->x[i]));
830:       PetscCall(VecDestroy(&bjac->y[i]));
831:       PetscCall(ISDestroy(&bjac->is[i]));
832:     }
833:   }
834:   PetscCall(PetscFree(jac->l_lens));
835:   PetscCall(PetscFree(jac->g_lens));
836:   PetscFunctionReturn(PETSC_SUCCESS);
837: }

839: static PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
840: {
841:   PC_BJacobi            *jac  = (PC_BJacobi *)pc->data;
842:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
843:   PetscInt               i;

845:   PetscFunctionBegin;
846:   PetscCall(PCReset_BJacobi_Multiblock(pc));
847:   if (bjac) {
848:     PetscCall(PetscFree2(bjac->x, bjac->y));
849:     PetscCall(PetscFree(bjac->starts));
850:     PetscCall(PetscFree(bjac->is));
851:   }
852:   PetscCall(PetscFree(jac->data));
853:   for (i = 0; i < jac->n_local; i++) PetscCall(KSPDestroy(&jac->ksp[i]));
854:   PetscCall(PetscFree(jac->ksp));
855:   PetscCall(PCDestroy_BJacobi(pc));
856:   PetscFunctionReturn(PETSC_SUCCESS);
857: }

859: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
860: {
861:   PC_BJacobi        *jac = (PC_BJacobi *)pc->data;
862:   PetscInt           i, n_local = jac->n_local;
863:   KSPConvergedReason reason;

865:   PetscFunctionBegin;
866:   for (i = 0; i < n_local; i++) {
867:     PetscCall(KSPSetUp(jac->ksp[i]));
868:     PetscCall(KSPGetConvergedReason(jac->ksp[i], &reason));
869:     if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
870:   }
871:   PetscFunctionReturn(PETSC_SUCCESS);
872: }

874: static PetscErrorCode PCApply_BJacobi_Multiblock(PC pc, Vec x, Vec y)
875: {
876:   PC_BJacobi            *jac = (PC_BJacobi *)pc->data;
877:   PetscInt               i, n_local = jac->n_local;
878:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
879:   PetscScalar           *yin;
880:   const PetscScalar     *xin;

882:   PetscFunctionBegin;
883:   PetscCall(VecGetArrayRead(x, &xin));
884:   PetscCall(VecGetArray(y, &yin));
885:   for (i = 0; i < n_local; i++) {
886:     /*
887:        To avoid copying the subvector from x into a workspace we instead
888:        make the workspace vector array point to the subpart of the array of
889:        the global vector.
890:     */
891:     PetscCall(VecPlaceArray(bjac->x[i], xin + bjac->starts[i]));
892:     PetscCall(VecPlaceArray(bjac->y[i], yin + bjac->starts[i]));

894:     PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
895:     PetscCall(KSPSolve(jac->ksp[i], bjac->x[i], bjac->y[i]));
896:     PetscCall(KSPCheckSolve(jac->ksp[i], pc, bjac->y[i]));
897:     PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));

899:     PetscCall(VecResetArray(bjac->x[i]));
900:     PetscCall(VecResetArray(bjac->y[i]));
901:   }
902:   PetscCall(VecRestoreArrayRead(x, &xin));
903:   PetscCall(VecRestoreArray(y, &yin));
904:   PetscFunctionReturn(PETSC_SUCCESS);
905: }

907: static PetscErrorCode PCApplySymmetricLeft_BJacobi_Multiblock(PC pc, Vec x, Vec y)
908: {
909:   PC_BJacobi            *jac = (PC_BJacobi *)pc->data;
910:   PetscInt               i, n_local = jac->n_local;
911:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
912:   PetscScalar           *yin;
913:   const PetscScalar     *xin;
914:   PC                     subpc;

916:   PetscFunctionBegin;
917:   PetscCall(VecGetArrayRead(x, &xin));
918:   PetscCall(VecGetArray(y, &yin));
919:   for (i = 0; i < n_local; i++) {
920:     /*
921:        To avoid copying the subvector from x into a workspace we instead
922:        make the workspace vector array point to the subpart of the array of
923:        the global vector.
924:     */
925:     PetscCall(VecPlaceArray(bjac->x[i], xin + bjac->starts[i]));
926:     PetscCall(VecPlaceArray(bjac->y[i], yin + bjac->starts[i]));

928:     PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
929:     /* apply the symmetric left portion of the inner PC operator */
930:     /* note this bypasses the inner KSP and its options completely */
931:     PetscCall(KSPGetPC(jac->ksp[i], &subpc));
932:     PetscCall(PCApplySymmetricLeft(subpc, bjac->x[i], bjac->y[i]));
933:     PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));

935:     PetscCall(VecResetArray(bjac->x[i]));
936:     PetscCall(VecResetArray(bjac->y[i]));
937:   }
938:   PetscCall(VecRestoreArrayRead(x, &xin));
939:   PetscCall(VecRestoreArray(y, &yin));
940:   PetscFunctionReturn(PETSC_SUCCESS);
941: }

943: static PetscErrorCode PCApplySymmetricRight_BJacobi_Multiblock(PC pc, Vec x, Vec y)
944: {
945:   PC_BJacobi            *jac = (PC_BJacobi *)pc->data;
946:   PetscInt               i, n_local = jac->n_local;
947:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
948:   PetscScalar           *yin;
949:   const PetscScalar     *xin;
950:   PC                     subpc;

952:   PetscFunctionBegin;
953:   PetscCall(VecGetArrayRead(x, &xin));
954:   PetscCall(VecGetArray(y, &yin));
955:   for (i = 0; i < n_local; i++) {
956:     /*
957:        To avoid copying the subvector from x into a workspace we instead
958:        make the workspace vector array point to the subpart of the array of
959:        the global vector.
960:     */
961:     PetscCall(VecPlaceArray(bjac->x[i], xin + bjac->starts[i]));
962:     PetscCall(VecPlaceArray(bjac->y[i], yin + bjac->starts[i]));

964:     PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
965:     /* apply the symmetric left portion of the inner PC operator */
966:     /* note this bypasses the inner KSP and its options completely */
967:     PetscCall(KSPGetPC(jac->ksp[i], &subpc));
968:     PetscCall(PCApplySymmetricRight(subpc, bjac->x[i], bjac->y[i]));
969:     PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));

971:     PetscCall(VecResetArray(bjac->x[i]));
972:     PetscCall(VecResetArray(bjac->y[i]));
973:   }
974:   PetscCall(VecRestoreArrayRead(x, &xin));
975:   PetscCall(VecRestoreArray(y, &yin));
976:   PetscFunctionReturn(PETSC_SUCCESS);
977: }

979: static PetscErrorCode PCApplyTranspose_BJacobi_Multiblock(PC pc, Vec x, Vec y)
980: {
981:   PC_BJacobi            *jac = (PC_BJacobi *)pc->data;
982:   PetscInt               i, n_local = jac->n_local;
983:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
984:   PetscScalar           *yin;
985:   const PetscScalar     *xin;

987:   PetscFunctionBegin;
988:   PetscCall(VecGetArrayRead(x, &xin));
989:   PetscCall(VecGetArray(y, &yin));
990:   for (i = 0; i < n_local; i++) {
991:     /*
992:        To avoid copying the subvector from x into a workspace we instead
993:        make the workspace vector array point to the subpart of the array of
994:        the global vector.
995:     */
996:     PetscCall(VecPlaceArray(bjac->x[i], xin + bjac->starts[i]));
997:     PetscCall(VecPlaceArray(bjac->y[i], yin + bjac->starts[i]));

999:     PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
1000:     PetscCall(KSPSolveTranspose(jac->ksp[i], bjac->x[i], bjac->y[i]));
1001:     PetscCall(KSPCheckSolve(jac->ksp[i], pc, bjac->y[i]));
1002:     PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));

1004:     PetscCall(VecResetArray(bjac->x[i]));
1005:     PetscCall(VecResetArray(bjac->y[i]));
1006:   }
1007:   PetscCall(VecRestoreArrayRead(x, &xin));
1008:   PetscCall(VecRestoreArray(y, &yin));
1009:   PetscFunctionReturn(PETSC_SUCCESS);
1010: }

1012: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC pc, Mat mat, Mat pmat)
1013: {
1014:   PC_BJacobi            *jac = (PC_BJacobi *)pc->data;
1015:   PetscInt               m, n_local, N, M, start, i;
1016:   const char            *prefix;
1017:   KSP                    ksp;
1018:   Vec                    x, y;
1019:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
1020:   PC                     subpc;
1021:   IS                     is;
1022:   MatReuse               scall;
1023:   VecType                vectype;
1024:   MatNullSpace          *nullsp_mat = NULL, *nullsp_pmat = NULL;

1026:   PetscFunctionBegin;
1027:   PetscCall(MatGetLocalSize(pc->pmat, &M, &N));

1029:   n_local = jac->n_local;

1031:   if (pc->useAmat) {
1032:     PetscBool same;
1033:     PetscCall(PetscObjectTypeCompare((PetscObject)mat, ((PetscObject)pmat)->type_name, &same));
1034:     PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Matrices not of same type");
1035:   }

1037:   if (!pc->setupcalled) {
1038:     PetscInt nestlevel;

1040:     scall = MAT_INITIAL_MATRIX;

1042:     if (!jac->ksp) {
1043:       pc->ops->reset               = PCReset_BJacobi_Multiblock;
1044:       pc->ops->destroy             = PCDestroy_BJacobi_Multiblock;
1045:       pc->ops->apply               = PCApply_BJacobi_Multiblock;
1046:       pc->ops->matapply            = NULL;
1047:       pc->ops->matapplytranspose   = NULL;
1048:       pc->ops->applysymmetricleft  = PCApplySymmetricLeft_BJacobi_Multiblock;
1049:       pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Multiblock;
1050:       pc->ops->applytranspose      = PCApplyTranspose_BJacobi_Multiblock;
1051:       pc->ops->setuponblocks       = PCSetUpOnBlocks_BJacobi_Multiblock;

1053:       PetscCall(PetscNew(&bjac));
1054:       PetscCall(PetscMalloc1(n_local, &jac->ksp));
1055:       PetscCall(PetscMalloc2(n_local, &bjac->x, n_local, &bjac->y));
1056:       PetscCall(PetscMalloc1(n_local, &bjac->starts));

1058:       jac->data = (void *)bjac;
1059:       PetscCall(PetscMalloc1(n_local, &bjac->is));

1061:       for (i = 0; i < n_local; i++) {
1062:         PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
1063:         PetscCall(PCGetKSPNestLevel(pc, &nestlevel));
1064:         PetscCall(KSPSetNestLevel(ksp, nestlevel + 1));
1065:         PetscCall(KSPSetErrorIfNotConverged(ksp, pc->erroriffailure));
1066:         PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1));
1067:         PetscCall(KSPSetType(ksp, KSPPREONLY));
1068:         PetscCall(KSPGetPC(ksp, &subpc));
1069:         PetscCall(PCGetOptionsPrefix(pc, &prefix));
1070:         PetscCall(KSPSetOptionsPrefix(ksp, prefix));
1071:         PetscCall(KSPAppendOptionsPrefix(ksp, "sub_"));

1073:         jac->ksp[i] = ksp;
1074:       }
1075:     } else {
1076:       bjac = (PC_BJacobi_Multiblock *)jac->data;
1077:     }

1079:     start = 0;
1080:     PetscCall(MatGetVecType(pmat, &vectype));
1081:     for (i = 0; i < n_local; i++) {
1082:       m = jac->l_lens[i];
1083:       /*
1084:       The reason we need to generate these vectors is to serve
1085:       as the right-hand side and solution vector for the solve on the
1086:       block. We do not need to allocate space for the vectors since
1087:       that is provided via VecPlaceArray() just before the call to
1088:       KSPSolve() on the block.

1090:       */
1091:       PetscCall(VecCreateSeq(PETSC_COMM_SELF, m, &x));
1092:       PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, m, NULL, &y));
1093:       PetscCall(VecSetType(x, vectype));
1094:       PetscCall(VecSetType(y, vectype));

1096:       bjac->x[i]      = x;
1097:       bjac->y[i]      = y;
1098:       bjac->starts[i] = start;

1100:       PetscCall(ISCreateStride(PETSC_COMM_SELF, m, start, 1, &is));
1101:       bjac->is[i] = is;

1103:       start += m;
1104:     }
1105:   } else {
1106:     bjac = (PC_BJacobi_Multiblock *)jac->data;
1107:     /*
1108:        Destroy the blocks from the previous iteration
1109:     */
1110:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1111:       PetscCall(MatGetNullSpaces(n_local, bjac->pmat, &nullsp_pmat));
1112:       PetscCall(MatDestroyMatrices(n_local, &bjac->pmat));
1113:       if (pc->useAmat) {
1114:         PetscCall(MatGetNullSpaces(n_local, bjac->mat, &nullsp_mat));
1115:         PetscCall(MatDestroyMatrices(n_local, &bjac->mat));
1116:       }
1117:       scall = MAT_INITIAL_MATRIX;
1118:     } else scall = MAT_REUSE_MATRIX;
1119:   }

1121:   PetscCall(MatCreateSubMatrices(pmat, n_local, bjac->is, bjac->is, scall, &bjac->pmat));
1122:   if (nullsp_pmat) PetscCall(MatRestoreNullSpaces(n_local, bjac->pmat, &nullsp_pmat));
1123:   if (pc->useAmat) {
1124:     PetscCall(MatCreateSubMatrices(mat, n_local, bjac->is, bjac->is, scall, &bjac->mat));
1125:     if (nullsp_mat) PetscCall(MatRestoreNullSpaces(n_local, bjac->mat, &nullsp_mat));
1126:   }
1127:   /* Return control to the user so that the submatrices can be modified (e.g., to apply
1128:      different boundary conditions for the submatrices than for the global problem) */
1129:   PetscCall(PCModifySubMatrices(pc, n_local, bjac->is, bjac->is, bjac->pmat, pc->modifysubmatricesP));

1131:   for (i = 0; i < n_local; i++) {
1132:     PetscCall(KSPGetOptionsPrefix(jac->ksp[i], &prefix));
1133:     if (pc->useAmat) {
1134:       PetscCall(KSPSetOperators(jac->ksp[i], bjac->mat[i], bjac->pmat[i]));
1135:       PetscCall(MatSetOptionsPrefix(bjac->mat[i], prefix));
1136:     } else {
1137:       PetscCall(KSPSetOperators(jac->ksp[i], bjac->pmat[i], bjac->pmat[i]));
1138:     }
1139:     PetscCall(MatSetOptionsPrefix(bjac->pmat[i], prefix));
1140:     if (pc->setfromoptionscalled) PetscCall(KSPSetFromOptions(jac->ksp[i]));
1141:   }
1142:   PetscFunctionReturn(PETSC_SUCCESS);
1143: }

1145: /*
1146:       These are for a single block with multiple processes
1147: */
1148: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiproc(PC pc)
1149: {
1150:   PC_BJacobi        *jac    = (PC_BJacobi *)pc->data;
1151:   KSP                subksp = jac->ksp[0];
1152:   KSPConvergedReason reason;

1154:   PetscFunctionBegin;
1155:   PetscCall(KSPSetUp(subksp));
1156:   PetscCall(KSPGetConvergedReason(subksp, &reason));
1157:   if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
1158:   PetscFunctionReturn(PETSC_SUCCESS);
1159: }

1161: static PetscErrorCode PCReset_BJacobi_Multiproc(PC pc)
1162: {
1163:   PC_BJacobi           *jac   = (PC_BJacobi *)pc->data;
1164:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;

1166:   PetscFunctionBegin;
1167:   PetscCall(VecDestroy(&mpjac->ysub));
1168:   PetscCall(VecDestroy(&mpjac->xsub));
1169:   PetscCall(MatDestroy(&mpjac->submats));
1170:   if (jac->ksp) PetscCall(KSPReset(jac->ksp[0]));
1171:   PetscFunctionReturn(PETSC_SUCCESS);
1172: }

1174: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
1175: {
1176:   PC_BJacobi           *jac   = (PC_BJacobi *)pc->data;
1177:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;

1179:   PetscFunctionBegin;
1180:   PetscCall(PCReset_BJacobi_Multiproc(pc));
1181:   PetscCall(KSPDestroy(&jac->ksp[0]));
1182:   PetscCall(PetscFree(jac->ksp));
1183:   PetscCall(PetscSubcommDestroy(&mpjac->psubcomm));

1185:   PetscCall(PetscFree(mpjac));
1186:   PetscCall(PCDestroy_BJacobi(pc));
1187:   PetscFunctionReturn(PETSC_SUCCESS);
1188: }

1190: static PetscErrorCode PCApply_BJacobi_Multiproc(PC pc, Vec x, Vec y)
1191: {
1192:   PC_BJacobi           *jac   = (PC_BJacobi *)pc->data;
1193:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
1194:   PetscScalar          *yarray;
1195:   const PetscScalar    *xarray;
1196:   KSPConvergedReason    reason;

1198:   PetscFunctionBegin;
1199:   /* place x's and y's local arrays into xsub and ysub */
1200:   PetscCall(VecGetArrayRead(x, &xarray));
1201:   PetscCall(VecGetArray(y, &yarray));
1202:   PetscCall(VecPlaceArray(mpjac->xsub, xarray));
1203:   PetscCall(VecPlaceArray(mpjac->ysub, yarray));

1205:   /* apply preconditioner on each matrix block */
1206:   PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], mpjac->xsub, mpjac->ysub, 0));
1207:   PetscCall(KSPSolve(jac->ksp[0], mpjac->xsub, mpjac->ysub));
1208:   PetscCall(KSPCheckSolve(jac->ksp[0], pc, mpjac->ysub));
1209:   PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], mpjac->xsub, mpjac->ysub, 0));
1210:   PetscCall(KSPGetConvergedReason(jac->ksp[0], &reason));
1211:   if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;

1213:   PetscCall(VecResetArray(mpjac->xsub));
1214:   PetscCall(VecResetArray(mpjac->ysub));
1215:   PetscCall(VecRestoreArrayRead(x, &xarray));
1216:   PetscCall(VecRestoreArray(y, &yarray));
1217:   PetscFunctionReturn(PETSC_SUCCESS);
1218: }

1220: static PetscErrorCode PCMatApply_BJacobi_Multiproc(PC pc, Mat X, Mat Y)
1221: {
1222:   PC_BJacobi        *jac = (PC_BJacobi *)pc->data;
1223:   KSPConvergedReason reason;
1224:   Mat                sX, sY;
1225:   const PetscScalar *x;
1226:   PetscScalar       *y;
1227:   PetscInt           m, N, lda, ldb;

1229:   PetscFunctionBegin;
1230:   /* apply preconditioner on each matrix block */
1231:   PetscCall(MatGetLocalSize(X, &m, NULL));
1232:   PetscCall(MatGetSize(X, NULL, &N));
1233:   PetscCall(MatDenseGetLDA(X, &lda));
1234:   PetscCall(MatDenseGetLDA(Y, &ldb));
1235:   PetscCall(MatDenseGetArrayRead(X, &x));
1236:   PetscCall(MatDenseGetArrayWrite(Y, &y));
1237:   PetscCall(MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]), m, PETSC_DECIDE, PETSC_DECIDE, N, (PetscScalar *)x, &sX));
1238:   PetscCall(MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]), m, PETSC_DECIDE, PETSC_DECIDE, N, y, &sY));
1239:   PetscCall(MatDenseSetLDA(sX, lda));
1240:   PetscCall(MatDenseSetLDA(sY, ldb));
1241:   PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], X, Y, 0));
1242:   PetscCall(KSPMatSolve(jac->ksp[0], sX, sY));
1243:   PetscCall(KSPCheckSolve(jac->ksp[0], pc, NULL));
1244:   PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], X, Y, 0));
1245:   PetscCall(MatDestroy(&sY));
1246:   PetscCall(MatDestroy(&sX));
1247:   PetscCall(MatDenseRestoreArrayWrite(Y, &y));
1248:   PetscCall(MatDenseRestoreArrayRead(X, &x));
1249:   PetscCall(KSPGetConvergedReason(jac->ksp[0], &reason));
1250:   if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
1251:   PetscFunctionReturn(PETSC_SUCCESS);
1252: }

1254: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1255: {
1256:   PC_BJacobi           *jac   = (PC_BJacobi *)pc->data;
1257:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
1258:   PetscInt              m, n;
1259:   MPI_Comm              comm, subcomm = 0;
1260:   const char           *prefix;
1261:   PetscBool             wasSetup = PETSC_TRUE;
1262:   VecType               vectype;

1264:   PetscFunctionBegin;
1265:   PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
1266:   PetscCheck(jac->n_local <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Only a single block in a subcommunicator is supported");
1267:   jac->n_local = 1; /* currently only a single block is supported for a subcommunicator */
1268:   if (!pc->setupcalled) {
1269:     PetscInt nestlevel;

1271:     wasSetup = PETSC_FALSE;
1272:     PetscCall(PetscNew(&mpjac));
1273:     jac->data = (void *)mpjac;

1275:     /* initialize datastructure mpjac */
1276:     if (!jac->psubcomm) {
1277:       /* Create default contiguous subcommunicatiors if user does not provide them */
1278:       PetscCall(PetscSubcommCreate(comm, &jac->psubcomm));
1279:       PetscCall(PetscSubcommSetNumber(jac->psubcomm, jac->n));
1280:       PetscCall(PetscSubcommSetType(jac->psubcomm, PETSC_SUBCOMM_CONTIGUOUS));
1281:     }
1282:     mpjac->psubcomm = jac->psubcomm;
1283:     subcomm         = PetscSubcommChild(mpjac->psubcomm);

1285:     /* Get matrix blocks of pmat */
1286:     PetscCall(MatGetMultiProcBlock(pc->pmat, subcomm, MAT_INITIAL_MATRIX, &mpjac->submats));

1288:     /* create a new PC that processors in each subcomm have copy of */
1289:     PetscCall(PetscMalloc1(1, &jac->ksp));
1290:     PetscCall(KSPCreate(subcomm, &jac->ksp[0]));
1291:     PetscCall(PCGetKSPNestLevel(pc, &nestlevel));
1292:     PetscCall(KSPSetNestLevel(jac->ksp[0], nestlevel + 1));
1293:     PetscCall(KSPSetErrorIfNotConverged(jac->ksp[0], pc->erroriffailure));
1294:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->ksp[0], (PetscObject)pc, 1));
1295:     PetscCall(KSPSetOperators(jac->ksp[0], mpjac->submats, mpjac->submats));
1296:     PetscCall(KSPGetPC(jac->ksp[0], &mpjac->pc));

1298:     PetscCall(PCGetOptionsPrefix(pc, &prefix));
1299:     PetscCall(KSPSetOptionsPrefix(jac->ksp[0], prefix));
1300:     PetscCall(KSPAppendOptionsPrefix(jac->ksp[0], "sub_"));
1301:     PetscCall(KSPGetOptionsPrefix(jac->ksp[0], &prefix));
1302:     PetscCall(MatSetOptionsPrefix(mpjac->submats, prefix));

1304:     /* create dummy vectors xsub and ysub */
1305:     PetscCall(MatGetLocalSize(mpjac->submats, &m, &n));
1306:     PetscCall(VecCreateMPIWithArray(subcomm, 1, n, PETSC_DECIDE, NULL, &mpjac->xsub));
1307:     PetscCall(VecCreateMPIWithArray(subcomm, 1, m, PETSC_DECIDE, NULL, &mpjac->ysub));
1308:     PetscCall(MatGetVecType(mpjac->submats, &vectype));
1309:     PetscCall(VecSetType(mpjac->xsub, vectype));
1310:     PetscCall(VecSetType(mpjac->ysub, vectype));

1312:     pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiproc;
1313:     pc->ops->reset         = PCReset_BJacobi_Multiproc;
1314:     pc->ops->destroy       = PCDestroy_BJacobi_Multiproc;
1315:     pc->ops->apply         = PCApply_BJacobi_Multiproc;
1316:     pc->ops->matapply      = PCMatApply_BJacobi_Multiproc;
1317:   } else { /* pc->setupcalled */
1318:     subcomm = PetscSubcommChild(mpjac->psubcomm);
1319:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1320:       /* destroy old matrix blocks, then get new matrix blocks */
1321:       if (mpjac->submats) PetscCall(MatDestroy(&mpjac->submats));
1322:       PetscCall(MatGetMultiProcBlock(pc->pmat, subcomm, MAT_INITIAL_MATRIX, &mpjac->submats));
1323:     } else {
1324:       PetscCall(MatGetMultiProcBlock(pc->pmat, subcomm, MAT_REUSE_MATRIX, &mpjac->submats));
1325:     }
1326:     PetscCall(KSPSetOperators(jac->ksp[0], mpjac->submats, mpjac->submats));
1327:   }

1329:   if (!wasSetup && pc->setfromoptionscalled) PetscCall(KSPSetFromOptions(jac->ksp[0]));
1330:   PetscFunctionReturn(PETSC_SUCCESS);
1331: }