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) PetscCall(PCSetUp_BJacobi_Singleblock(pc, mat, pmat));
    127:   else PetscCall(PCSetUp_BJacobi_Multiblock(pc, mat, pmat));
    128:   PetscFunctionReturn(PETSC_SUCCESS);
    129: }

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

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

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

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

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

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

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

    255: static PetscErrorCode PCBJacobiGetSubKSP_BJacobi(PC pc, PetscInt *n_local, PetscInt *first_local, KSP **ksp)
    256: {
    257:   PC_BJacobi *jac = (PC_BJacobi *)pc->data;

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

    262:   if (n_local) *n_local = jac->n_local;
    263:   if (first_local) *first_local = jac->first_local;
    264:   if (ksp) *ksp = jac->ksp;
    265:   PetscFunctionReturn(PETSC_SUCCESS);
    266: }

    268: static PetscErrorCode PCBJacobiSetTotalBlocks_BJacobi(PC pc, PetscInt blocks, const PetscInt *lens)
    269: {
    270:   PC_BJacobi *jac = (PC_BJacobi *)pc->data;

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

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

    287:   PetscFunctionBegin;
    288:   *blocks = jac->n;
    289:   if (lens) *lens = jac->g_lens;
    290:   PetscFunctionReturn(PETSC_SUCCESS);
    291: }

    293: static PetscErrorCode PCBJacobiSetLocalBlocks_BJacobi(PC pc, PetscInt blocks, const PetscInt lens[])
    294: {
    295:   PC_BJacobi *jac;

    297:   PetscFunctionBegin;
    298:   jac = (PC_BJacobi *)pc->data;

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

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

    313:   PetscFunctionBegin;
    314:   *blocks = jac->n_local;
    315:   if (lens) *lens = jac->l_lens;
    316:   PetscFunctionReturn(PETSC_SUCCESS);
    317: }

    319: /*@C
    320:   PCBJacobiGetSubKSP - Gets the local `KSP` contexts for all blocks on
    321:   this processor.

    323:   Not Collective

    325:   Input Parameter:
    326: . pc - the preconditioner context

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

    333:   Level: advanced

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

    338:   Currently for some matrix implementations only 1 block per processor
    339:   is supported.

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

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

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

    356: /*@
    357:   PCBJacobiSetTotalBlocks - Sets the global number of blocks for the block
    358:   Jacobi preconditioner.

    360:   Collective

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

    367:   Options Database Key:
    368: . -pc_bjacobi_blocks blocks - Sets the number of global blocks

    370:   Level: intermediate

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

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

    387: /*@C
    388:   PCBJacobiGetTotalBlocks - Gets the global number of blocks for the block
    389:   Jacobi, `PCBJACOBI`, preconditioner.

    391:   Not Collective

    393:   Input Parameter:
    394: . pc - the preconditioner context

    396:   Output Parameters:
    397: + blocks - the number of blocks
    398: - lens   - integer array containing the size of each block

    400:   Level: intermediate

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

    413: /*@
    414:   PCBJacobiSetLocalBlocks - Sets the local number of blocks for the block
    415:   Jacobi, `PCBJACOBI`,  preconditioner.

    417:   Not Collective

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

    424:   Options Database Key:
    425: . -pc_bjacobi_local_blocks blocks - Sets the number of local blocks

    427:   Level: intermediate

    429:   Note:
    430:   Currently only a limited number of blocking configurations are supported.

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

    443: /*@C
    444:   PCBJacobiGetLocalBlocks - Gets the local number of blocks for the block
    445:   Jacobi, `PCBJACOBI`, preconditioner.

    447:   Not Collective

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

    454:   Level: intermediate

    456:   Note:
    457:   Currently only a limited number of blocking configurations are supported.

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

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

    473:    Options Database Keys:
    474: +  -pc_use_amat (true|false) - use `Amat` to apply block of operator in inner Krylov method
    475: -  -pc_bjacobi_blocks n      - use n total blocks

    477:    Level: beginner

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

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

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

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

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

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

    497: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
    498:           `PCASM`, `PCSetUseAmat()`, `PCGetUseAmat()`, `PCBJacobiGetSubKSP()`, `PCBJacobiSetTotalBlocks()`,
    499:           `PCBJacobiSetLocalBlocks()`, `PCSetModifySubMatrices()`, `PCJACOBI`, `PCVPBJACOBI`, `PCPBJACOBI`
    500: M*/

    502: PETSC_EXTERN PetscErrorCode PCCreate_BJacobi(PC pc)
    503: {
    504:   PetscMPIInt rank;
    505:   PC_BJacobi *jac;

    507:   PetscFunctionBegin;
    508:   PetscCall(PetscNew(&jac));
    509:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));

    511:   pc->ops->apply             = NULL;
    512:   pc->ops->matapply          = NULL;
    513:   pc->ops->applytranspose    = NULL;
    514:   pc->ops->matapplytranspose = NULL;
    515:   pc->ops->setup             = PCSetUp_BJacobi;
    516:   pc->ops->destroy           = PCDestroy_BJacobi;
    517:   pc->ops->setfromoptions    = PCSetFromOptions_BJacobi;
    518:   pc->ops->view              = PCView_BJacobi;
    519:   pc->ops->applyrichardson   = NULL;

    521:   pc->data         = (void *)jac;
    522:   jac->n           = -1;
    523:   jac->n_local     = -1;
    524:   jac->first_local = rank;
    525:   jac->ksp         = NULL;
    526:   jac->g_lens      = NULL;
    527:   jac->l_lens      = NULL;
    528:   jac->psubcomm    = NULL;

    530:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetSubKSP_C", PCBJacobiGetSubKSP_BJacobi));
    531:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetTotalBlocks_C", PCBJacobiSetTotalBlocks_BJacobi));
    532:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetTotalBlocks_C", PCBJacobiGetTotalBlocks_BJacobi));
    533:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetLocalBlocks_C", PCBJacobiSetLocalBlocks_BJacobi));
    534:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetLocalBlocks_C", PCBJacobiGetLocalBlocks_BJacobi));
    535:   PetscFunctionReturn(PETSC_SUCCESS);
    536: }

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

    546:   PetscFunctionBegin;
    547:   PetscCall(KSPReset(jac->ksp[0]));
    548:   PetscCall(VecDestroy(&bjac->x));
    549:   PetscCall(VecDestroy(&bjac->y));
    550:   PetscFunctionReturn(PETSC_SUCCESS);
    551: }

    553: static PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
    554: {
    555:   PC_BJacobi             *jac  = (PC_BJacobi *)pc->data;
    556:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;

    558:   PetscFunctionBegin;
    559:   PetscCall(PCReset_BJacobi_Singleblock(pc));
    560:   PetscCall(KSPDestroy(&jac->ksp[0]));
    561:   PetscCall(PetscFree(jac->ksp));
    562:   PetscCall(PetscFree(bjac));
    563:   PetscCall(PCDestroy_BJacobi(pc));
    564:   PetscFunctionReturn(PETSC_SUCCESS);
    565: }

    567: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
    568: {
    569:   PC_BJacobi        *jac    = (PC_BJacobi *)pc->data;
    570:   KSP                subksp = jac->ksp[0];
    571:   KSPConvergedReason reason;

    573:   PetscFunctionBegin;
    574:   PetscCall(KSPSetUp(subksp));
    575:   PetscCall(KSPGetConvergedReason(subksp, &reason));
    576:   if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
    577:   PetscFunctionReturn(PETSC_SUCCESS);
    578: }

    580: static PetscErrorCode PCApply_BJacobi_Singleblock(PC pc, Vec x, Vec y)
    581: {
    582:   PC_BJacobi             *jac  = (PC_BJacobi *)pc->data;
    583:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;

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

    601: static PetscErrorCode PCMatApply_BJacobi_Singleblock_Private(PC pc, Mat X, Mat Y, PetscBool transpose)
    602: {
    603:   PC_BJacobi *jac = (PC_BJacobi *)pc->data;
    604:   Mat         sX, sY;

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

    625: static PetscErrorCode PCMatApply_BJacobi_Singleblock(PC pc, Mat X, Mat Y)
    626: {
    627:   PetscFunctionBegin;
    628:   PetscCall(PCMatApply_BJacobi_Singleblock_Private(pc, X, Y, PETSC_FALSE));
    629:   PetscFunctionReturn(PETSC_SUCCESS);
    630: }

    632: static PetscErrorCode PCMatApplyTranspose_BJacobi_Singleblock(PC pc, Mat X, Mat Y)
    633: {
    634:   PetscFunctionBegin;
    635:   PetscCall(PCMatApply_BJacobi_Singleblock_Private(pc, X, Y, PETSC_TRUE));
    636:   PetscFunctionReturn(PETSC_SUCCESS);
    637: }

    639: static PetscErrorCode PCApplySymmetricLeft_BJacobi_Singleblock(PC pc, Vec x, Vec y)
    640: {
    641:   PC_BJacobi             *jac  = (PC_BJacobi *)pc->data;
    642:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
    643:   PetscScalar            *y_array;
    644:   const PetscScalar      *x_array;
    645:   PC                      subpc;

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

    669: static PetscErrorCode PCApplySymmetricRight_BJacobi_Singleblock(PC pc, Vec x, Vec y)
    670: {
    671:   PC_BJacobi             *jac  = (PC_BJacobi *)pc->data;
    672:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
    673:   PetscScalar            *y_array;
    674:   const PetscScalar      *x_array;
    675:   PC                      subpc;

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

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

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

    695:   PetscCall(VecResetArray(bjac->x));
    696:   PetscCall(VecResetArray(bjac->y));
    697:   PetscCall(VecRestoreArrayRead(x, &x_array));
    698:   PetscCall(VecRestoreArray(y, &y_array));
    699:   PetscFunctionReturn(PETSC_SUCCESS);
    700: }

    702: static PetscErrorCode PCApplyTranspose_BJacobi_Singleblock(PC pc, Vec x, Vec y)
    703: {
    704:   PC_BJacobi             *jac  = (PC_BJacobi *)pc->data;
    705:   PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
    706:   PetscScalar            *y_array;
    707:   const PetscScalar      *x_array;

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

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

    741:   PetscFunctionBegin;
    742:   if (!pc->setupcalled) {
    743:     if (!jac->ksp) {
    744:       PetscInt nestlevel;

    746:       wasSetup = PETSC_FALSE;

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

    758:       pc->ops->reset               = PCReset_BJacobi_Singleblock;
    759:       pc->ops->destroy             = PCDestroy_BJacobi_Singleblock;
    760:       pc->ops->apply               = PCApply_BJacobi_Singleblock;
    761:       pc->ops->matapply            = PCMatApply_BJacobi_Singleblock;
    762:       pc->ops->matapplytranspose   = PCMatApplyTranspose_BJacobi_Singleblock;
    763:       pc->ops->applysymmetricleft  = PCApplySymmetricLeft_BJacobi_Singleblock;
    764:       pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Singleblock;
    765:       pc->ops->applytranspose      = PCApplyTranspose_BJacobi_Singleblock;
    766:       pc->ops->setuponblocks       = PCSetUpOnBlocks_BJacobi_Singleblock;

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

    771:       PetscCall(PetscNew(&bjac));
    772:       jac->data = (void *)bjac;
    773:     } else {
    774:       ksp  = jac->ksp[0];
    775:       bjac = (PC_BJacobi_Singleblock *)jac->data;
    776:     }

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

    810: static PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
    811: {
    812:   PC_BJacobi            *jac  = (PC_BJacobi *)pc->data;
    813:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
    814:   PetscInt               i;

    816:   PetscFunctionBegin;
    817:   if (bjac && bjac->pmat) {
    818:     PetscCall(MatDestroyMatrices(jac->n_local, &bjac->pmat));
    819:     if (pc->useAmat) PetscCall(MatDestroyMatrices(jac->n_local, &bjac->mat));
    820:   }

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

    835: static PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
    836: {
    837:   PC_BJacobi            *jac  = (PC_BJacobi *)pc->data;
    838:   PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
    839:   PetscInt               i;

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

    855: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
    856: {
    857:   PC_BJacobi        *jac = (PC_BJacobi *)pc->data;
    858:   PetscInt           i, n_local = jac->n_local;
    859:   KSPConvergedReason reason;

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

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

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

    890:     PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
    891:     PetscCall(KSPSolve(jac->ksp[i], bjac->x[i], bjac->y[i]));
    892:     PetscCall(KSPCheckSolve(jac->ksp[i], pc, bjac->y[i]));
    893:     PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));

    895:     PetscCall(VecResetArray(bjac->x[i]));
    896:     PetscCall(VecResetArray(bjac->y[i]));
    897:   }
    898:   PetscCall(VecRestoreArrayRead(x, &xin));
    899:   PetscCall(VecRestoreArray(y, &yin));
    900:   PetscFunctionReturn(PETSC_SUCCESS);
    901: }

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

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

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

    931:     PetscCall(VecResetArray(bjac->x[i]));
    932:     PetscCall(VecResetArray(bjac->y[i]));
    933:   }
    934:   PetscCall(VecRestoreArrayRead(x, &xin));
    935:   PetscCall(VecRestoreArray(y, &yin));
    936:   PetscFunctionReturn(PETSC_SUCCESS);
    937: }

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

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

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

    967:     PetscCall(VecResetArray(bjac->x[i]));
    968:     PetscCall(VecResetArray(bjac->y[i]));
    969:   }
    970:   PetscCall(VecRestoreArrayRead(x, &xin));
    971:   PetscCall(VecRestoreArray(y, &yin));
    972:   PetscFunctionReturn(PETSC_SUCCESS);
    973: }

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

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

    995:     PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
    996:     PetscCall(KSPSolveTranspose(jac->ksp[i], bjac->x[i], bjac->y[i]));
    997:     PetscCall(KSPCheckSolve(jac->ksp[i], pc, bjac->y[i]));
    998:     PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));

    1000:     PetscCall(VecResetArray(bjac->x[i]));
    1001:     PetscCall(VecResetArray(bjac->y[i]));
    1002:   }
    1003:   PetscCall(VecRestoreArrayRead(x, &xin));
    1004:   PetscCall(VecRestoreArray(y, &yin));
    1005:   PetscFunctionReturn(PETSC_SUCCESS);
    1006: }

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

    1022:   PetscFunctionBegin;
    1023:   PetscCall(MatGetLocalSize(pc->pmat, &M, &N));

    1025:   n_local = jac->n_local;

    1027:   if (pc->useAmat) {
    1028:     PetscBool same;
    1029:     PetscCall(PetscObjectTypeCompare((PetscObject)mat, ((PetscObject)pmat)->type_name, &same));
    1030:     PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Matrices not of same type");
    1031:   }

    1033:   if (!pc->setupcalled) {
    1034:     PetscInt nestlevel;

    1036:     scall = MAT_INITIAL_MATRIX;

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

    1049:       PetscCall(PetscNew(&bjac));
    1050:       PetscCall(PetscMalloc1(n_local, &jac->ksp));
    1051:       PetscCall(PetscMalloc2(n_local, &bjac->x, n_local, &bjac->y));
    1052:       PetscCall(PetscMalloc1(n_local, &bjac->starts));

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

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

    1069:         jac->ksp[i] = ksp;
    1070:       }
    1071:     } else {
    1072:       bjac = (PC_BJacobi_Multiblock *)jac->data;
    1073:     }

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

    1086:       */
    1087:       PetscCall(VecCreateSeq(PETSC_COMM_SELF, m, &x));
    1088:       PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, m, NULL, &y));
    1089:       PetscCall(VecSetType(x, vectype));
    1090:       PetscCall(VecSetType(y, vectype));

    1092:       bjac->x[i]      = x;
    1093:       bjac->y[i]      = y;
    1094:       bjac->starts[i] = start;

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

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

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

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

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

    1150:   PetscFunctionBegin;
    1151:   PetscCall(KSPSetUp(subksp));
    1152:   PetscCall(KSPGetConvergedReason(subksp, &reason));
    1153:   if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
    1154:   PetscFunctionReturn(PETSC_SUCCESS);
    1155: }

    1157: static PetscErrorCode PCReset_BJacobi_Multiproc(PC pc)
    1158: {
    1159:   PC_BJacobi           *jac   = (PC_BJacobi *)pc->data;
    1160:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;

    1162:   PetscFunctionBegin;
    1163:   PetscCall(VecDestroy(&mpjac->ysub));
    1164:   PetscCall(VecDestroy(&mpjac->xsub));
    1165:   PetscCall(MatDestroy(&mpjac->submats));
    1166:   if (jac->ksp) PetscCall(KSPReset(jac->ksp[0]));
    1167:   PetscFunctionReturn(PETSC_SUCCESS);
    1168: }

    1170: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
    1171: {
    1172:   PC_BJacobi           *jac   = (PC_BJacobi *)pc->data;
    1173:   PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;

    1175:   PetscFunctionBegin;
    1176:   PetscCall(PCReset_BJacobi_Multiproc(pc));
    1177:   PetscCall(KSPDestroy(&jac->ksp[0]));
    1178:   PetscCall(PetscFree(jac->ksp));
    1179:   PetscCall(PetscSubcommDestroy(&mpjac->psubcomm));

    1181:   PetscCall(PetscFree(mpjac));
    1182:   PetscCall(PCDestroy_BJacobi(pc));
    1183:   PetscFunctionReturn(PETSC_SUCCESS);
    1184: }

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

    1194:   PetscFunctionBegin;
    1195:   /* place x's and y's local arrays into xsub and ysub */
    1196:   PetscCall(VecGetArrayRead(x, &xarray));
    1197:   PetscCall(VecGetArray(y, &yarray));
    1198:   PetscCall(VecPlaceArray(mpjac->xsub, xarray));
    1199:   PetscCall(VecPlaceArray(mpjac->ysub, yarray));

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

    1209:   PetscCall(VecResetArray(mpjac->xsub));
    1210:   PetscCall(VecResetArray(mpjac->ysub));
    1211:   PetscCall(VecRestoreArrayRead(x, &xarray));
    1212:   PetscCall(VecRestoreArray(y, &yarray));
    1213:   PetscFunctionReturn(PETSC_SUCCESS);
    1214: }

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

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

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

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

    1267:     wasSetup = PETSC_FALSE;
    1268:     PetscCall(PetscNew(&mpjac));
    1269:     jac->data = (void *)mpjac;

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

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

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

    1294:     PetscCall(PCGetOptionsPrefix(pc, &prefix));
    1295:     PetscCall(KSPSetOptionsPrefix(jac->ksp[0], prefix));
    1296:     PetscCall(KSPAppendOptionsPrefix(jac->ksp[0], "sub_"));
    1297:     PetscCall(KSPGetOptionsPrefix(jac->ksp[0], &prefix));
    1298:     PetscCall(MatSetOptionsPrefix(mpjac->submats, prefix));

    1300:     PetscCall(MatGetLocalSize(mpjac->submats, &m, &n));
    1301:     PetscCall(VecCreateMPIWithArray(subcomm, 1, n, PETSC_DECIDE, NULL, &mpjac->xsub));
    1302:     PetscCall(VecCreateMPIWithArray(subcomm, 1, m, PETSC_DECIDE, NULL, &mpjac->ysub));
    1303:     PetscCall(MatGetVecType(mpjac->submats, &vectype));
    1304:     PetscCall(VecSetType(mpjac->xsub, vectype));
    1305:     PetscCall(VecSetType(mpjac->ysub, vectype));

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

    1324:   if (!wasSetup && pc->setfromoptionscalled) PetscCall(KSPSetFromOptions(jac->ksp[0]));
    1325:   PetscFunctionReturn(PETSC_SUCCESS);
    1326: }