Actual source code: fieldsplit.c

      1: #include <petsc/private/pcimpl.h>
      2: #include <petsc/private/kspimpl.h>
      3: #include <petsc/private/matimpl.h>
      4: #include <petscdm.h>
      5: #include <petscdevice.h>
      6: #if PetscDefined(HAVE_CUDA)
      7: #include <petscdevice_cuda.h>
      8: #endif
      9: #if PetscDefined(HAVE_HIP)
     10: #include <petscdevice_hip.h>
     11: #endif

     13: const char *const PCFieldSplitSchurPreTypes[]  = {"SELF", "SELFP", "A11", "USER", "FULL", "PCFieldSplitSchurPreType", "PC_FIELDSPLIT_SCHUR_PRE_", NULL};
     14: const char *const PCFieldSplitSchurFactTypes[] = {"DIAG", "LOWER", "UPPER", "FULL", "PCFieldSplitSchurFactType", "PC_FIELDSPLIT_SCHUR_FACT_", NULL};

     16: PetscLogEvent KSP_Solve_FS_0, KSP_Solve_FS_1, KSP_Solve_FS_S, KSP_Solve_FS_U, KSP_Solve_FS_L, KSP_Solve_FS_2, KSP_Solve_FS_3, KSP_Solve_FS_4;

     18: typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
     19: struct _PC_FieldSplitLink {
     20:   KSP               ksp;
     21:   Vec               x, y, z;
     22:   Mat               X, Y, Z;
     23:   char             *splitname;
     24:   PetscInt          nfields;
     25:   PetscInt         *fields, *fields_col;
     26:   VecScatter        sctx;
     27:   IS                is, is_col;
     28:   PC_FieldSplitLink next, previous;
     29:   PetscLogEvent     event;

     31:   /* Used only when setting coordinates with PCSetCoordinates */
     32:   PetscInt   dim;
     33:   PetscInt   ndofs;
     34:   PetscReal *coords;
     35: };

     37: typedef struct {
     38:   PCCompositeType type;
     39:   PetscBool       defaultsplit; /* Flag for a system with a set of 'k' scalar fields with the same layout (and bs = k) */
     40:   PetscBool       splitdefined; /* Flag is set after the splits have been defined, to prevent more splits from being added */
     41:   PetscInt        bs;           /* Block size for IS and Mat structures */
     42:   PetscInt        nsplits;      /* Number of field divisions defined */
     43:   Vec            *x, *y, w1, w2;
     44:   Mat            *mat;    /* The diagonal block for each split */
     45:   Mat            *pmat;   /* The preconditioning diagonal block for each split */
     46:   Mat            *Afield; /* The rows of the matrix associated with each split */
     47:   PetscBool       issetup;

     49:   /* Only used when Schur complement preconditioning is used */
     50:   Mat                       B;          /* The (0,1) block */
     51:   Mat                       C;          /* The (1,0) block */
     52:   Mat                       schur;      /* The Schur complement S = A11 - A10 A00^{-1} A01, the KSP here, kspinner, is H_1 in [El08] */
     53:   Mat                       schurp;     /* Assembled approximation to S built by MatSchurComplement to be used as a matrix for constructing the preconditioner when solving with S */
     54:   Mat                       schur_user; /* User-provided matrix for constructing the preconditioner for the Schur complement */
     55:   PCFieldSplitSchurPreType  schurpre;   /* Determines which matrix is used for the Schur complement */
     56:   PCFieldSplitSchurFactType schurfactorization;
     57:   KSP                       kspschur;   /* The solver for S */
     58:   KSP                       kspupper;   /* The solver for A in the upper diagonal part of the factorization (H_2 in [El08]) */
     59:   PetscScalar               schurscale; /* Scaling factor for the Schur complement solution with DIAG factorization */

     61:   /* Only used when Golub-Kahan bidiagonalization preconditioning is used */
     62:   Mat          H;           /* The modified matrix H = A00 + nu*A01*A01'              */
     63:   PetscReal    gkbtol;      /* Stopping tolerance for lower bound estimate            */
     64:   PetscInt     gkbdelay;    /* The delay window for the stopping criterion            */
     65:   PetscReal    gkbnu;       /* Parameter for augmented Lagrangian H = A + nu*A01*A01' */
     66:   PetscInt     gkbmaxit;    /* Maximum number of iterations for outer loop            */
     67:   PetscBool    gkbmonitor;  /* Monitor for gkb iterations and the lower bound error   */
     68:   PetscViewer  gkbviewer;   /* Viewer context for gkbmonitor                          */
     69:   Vec          u, v, d, Hu; /* Work vectors for the GKB algorithm                     */
     70:   PetscScalar *vecz;        /* Contains intermediate values, eg for lower bound       */

     72:   PC_FieldSplitLink head;
     73:   PetscBool         isrestrict;       /* indicates PCFieldSplitRestrictIS() has been last called on this object, hack */
     74:   PetscBool         suboptionsset;    /* Indicates that the KSPSetFromOptions() has been called on the sub-KSPs */
     75:   PetscBool         dm_splits;        /* Whether to use DMCreateFieldDecomposition() whenever possible */
     76:   PetscBool         diag_use_amat;    /* Whether to extract diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
     77:   PetscBool         offdiag_use_amat; /* Whether to extract off-diagonal matrix blocks from Amat, rather than Pmat (weaker than -pc_use_amat) */
     78:   PetscBool         detect;           /* Whether to form 2-way split by finding zero diagonal entries */
     79:   PetscBool         coordinates_set;  /* Whether PCSetCoordinates has been called */
     80: } PC_FieldSplit;

     82: /*
     83:     Note:
     84:     there is no particular reason that pmat, x, and y are stored as arrays in PC_FieldSplit instead of
     85:    inside PC_FieldSplitLink, just historical. If you want to be able to add new fields after already using the
     86:    PC you could change this.
     87: */

     89: /* This helper is so that setting a user-provided matrix is orthogonal to choosing to use it.  This way the
     90: * application-provided FormJacobian can provide this matrix without interfering with the user's (command-line) choices. */
     91: static Mat FieldSplitSchurPre(PC_FieldSplit *jac)
     92: {
     93:   switch (jac->schurpre) {
     94:   case PC_FIELDSPLIT_SCHUR_PRE_SELF:
     95:     return jac->schur;
     96:   case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
     97:     return jac->schurp;
     98:   case PC_FIELDSPLIT_SCHUR_PRE_A11:
     99:     return jac->pmat[1];
    100:   case PC_FIELDSPLIT_SCHUR_PRE_FULL: /* We calculate this and store it in schur_user */
    101:   case PC_FIELDSPLIT_SCHUR_PRE_USER: /* Use a user-provided matrix if it is given, otherwise diagonal block */
    102:   default:
    103:     return jac->schur_user ? jac->schur_user : jac->pmat[1];
    104:   }
    105: }

    107: #include <petscdraw.h>
    108: static PetscErrorCode PCView_FieldSplit(PC pc, PetscViewer viewer)
    109: {
    110:   PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
    111:   PetscBool         isascii, isdraw;
    112:   PetscInt          i, j;
    113:   PC_FieldSplitLink ilink = jac->head;

    115:   PetscFunctionBegin;
    116:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    117:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
    118:   if (isascii) {
    119:     if (jac->bs > 0) {
    120:       PetscCall(PetscViewerASCIIPrintf(viewer, "  FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits, jac->bs));
    121:     } else {
    122:       PetscCall(PetscViewerASCIIPrintf(viewer, "  FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits));
    123:     }
    124:     if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for blocks\n"));
    125:     if (jac->diag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for diagonal blocks\n"));
    126:     if (jac->offdiag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for off-diagonal blocks\n"));
    127:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Solver info for each split is in the following KSP objects:\n"));
    128:     for (i = 0; i < jac->nsplits; i++) {
    129:       if (ilink->fields) {
    130:         PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Fields ", i));
    131:         PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
    132:         for (j = 0; j < ilink->nfields; j++) {
    133:           if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ","));
    134:           PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]));
    135:         }
    136:         PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
    137:         PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
    138:       } else {
    139:         PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Defined by IS\n", i));
    140:       }
    141:       PetscCall(KSPView(ilink->ksp, viewer));
    142:       ilink = ilink->next;
    143:     }
    144:   }

    146:   if (isdraw) {
    147:     PetscDraw draw;
    148:     PetscReal x, y, w, wd;

    150:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
    151:     PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
    152:     w  = 2 * PetscMin(1.0 - x, x);
    153:     wd = w / (jac->nsplits + 1);
    154:     x  = x - wd * (jac->nsplits - 1) / 2.0;
    155:     for (i = 0; i < jac->nsplits; i++) {
    156:       PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
    157:       PetscCall(KSPView(ilink->ksp, viewer));
    158:       PetscCall(PetscDrawPopCurrentPoint(draw));
    159:       x += wd;
    160:       ilink = ilink->next;
    161:     }
    162:   }
    163:   PetscFunctionReturn(PETSC_SUCCESS);
    164: }

    166: static PetscErrorCode PCView_FieldSplit_Schur(PC pc, PetscViewer viewer)
    167: {
    168:   PC_FieldSplit             *jac = (PC_FieldSplit *)pc->data;
    169:   PetscBool                  isascii, isdraw;
    170:   PetscInt                   i, j;
    171:   PC_FieldSplitLink          ilink = jac->head;
    172:   MatSchurComplementAinvType atype;

    174:   PetscFunctionBegin;
    175:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    176:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
    177:   if (isascii) {
    178:     if (jac->bs > 0) {
    179:       PetscCall(PetscViewerASCIIPrintf(viewer, "  FieldSplit with Schur preconditioner, blocksize = %" PetscInt_FMT ", factorization %s\n", jac->bs, PCFieldSplitSchurFactTypes[jac->schurfactorization]));
    180:     } else {
    181:       PetscCall(PetscViewerASCIIPrintf(viewer, "  FieldSplit with Schur preconditioner, factorization %s\n", PCFieldSplitSchurFactTypes[jac->schurfactorization]));
    182:     }
    183:     if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for blocks\n"));
    184:     switch (jac->schurpre) {
    185:     case PC_FIELDSPLIT_SCHUR_PRE_SELF:
    186:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from S itself\n"));
    187:       break;
    188:     case PC_FIELDSPLIT_SCHUR_PRE_SELFP:
    189:       if (jac->schur) {
    190:         PetscCall(MatSchurComplementGetAinvType(jac->schur, &atype));
    191:         PetscCall(PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from Sp, an assembled approximation to S, which uses A00's %sinverse\n", atype == MAT_SCHUR_COMPLEMENT_AINV_DIAG ? "diagonal's " : (atype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG ? "block diagonal's " : (atype == MAT_SCHUR_COMPLEMENT_AINV_FULL ? "full " : "lumped diagonal's "))));
    192:       }
    193:       break;
    194:     case PC_FIELDSPLIT_SCHUR_PRE_A11:
    195:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from A11\n"));
    196:       break;
    197:     case PC_FIELDSPLIT_SCHUR_PRE_FULL:
    198:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from the exact Schur complement\n"));
    199:       break;
    200:     case PC_FIELDSPLIT_SCHUR_PRE_USER:
    201:       if (jac->schur_user) {
    202:         PetscCall(PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from user provided matrix\n"));
    203:       } else {
    204:         PetscCall(PetscViewerASCIIPrintf(viewer, "  Preconditioner for the Schur complement formed from A11\n"));
    205:       }
    206:       break;
    207:     default:
    208:       SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Invalid Schur preconditioning type: %d", jac->schurpre);
    209:     }
    210:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Split info:\n"));
    211:     PetscCall(PetscViewerASCIIPushTab(viewer));
    212:     for (i = 0; i < jac->nsplits; i++) {
    213:       if (ilink->fields) {
    214:         PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Fields ", i));
    215:         PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
    216:         for (j = 0; j < ilink->nfields; j++) {
    217:           if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ","));
    218:           PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]));
    219:         }
    220:         PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
    221:         PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
    222:       } else {
    223:         PetscCall(PetscViewerASCIIPrintf(viewer, "Split number %" PetscInt_FMT " Defined by IS\n", i));
    224:       }
    225:       ilink = ilink->next;
    226:     }
    227:     PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for A00 block\n"));
    228:     PetscCall(PetscViewerASCIIPushTab(viewer));
    229:     if (jac->head) PetscCall(KSPView(jac->head->ksp, viewer));
    230:     else PetscCall(PetscViewerASCIIPrintf(viewer, "  not yet available\n"));
    231:     PetscCall(PetscViewerASCIIPopTab(viewer));
    232:     if (jac->head && jac->kspupper != jac->head->ksp) {
    233:       PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for upper A00 in upper triangular factor\n"));
    234:       PetscCall(PetscViewerASCIIPushTab(viewer));
    235:       if (jac->kspupper) PetscCall(KSPView(jac->kspupper, viewer));
    236:       else PetscCall(PetscViewerASCIIPrintf(viewer, "  not yet available\n"));
    237:       PetscCall(PetscViewerASCIIPopTab(viewer));
    238:     }
    239:     PetscCall(PetscViewerASCIIPrintf(viewer, "KSP solver for S = A11 - A10 inv(A00) A01\n"));
    240:     PetscCall(PetscViewerASCIIPushTab(viewer));
    241:     if (jac->kspschur) {
    242:       PetscCall(KSPView(jac->kspschur, viewer));
    243:     } else {
    244:       PetscCall(PetscViewerASCIIPrintf(viewer, "  not yet available\n"));
    245:     }
    246:     PetscCall(PetscViewerASCIIPopTab(viewer));
    247:     PetscCall(PetscViewerASCIIPopTab(viewer));
    248:   } else if (isdraw && jac->head) {
    249:     PetscDraw draw;
    250:     PetscReal x, y, w, wd, h;
    251:     PetscInt  cnt = 2;
    252:     char      str[32];

    254:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
    255:     PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
    256:     if (jac->kspupper != jac->head->ksp) cnt++;
    257:     w  = 2 * PetscMin(1.0 - x, x);
    258:     wd = w / (cnt + 1);

    260:     PetscCall(PetscSNPrintf(str, 32, "Schur fact. %s", PCFieldSplitSchurFactTypes[jac->schurfactorization]));
    261:     PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
    262:     y -= h;
    263:     if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_USER && !jac->schur_user) {
    264:       PetscCall(PetscSNPrintf(str, 32, "Prec. for Schur from %s", PCFieldSplitSchurPreTypes[PC_FIELDSPLIT_SCHUR_PRE_A11]));
    265:     } else {
    266:       PetscCall(PetscSNPrintf(str, 32, "Prec. for Schur from %s", PCFieldSplitSchurPreTypes[jac->schurpre]));
    267:     }
    268:     PetscCall(PetscDrawStringBoxed(draw, x + wd * (cnt - 1) / 2.0, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
    269:     y -= h;
    270:     x = x - wd * (cnt - 1) / 2.0;

    272:     PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
    273:     PetscCall(KSPView(jac->head->ksp, viewer));
    274:     PetscCall(PetscDrawPopCurrentPoint(draw));
    275:     if (jac->kspupper != jac->head->ksp) {
    276:       x += wd;
    277:       PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
    278:       PetscCall(KSPView(jac->kspupper, viewer));
    279:       PetscCall(PetscDrawPopCurrentPoint(draw));
    280:     }
    281:     x += wd;
    282:     PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
    283:     PetscCall(KSPView(jac->kspschur, viewer));
    284:     PetscCall(PetscDrawPopCurrentPoint(draw));
    285:   }
    286:   PetscFunctionReturn(PETSC_SUCCESS);
    287: }

    289: static PetscErrorCode PCView_FieldSplit_GKB(PC pc, PetscViewer viewer)
    290: {
    291:   PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
    292:   PetscBool         isascii, isdraw;
    293:   PetscInt          i, j;
    294:   PC_FieldSplitLink ilink = jac->head;

    296:   PetscFunctionBegin;
    297:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    298:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
    299:   if (isascii) {
    300:     if (jac->bs > 0) {
    301:       PetscCall(PetscViewerASCIIPrintf(viewer, "  FieldSplit with %s composition: total splits = %" PetscInt_FMT ", blocksize = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits, jac->bs));
    302:     } else {
    303:       PetscCall(PetscViewerASCIIPrintf(viewer, "  FieldSplit with %s composition: total splits = %" PetscInt_FMT "\n", PCCompositeTypes[jac->type], jac->nsplits));
    304:     }
    305:     if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for blocks\n"));
    306:     if (jac->diag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for diagonal blocks\n"));
    307:     if (jac->offdiag_use_amat) PetscCall(PetscViewerASCIIPrintf(viewer, "  using Amat (not Pmat) as operator for off-diagonal blocks\n"));

    309:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Stopping tolerance=%.1e, delay in error estimate=%" PetscInt_FMT ", maximum iterations=%" PetscInt_FMT "\n", (double)jac->gkbtol, jac->gkbdelay, jac->gkbmaxit));
    310:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Solver info for H = A00 + nu*A01*A01' matrix:\n"));
    311:     PetscCall(PetscViewerASCIIPushTab(viewer));

    313:     if (ilink->fields) {
    314:       PetscCall(PetscViewerASCIIPrintf(viewer, "Split number 0 Fields "));
    315:       PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
    316:       for (j = 0; j < ilink->nfields; j++) {
    317:         if (j > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ","));
    318:         PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, ilink->fields[j]));
    319:       }
    320:       PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
    321:       PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
    322:     } else {
    323:       PetscCall(PetscViewerASCIIPrintf(viewer, "Split number 0 Defined by IS\n"));
    324:     }
    325:     PetscCall(KSPView(ilink->ksp, viewer));

    327:     PetscCall(PetscViewerASCIIPopTab(viewer));
    328:   }

    330:   if (isdraw) {
    331:     PetscDraw draw;
    332:     PetscReal x, y, w, wd;

    334:     PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
    335:     PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
    336:     w  = 2 * PetscMin(1.0 - x, x);
    337:     wd = w / (jac->nsplits + 1);
    338:     x  = x - wd * (jac->nsplits - 1) / 2.0;
    339:     for (i = 0; i < jac->nsplits; i++) {
    340:       PetscCall(PetscDrawPushCurrentPoint(draw, x, y));
    341:       PetscCall(KSPView(ilink->ksp, viewer));
    342:       PetscCall(PetscDrawPopCurrentPoint(draw));
    343:       x += wd;
    344:       ilink = ilink->next;
    345:     }
    346:   }
    347:   PetscFunctionReturn(PETSC_SUCCESS);
    348: }

    350: /* Precondition: jac->bs is set to a meaningful value or MATNEST */
    351: static PetscErrorCode PCFieldSplitSetRuntimeSplits_Private(PC pc)
    352: {
    353:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
    354:   PetscInt       bs, i, nfields, *ifields, nfields_col, *ifields_col;
    355:   PetscBool      flg, flg_col, mnest;
    356:   char           optionname[128], splitname[8], optionname_col[128];

    358:   PetscFunctionBegin;
    359:   PetscCall(PetscObjectTypeCompare((PetscObject)pc->mat, MATNEST, &mnest));
    360:   if (mnest) PetscCall(MatNestGetSize(pc->pmat, &bs, NULL));
    361:   else bs = jac->bs;
    362:   PetscCall(PetscMalloc2(bs, &ifields, bs, &ifields_col));
    363:   for (i = 0, flg = PETSC_TRUE;; i++) {
    364:     PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i));
    365:     PetscCall(PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i));
    366:     PetscCall(PetscSNPrintf(optionname_col, sizeof(optionname_col), "-pc_fieldsplit_%" PetscInt_FMT "_fields_col", i));
    367:     nfields     = bs;
    368:     nfields_col = bs;
    369:     PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg));
    370:     PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname_col, ifields_col, &nfields_col, &flg_col));
    371:     if (!flg) break;
    372:     else if (flg && !flg_col) {
    373:       PetscCheck(nfields, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot list zero fields");
    374:       PetscCall(PCFieldSplitSetFields(pc, splitname, nfields, ifields, ifields));
    375:     } else {
    376:       PetscCheck(nfields && nfields_col, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot list zero fields");
    377:       PetscCheck(nfields == nfields_col, PETSC_COMM_SELF, PETSC_ERR_USER, "Number of row and column fields must match");
    378:       PetscCall(PCFieldSplitSetFields(pc, splitname, nfields, ifields, ifields_col));
    379:     }
    380:   }
    381:   if (i > 0) {
    382:     /* Makes command-line setting of splits take precedence over setting them in code.
    383:        Otherwise subsequent calls to PCFieldSplitSetIS() or PCFieldSplitSetFields() would
    384:        create new splits, which would probably not be what the user wanted. */
    385:     jac->splitdefined = PETSC_TRUE;
    386:   }
    387:   PetscCall(PetscFree2(ifields, ifields_col));
    388:   PetscFunctionReturn(PETSC_SUCCESS);
    389: }

    391: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
    392: {
    393:   PC_FieldSplit    *jac                = (PC_FieldSplit *)pc->data;
    394:   PC_FieldSplitLink ilink              = jac->head;
    395:   PetscBool         fieldsplit_default = PETSC_FALSE, coupling = PETSC_FALSE;
    396:   PetscInt          i;

    398:   PetscFunctionBegin;
    399:   /*
    400:    Kinda messy, but at least this now uses DMCreateFieldDecomposition().
    401:    Should probably be rewritten.
    402:    */
    403:   if (!ilink) {
    404:     PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_detect_coupling", &coupling, NULL));
    405:     if (pc->dm && jac->dm_splits && !jac->detect && !coupling) {
    406:       PetscInt  numFields, f, i, j;
    407:       char    **fieldNames;
    408:       IS       *fields;
    409:       DM       *dms;
    410:       DM        subdm[128];
    411:       PetscBool flg;

    413:       PetscCall(DMCreateFieldDecomposition(pc->dm, &numFields, &fieldNames, &fields, &dms));
    414:       /* Allow the user to prescribe the splits */
    415:       for (i = 0, flg = PETSC_TRUE;; i++) {
    416:         PetscInt ifields[128];
    417:         IS       compField;
    418:         char     optionname[128], splitname[8];
    419:         PetscInt nfields = numFields;

    421:         PetscCall(PetscSNPrintf(optionname, sizeof(optionname), "-pc_fieldsplit_%" PetscInt_FMT "_fields", i));
    422:         PetscCall(PetscOptionsGetIntArray(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, optionname, ifields, &nfields, &flg));
    423:         if (!flg) break;
    424:         PetscCheck(numFields <= 128, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot currently support %" PetscInt_FMT " > 128 fields", numFields);
    425:         PetscCall(DMCreateSubDM(pc->dm, nfields, ifields, &compField, &subdm[i]));
    426:         if (nfields == 1) PetscCall(PCFieldSplitSetIS(pc, fieldNames[ifields[0]], compField));
    427:         else {
    428:           PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i));
    429:           PetscCall(PCFieldSplitSetIS(pc, splitname, compField));
    430:         }
    431:         PetscCall(ISDestroy(&compField));
    432:         for (j = 0; j < nfields; ++j) {
    433:           f = ifields[j];
    434:           PetscCall(PetscFree(fieldNames[f]));
    435:           PetscCall(ISDestroy(&fields[f]));
    436:         }
    437:       }
    438:       if (i == 0) {
    439:         for (f = 0; f < numFields; ++f) {
    440:           PetscCall(PCFieldSplitSetIS(pc, fieldNames[f], fields[f]));
    441:           PetscCall(PetscFree(fieldNames[f]));
    442:           PetscCall(ISDestroy(&fields[f]));
    443:         }
    444:       } else {
    445:         for (j = 0; j < numFields; j++) PetscCall(DMDestroy(dms + j));
    446:         PetscCall(PetscFree(dms));
    447:         PetscCall(PetscMalloc1(i, &dms));
    448:         for (j = 0; j < i; ++j) dms[j] = subdm[j];
    449:       }
    450:       PetscCall(PetscFree(fieldNames));
    451:       PetscCall(PetscFree(fields));
    452:       if (dms) {
    453:         PetscCall(PetscInfo(pc, "Setting up physics based fieldsplit preconditioner using the embedded DM\n"));
    454:         for (ilink = jac->head, i = 0; ilink; ilink = ilink->next, ++i) {
    455:           const char *prefix;
    456:           PetscCall(PetscObjectGetOptionsPrefix((PetscObject)ilink->ksp, &prefix));
    457:           PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dms[i], prefix));
    458:           PetscCall(KSPSetDM(ilink->ksp, dms[i]));
    459:           PetscCall(KSPSetDMActive(ilink->ksp, KSP_DMACTIVE_ALL, PETSC_FALSE));
    460:           PetscCall(PetscObjectIncrementTabLevel((PetscObject)dms[i], (PetscObject)ilink->ksp, 0));
    461:           PetscCall(DMDestroy(&dms[i]));
    462:         }
    463:         PetscCall(PetscFree(dms));
    464:       }
    465:     } else {
    466:       if (jac->bs <= 0) {
    467:         if (pc->pmat) PetscCall(MatGetBlockSize(pc->pmat, &jac->bs));
    468:         else jac->bs = 1;
    469:       }

    471:       if (jac->detect) {
    472:         IS       zerodiags, rest;
    473:         PetscInt nmin, nmax;

    475:         PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax));
    476:         if (jac->diag_use_amat) {
    477:           PetscCall(MatFindZeroDiagonals(pc->mat, &zerodiags));
    478:         } else {
    479:           PetscCall(MatFindZeroDiagonals(pc->pmat, &zerodiags));
    480:         }
    481:         PetscCall(ISComplement(zerodiags, nmin, nmax, &rest));
    482:         PetscCall(PCFieldSplitSetIS(pc, "0", rest));
    483:         PetscCall(PCFieldSplitSetIS(pc, "1", zerodiags));
    484:         PetscCall(ISDestroy(&zerodiags));
    485:         PetscCall(ISDestroy(&rest));
    486:       } else if (coupling) {
    487:         IS       coupling, rest;
    488:         PetscInt nmin, nmax;

    490:         PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax));
    491:         if (jac->offdiag_use_amat) {
    492:           PetscCall(MatFindOffBlockDiagonalEntries(pc->mat, &coupling));
    493:         } else {
    494:           PetscCall(MatFindOffBlockDiagonalEntries(pc->pmat, &coupling));
    495:         }
    496:         PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc->mat), nmax - nmin, nmin, 1, &rest));
    497:         PetscCall(ISSetIdentity(rest));
    498:         PetscCall(PCFieldSplitSetIS(pc, "0", rest));
    499:         PetscCall(PCFieldSplitSetIS(pc, "1", coupling));
    500:         PetscCall(ISDestroy(&coupling));
    501:         PetscCall(ISDestroy(&rest));
    502:       } else {
    503:         PetscCall(PetscOptionsGetBool(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_default", &fieldsplit_default, NULL));
    504:         if (!fieldsplit_default) {
    505:           /* Allow user to set fields from command line,  if bs was known at the time of PCSetFromOptions_FieldSplit()
    506:            then it is set there. This is not ideal because we should only have options set in XXSetFromOptions(). */
    507:           PetscCall(PCFieldSplitSetRuntimeSplits_Private(pc));
    508:           if (jac->splitdefined) PetscCall(PetscInfo(pc, "Splits defined using the options database\n"));
    509:         }
    510:         if ((fieldsplit_default || !jac->splitdefined) && !jac->isrestrict) {
    511:           Mat       M = pc->pmat;
    512:           PetscBool isnest;
    513:           PetscInt  nf;

    515:           PetscCall(PetscInfo(pc, "Using default splitting of fields\n"));
    516:           PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &isnest));
    517:           if (!isnest) {
    518:             M = pc->mat;
    519:             PetscCall(PetscObjectTypeCompare((PetscObject)pc->mat, MATNEST, &isnest));
    520:           }
    521:           if (!isnest) nf = jac->bs;
    522:           else PetscCall(MatNestGetSize(M, &nf, NULL));
    523:           for (i = 0; i < nf; i++) {
    524:             char splitname[8];

    526:             PetscCall(PetscSNPrintf(splitname, sizeof(splitname), "%" PetscInt_FMT, i));
    527:             PetscCall(PCFieldSplitSetFields(pc, splitname, 1, &i, &i));
    528:           }
    529:           jac->defaultsplit = PETSC_TRUE;
    530:         }
    531:       }
    532:     }
    533:   } else if (jac->nsplits == 1) {
    534:     IS       is2;
    535:     PetscInt nmin, nmax;

    537:     PetscCheck(ilink->is, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Must provide at least two sets of fields to PCFieldSplit()");
    538:     PetscCall(MatGetOwnershipRange(pc->mat, &nmin, &nmax));
    539:     PetscCall(ISComplement(ilink->is, nmin, nmax, &is2));
    540:     PetscCall(PCFieldSplitSetIS(pc, "1", is2));
    541:     PetscCall(ISDestroy(&is2));
    542:   }

    544:   PetscCheck(jac->nsplits >= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Unhandled case, must have at least two fields, not %" PetscInt_FMT, jac->nsplits);
    545:   PetscFunctionReturn(PETSC_SUCCESS);
    546: }

    548: static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A, Mat B, Mat C, Mat *H, PetscReal gkbnu)
    549: {
    550:   Mat       BT, T;
    551:   PetscReal nrmT, nrmB;

    553:   PetscFunctionBegin;
    554:   PetscCall(MatHermitianTranspose(C, MAT_INITIAL_MATRIX, &T)); /* Test if augmented matrix is symmetric */
    555:   PetscCall(MatAXPY(T, -1.0, B, DIFFERENT_NONZERO_PATTERN));
    556:   PetscCall(MatNorm(T, NORM_1, &nrmT));
    557:   PetscCall(MatNorm(B, NORM_1, &nrmB));
    558:   PetscCheck(nrmB <= 0 || nrmT / nrmB < PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Matrix is not symmetric/Hermitian, GKB is not applicable.");

    560:   /* Compute augmented Lagrangian matrix H = A00 + nu*A01*A01'. This corresponds to */
    561:   /* setting N := 1/nu*I in [Ar13].                                                 */
    562:   PetscCall(MatHermitianTranspose(B, MAT_INITIAL_MATRIX, &BT));
    563:   PetscCall(MatMatMult(B, BT, MAT_INITIAL_MATRIX, PETSC_CURRENT, H)); /* H = A01*A01'          */
    564:   PetscCall(MatAYPX(*H, gkbnu, A, DIFFERENT_NONZERO_PATTERN));        /* H = A00 + nu*A01*A01' */

    566:   PetscCall(MatDestroy(&BT));
    567:   PetscCall(MatDestroy(&T));
    568:   PetscFunctionReturn(PETSC_SUCCESS);
    569: }

    571: PETSC_EXTERN PetscErrorCode PetscOptionsFindPairPrefix_Private(PetscOptions, const char pre[], const char name[], const char *option[], const char *value[], PetscBool *flg);

    573: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
    574: {
    575:   PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
    576:   PC_FieldSplitLink ilink;
    577:   PetscInt          i, nsplit;
    578:   PetscBool         matnest = PETSC_FALSE;

    580:   PetscFunctionBegin;
    581:   pc->failedreason = PC_NOERROR;
    582:   PetscCall(PCFieldSplitSetDefaults(pc));
    583:   nsplit = jac->nsplits;
    584:   ilink  = jac->head;
    585:   if (pc->pmat) PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &matnest));

    587:   /* get the matrices for each split */
    588:   if (!jac->issetup) {
    589:     PetscInt rstart, rend, nslots, bs;

    591:     jac->issetup = PETSC_TRUE;

    593:     /* This is done here instead of in PCFieldSplitSetFields() because may not have matrix at that point */
    594:     if (jac->defaultsplit || !ilink->is) {
    595:       if (jac->bs <= 0) jac->bs = nsplit;
    596:     }

    598:     /*  MatCreateSubMatrix() for [S]BAIJ matrices can only work if the indices include entire blocks of the matrix */
    599:     PetscCall(MatGetBlockSize(pc->pmat, &bs));
    600:     if (bs > 1 && (jac->bs <= bs || jac->bs % bs)) {
    601:       PetscBool blk;

    603:       PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &blk, MATBAIJ, MATSBAIJ, MATSEQBAIJ, MATSEQSBAIJ, MATMPIBAIJ, MATMPISBAIJ, NULL));
    604:       PetscCheck(!blk, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Cannot use MATBAIJ with PCFIELDSPLIT and currently set matrix and PC blocksizes");
    605:     }

    607:     if (!matnest) { /* use the matrix blocksize and stride IS to determine the index sets that define the submatrices */
    608:       bs = jac->bs;
    609:       PetscCall(MatGetOwnershipRange(pc->pmat, &rstart, &rend));
    610:       nslots = (rend - rstart) / bs;
    611:       for (i = 0; i < nsplit; i++) {
    612:         if (jac->defaultsplit) {
    613:           PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + i, nsplit, &ilink->is));
    614:           PetscCall(PetscObjectReference((PetscObject)ilink->is));
    615:           ilink->is_col = ilink->is;
    616:         } else if (!ilink->is) {
    617:           PetscBool same_fields = PETSC_TRUE;

    619:           for (PetscInt k = 0; k < ilink->nfields; k++) {
    620:             if (ilink->fields[k] != ilink->fields_col[k]) same_fields = PETSC_FALSE;
    621:           }

    623:           if (ilink->nfields > 1) {
    624:             PetscInt *ii, *jj, j, k, nfields = ilink->nfields, *fields = ilink->fields, *fields_col = ilink->fields_col;

    626:             PetscCall(PetscMalloc1(ilink->nfields * nslots, &ii));
    627:             if (!same_fields) PetscCall(PetscMalloc1(ilink->nfields * nslots, &jj));
    628:             for (j = 0; j < nslots; j++) {
    629:               for (k = 0; k < nfields; k++) {
    630:                 ii[nfields * j + k] = rstart + bs * j + fields[k];
    631:                 if (!same_fields) jj[nfields * j + k] = rstart + bs * j + fields_col[k];
    632:               }
    633:             }
    634:             PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, ii, PETSC_OWN_POINTER, &ilink->is));
    635:             if (!same_fields) PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, jj, PETSC_OWN_POINTER, &ilink->is_col));
    636:             else {
    637:               PetscCall(PetscObjectReference((PetscObject)ilink->is));
    638:               ilink->is_col = ilink->is;
    639:             }
    640:             PetscCall(ISSetBlockSize(ilink->is, nfields));
    641:             PetscCall(ISSetBlockSize(ilink->is_col, nfields));
    642:           } else {
    643:             PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields[0], bs, &ilink->is));
    644:             if (!same_fields) PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields_col[0], bs, &ilink->is_col));
    645:             else {
    646:               PetscCall(PetscObjectReference((PetscObject)ilink->is));
    647:               ilink->is_col = ilink->is;
    648:             }
    649:           }
    650:         }
    651:         ilink = ilink->next;
    652:       }
    653:     } else { /* use the IS that define the MATNEST to determine the index sets that define the submatrices */
    654:       IS      *rowis, *colis, *ises = NULL;
    655:       PetscInt mis, nis;

    657:       PetscCall(MatNestGetSize(pc->pmat, &mis, &nis));
    658:       PetscCall(PetscMalloc2(mis, &rowis, nis, &colis));
    659:       PetscCall(MatNestGetISs(pc->pmat, rowis, colis));
    660:       if (!jac->defaultsplit) PetscCall(PetscMalloc1(mis, &ises));

    662:       for (i = 0; i < nsplit; i++) {
    663:         if (jac->defaultsplit) {
    664:           PetscCall(ISDuplicate(rowis[i], &ilink->is));
    665:           PetscCall(PetscObjectReference((PetscObject)ilink->is));
    666:           ilink->is_col = ilink->is;
    667:         } else if (!ilink->is) {
    668:           if (ilink->nfields > 1) {
    669:             for (PetscInt j = 0; j < ilink->nfields; j++) ises[j] = rowis[ilink->fields[j]];
    670:             PetscCall(ISConcatenate(PetscObjectComm((PetscObject)pc), ilink->nfields, ises, &ilink->is));
    671:           } else {
    672:             PetscCall(ISDuplicate(rowis[ilink->fields[0]], &ilink->is));
    673:           }
    674:           PetscCall(PetscObjectReference((PetscObject)ilink->is));
    675:           ilink->is_col = ilink->is;
    676:         }
    677:         ilink = ilink->next;
    678:       }
    679:       PetscCall(PetscFree2(rowis, colis));
    680:       PetscCall(PetscFree(ises));
    681:     }
    682:   }

    684:   ilink = jac->head;
    685:   if (!jac->pmat) {
    686:     Vec xtmp;

    688:     PetscCall(MatCreateVecs(pc->pmat, &xtmp, NULL));
    689:     PetscCall(PetscMalloc1(nsplit, &jac->pmat));
    690:     PetscCall(PetscMalloc2(nsplit, &jac->x, nsplit, &jac->y));
    691:     for (i = 0; i < nsplit; i++) {
    692:       MatNullSpace sp;

    694:       /* Check for matrix attached to IS */
    695:       PetscCall(PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&jac->pmat[i]));
    696:       if (jac->pmat[i]) {
    697:         PetscCall(PetscObjectReference((PetscObject)jac->pmat[i]));
    698:         if (jac->type == PC_COMPOSITE_SCHUR) {
    699:           jac->schur_user = jac->pmat[i];

    701:           PetscCall(PetscObjectReference((PetscObject)jac->schur_user));
    702:         }
    703:       } else {
    704:         const char *prefix;
    705:         PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->pmat[i]));
    706:         PetscCall(MatGetOptionsPrefix(jac->pmat[i], &prefix));
    707:         if (!prefix) {
    708:           PetscCall(KSPGetOptionsPrefix(ilink->ksp, &prefix));
    709:           PetscCall(MatSetOptionsPrefix(jac->pmat[i], prefix));
    710:         }
    711:         PetscCall(MatSetFromOptions(jac->pmat[i]));
    712:         PetscCall(MatViewFromOptions(jac->pmat[i], NULL, "-mat_view"));
    713:       }
    714:       /* create work vectors for each split */
    715:       PetscCall(MatCreateVecs(jac->pmat[i], &jac->x[i], &jac->y[i]));
    716:       ilink->x = jac->x[i];
    717:       ilink->y = jac->y[i];
    718:       ilink->z = NULL;
    719:       /* compute scatter contexts needed by multiplicative versions and non-default splits */
    720:       PetscCall(VecScatterCreate(xtmp, ilink->is, jac->x[i], NULL, &ilink->sctx));
    721:       PetscCall(PetscObjectQuery((PetscObject)ilink->is, "nearnullspace", (PetscObject *)&sp));
    722:       if (sp) PetscCall(MatSetNearNullSpace(jac->pmat[i], sp));
    723:       ilink = ilink->next;
    724:     }
    725:     PetscCall(VecDestroy(&xtmp));
    726:   } else {
    727:     MatReuse      scall;
    728:     MatNullSpace *nullsp = NULL;

    730:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
    731:       PetscCall(MatGetNullSpaces(nsplit, jac->pmat, &nullsp));
    732:       for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->pmat[i]));
    733:       scall = MAT_INITIAL_MATRIX;
    734:     } else scall = MAT_REUSE_MATRIX;

    736:     for (i = 0; i < nsplit; i++) {
    737:       Mat pmat;

    739:       /* Check for matrix attached to IS */
    740:       PetscCall(PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&pmat));
    741:       if (!pmat) PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, scall, &jac->pmat[i]));
    742:       ilink = ilink->next;
    743:     }
    744:     if (nullsp) PetscCall(MatRestoreNullSpaces(nsplit, jac->pmat, &nullsp));
    745:   }
    746:   if (jac->diag_use_amat) {
    747:     ilink = jac->head;
    748:     if (!jac->mat) {
    749:       PetscCall(PetscMalloc1(nsplit, &jac->mat));
    750:       for (i = 0; i < nsplit; i++) {
    751:         PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->mat[i]));
    752:         ilink = ilink->next;
    753:       }
    754:     } else {
    755:       MatReuse      scall;
    756:       MatNullSpace *nullsp = NULL;

    758:       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
    759:         PetscCall(MatGetNullSpaces(nsplit, jac->mat, &nullsp));
    760:         for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->mat[i]));
    761:         scall = MAT_INITIAL_MATRIX;
    762:       } else scall = MAT_REUSE_MATRIX;

    764:       for (i = 0; i < nsplit; i++) {
    765:         PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, scall, &jac->mat[i]));
    766:         ilink = ilink->next;
    767:       }
    768:       if (nullsp) PetscCall(MatRestoreNullSpaces(nsplit, jac->mat, &nullsp));
    769:     }
    770:   } else {
    771:     jac->mat = jac->pmat;
    772:   }

    774:   /* Check for null space attached to IS */
    775:   ilink = jac->head;
    776:   for (i = 0; i < nsplit; i++) {
    777:     MatNullSpace sp;

    779:     PetscCall(PetscObjectQuery((PetscObject)ilink->is, "nullspace", (PetscObject *)&sp));
    780:     if (sp) PetscCall(MatSetNullSpace(jac->mat[i], sp));
    781:     ilink = ilink->next;
    782:   }

    784:   if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) {
    785:     /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
    786:     /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
    787:     ilink = jac->head;
    788:     if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
    789:       /* special case need where Afield[0] is not needed and only certain columns of Afield[1] are needed since update is only on those rows of the solution */
    790:       if (!jac->Afield) {
    791:         PetscCall(PetscCalloc1(nsplit, &jac->Afield));
    792:         if (jac->offdiag_use_amat) {
    793:           PetscCall(MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1]));
    794:         } else {
    795:           PetscCall(MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1]));
    796:         }
    797:       } else {
    798:         MatReuse scall;

    800:         if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
    801:           PetscCall(MatDestroy(&jac->Afield[1]));
    802:           scall = MAT_INITIAL_MATRIX;
    803:         } else scall = MAT_REUSE_MATRIX;

    805:         if (jac->offdiag_use_amat) {
    806:           PetscCall(MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, scall, &jac->Afield[1]));
    807:         } else {
    808:           PetscCall(MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, scall, &jac->Afield[1]));
    809:         }
    810:       }
    811:     } else {
    812:       if (!jac->Afield) {
    813:         PetscCall(PetscMalloc1(nsplit, &jac->Afield));
    814:         for (i = 0; i < nsplit; i++) {
    815:           if (jac->offdiag_use_amat) {
    816:             PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i]));
    817:           } else {
    818:             PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i]));
    819:           }
    820:           ilink = ilink->next;
    821:         }
    822:       } else {
    823:         MatReuse scall;
    824:         if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
    825:           for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->Afield[i]));
    826:           scall = MAT_INITIAL_MATRIX;
    827:         } else scall = MAT_REUSE_MATRIX;

    829:         for (i = 0; i < nsplit; i++) {
    830:           if (jac->offdiag_use_amat) {
    831:             PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, NULL, scall, &jac->Afield[i]));
    832:           } else {
    833:             PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, NULL, scall, &jac->Afield[i]));
    834:           }
    835:           ilink = ilink->next;
    836:         }
    837:       }
    838:     }
    839:   }

    841:   if (jac->type == PC_COMPOSITE_SCHUR) {
    842:     PetscBool   isset, isspd = PETSC_FALSE, issym = PETSC_FALSE, flg;
    843:     char        lscname[256];
    844:     PetscObject LSC_L;

    846:     PetscCheck(nsplit == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "To use Schur complement preconditioner you must have exactly 2 fields");

    848:     /* If pc->mat is SPD, don't scale by -1 the Schur complement */
    849:     PetscCall(MatIsSPDKnown(pc->pmat, &isset, &isspd));
    850:     if (jac->schurscale == (PetscScalar)-1.0) jac->schurscale = (isset && isspd) ? 1.0 : -1.0;
    851:     PetscCall(MatIsSymmetricKnown(pc->pmat, &isset, &issym));

    853:     PetscCall(PetscObjectTypeCompareAny(jac->offdiag_use_amat ? (PetscObject)pc->mat : (PetscObject)pc->pmat, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));

    855:     if (jac->schur) {
    856:       KSP      kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
    857:       MatReuse scall;

    859:       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
    860:         scall = MAT_INITIAL_MATRIX;
    861:         PetscCall(MatDestroy(&jac->B));
    862:         PetscCall(MatDestroy(&jac->C));
    863:       } else scall = MAT_REUSE_MATRIX;

    865:       PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner));
    866:       ilink = jac->head;
    867:       PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->is, ilink->next->is, scall, &jac->B));
    868:       if (!flg) PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->next->is, ilink->is, scall, &jac->C));
    869:       else {
    870:         PetscCall(MatIsHermitianKnown(jac->offdiag_use_amat ? pc->mat : pc->pmat, &isset, &flg));
    871:         if (isset && flg) PetscCall(MatCreateHermitianTranspose(jac->B, &jac->C));
    872:         else PetscCall(MatCreateTranspose(jac->B, &jac->C));
    873:       }
    874:       ilink = ilink->next;
    875:       PetscCall(MatSchurComplementUpdateSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1]));
    876:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
    877:         PetscCall(MatDestroy(&jac->schurp));
    878:         PetscCall(MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp));
    879:       } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL && jac->kspupper != jac->head->ksp) {
    880:         PetscCall(MatDestroy(&jac->schur_user));
    881:         PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user));
    882:       }
    883:       if (kspA != kspInner) PetscCall(KSPSetOperators(kspA, jac->mat[0], jac->pmat[0]));
    884:       if (kspUpper != kspA) PetscCall(KSPSetOperators(kspUpper, jac->mat[0], jac->pmat[0]));
    885:       PetscCall(KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac)));
    886:     } else {
    887:       const char  *Dprefix;
    888:       char         schurprefix[256], schurmatprefix[256];
    889:       char         schurtestoption[256];
    890:       MatNullSpace sp;
    891:       KSP          kspt;

    893:       /* extract the A01 and A10 matrices */
    894:       ilink = jac->head;
    895:       PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->is, ilink->next->is, MAT_INITIAL_MATRIX, &jac->B));
    896:       if (!flg) PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->C));
    897:       else {
    898:         PetscCall(MatIsHermitianKnown(jac->offdiag_use_amat ? pc->mat : pc->pmat, &isset, &flg));
    899:         if (isset && flg) PetscCall(MatCreateHermitianTranspose(jac->B, &jac->C));
    900:         else PetscCall(MatCreateTranspose(jac->B, &jac->C));
    901:       }
    902:       ilink = ilink->next;
    903:       /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
    904:       PetscCall(MatCreate(((PetscObject)jac->mat[0])->comm, &jac->schur));
    905:       PetscCall(MatSetType(jac->schur, MATSCHURCOMPLEMENT));
    906:       PetscCall(MatSchurComplementSetSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1]));
    907:       PetscCall(PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
    908:       PetscCall(MatSetOptionsPrefix(jac->schur, schurmatprefix));
    909:       PetscCall(MatSchurComplementGetKSP(jac->schur, &kspt));
    910:       PetscCall(KSPSetOptionsPrefix(kspt, schurmatprefix));

    912:       /* Note: this is not true in general */
    913:       PetscCall(MatGetNullSpace(jac->mat[1], &sp));
    914:       if (sp) PetscCall(MatSetNullSpace(jac->schur, sp));

    916:       PetscCall(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname));
    917:       PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, NULL, &flg));
    918:       if (flg) {
    919:         DM  dmInner;
    920:         KSP kspInner;
    921:         PC  pcInner;

    923:         PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner));
    924:         PetscCall(KSPReset(kspInner));
    925:         PetscCall(KSPSetOperators(kspInner, jac->mat[0], jac->pmat[0]));
    926:         PetscCall(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
    927:         /* Indent this deeper to emphasize the "inner" nature of this solver. */
    928:         PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject)pc, 2));
    929:         PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject)pc, 2));
    930:         PetscCall(KSPSetOptionsPrefix(kspInner, schurprefix));

    932:         /* Set DM for new solver */
    933:         PetscCall(KSPGetDM(jac->head->ksp, &dmInner));
    934:         PetscCall(KSPSetDM(kspInner, dmInner));
    935:         PetscCall(KSPSetDMActive(kspInner, KSP_DMACTIVE_ALL, PETSC_FALSE));

    937:         /* Defaults to PCKSP as preconditioner */
    938:         PetscCall(KSPGetPC(kspInner, &pcInner));
    939:         PetscCall(PCSetType(pcInner, PCKSP));
    940:         PetscCall(PCKSPSetKSP(pcInner, jac->head->ksp));
    941:       } else {
    942:         /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
    943:           * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
    944:           * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
    945:           * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
    946:           * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
    947:           * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
    948:         PetscCall(KSPSetType(jac->head->ksp, KSPGMRES));
    949:         PetscCall(MatSchurComplementSetKSP(jac->schur, jac->head->ksp));
    950:       }
    951:       PetscCall(KSPSetOperators(jac->head->ksp, jac->mat[0], jac->pmat[0]));
    952:       PetscCall(KSPSetFromOptions(jac->head->ksp));
    953:       PetscCall(MatSetFromOptions(jac->schur));

    955:       PetscCall(PetscObjectTypeCompare((PetscObject)jac->schur, MATSCHURCOMPLEMENT, &flg));
    956:       if (flg) { /* Need to do this otherwise PCSetUp_KSP will overwrite the amat of jac->head->ksp */
    957:         KSP kspInner;
    958:         PC  pcInner;

    960:         PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner));
    961:         PetscCall(KSPGetPC(kspInner, &pcInner));
    962:         PetscCall(PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg));
    963:         if (flg) {
    964:           KSP ksp;

    966:           PetscCall(PCKSPGetKSP(pcInner, &ksp));
    967:           if (ksp == jac->head->ksp) PetscCall(PCSetUseAmat(pcInner, PETSC_TRUE));
    968:         }
    969:       }
    970:       PetscCall(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname));
    971:       PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, NULL, &flg));
    972:       if (flg) {
    973:         DM dmInner;

    975:         PetscCall(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
    976:         PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper));
    977:         PetscCall(KSPSetNestLevel(jac->kspupper, pc->kspnestlevel));
    978:         PetscCall(KSPSetErrorIfNotConverged(jac->kspupper, pc->erroriffailure));
    979:         PetscCall(KSPSetOptionsPrefix(jac->kspupper, schurprefix));
    980:         PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject)pc, 1));
    981:         PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject)pc, 1));
    982:         PetscCall(KSPGetDM(jac->head->ksp, &dmInner));
    983:         PetscCall(KSPSetDM(jac->kspupper, dmInner));
    984:         PetscCall(KSPSetDMActive(jac->kspupper, KSP_DMACTIVE_ALL, PETSC_FALSE));
    985:         PetscCall(KSPSetFromOptions(jac->kspupper));
    986:         PetscCall(KSPSetOperators(jac->kspupper, jac->mat[0], jac->pmat[0]));
    987:         PetscCall(VecDuplicate(jac->head->x, &jac->head->z));
    988:       } else {
    989:         jac->kspupper = jac->head->ksp;
    990:         PetscCall(PetscObjectReference((PetscObject)jac->head->ksp));
    991:       }

    993:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) PetscCall(MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp));
    994:       PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspschur));
    995:       PetscCall(KSPSetNestLevel(jac->kspschur, pc->kspnestlevel));
    996:       PetscCall(KSPSetErrorIfNotConverged(jac->kspschur, pc->erroriffailure));
    997:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspschur, (PetscObject)pc, 1));
    998:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
    999:         PC pcschur;
    1000:         PetscCall(KSPGetPC(jac->kspschur, &pcschur));
    1001:         PetscCall(PCSetType(pcschur, PCNONE));
    1002:         /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
    1003:       } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
    1004:         if (jac->schurfactorization != PC_FIELDSPLIT_SCHUR_FACT_FULL || jac->kspupper != jac->head->ksp) PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user));
    1005:       }
    1006:       PetscCall(KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac)));
    1007:       PetscCall(KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix));
    1008:       PetscCall(KSPSetOptionsPrefix(jac->kspschur, Dprefix));
    1009:       /* propagate DM */
    1010:       {
    1011:         DM sdm;
    1012:         PetscCall(KSPGetDM(jac->head->next->ksp, &sdm));
    1013:         if (sdm) {
    1014:           PetscCall(KSPSetDM(jac->kspschur, sdm));
    1015:           PetscCall(KSPSetDMActive(jac->kspschur, KSP_DMACTIVE_ALL, PETSC_FALSE));
    1016:         }
    1017:       }
    1018:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
    1019:       /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
    1020:       PetscCall(KSPSetFromOptions(jac->kspschur));
    1021:     }
    1022:     PetscCall(MatAssemblyBegin(jac->schur, MAT_FINAL_ASSEMBLY));
    1023:     PetscCall(MatAssemblyEnd(jac->schur, MAT_FINAL_ASSEMBLY));
    1024:     if (issym) PetscCall(MatSetOption(jac->schur, MAT_SYMMETRIC, PETSC_TRUE));
    1025:     if (isspd) PetscCall(MatSetOption(jac->schur, MAT_SPD, PETSC_TRUE));

    1027:     /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
    1028:     PetscCall(PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_L", ilink->splitname));
    1029:     PetscCall(PetscObjectQuery((PetscObject)pc->mat, lscname, &LSC_L));
    1030:     if (!LSC_L) PetscCall(PetscObjectQuery((PetscObject)pc->pmat, lscname, &LSC_L));
    1031:     if (LSC_L) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "LSC_L", LSC_L));
    1032:     PetscCall(PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_Lp", ilink->splitname));
    1033:     PetscCall(PetscObjectQuery((PetscObject)pc->pmat, lscname, &LSC_L));
    1034:     if (!LSC_L) PetscCall(PetscObjectQuery((PetscObject)pc->mat, lscname, &LSC_L));
    1035:     if (LSC_L) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "LSC_Lp", LSC_L));
    1036:   } else if (jac->type == PC_COMPOSITE_GKB) {
    1037:     PetscCheck(nsplit == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "To use GKB preconditioner you must have exactly 2 fields");
    1038:     ilink = jac->head;
    1039:     PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->is, ilink->next->is, MAT_INITIAL_MATRIX, &jac->B));
    1040:     /* Create work vectors for GKB algorithm */
    1041:     PetscCall(VecDuplicate(ilink->x, &jac->u));
    1042:     PetscCall(VecDuplicate(ilink->x, &jac->Hu));
    1043:     PetscCall(VecDuplicate(ilink->x, &jac->w2));
    1044:     PetscCall(MatCreateSubMatrix(jac->offdiag_use_amat ? pc->mat : pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->C));
    1045:     ilink = ilink->next;
    1046:     /* Create work vectors for GKB algorithm */
    1047:     PetscCall(VecDuplicate(ilink->x, &jac->v));
    1048:     PetscCall(VecDuplicate(ilink->x, &jac->d));
    1049:     PetscCall(VecDuplicate(ilink->x, &jac->w1));
    1050:     PetscCall(MatGolubKahanComputeExplicitOperator(jac->mat[0], jac->B, jac->C, &jac->H, jac->gkbnu));
    1051:     PetscCall(PetscCalloc1(jac->gkbdelay, &jac->vecz));

    1053:     ilink = jac->head;
    1054:     PetscCall(KSPSetOperators(ilink->ksp, jac->H, jac->H));
    1055:     if (!jac->suboptionsset) PetscCall(KSPSetFromOptions(ilink->ksp));
    1056:     /* Create gkb_monitor context */
    1057:     if (jac->gkbmonitor) {
    1058:       PetscInt tablevel;
    1059:       PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &jac->gkbviewer));
    1060:       PetscCall(PetscViewerSetType(jac->gkbviewer, PETSCVIEWERASCII));
    1061:       PetscCall(PetscObjectGetTabLevel((PetscObject)ilink->ksp, &tablevel));
    1062:       PetscCall(PetscViewerASCIISetTab(jac->gkbviewer, tablevel));
    1063:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)ilink->ksp, 1));
    1064:     }
    1065:   } else {
    1066:     /* set up the individual splits' PCs */
    1067:     i     = 0;
    1068:     ilink = jac->head;
    1069:     while (ilink) {
    1070:       PetscCall(KSPSetOperators(ilink->ksp, jac->mat[i], jac->pmat[i]));
    1071:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
    1072:       if (!jac->suboptionsset) PetscCall(KSPSetFromOptions(ilink->ksp));
    1073:       i++;
    1074:       ilink = ilink->next;
    1075:     }
    1076:   }

    1078:   /* Set coordinates to the sub PC objects whenever these are set */
    1079:   if (jac->coordinates_set) {
    1080:     PC pc_coords;
    1081:     if (jac->type == PC_COMPOSITE_SCHUR) {
    1082:       // Head is first block.
    1083:       PetscCall(KSPGetPC(jac->head->ksp, &pc_coords));
    1084:       PetscCall(PCSetCoordinates(pc_coords, jac->head->dim, jac->head->ndofs, jac->head->coords));
    1085:       // Second one is Schur block, but its KSP object is in kspschur.
    1086:       PetscCall(KSPGetPC(jac->kspschur, &pc_coords));
    1087:       PetscCall(PCSetCoordinates(pc_coords, jac->head->next->dim, jac->head->next->ndofs, jac->head->next->coords));
    1088:     } else if (jac->type == PC_COMPOSITE_GKB) {
    1089:       PetscCall(PetscInfo(pc, "Warning: Setting coordinates does nothing for the GKB Fieldpslit preconditioner\n"));
    1090:     } else {
    1091:       ilink = jac->head;
    1092:       while (ilink) {
    1093:         PetscCall(KSPGetPC(ilink->ksp, &pc_coords));
    1094:         PetscCall(PCSetCoordinates(pc_coords, ilink->dim, ilink->ndofs, ilink->coords));
    1095:         ilink = ilink->next;
    1096:       }
    1097:     }
    1098:   }

    1100:   jac->suboptionsset = PETSC_TRUE;
    1101:   PetscFunctionReturn(PETSC_SUCCESS);
    1102: }

    1104: static PetscErrorCode PCSetUpOnBlocks_FieldSplit_Schur(PC pc)
    1105: {
    1106:   PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
    1107:   PC_FieldSplitLink ilinkA = jac->head;
    1108:   KSP               kspA = ilinkA->ksp, kspUpper = jac->kspupper;

    1110:   PetscFunctionBegin;
    1111:   if (jac->schurfactorization == PC_FIELDSPLIT_SCHUR_FACT_FULL && kspUpper != kspA) {
    1112:     PetscCall(KSPSetUp(kspUpper));
    1113:     PetscCall(KSPSetUpOnBlocks(kspUpper));
    1114:   }
    1115:   PetscCall(KSPSetUp(kspA));
    1116:   PetscCall(KSPSetUpOnBlocks(kspA));
    1117:   if (jac->schurpre != PC_FIELDSPLIT_SCHUR_PRE_FULL) {
    1118:     PetscCall(KSPSetUp(jac->kspschur));
    1119:     PetscCall(KSPSetUpOnBlocks(jac->kspschur));
    1120:   } else if (kspUpper == kspA && jac->schurfactorization == PC_FIELDSPLIT_SCHUR_FACT_FULL) {
    1121:     Mat          A;
    1122:     PetscInt     m, M, N;
    1123:     VecType      vtype;
    1124:     PetscMemType mtype;
    1125:     PetscScalar *array;

    1127:     PetscCall(MatGetSize(jac->B, &M, &N));
    1128:     PetscCall(MatGetLocalSize(jac->B, &m, NULL));
    1129:     PetscCall(MatGetVecType(jac->B, &vtype));
    1130:     PetscCall(VecGetArrayAndMemType(ilinkA->x, &array, &mtype));
    1131:     PetscCall(VecRestoreArrayAndMemType(ilinkA->x, &array));
    1132:     PetscCall(PetscObjectQuery((PetscObject)jac->schur, "AinvB", (PetscObject *)&A));
    1133:     if (A) {
    1134:       PetscInt P;

    1136:       PetscCall(MatGetSize(A, NULL, &P));
    1137:       if (P < N + 1) { // need to recreate AinvB, otherwise, the Schur complement won't be updated
    1138:         PetscCall(PetscObjectCompose((PetscObject)jac->schur, "AinvB", NULL));
    1139:         A = NULL;
    1140:       }
    1141:     }
    1142:     if (!A) {
    1143:       if (PetscMemTypeHost(mtype) || (!PetscDefined(HAVE_CUDA) && !PetscDefined(HAVE_HIP))) PetscCall(PetscMalloc1(m * (N + 1), &array));
    1144: #if PetscDefined(HAVE_CUDA)
    1145:       else if (PetscMemTypeCUDA(mtype)) PetscCallCUDA(cudaMalloc((void **)&array, sizeof(PetscScalar) * m * (N + 1)));
    1146: #endif
    1147: #if PetscDefined(HAVE_HIP)
    1148:       else if (PetscMemTypeHIP(mtype)) PetscCallHIP(hipMalloc((void **)&array, sizeof(PetscScalar) * m * (N + 1)));
    1149: #endif
    1150:       PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)jac->schur), vtype, m, PETSC_DECIDE, M, N + 1, PETSC_DECIDE, array, &A)); // number of columns of the Schur complement plus one
    1151:       PetscCall(PetscObjectCompose((PetscObject)jac->schur, "AinvB", (PetscObject)A));
    1152:       PetscCall(MatDestroy(&A));
    1153:     }
    1154:   }
    1155:   PetscFunctionReturn(PETSC_SUCCESS);
    1156: }

    1158: static PetscErrorCode PCSetUpOnBlocks_FieldSplit(PC pc)
    1159: {
    1160:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    1161:   PC_FieldSplitLink ilink = jac->head;

    1163:   PetscFunctionBegin;
    1164:   while (ilink) {
    1165:     PetscCall(KSPSetUp(ilink->ksp));
    1166:     PetscCall(KSPSetUpOnBlocks(ilink->ksp));
    1167:     ilink = ilink->next;
    1168:   }
    1169:   PetscFunctionReturn(PETSC_SUCCESS);
    1170: }

    1172: static PetscErrorCode PCSetUpOnBlocks_FieldSplit_GKB(PC pc)
    1173: {
    1174:   PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
    1175:   PC_FieldSplitLink ilinkA = jac->head;
    1176:   KSP               ksp    = ilinkA->ksp;

    1178:   PetscFunctionBegin;
    1179:   PetscCall(KSPSetUp(ksp));
    1180:   PetscCall(KSPSetUpOnBlocks(ksp));
    1181:   PetscFunctionReturn(PETSC_SUCCESS);
    1182: }

    1184: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc, Vec x, Vec y)
    1185: {
    1186:   PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
    1187:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
    1188:   KSP               kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
    1189:   Mat               AinvB = NULL;
    1190:   PetscInt          N, P;

    1192:   PetscFunctionBegin;
    1193:   switch (jac->schurfactorization) {
    1194:   case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
    1195:     /* [A00 0; 0 -S], positive definite, suitable for MINRES */
    1196:     PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1197:     PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1198:     PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1199:     PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1200:     PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
    1201:     PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1202:     PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1203:     PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1204:     PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1205:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1206:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1207:     PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
    1208:     PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    1209:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1210:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1211:     PetscCall(VecScale(ilinkD->y, jac->schurscale));
    1212:     PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1213:     PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1214:     PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1215:     break;
    1216:   case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
    1217:     /* [A00 0; A10 S], suitable for left preconditioning */
    1218:     PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1219:     PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1220:     PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1221:     PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
    1222:     PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1223:     PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1224:     PetscCall(MatMult(jac->C, ilinkA->y, ilinkD->x));
    1225:     PetscCall(VecScale(ilinkD->x, -1.));
    1226:     PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    1227:     PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1228:     PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    1229:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1230:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1231:     PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
    1232:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1233:     PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    1234:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1235:     PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1236:     PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1237:     PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1238:     break;
    1239:   case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
    1240:     /* [A00 A01; 0 S], suitable for right preconditioning */
    1241:     PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1242:     PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1243:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1244:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1245:     PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
    1246:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1247:     PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    1248:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1249:     PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->x));
    1250:     PetscCall(VecScale(ilinkA->x, -1.));
    1251:     PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
    1252:     PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1253:     PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
    1254:     PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1255:     PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
    1256:     PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1257:     PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1258:     PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1259:     PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1260:     PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1261:     break;
    1262:   case PC_FIELDSPLIT_SCHUR_FACT_FULL:
    1263:     /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */
    1264:     PetscCall(MatGetSize(jac->B, NULL, &P));
    1265:     N = P;
    1266:     PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1267:     PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1268:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL));
    1269:     if (kspUpper == kspA) {
    1270:       PetscCall(PetscObjectQuery((PetscObject)jac->schur, "AinvB", (PetscObject *)&AinvB));
    1271:       if (AinvB) {
    1272:         PetscCall(MatGetSize(AinvB, NULL, &N));
    1273:         if (N > P) { // first time PCApply_FieldSplit_Schur() is called
    1274:           PetscMemType mtype;
    1275:           Vec          c = NULL;
    1276:           PetscScalar *array;
    1277:           PetscInt     m, M;

    1279:           PetscCall(MatGetSize(jac->B, &M, NULL));
    1280:           PetscCall(MatGetLocalSize(jac->B, &m, NULL));
    1281:           PetscCall(MatDenseGetArrayAndMemType(AinvB, &array, &mtype));
    1282:           if (PetscMemTypeHost(mtype) || (!PetscDefined(HAVE_CUDA) && !PetscDefined(HAVE_HIP))) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)jac->schur), 1, m, M, array + m * P, &c));
    1283: #if PetscDefined(HAVE_CUDA)
    1284:           else if (PetscMemTypeCUDA(mtype)) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)jac->schur), 1, m, M, array + m * P, &c));
    1285: #endif
    1286: #if PetscDefined(HAVE_HIP)
    1287:           else if (PetscMemTypeHIP(mtype)) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)jac->schur), 1, m, M, array + m * P, &c));
    1288: #endif
    1289:           PetscCall(MatDenseRestoreArrayAndMemType(AinvB, &array));
    1290:           PetscCall(VecCopy(ilinkA->x, c));
    1291:           PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user));
    1292:           PetscCall(KSPSetOperators(jac->kspschur, jac->schur, jac->schur_user));
    1293:           PetscCall(VecCopy(c, ilinkA->y)); // retrieve the solution as the last column of the composed Mat
    1294:           PetscCall(VecDestroy(&c));
    1295:         }
    1296:       }
    1297:     }
    1298:     if (N == P) PetscCall(KSPSolve(kspLower, ilinkA->x, ilinkA->y));
    1299:     PetscCall(KSPCheckSolve(kspLower, pc, ilinkA->y));
    1300:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL));
    1301:     PetscCall(MatMult(jac->C, ilinkA->y, ilinkD->x));
    1302:     PetscCall(VecScale(ilinkD->x, -1.0));
    1303:     PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    1304:     PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));

    1306:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1307:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1308:     PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
    1309:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1310:     PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    1311:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1312:     PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));

    1314:     if (kspUpper == kspA) {
    1315:       if (!AinvB) {
    1316:         PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->y));
    1317:         PetscCall(VecAXPY(ilinkA->x, -1.0, ilinkA->y));
    1318:         PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1319:         PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
    1320:         PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1321:         PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1322:       } else PetscCall(MatMultAdd(AinvB, ilinkD->y, ilinkA->y, ilinkA->y));
    1323:     } else {
    1324:       PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1325:       PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
    1326:       PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1327:       PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->x));
    1328:       PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL));
    1329:       PetscCall(KSPSolve(kspUpper, ilinkA->x, ilinkA->z));
    1330:       PetscCall(KSPCheckSolve(kspUpper, pc, ilinkA->z));
    1331:       PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL));
    1332:       PetscCall(VecAXPY(ilinkA->y, -1.0, ilinkA->z));
    1333:     }
    1334:     PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1335:     PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1336:     PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1337:   }
    1338:   PetscFunctionReturn(PETSC_SUCCESS);
    1339: }

    1341: /*
    1342:   PCFieldSplitCreateWorkMats_Private - Allocate per-field dense work matrices for multi-RHS

    1344:   Input Parameters:
    1345: + pc - the PC context
    1346: - X  - matrix to copy column-layout from

    1348:   Notes:
    1349:   If matrices already exist with correct column count, they are reused.
    1350:   If column count changed, old matrices are destroyed and new ones created.
    1351: */
    1352: static PetscErrorCode PCFieldSplitCreateWorkMats_Private(PC pc, Mat X)
    1353: {
    1354:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    1355:   PC_FieldSplitLink ilink = jac->head;
    1356:   PetscInt          mx, Mx, my, My, N;

    1358:   PetscFunctionBegin;
    1359:   while (ilink) {
    1360:     /* check if reallocation needed (previous allocation with wrong column count) */
    1361:     if (ilink->X) {
    1362:       PetscCall(MatGetSize(ilink->X, NULL, &N));
    1363:       if (N != X->cmap->N) {
    1364:         PetscCall(MatDestroy(&ilink->X));
    1365:         PetscCall(MatDestroy(&ilink->Y));
    1366:         PetscCall(MatDestroy(&ilink->Z));
    1367:       }
    1368:     }
    1369:     /* create if needed */
    1370:     if (!ilink->X) {
    1371:       VecType xtype, ytype;

    1373:       PetscCall(VecGetType(ilink->x, &xtype));
    1374:       PetscCall(VecGetType(ilink->y, &ytype));
    1375:       PetscCall(VecGetLocalSize(ilink->x, &mx));
    1376:       PetscCall(VecGetSize(ilink->x, &Mx));
    1377:       PetscCall(VecGetLocalSize(ilink->y, &my));
    1378:       PetscCall(VecGetSize(ilink->y, &My));
    1379:       /* use default lda */
    1380:       PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)pc), xtype, mx, X->cmap->n, Mx, X->cmap->N, PETSC_DECIDE, NULL, &ilink->X));
    1381:       PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)pc), ytype, my, X->cmap->n, My, X->cmap->N, PETSC_DECIDE, NULL, &ilink->Y));
    1382:     }
    1383:     ilink = ilink->next;
    1384:   }
    1385:   PetscFunctionReturn(PETSC_SUCCESS);
    1386: }

    1388: static PetscErrorCode PCMatApply_FieldSplit_Schur(PC pc, Mat X, Mat Y)
    1389: {
    1390:   PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
    1391:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
    1392:   KSP               kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
    1393:   Mat               AinvB = NULL;
    1394:   PetscInt          N, P;

    1396:   PetscFunctionBegin;
    1397:   /* create working matrices with the correct number of columns */
    1398:   PetscCall(PCFieldSplitCreateWorkMats_Private(pc, X));
    1399:   switch (jac->schurfactorization) {
    1400:   case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
    1401:     /* [A00 0; 0 -S], positive definite, suitable for MINRES */
    1402:     PetscCall(MatDenseScatter_Private(ilinkA->sctx, X, ilinkA->X, INSERT_VALUES, SCATTER_FORWARD));
    1403:     PetscCall(MatDenseScatter_Private(ilinkD->sctx, X, ilinkD->X, INSERT_VALUES, SCATTER_FORWARD));
    1404:     PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
    1405:     PetscCall(KSPMatSolve(kspA, ilinkA->X, ilinkA->Y));
    1406:     PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
    1407:     PetscCall(MatDenseScatter_Private(ilinkA->sctx, ilinkA->Y, Y, INSERT_VALUES, SCATTER_REVERSE));
    1408:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL));
    1409:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1410:     PetscCall(KSPMatSolve(jac->kspschur, ilinkD->X, ilinkD->Y));
    1411:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1412:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL));
    1413:     PetscCall(MatScale(ilinkD->Y, jac->schurscale));
    1414:     PetscCall(MatDenseScatter_Private(ilinkD->sctx, ilinkD->Y, Y, INSERT_VALUES, SCATTER_REVERSE));
    1415:     break;
    1416:   case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
    1417:     /* [A00 0; A10 S], suitable for left preconditioning */
    1418:     PetscCall(MatDenseScatter_Private(ilinkA->sctx, X, ilinkA->X, INSERT_VALUES, SCATTER_FORWARD));
    1419:     PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
    1420:     PetscCall(KSPMatSolve(kspA, ilinkA->X, ilinkA->Y));
    1421:     PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
    1422:     PetscCall(MatMatMult(jac->C, ilinkA->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkD->X));
    1423:     PetscCall(MatScale(ilinkD->X, -1.0));
    1424:     PetscCall(MatDenseScatter_Private(ilinkD->sctx, X, ilinkD->X, ADD_VALUES, SCATTER_FORWARD));
    1425:     PetscCall(MatDenseScatter_Private(ilinkA->sctx, ilinkA->Y, Y, INSERT_VALUES, SCATTER_REVERSE));
    1426:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL));
    1427:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1428:     PetscCall(KSPMatSolve(jac->kspschur, ilinkD->X, ilinkD->Y));
    1429:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1430:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL));
    1431:     PetscCall(MatDenseScatter_Private(ilinkD->sctx, ilinkD->Y, Y, INSERT_VALUES, SCATTER_REVERSE));
    1432:     break;
    1433:   case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
    1434:     /* [A00 A01; 0 S], suitable for right preconditioning */
    1435:     PetscCall(MatDenseScatter_Private(ilinkD->sctx, X, ilinkD->X, INSERT_VALUES, SCATTER_FORWARD));
    1436:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL));
    1437:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1438:     PetscCall(KSPMatSolve(jac->kspschur, ilinkD->X, ilinkD->Y));
    1439:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1440:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL));
    1441:     PetscCall(MatMatMult(jac->B, ilinkD->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkA->X));
    1442:     PetscCall(MatScale(ilinkA->X, -1.0));
    1443:     PetscCall(MatDenseScatter_Private(ilinkA->sctx, X, ilinkA->X, ADD_VALUES, SCATTER_FORWARD));
    1444:     PetscCall(MatDenseScatter_Private(ilinkD->sctx, ilinkD->Y, Y, INSERT_VALUES, SCATTER_REVERSE));
    1445:     PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
    1446:     PetscCall(KSPMatSolve(kspA, ilinkA->X, ilinkA->Y));
    1447:     PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
    1448:     PetscCall(MatDenseScatter_Private(ilinkA->sctx, ilinkA->Y, Y, INSERT_VALUES, SCATTER_REVERSE));
    1449:     break;
    1450:   case PC_FIELDSPLIT_SCHUR_FACT_FULL:
    1451:     /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */
    1452:     PetscCall(MatGetSize(jac->B, NULL, &P));
    1453:     N = P;
    1454:     PetscCall(MatDenseScatter_Private(ilinkA->sctx, X, ilinkA->X, INSERT_VALUES, SCATTER_FORWARD));
    1455:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->X, ilinkA->Y, NULL));
    1456:     if (kspUpper == kspA) {
    1457:       PetscCall(PetscObjectQuery((PetscObject)jac->schur, "AinvB", (PetscObject *)&AinvB));
    1458:       if (AinvB) {
    1459:         PetscCall(MatGetSize(AinvB, NULL, &N));
    1460:         if (N > P) { // first time PCApply_FieldSplit_Schur() is called
    1461:           PetscMemType mtype;
    1462:           Mat          C = NULL;
    1463:           PetscScalar *array;
    1464:           PetscInt     m, M, q, Q, p;

    1466:           PetscCall(MatGetSize(jac->B, &M, NULL));
    1467:           PetscCall(MatGetLocalSize(jac->B, &m, NULL));
    1468:           PetscCall(MatGetSize(X, NULL, &Q));
    1469:           PetscCall(MatGetLocalSize(X, NULL, &q));
    1470:           PetscCall(MatDenseGetArrayAndMemType(AinvB, &array, &mtype));
    1471:           if (N != P + Q) {
    1472:             Mat replace;

    1474:             PetscCall(MatGetLocalSize(jac->B, NULL, &p));
    1475:             if (PetscMemTypeCUDA(mtype)) {
    1476: #if PetscDefined(HAVE_CUDA)
    1477:               PetscCallCUDA(cudaFree(array));
    1478:               PetscCallCUDA(cudaMalloc((void **)&array, sizeof(PetscScalar) * m * (P + Q)));
    1479: #endif
    1480:             } else if (PetscMemTypeHIP(mtype)) {
    1481: #if PetscDefined(HAVE_HIP)
    1482:               PetscCallHIP(hipFree(array));
    1483:               PetscCallHIP(hipMalloc((void **)&array, sizeof(PetscScalar) * m * (P + Q)));
    1484: #endif
    1485:             } else {
    1486:               PetscCheck(PetscMemTypeHost(mtype), PetscObjectComm((PetscObject)jac->schur), PETSC_ERR_SUP, "PetscMemType should be either PETSC_MEMTYPE_HOST, PETSC_MEMTYPE_CUDA, or PETSC_MEMTYPE_HIP");
    1487:               PetscCall(PetscFree(array));
    1488:               PetscCall(PetscMalloc1(m * (P + Q), &array));
    1489:             }
    1490:             PetscCall(MatCreateDenseWithMemType(PetscObjectComm((PetscObject)jac->schur), mtype, m, PETSC_DECIDE, M, P + Q, PETSC_DECIDE, array, &replace));
    1491:             PetscCall(MatHeaderReplace(AinvB, &replace));
    1492:           }
    1493:           PetscCall(MatCreateDenseWithMemType(PetscObjectComm((PetscObject)jac->schur), mtype, m, q, M, Q, PETSC_DECIDE, array + m * P, &C));
    1494:           PetscCall(MatDenseRestoreArrayAndMemType(AinvB, &array));
    1495:           PetscCall(MatCopy(ilinkA->X, C, SAME_NONZERO_PATTERN));
    1496:           PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user));
    1497:           PetscCall(KSPSetOperators(jac->kspschur, jac->schur, jac->schur_user));
    1498:           PetscCall(MatCopy(C, ilinkA->Y, SAME_NONZERO_PATTERN)); // retrieve solutions as last columns of the composed Mat
    1499:           PetscCall(MatDestroy(&C));
    1500:         }
    1501:       }
    1502:     }
    1503:     if (N == P) PetscCall(KSPMatSolve(kspLower, ilinkA->X, ilinkA->Y));
    1504:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->X, ilinkA->Y, NULL));
    1505:     PetscCall(MatMatMult(jac->C, ilinkA->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkD->X));
    1506:     PetscCall(MatScale(ilinkD->X, -1.0));
    1507:     PetscCall(MatDenseScatter_Private(ilinkD->sctx, X, ilinkD->X, ADD_VALUES, SCATTER_FORWARD));

    1509:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL));
    1510:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1511:     PetscCall(KSPMatSolve(jac->kspschur, ilinkD->X, ilinkD->Y));
    1512:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1513:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->X, ilinkD->Y, NULL));
    1514:     PetscCall(MatDenseScatter_Private(ilinkD->sctx, ilinkD->Y, Y, INSERT_VALUES, SCATTER_REVERSE));

    1516:     if (kspUpper == kspA) {
    1517:       if (!AinvB) {
    1518:         PetscCall(MatMatMult(jac->B, ilinkD->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkA->Y));
    1519:         PetscCall(MatAXPY(ilinkA->X, -1.0, ilinkA->Y, SAME_NONZERO_PATTERN));
    1520:         PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
    1521:         PetscCall(KSPMatSolve(kspA, ilinkA->X, ilinkA->Y));
    1522:         PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
    1523:       } else {
    1524:         PetscCall(MatMatMult(AinvB, ilinkD->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkA->X));
    1525:         PetscCall(MatAXPY(ilinkA->Y, 1.0, ilinkA->X, SAME_NONZERO_PATTERN));
    1526:       }
    1527:     } else {
    1528:       PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->X, ilinkA->Y, NULL));
    1529:       PetscCall(KSPMatSolve(kspA, ilinkA->X, ilinkA->Y));
    1530:       PetscCall(MatMatMult(jac->B, ilinkD->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilinkA->X));
    1531:       if (!ilinkA->Z) PetscCall(MatDuplicate(ilinkA->X, MAT_DO_NOT_COPY_VALUES, &ilinkA->Z));
    1532:       PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->X, ilinkA->Z, NULL));
    1533:       PetscCall(KSPMatSolve(kspUpper, ilinkA->X, ilinkA->Z));
    1534:       PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->X, ilinkA->Z, NULL));
    1535:       PetscCall(MatAXPY(ilinkA->Y, -1.0, ilinkA->Z, SAME_NONZERO_PATTERN));
    1536:     }
    1537:     PetscCall(MatDenseScatter_Private(ilinkA->sctx, ilinkA->Y, Y, INSERT_VALUES, SCATTER_REVERSE));
    1538:   }
    1539:   PetscFunctionReturn(PETSC_SUCCESS);
    1540: }

    1542: static PetscErrorCode PCApplyTranspose_FieldSplit_Schur(PC pc, Vec x, Vec y)
    1543: {
    1544:   PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
    1545:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
    1546:   KSP               kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;

    1548:   PetscFunctionBegin;
    1549:   switch (jac->schurfactorization) {
    1550:   case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
    1551:     /* [A00 0; 0 -S], positive definite, suitable for MINRES */
    1552:     PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1553:     PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1554:     PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1555:     PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1556:     PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
    1557:     PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1558:     PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1559:     PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1560:     PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1561:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1562:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1563:     PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
    1564:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1565:     PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    1566:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1567:     PetscCall(VecScale(ilinkD->y, jac->schurscale));
    1568:     PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1569:     PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1570:     PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1571:     break;
    1572:   case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
    1573:     PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1574:     PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1575:     PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1576:     PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
    1577:     PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1578:     PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1579:     PetscCall(MatMultTranspose(jac->B, ilinkA->y, ilinkD->x));
    1580:     PetscCall(VecScale(ilinkD->x, -1.));
    1581:     PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    1582:     PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1583:     PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    1584:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1585:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1586:     PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
    1587:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1588:     PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    1589:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1590:     PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1591:     PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1592:     PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1593:     break;
    1594:   case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
    1595:     PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1596:     PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1597:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1598:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1599:     PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
    1600:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1601:     PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    1602:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1603:     PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->x));
    1604:     PetscCall(VecScale(ilinkA->x, -1.));
    1605:     PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
    1606:     PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1607:     PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
    1608:     PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1609:     PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
    1610:     PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1611:     PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1612:     PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1613:     PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1614:     PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1615:     break;
    1616:   case PC_FIELDSPLIT_SCHUR_FACT_FULL:
    1617:     PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1618:     PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1619:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->y, NULL));
    1620:     PetscCall(KSPSolveTranspose(kspUpper, ilinkA->x, ilinkA->y));
    1621:     PetscCall(KSPCheckSolve(kspUpper, pc, ilinkA->y));
    1622:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->y, NULL));
    1623:     PetscCall(MatMultTranspose(jac->B, ilinkA->y, ilinkD->x));
    1624:     PetscCall(VecScale(ilinkD->x, -1.0));
    1625:     PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    1626:     PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));

    1628:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1629:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1630:     PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
    1631:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1632:     PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    1633:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1634:     PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));

    1636:     if (kspLower == kspA) {
    1637:       PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->y));
    1638:       PetscCall(VecAXPY(ilinkA->x, -1.0, ilinkA->y));
    1639:       PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1640:       PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
    1641:       PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1642:       PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1643:     } else {
    1644:       PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1645:       PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
    1646:       PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1647:       PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->x));
    1648:       PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->z, NULL));
    1649:       PetscCall(KSPSolveTranspose(kspLower, ilinkA->x, ilinkA->z));
    1650:       PetscCall(KSPCheckSolve(kspLower, pc, ilinkA->z));
    1651:       PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->z, NULL));
    1652:       PetscCall(VecAXPY(ilinkA->y, -1.0, ilinkA->z));
    1653:     }
    1654:     PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1655:     PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1656:     PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1657:   }
    1658:   PetscFunctionReturn(PETSC_SUCCESS);
    1659: }

    1661: #define FieldSplitSplitSolveAdd(ilink, xx, yy) \
    1662:   ((PetscErrorCode)(VecScatterBegin(ilink->sctx, xx, ilink->x, INSERT_VALUES, SCATTER_FORWARD) || VecScatterEnd(ilink->sctx, xx, ilink->x, INSERT_VALUES, SCATTER_FORWARD) || PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL) || \
    1663:                     KSPSolve(ilink->ksp, ilink->x, ilink->y) || KSPCheckSolve(ilink->ksp, pc, ilink->y) || PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL) || VecScatterBegin(ilink->sctx, ilink->y, yy, ADD_VALUES, SCATTER_REVERSE) || \
    1664:                     VecScatterEnd(ilink->sctx, ilink->y, yy, ADD_VALUES, SCATTER_REVERSE)))

    1666: static PetscErrorCode PCApply_FieldSplit(PC pc, Vec x, Vec y)
    1667: {
    1668:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    1669:   PC_FieldSplitLink ilink = jac->head;
    1670:   PetscInt          cnt, bs;

    1672:   PetscFunctionBegin;
    1673:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
    1674:     PetscBool matnest;

    1676:     PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &matnest));
    1677:     if (jac->defaultsplit && !matnest) {
    1678:       PetscCall(VecGetBlockSize(x, &bs));
    1679:       PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of x vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs);
    1680:       PetscCall(VecGetBlockSize(y, &bs));
    1681:       PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of y vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs);
    1682:       PetscCall(VecStrideGatherAll(x, jac->x, INSERT_VALUES));
    1683:       while (ilink) {
    1684:         PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1685:         PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
    1686:         PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
    1687:         PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1688:         ilink = ilink->next;
    1689:       }
    1690:       PetscCall(VecStrideScatterAll(jac->y, y, INSERT_VALUES));
    1691:     } else {
    1692:       PetscCall(VecSet(y, 0.0));
    1693:       while (ilink) {
    1694:         PetscCall(FieldSplitSplitSolveAdd(ilink, x, y));
    1695:         ilink = ilink->next;
    1696:       }
    1697:     }
    1698:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
    1699:     PetscCall(VecSet(y, 0.0));
    1700:     /* solve on first block for first block variables */
    1701:     PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD));
    1702:     PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD));
    1703:     PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1704:     PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
    1705:     PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
    1706:     PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1707:     PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    1708:     PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));

    1710:     /* compute the residual only onto second block variables using first block variables */
    1711:     PetscCall(MatMult(jac->Afield[1], ilink->y, ilink->next->x));
    1712:     ilink = ilink->next;
    1713:     PetscCall(VecScale(ilink->x, -1.0));
    1714:     PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
    1715:     PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));

    1717:     /* solve on second block variables */
    1718:     PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1719:     PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
    1720:     PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
    1721:     PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1722:     PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    1723:     PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    1724:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
    1725:     if (!jac->w1) {
    1726:       PetscCall(VecDuplicate(x, &jac->w1));
    1727:       PetscCall(VecDuplicate(x, &jac->w2));
    1728:     }
    1729:     PetscCall(VecSet(y, 0.0));
    1730:     PetscCall(FieldSplitSplitSolveAdd(ilink, x, y));
    1731:     cnt = 1;
    1732:     while (ilink->next) {
    1733:       ilink = ilink->next;
    1734:       /* compute the residual only over the part of the vector needed */
    1735:       PetscCall(MatMult(jac->Afield[cnt++], y, ilink->x));
    1736:       PetscCall(VecScale(ilink->x, -1.0));
    1737:       PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
    1738:       PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
    1739:       PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1740:       PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
    1741:       PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
    1742:       PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1743:       PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    1744:       PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    1745:     }
    1746:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
    1747:       cnt -= 2;
    1748:       while (ilink->previous) {
    1749:         ilink = ilink->previous;
    1750:         /* compute the residual only over the part of the vector needed */
    1751:         PetscCall(MatMult(jac->Afield[cnt--], y, ilink->x));
    1752:         PetscCall(VecScale(ilink->x, -1.0));
    1753:         PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
    1754:         PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
    1755:         PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1756:         PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
    1757:         PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
    1758:         PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1759:         PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    1760:         PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    1761:       }
    1762:     }
    1763:   } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Unsupported or unknown composition %d", (int)jac->type);
    1764:   PetscFunctionReturn(PETSC_SUCCESS);
    1765: }

    1767: static PetscErrorCode PCMatApply_FieldSplit(PC pc, Mat X, Mat Y)
    1768: {
    1769:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    1770:   PC_FieldSplitLink ilink = jac->head;
    1771:   PetscInt          cnt;

    1773:   PetscFunctionBegin;
    1774:   /* create working matrices with the correct number of columns */
    1775:   PetscCall(PCFieldSplitCreateWorkMats_Private(pc, X));
    1776:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
    1777:     PetscCall(MatZeroEntries(Y));
    1778:     while (ilink) {
    1779:       PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, INSERT_VALUES, SCATTER_FORWARD));
    1780:       PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    1781:       PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
    1782:       PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    1783:       PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));
    1784:       ilink = ilink->next;
    1785:     }
    1786:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
    1787:     PetscCall(MatZeroEntries(Y));
    1788:     PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, INSERT_VALUES, SCATTER_FORWARD));
    1789:     PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    1790:     PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
    1791:     PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    1792:     PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));

    1794:     /* compute the residual only onto second block variables using first block variables */
    1795:     PetscCall(MatMatMult(jac->Afield[1], ilink->Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilink->next->X));
    1796:     ilink = ilink->next;
    1797:     PetscCall(MatScale(ilink->X, -1.0));
    1798:     PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, ADD_VALUES, SCATTER_FORWARD));

    1800:     /* solve on second block variables */
    1801:     PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    1802:     PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
    1803:     PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    1804:     PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));
    1805:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
    1806:     /* general multiplicative with any number of splits */
    1807:     PetscCall(MatZeroEntries(Y));
    1808:     /* first split */
    1809:     PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, INSERT_VALUES, SCATTER_FORWARD));
    1810:     PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    1811:     PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
    1812:     PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    1813:     PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));
    1814:     cnt = 1;
    1815:     /* forward sweep */
    1816:     while (ilink->next) {
    1817:       ilink = ilink->next;
    1818:       /* compute the residual only over the part of the vector needed */
    1819:       PetscCall(MatMatMult(jac->Afield[cnt++], Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilink->X));
    1820:       PetscCall(MatScale(ilink->X, -1.0));
    1821:       PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, ADD_VALUES, SCATTER_FORWARD));
    1822:       PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    1823:       PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
    1824:       PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    1825:       PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));
    1826:     }
    1827:     /* backward sweep for symmetric multiplicative */
    1828:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
    1829:       cnt -= 2;
    1830:       while (ilink->previous) {
    1831:         ilink = ilink->previous;
    1832:         /* compute the residual only over the part of the vector needed */
    1833:         PetscCall(MatMatMult(jac->Afield[cnt--], Y, MAT_REUSE_MATRIX, PETSC_DETERMINE, &ilink->X));
    1834:         PetscCall(MatScale(ilink->X, -1.0));
    1835:         PetscCall(MatDenseScatter_Private(ilink->sctx, X, ilink->X, ADD_VALUES, SCATTER_FORWARD));
    1836:         PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    1837:         PetscCall(KSPMatSolve(ilink->ksp, ilink->X, ilink->Y));
    1838:         PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->X, ilink->Y, NULL));
    1839:         PetscCall(MatDenseScatter_Private(ilink->sctx, ilink->Y, Y, ADD_VALUES, SCATTER_REVERSE));
    1840:       }
    1841:     }
    1842:   } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PCMatApply() not implemented for this fieldsplit type");
    1843:   PetscFunctionReturn(PETSC_SUCCESS);
    1844: }

    1846: static PetscErrorCode PCApply_FieldSplit_GKB(PC pc, Vec x, Vec y)
    1847: {
    1848:   PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
    1849:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
    1850:   KSP               ksp = ilinkA->ksp;
    1851:   Vec               u, v, Hu, d, work1, work2;
    1852:   PetscScalar       alpha, z, nrmz2, *vecz;
    1853:   PetscReal         lowbnd, nu, beta;
    1854:   PetscInt          j, iterGKB;

    1856:   PetscFunctionBegin;
    1857:   PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1858:   PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1859:   PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1860:   PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));

    1862:   u     = jac->u;
    1863:   v     = jac->v;
    1864:   Hu    = jac->Hu;
    1865:   d     = jac->d;
    1866:   work1 = jac->w1;
    1867:   work2 = jac->w2;
    1868:   vecz  = jac->vecz;

    1870:   /* Change RHS to comply with matrix regularization H = A + nu*B*B' */
    1871:   /* Add q = q + nu*B*b */
    1872:   if (jac->gkbnu) {
    1873:     nu = jac->gkbnu;
    1874:     PetscCall(VecScale(ilinkD->x, jac->gkbnu));
    1875:     PetscCall(MatMultAdd(jac->B, ilinkD->x, ilinkA->x, ilinkA->x)); /* q = q + nu*B*b */
    1876:   } else {
    1877:     /* Situation when no augmented Lagrangian is used. Then we set inner  */
    1878:     /* matrix N = I in [Ar13], and thus nu = 1.                           */
    1879:     nu = 1;
    1880:   }

    1882:   /* Transform rhs from [q,tilde{b}] to [0,b] */
    1883:   PetscCall(PetscLogEventBegin(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL));
    1884:   PetscCall(KSPSolve(ksp, ilinkA->x, ilinkA->y));
    1885:   PetscCall(KSPCheckSolve(ksp, pc, ilinkA->y));
    1886:   PetscCall(PetscLogEventEnd(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL));
    1887:   PetscCall(MatMultHermitianTranspose(jac->B, ilinkA->y, work1));
    1888:   PetscCall(VecAXPBY(work1, 1.0 / nu, -1.0, ilinkD->x)); /* c = b - B'*x        */

    1890:   /* First step of algorithm */
    1891:   PetscCall(VecNorm(work1, NORM_2, &beta)); /* beta = sqrt(nu*c'*c)*/
    1892:   KSPCheckDot(ksp, beta);
    1893:   beta = PetscSqrtReal(nu) * beta;
    1894:   PetscCall(VecAXPBY(v, nu / beta, 0.0, work1)); /* v = nu/beta *c      */
    1895:   PetscCall(MatMult(jac->B, v, work2));          /* u = H^{-1}*B*v      */
    1896:   PetscCall(PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL));
    1897:   PetscCall(KSPSolve(ksp, work2, u));
    1898:   PetscCall(KSPCheckSolve(ksp, pc, u));
    1899:   PetscCall(PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL));
    1900:   PetscCall(MatMult(jac->H, u, Hu)); /* alpha = u'*H*u      */
    1901:   PetscCall(VecDot(Hu, u, &alpha));
    1902:   KSPCheckDot(ksp, alpha);
    1903:   PetscCheck(PetscRealPart(alpha) > 0.0, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "GKB preconditioner diverged, H is not positive definite");
    1904:   alpha = PetscSqrtReal(PetscAbsScalar(alpha));
    1905:   PetscCall(VecScale(u, 1.0 / alpha));
    1906:   PetscCall(VecAXPBY(d, 1.0 / alpha, 0.0, v)); /* v = nu/beta *c      */

    1908:   z       = beta / alpha;
    1909:   vecz[1] = z;

    1911:   /* Computation of first iterate x(1) and p(1) */
    1912:   PetscCall(VecAXPY(ilinkA->y, z, u));
    1913:   PetscCall(VecCopy(d, ilinkD->y));
    1914:   PetscCall(VecScale(ilinkD->y, -z));

    1916:   iterGKB = 1;
    1917:   lowbnd  = 2 * jac->gkbtol;
    1918:   if (jac->gkbmonitor) PetscCall(PetscViewerASCIIPrintf(jac->gkbviewer, "%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n", iterGKB, (double)lowbnd));

    1920:   while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) {
    1921:     iterGKB += 1;
    1922:     PetscCall(MatMultHermitianTranspose(jac->B, u, work1)); /* v <- nu*(B'*u-alpha/nu*v) */
    1923:     PetscCall(VecAXPBY(v, nu, -alpha, work1));
    1924:     PetscCall(VecNorm(v, NORM_2, &beta)); /* beta = sqrt(nu)*v'*v      */
    1925:     beta = beta / PetscSqrtReal(nu);
    1926:     PetscCall(VecScale(v, 1.0 / beta));
    1927:     PetscCall(MatMult(jac->B, v, work2)); /* u <- H^{-1}*(B*v-beta*H*u) */
    1928:     PetscCall(MatMult(jac->H, u, Hu));
    1929:     PetscCall(VecAXPY(work2, -beta, Hu));
    1930:     PetscCall(PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL));
    1931:     PetscCall(KSPSolve(ksp, work2, u));
    1932:     PetscCall(KSPCheckSolve(ksp, pc, u));
    1933:     PetscCall(PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL));
    1934:     PetscCall(MatMult(jac->H, u, Hu)); /* alpha = u'*H*u            */
    1935:     PetscCall(VecDot(Hu, u, &alpha));
    1936:     KSPCheckDot(ksp, alpha);
    1937:     PetscCheck(PetscRealPart(alpha) > 0.0, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "GKB preconditioner diverged, H is not positive definite");
    1938:     alpha = PetscSqrtReal(PetscAbsScalar(alpha));
    1939:     PetscCall(VecScale(u, 1.0 / alpha));

    1941:     z       = -beta / alpha * z; /* z <- beta/alpha*z     */
    1942:     vecz[0] = z;

    1944:     /* Computation of new iterate x(i+1) and p(i+1) */
    1945:     PetscCall(VecAXPBY(d, 1.0 / alpha, -beta / alpha, v)); /* d = (v-beta*d)/alpha */
    1946:     PetscCall(VecAXPY(ilinkA->y, z, u));                   /* r = r + z*u          */
    1947:     PetscCall(VecAXPY(ilinkD->y, -z, d));                  /* p = p - z*d          */
    1948:     PetscCall(MatMult(jac->H, ilinkA->y, Hu));             /* ||u||_H = u'*H*u     */
    1949:     PetscCall(VecDot(Hu, ilinkA->y, &nrmz2));

    1951:     /* Compute Lower Bound estimate */
    1952:     if (iterGKB > jac->gkbdelay) {
    1953:       lowbnd = 0.0;
    1954:       for (j = 0; j < jac->gkbdelay; j++) lowbnd += PetscAbsScalar(vecz[j] * vecz[j]);
    1955:       lowbnd = PetscSqrtReal(lowbnd / PetscAbsScalar(nrmz2));
    1956:     }

    1958:     for (j = 0; j < jac->gkbdelay - 1; j++) vecz[jac->gkbdelay - j - 1] = vecz[jac->gkbdelay - j - 2];
    1959:     if (jac->gkbmonitor) PetscCall(PetscViewerASCIIPrintf(jac->gkbviewer, "%3" PetscInt_FMT " GKB Lower bound estimate %14.12e\n", iterGKB, (double)lowbnd));
    1960:   }

    1962:   PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1963:   PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1964:   PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1965:   PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1966:   PetscFunctionReturn(PETSC_SUCCESS);
    1967: }

    1969: #define FieldSplitSplitSolveAddTranspose(ilink, xx, yy) \
    1970:   ((PetscErrorCode)(VecScatterBegin(ilink->sctx, xx, ilink->y, INSERT_VALUES, SCATTER_FORWARD) || VecScatterEnd(ilink->sctx, xx, ilink->y, INSERT_VALUES, SCATTER_FORWARD) || PetscLogEventBegin(ilink->event, ilink->ksp, ilink->y, ilink->x, NULL) || \
    1971:                     KSPSolveTranspose(ilink->ksp, ilink->y, ilink->x) || KSPCheckSolve(ilink->ksp, pc, ilink->x) || PetscLogEventEnd(ilink->event, ilink->ksp, ilink->y, ilink->x, NULL) || VecScatterBegin(ilink->sctx, ilink->x, yy, ADD_VALUES, SCATTER_REVERSE) || \
    1972:                     VecScatterEnd(ilink->sctx, ilink->x, yy, ADD_VALUES, SCATTER_REVERSE)))

    1974: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc, Vec x, Vec y)
    1975: {
    1976:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    1977:   PC_FieldSplitLink ilink = jac->head;
    1978:   PetscInt          bs;

    1980:   PetscFunctionBegin;
    1981:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
    1982:     PetscBool matnest;

    1984:     PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &matnest));
    1985:     if (jac->defaultsplit && !matnest) {
    1986:       PetscCall(VecGetBlockSize(x, &bs));
    1987:       PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of x vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs);
    1988:       PetscCall(VecGetBlockSize(y, &bs));
    1989:       PetscCheck(jac->bs <= 0 || bs == jac->bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Blocksize of y vector %" PetscInt_FMT " does not match fieldsplit blocksize %" PetscInt_FMT, bs, jac->bs);
    1990:       PetscCall(VecStrideGatherAll(x, jac->x, INSERT_VALUES));
    1991:       while (ilink) {
    1992:         PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1993:         PetscCall(KSPSolveTranspose(ilink->ksp, ilink->x, ilink->y));
    1994:         PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
    1995:         PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1996:         ilink = ilink->next;
    1997:       }
    1998:       PetscCall(VecStrideScatterAll(jac->y, y, INSERT_VALUES));
    1999:     } else {
    2000:       PetscCall(VecSet(y, 0.0));
    2001:       while (ilink) {
    2002:         PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y));
    2003:         ilink = ilink->next;
    2004:       }
    2005:     }
    2006:   } else {
    2007:     if (!jac->w1) {
    2008:       PetscCall(VecDuplicate(x, &jac->w1));
    2009:       PetscCall(VecDuplicate(x, &jac->w2));
    2010:     }
    2011:     PetscCall(VecSet(y, 0.0));
    2012:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
    2013:       PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y));
    2014:       while (ilink->next) {
    2015:         ilink = ilink->next;
    2016:         PetscCall(MatMultTranspose(pc->mat, y, jac->w1));
    2017:         PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x));
    2018:         PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y));
    2019:       }
    2020:       while (ilink->previous) {
    2021:         ilink = ilink->previous;
    2022:         PetscCall(MatMultTranspose(pc->mat, y, jac->w1));
    2023:         PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x));
    2024:         PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y));
    2025:       }
    2026:     } else {
    2027:       while (ilink->next) { /* get to last entry in linked list */
    2028:         ilink = ilink->next;
    2029:       }
    2030:       PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y));
    2031:       while (ilink->previous) {
    2032:         ilink = ilink->previous;
    2033:         PetscCall(MatMultTranspose(pc->mat, y, jac->w1));
    2034:         PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x));
    2035:         PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y));
    2036:       }
    2037:     }
    2038:   }
    2039:   PetscFunctionReturn(PETSC_SUCCESS);
    2040: }

    2042: static PetscErrorCode PCReset_FieldSplit(PC pc)
    2043: {
    2044:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    2045:   PC_FieldSplitLink ilink = jac->head, next;

    2047:   PetscFunctionBegin;
    2048:   while (ilink) {
    2049:     PetscCall(KSPDestroy(&ilink->ksp));
    2050:     PetscCall(VecDestroy(&ilink->x));
    2051:     PetscCall(VecDestroy(&ilink->y));
    2052:     PetscCall(VecDestroy(&ilink->z));
    2053:     PetscCall(MatDestroy(&ilink->X));
    2054:     PetscCall(MatDestroy(&ilink->Y));
    2055:     PetscCall(MatDestroy(&ilink->Z));
    2056:     PetscCall(VecScatterDestroy(&ilink->sctx));
    2057:     PetscCall(ISDestroy(&ilink->is));
    2058:     PetscCall(ISDestroy(&ilink->is_col));
    2059:     PetscCall(PetscFree(ilink->splitname));
    2060:     PetscCall(PetscFree(ilink->fields));
    2061:     PetscCall(PetscFree(ilink->fields_col));
    2062:     next = ilink->next;
    2063:     PetscCall(PetscFree(ilink));
    2064:     ilink = next;
    2065:   }
    2066:   jac->head = NULL;
    2067:   PetscCall(PetscFree2(jac->x, jac->y));
    2068:   if (jac->mat && jac->mat != jac->pmat) {
    2069:     PetscCall(MatDestroyMatrices(jac->nsplits, &jac->mat));
    2070:   } else if (jac->mat) {
    2071:     jac->mat = NULL;
    2072:   }
    2073:   if (jac->pmat) PetscCall(MatDestroyMatrices(jac->nsplits, &jac->pmat));
    2074:   if (jac->Afield) PetscCall(MatDestroyMatrices(jac->nsplits, &jac->Afield));
    2075:   jac->nsplits = 0;
    2076:   PetscCall(VecDestroy(&jac->w1));
    2077:   PetscCall(VecDestroy(&jac->w2));
    2078:   if (jac->schur) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "AinvB", NULL));
    2079:   PetscCall(MatDestroy(&jac->schur));
    2080:   PetscCall(MatDestroy(&jac->schurp));
    2081:   PetscCall(MatDestroy(&jac->schur_user));
    2082:   PetscCall(KSPDestroy(&jac->kspschur));
    2083:   PetscCall(KSPDestroy(&jac->kspupper));
    2084:   PetscCall(MatDestroy(&jac->B));
    2085:   PetscCall(MatDestroy(&jac->C));
    2086:   PetscCall(MatDestroy(&jac->H));
    2087:   PetscCall(VecDestroy(&jac->u));
    2088:   PetscCall(VecDestroy(&jac->v));
    2089:   PetscCall(VecDestroy(&jac->Hu));
    2090:   PetscCall(VecDestroy(&jac->d));
    2091:   PetscCall(PetscFree(jac->vecz));
    2092:   PetscCall(PetscViewerDestroy(&jac->gkbviewer));
    2093:   jac->isrestrict = PETSC_FALSE;
    2094:   PetscFunctionReturn(PETSC_SUCCESS);
    2095: }

    2097: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
    2098: {
    2099:   PetscFunctionBegin;
    2100:   PetscCall(PCReset_FieldSplit(pc));
    2101:   PetscCall(PetscFree(pc->data));
    2102:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
    2103:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", NULL));
    2104:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", NULL));
    2105:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", NULL));
    2106:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", NULL));
    2107:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", NULL));
    2108:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", NULL));
    2109:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL));
    2110:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL));
    2111:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL));
    2112:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL));
    2113:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL));
    2114:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL));
    2115:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL));
    2116:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL));
    2117:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL));
    2118:   PetscFunctionReturn(PETSC_SUCCESS);
    2119: }

    2121: static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc, PetscOptionItems PetscOptionsObject)
    2122: {
    2123:   PetscInt        bs;
    2124:   PetscBool       flg;
    2125:   PC_FieldSplit  *jac = (PC_FieldSplit *)pc->data;
    2126:   PCCompositeType ctype;

    2128:   PetscFunctionBegin;
    2129:   PetscOptionsHeadBegin(PetscOptionsObject, "FieldSplit options");
    2130:   PetscCall(PetscOptionsBool("-pc_fieldsplit_dm_splits", "Whether to use DMCreateFieldDecomposition() for splits", "PCFieldSplitSetDMSplits", jac->dm_splits, &jac->dm_splits, NULL));
    2131:   PetscCall(PetscOptionsInt("-pc_fieldsplit_block_size", "Blocksize that defines number of fields", "PCFieldSplitSetBlockSize", jac->bs, &bs, &flg));
    2132:   if (flg) PetscCall(PCFieldSplitSetBlockSize(pc, bs));
    2133:   jac->diag_use_amat = pc->useAmat;
    2134:   PetscCall(PetscOptionsBool("-pc_fieldsplit_diag_use_amat", "Use Amat (not Pmat) to extract diagonal fieldsplit blocks", "PCFieldSplitSetDiagUseAmat", jac->diag_use_amat, &jac->diag_use_amat, NULL));
    2135:   jac->offdiag_use_amat = pc->useAmat;
    2136:   PetscCall(PetscOptionsBool("-pc_fieldsplit_off_diag_use_amat", "Use Amat (not Pmat) to extract off-diagonal fieldsplit blocks", "PCFieldSplitSetOffDiagUseAmat", jac->offdiag_use_amat, &jac->offdiag_use_amat, NULL));
    2137:   PetscCall(PetscOptionsBool("-pc_fieldsplit_detect_saddle_point", "Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint", jac->detect, &jac->detect, NULL));
    2138:   PetscCall(PCFieldSplitSetDetectSaddlePoint(pc, jac->detect)); /* Sets split type and Schur PC type */
    2139:   PetscCall(PetscOptionsEnum("-pc_fieldsplit_type", "Type of composition", "PCFieldSplitSetType", PCCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&ctype, &flg));
    2140:   if (flg) PetscCall(PCFieldSplitSetType(pc, ctype));
    2141:   /* Only setup fields once */
    2142:   if (jac->bs > 0 && jac->nsplits == 0) {
    2143:     /* only allow user to set fields from command line.
    2144:        otherwise user can set them in PCFieldSplitSetDefaults() */
    2145:     PetscCall(PCFieldSplitSetRuntimeSplits_Private(pc));
    2146:     if (jac->splitdefined) PetscCall(PetscInfo(pc, "Splits defined using the options database\n"));
    2147:   }
    2148:   if (jac->type == PC_COMPOSITE_SCHUR) {
    2149:     PetscCall(PetscOptionsGetEnum(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_schur_factorization_type", PCFieldSplitSchurFactTypes, (PetscEnum *)&jac->schurfactorization, &flg));
    2150:     if (flg) PetscCall(PetscInfo(pc, "Deprecated use of -pc_fieldsplit_schur_factorization_type\n"));
    2151:     PetscCall(PetscOptionsEnum("-pc_fieldsplit_schur_fact_type", "Which off-diagonal parts of the block factorization to use", "PCFieldSplitSetSchurFactType", PCFieldSplitSchurFactTypes, (PetscEnum)jac->schurfactorization, (PetscEnum *)&jac->schurfactorization, NULL));
    2152:     PetscCall(PetscOptionsEnum("-pc_fieldsplit_schur_precondition", "How to build preconditioner for Schur complement", "PCFieldSplitSetSchurPre", PCFieldSplitSchurPreTypes, (PetscEnum)jac->schurpre, (PetscEnum *)&jac->schurpre, NULL));
    2153:     PetscCall(PetscOptionsScalar("-pc_fieldsplit_schur_scale", "Scale Schur complement", "PCFieldSplitSetSchurScale", jac->schurscale, &jac->schurscale, NULL));
    2154:   } else if (jac->type == PC_COMPOSITE_GKB) {
    2155:     PetscCall(PetscOptionsReal("-pc_fieldsplit_gkb_tol", "The tolerance for the lower bound stopping criterion", "PCFieldSplitSetGKBTol", jac->gkbtol, &jac->gkbtol, NULL));
    2156:     PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_delay", "The delay value for lower bound criterion", "PCFieldSplitSetGKBDelay", jac->gkbdelay, &jac->gkbdelay, NULL));
    2157:     PetscCall(PetscOptionsBoundedReal("-pc_fieldsplit_gkb_nu", "Parameter in augmented Lagrangian approach", "PCFieldSplitSetGKBNu", jac->gkbnu, &jac->gkbnu, NULL, 0.0));
    2158:     PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_maxit", "Maximum allowed number of iterations", "PCFieldSplitSetGKBMaxit", jac->gkbmaxit, &jac->gkbmaxit, NULL));
    2159:     PetscCall(PetscOptionsBool("-pc_fieldsplit_gkb_monitor", "Prints number of GKB iterations and error", "PCFieldSplitGKB", jac->gkbmonitor, &jac->gkbmonitor, NULL));
    2160:   }
    2161:   /*
    2162:     In the initial call to this routine the sub-solver data structures do not exist so we cannot call KSPSetFromOptions() on them yet.
    2163:     But after the initial setup of ALL the layers of sub-solvers is completed we do want to call KSPSetFromOptions() on the sub-solvers every time it
    2164:     is called on the outer solver in case changes were made in the options database

    2166:     But even after PCSetUp_FieldSplit() is called all the options inside the inner levels of sub-solvers may still not have been set thus we only call the KSPSetFromOptions()
    2167:     if we know that the entire stack of sub-solvers below this have been complete instantiated, we check this by seeing if any solver iterations are complete.
    2168:     Without this extra check test p2p1fetidp_olof_full and others fail with incorrect matrix types.

    2170:     There could be a negative side effect of calling the KSPSetFromOptions() below.

    2172:     If one captured the PetscObjectState of the options database one could skip these calls if the database has not changed from the previous call
    2173:   */
    2174:   if (jac->issetup) {
    2175:     PC_FieldSplitLink ilink = jac->head;
    2176:     if (jac->type == PC_COMPOSITE_SCHUR) {
    2177:       if (jac->kspupper && jac->kspupper->totalits > 0) PetscCall(KSPSetFromOptions(jac->kspupper));
    2178:       if (jac->kspschur && jac->kspschur->totalits > 0) PetscCall(KSPSetFromOptions(jac->kspschur));
    2179:     }
    2180:     while (ilink) {
    2181:       if (ilink->ksp->totalits > 0) PetscCall(KSPSetFromOptions(ilink->ksp));
    2182:       ilink = ilink->next;
    2183:     }
    2184:   }
    2185:   PetscOptionsHeadEnd();
    2186:   PetscFunctionReturn(PETSC_SUCCESS);
    2187: }

    2189: static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc, const char splitname[], PetscInt n, const PetscInt *fields, const PetscInt *fields_col)
    2190: {
    2191:   PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
    2192:   PC_FieldSplitLink ilink, next = jac->head;
    2193:   char              prefix[128];
    2194:   PetscInt          i;
    2195:   PetscLogEvent     nse;

    2197:   PetscFunctionBegin;
    2198:   if (jac->splitdefined) {
    2199:     PetscCall(PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname));
    2200:     PetscFunctionReturn(PETSC_SUCCESS);
    2201:   }
    2202:   for (i = 0; i < n; i++) PetscCheck(fields[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative field %" PetscInt_FMT " requested", fields[i]);
    2203:   PetscCall(PetscNew(&ilink));
    2204:   if (splitname) {
    2205:     PetscCall(PetscStrallocpy(splitname, &ilink->splitname));
    2206:   } else {
    2207:     PetscCall(PetscMalloc1(3, &ilink->splitname));
    2208:     PetscCall(PetscSNPrintf(ilink->splitname, 2, "%" PetscInt_FMT, jac->nsplits));
    2209:   }
    2210:   PetscCall(PetscMPIIntCast(jac->nsplits, &nse));
    2211:   ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + nse : KSP_Solve_FS_0 + 4; /* Splits greater than 4 logged in 4th split */
    2212:   PetscCall(PetscMalloc1(n, &ilink->fields));
    2213:   PetscCall(PetscArraycpy(ilink->fields, fields, n));
    2214:   PetscCall(PetscMalloc1(n, &ilink->fields_col));
    2215:   PetscCall(PetscArraycpy(ilink->fields_col, fields_col, n));

    2217:   ilink->nfields = n;
    2218:   ilink->next    = NULL;
    2219:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp));
    2220:   PetscCall(KSPSetNestLevel(ilink->ksp, pc->kspnestlevel));
    2221:   PetscCall(KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure));
    2222:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1));
    2223:   PetscCall(KSPSetType(ilink->ksp, KSPPREONLY));

    2225:   PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
    2226:   PetscCall(KSPSetOptionsPrefix(ilink->ksp, prefix));

    2228:   if (!next) {
    2229:     jac->head       = ilink;
    2230:     ilink->previous = NULL;
    2231:   } else {
    2232:     while (next->next) next = next->next;
    2233:     next->next      = ilink;
    2234:     ilink->previous = next;
    2235:   }
    2236:   jac->nsplits++;
    2237:   PetscFunctionReturn(PETSC_SUCCESS);
    2238: }

    2240: static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp)
    2241: {
    2242:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    2244:   PetscFunctionBegin;
    2245:   *subksp = NULL;
    2246:   if (n) *n = 0;
    2247:   if (jac->type == PC_COMPOSITE_SCHUR) {
    2248:     PetscInt nn;

    2250:     PetscCheck(jac->schur, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()");
    2251:     PetscCheck(jac->nsplits == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Unexpected number of splits %" PetscInt_FMT " != 2", jac->nsplits);
    2252:     nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0);
    2253:     PetscCall(PetscMalloc1(nn, subksp));
    2254:     (*subksp)[0] = jac->head->ksp;
    2255:     (*subksp)[1] = jac->kspschur;
    2256:     if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper;
    2257:     if (n) *n = nn;
    2258:   }
    2259:   PetscFunctionReturn(PETSC_SUCCESS);
    2260: }

    2262: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc, PetscInt *n, KSP **subksp)
    2263: {
    2264:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    2266:   PetscFunctionBegin;
    2267:   PetscCheck(jac->schur, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()");
    2268:   PetscCall(PetscMalloc1(jac->nsplits, subksp));
    2269:   PetscCall(MatSchurComplementGetKSP(jac->schur, *subksp));

    2271:   (*subksp)[1] = jac->kspschur;
    2272:   if (n) *n = jac->nsplits;
    2273:   PetscFunctionReturn(PETSC_SUCCESS);
    2274: }

    2276: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp)
    2277: {
    2278:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    2279:   PetscInt          cnt   = 0;
    2280:   PC_FieldSplitLink ilink = jac->head;

    2282:   PetscFunctionBegin;
    2283:   PetscCall(PetscMalloc1(jac->nsplits, subksp));
    2284:   while (ilink) {
    2285:     (*subksp)[cnt++] = ilink->ksp;
    2286:     ilink            = ilink->next;
    2287:   }
    2288:   PetscCheck(cnt == jac->nsplits, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Corrupt PCFIELDSPLIT object: number of splits in linked list %" PetscInt_FMT " does not match number in object %" PetscInt_FMT, cnt, jac->nsplits);
    2289:   if (n) *n = jac->nsplits;
    2290:   PetscFunctionReturn(PETSC_SUCCESS);
    2291: }

    2293: /*@
    2294:   PCFieldSplitRestrictIS - Restricts the fieldsplit `IS`s to be within a given `IS`.

    2296:   Input Parameters:
    2297: + pc  - the preconditioner context
    2298: - isy - the index set that defines the indices to which the fieldsplit is to be restricted

    2300:   Level: advanced

    2302:   Developer Notes:
    2303:   It seems the resulting `IS`s will not cover the entire space, so
    2304:   how can they define a convergent preconditioner? Needs explaining.

    2306: .seealso: [](sec_block_matrices), `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
    2307: @*/
    2308: PetscErrorCode PCFieldSplitRestrictIS(PC pc, IS isy)
    2309: {
    2310:   PetscFunctionBegin;
    2313:   PetscTryMethod(pc, "PCFieldSplitRestrictIS_C", (PC, IS), (pc, isy));
    2314:   PetscFunctionReturn(PETSC_SUCCESS);
    2315: }

    2317: static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
    2318: {
    2319:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    2320:   PC_FieldSplitLink ilink = jac->head, next;
    2321:   PetscInt          localsize, size, sizez, i;
    2322:   const PetscInt   *ind, *indz;
    2323:   PetscInt         *indc, *indcz;
    2324:   PetscBool         flg;

    2326:   PetscFunctionBegin;
    2327:   PetscCall(ISGetLocalSize(isy, &localsize));
    2328:   PetscCallMPI(MPI_Scan(&localsize, &size, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)isy)));
    2329:   size -= localsize;
    2330:   while (ilink) {
    2331:     IS isrl, isr;
    2332:     PC subpc;
    2333:     PetscCall(ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl));
    2334:     PetscCall(ISGetLocalSize(isrl, &localsize));
    2335:     PetscCall(PetscMalloc1(localsize, &indc));
    2336:     PetscCall(ISGetIndices(isrl, &ind));
    2337:     PetscCall(PetscArraycpy(indc, ind, localsize));
    2338:     PetscCall(ISRestoreIndices(isrl, &ind));
    2339:     PetscCall(ISDestroy(&isrl));
    2340:     for (i = 0; i < localsize; i++) *(indc + i) += size;
    2341:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)isy), localsize, indc, PETSC_OWN_POINTER, &isr));
    2342:     PetscCall(PetscObjectReference((PetscObject)isr));
    2343:     PetscCall(ISDestroy(&ilink->is));
    2344:     ilink->is = isr;
    2345:     PetscCall(PetscObjectReference((PetscObject)isr));
    2346:     PetscCall(ISDestroy(&ilink->is_col));
    2347:     ilink->is_col = isr;
    2348:     PetscCall(ISDestroy(&isr));
    2349:     PetscCall(KSPGetPC(ilink->ksp, &subpc));
    2350:     PetscCall(PetscObjectTypeCompare((PetscObject)subpc, PCFIELDSPLIT, &flg));
    2351:     if (flg) {
    2352:       IS       iszl, isz;
    2353:       MPI_Comm comm;
    2354:       PetscCall(ISGetLocalSize(ilink->is, &localsize));
    2355:       comm = PetscObjectComm((PetscObject)ilink->is);
    2356:       PetscCall(ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl));
    2357:       PetscCallMPI(MPI_Scan(&localsize, &sizez, 1, MPIU_INT, MPI_SUM, comm));
    2358:       sizez -= localsize;
    2359:       PetscCall(ISGetLocalSize(iszl, &localsize));
    2360:       PetscCall(PetscMalloc1(localsize, &indcz));
    2361:       PetscCall(ISGetIndices(iszl, &indz));
    2362:       PetscCall(PetscArraycpy(indcz, indz, localsize));
    2363:       PetscCall(ISRestoreIndices(iszl, &indz));
    2364:       PetscCall(ISDestroy(&iszl));
    2365:       for (i = 0; i < localsize; i++) *(indcz + i) += sizez;
    2366:       PetscCall(ISCreateGeneral(comm, localsize, indcz, PETSC_OWN_POINTER, &isz));
    2367:       PetscCall(PCFieldSplitRestrictIS(subpc, isz));
    2368:       PetscCall(ISDestroy(&isz));
    2369:     }
    2370:     next  = ilink->next;
    2371:     ilink = next;
    2372:   }
    2373:   jac->isrestrict = PETSC_TRUE;
    2374:   PetscFunctionReturn(PETSC_SUCCESS);
    2375: }

    2377: static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc, const char splitname[], IS is)
    2378: {
    2379:   PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
    2380:   PC_FieldSplitLink ilink, next = jac->head;
    2381:   char              prefix[128];
    2382:   PetscLogEvent     nse;

    2384:   PetscFunctionBegin;
    2385:   if (jac->splitdefined) {
    2386:     PetscCall(PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname));
    2387:     PetscFunctionReturn(PETSC_SUCCESS);
    2388:   }
    2389:   PetscCall(PetscNew(&ilink));
    2390:   if (splitname) {
    2391:     PetscCall(PetscStrallocpy(splitname, &ilink->splitname));
    2392:   } else {
    2393:     PetscCall(PetscMalloc1(8, &ilink->splitname));
    2394:     PetscCall(PetscSNPrintf(ilink->splitname, 7, "%" PetscInt_FMT, jac->nsplits));
    2395:   }
    2396:   PetscCall(PetscMPIIntCast(jac->nsplits, &nse));
    2397:   ilink->event = jac->nsplits < 5 ? KSP_Solve_FS_0 + nse : KSP_Solve_FS_0 + 4; /* Splits greater than 4 logged in 4th split */
    2398:   PetscCall(PetscObjectReference((PetscObject)is));
    2399:   PetscCall(ISDestroy(&ilink->is));
    2400:   ilink->is = is;
    2401:   PetscCall(PetscObjectReference((PetscObject)is));
    2402:   PetscCall(ISDestroy(&ilink->is_col));
    2403:   ilink->is_col = is;
    2404:   ilink->next   = NULL;
    2405:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp));
    2406:   PetscCall(KSPSetNestLevel(ilink->ksp, pc->kspnestlevel));
    2407:   PetscCall(KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure));
    2408:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1));
    2409:   PetscCall(KSPSetType(ilink->ksp, KSPPREONLY));

    2411:   PetscCall(PetscSNPrintf(prefix, sizeof(prefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
    2412:   PetscCall(KSPSetOptionsPrefix(ilink->ksp, prefix));

    2414:   if (!next) {
    2415:     jac->head       = ilink;
    2416:     ilink->previous = NULL;
    2417:   } else {
    2418:     while (next->next) next = next->next;
    2419:     next->next      = ilink;
    2420:     ilink->previous = next;
    2421:   }
    2422:   jac->nsplits++;
    2423:   PetscFunctionReturn(PETSC_SUCCESS);
    2424: }

    2426: /*@
    2427:   PCFieldSplitSetFields - Sets the fields that define one particular split in `PCFIELDSPLIT`

    2429:   Logically Collective

    2431:   Input Parameters:
    2432: + pc         - the preconditioner context
    2433: . splitname  - name of this split, if `NULL` the number of the split is used
    2434: . n          - the number of fields in this split
    2435: . fields     - the fields in this split
    2436: - fields_col - generally the same as `fields`, if it does not match `fields` then the submatrix that is solved for this set of fields comes from an off-diagonal block
    2437:                of the matrix and `fields_col` provides the column indices for that block

    2439:   Options Database Key:
    2440: . -pc_fieldsplit_%d_fields a,b,... - indicates the fields to be used in the `%d`'th split

    2442:   Level: intermediate

    2444:   Notes:
    2445:   Use `PCFieldSplitSetIS()` to set a  general set of indices as a split.

    2447:   If the matrix used to construct the preconditioner is `MATNEST` then field i refers to the `is_row[i]` `IS` passed to `MatCreateNest()`.

    2449:   If the matrix used to construct the preconditioner is not `MATNEST` then
    2450:   `PCFieldSplitSetFields()` is for defining fields as strided blocks (based on the block size provided to the matrix with `MatSetBlockSize()` or
    2451:   to the `PC` with `PCFieldSplitSetBlockSize()`). For example, if the block
    2452:   size is three then one can define a split as 0, or 1 or 2 or 0,1 or 0,2 or 1,2 which mean
    2453:   0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x23x56x8.. x12x45x78x....
    2454:   where the numbered entries indicate what is in the split.

    2456:   This function is called once per split (it creates a new split each time).  Solve options
    2457:   for this split will be available under the prefix `-fieldsplit_SPLITNAME_`.

    2459:   `PCFieldSplitSetIS()` does not support having a `fields_col` different from `fields`

    2461:   Developer Notes:
    2462:   This routine does not actually create the `IS` representing the split, that is delayed
    2463:   until `PCSetUp_FieldSplit()`, because information about the vector/matrix layouts may not be
    2464:   available when this routine is called.

    2466: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetIS()`, `PCFieldSplitRestrictIS()`,
    2467:           `MatSetBlockSize()`, `MatCreateNest()`
    2468: @*/
    2469: PetscErrorCode PCFieldSplitSetFields(PC pc, const char splitname[], PetscInt n, const PetscInt fields[], const PetscInt fields_col[])
    2470: {
    2471:   PetscFunctionBegin;
    2473:   PetscAssertPointer(splitname, 2);
    2474:   PetscCheck(n >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Provided number of fields %" PetscInt_FMT " in split \"%s\" not positive", n, splitname);
    2475:   PetscAssertPointer(fields, 4);
    2476:   PetscTryMethod(pc, "PCFieldSplitSetFields_C", (PC, const char[], PetscInt, const PetscInt *, const PetscInt *), (pc, splitname, n, fields, fields_col));
    2477:   PetscFunctionReturn(PETSC_SUCCESS);
    2478: }

    2480: /*@
    2481:   PCFieldSplitSetDiagUseAmat - set flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) to build
    2482:   the sub-matrices associated with each split. Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators.

    2484:   Logically Collective

    2486:   Input Parameters:
    2487: + pc  - the preconditioner object
    2488: - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

    2490:   Options Database Key:
    2491: . -pc_fieldsplit_diag_use_amat - use the Amat to provide the diagonal blocks

    2493:   Level: intermediate

    2495: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetDiagUseAmat()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFIELDSPLIT`
    2496: @*/
    2497: PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc, PetscBool flg)
    2498: {
    2499:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
    2500:   PetscBool      isfs;

    2502:   PetscFunctionBegin;
    2504:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
    2505:   PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
    2506:   jac->diag_use_amat = flg;
    2507:   PetscFunctionReturn(PETSC_SUCCESS);
    2508: }

    2510: /*@
    2511:   PCFieldSplitGetDiagUseAmat - get the flag indicating whether to extract diagonal blocks from Amat (rather than Pmat) to build
    2512:   the sub-matrices associated with each split.  Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators.

    2514:   Logically Collective

    2516:   Input Parameter:
    2517: . pc - the preconditioner object

    2519:   Output Parameter:
    2520: . flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

    2522:   Level: intermediate

    2524: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetDiagUseAmat()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFIELDSPLIT`
    2525: @*/
    2526: PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc, PetscBool *flg)
    2527: {
    2528:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
    2529:   PetscBool      isfs;

    2531:   PetscFunctionBegin;
    2533:   PetscAssertPointer(flg, 2);
    2534:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
    2535:   PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
    2536:   *flg = jac->diag_use_amat;
    2537:   PetscFunctionReturn(PETSC_SUCCESS);
    2538: }

    2540: /*@
    2541:   PCFieldSplitSetOffDiagUseAmat - set flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) to build
    2542:   the sub-matrices associated with each split.  Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators.

    2544:   Logically Collective

    2546:   Input Parameters:
    2547: + pc  - the preconditioner object
    2548: - flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from

    2550:   Options Database Key:
    2551: . -pc_fieldsplit_off_diag_use_amat (true|false) - use the Amat to extract the off-diagonal blocks

    2553:   Level: intermediate

    2555: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFieldSplitSetDiagUseAmat()`, `PCFIELDSPLIT`
    2556: @*/
    2557: PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc, PetscBool flg)
    2558: {
    2559:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
    2560:   PetscBool      isfs;

    2562:   PetscFunctionBegin;
    2564:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
    2565:   PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
    2566:   jac->offdiag_use_amat = flg;
    2567:   PetscFunctionReturn(PETSC_SUCCESS);
    2568: }

    2570: /*@
    2571:   PCFieldSplitGetOffDiagUseAmat - get the flag indicating whether to extract off-diagonal blocks from Amat (rather than Pmat) to build
    2572:   the sub-matrices associated with each split.  Where `KSPSetOperators`(ksp,Amat,Pmat) was used to supply the operators.

    2574:   Logically Collective

    2576:   Input Parameter:
    2577: . pc - the preconditioner object

    2579:   Output Parameter:
    2580: . flg - boolean flag indicating whether or not to use Amat to extract the off-diagonal blocks from

    2582:   Level: intermediate

    2584: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFieldSplitGetDiagUseAmat()`, `PCFIELDSPLIT`
    2585: @*/
    2586: PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc, PetscBool *flg)
    2587: {
    2588:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
    2589:   PetscBool      isfs;

    2591:   PetscFunctionBegin;
    2593:   PetscAssertPointer(flg, 2);
    2594:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
    2595:   PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
    2596:   *flg = jac->offdiag_use_amat;
    2597:   PetscFunctionReturn(PETSC_SUCCESS);
    2598: }

    2600: /*@
    2601:   PCFieldSplitSetIS - Sets the exact elements for a split in a `PCFIELDSPLIT`

    2603:   Logically Collective

    2605:   Input Parameters:
    2606: + pc        - the preconditioner context
    2607: . splitname - name of this split, if `NULL` the number of the split is used
    2608: - is        - the index set that defines the elements in this split

    2610:   Level: intermediate

    2612:   Notes:
    2613:   Use `PCFieldSplitSetFields()`, for splits defined by strided `IS` based on the matrix block size or the `is_rows[]` passed into `MATNEST`

    2615:   This function is called once per split (it creates a new split each time).  Solve options
    2616:   for this split will be available under the prefix -fieldsplit_SPLITNAME_.

    2618: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetFields()`
    2619: @*/
    2620: PetscErrorCode PCFieldSplitSetIS(PC pc, const char splitname[], IS is)
    2621: {
    2622:   PetscFunctionBegin;
    2624:   if (splitname) PetscAssertPointer(splitname, 2);
    2626:   PetscTryMethod(pc, "PCFieldSplitSetIS_C", (PC, const char[], IS), (pc, splitname, is));
    2627:   PetscFunctionReturn(PETSC_SUCCESS);
    2628: }

    2630: /*@
    2631:   PCFieldSplitGetIS - Retrieves the elements for a split as an `IS`

    2633:   Logically Collective

    2635:   Input Parameters:
    2636: + pc        - the preconditioner context
    2637: - splitname - name of this split

    2639:   Output Parameter:
    2640: . is - the index set that defines the elements in this split, or `NULL` if the split is not found

    2642:   Level: intermediate

    2644: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetIS()`, `PCFieldSplitGetISByIndex()`
    2645: @*/
    2646: PetscErrorCode PCFieldSplitGetIS(PC pc, const char splitname[], IS *is)
    2647: {
    2648:   PetscFunctionBegin;
    2650:   PetscAssertPointer(splitname, 2);
    2651:   PetscAssertPointer(is, 3);
    2652:   {
    2653:     PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    2654:     PC_FieldSplitLink ilink = jac->head;
    2655:     PetscBool         found;

    2657:     *is = NULL;
    2658:     while (ilink) {
    2659:       PetscCall(PetscStrcmp(ilink->splitname, splitname, &found));
    2660:       if (found) {
    2661:         *is = ilink->is;
    2662:         break;
    2663:       }
    2664:       ilink = ilink->next;
    2665:     }
    2666:   }
    2667:   PetscFunctionReturn(PETSC_SUCCESS);
    2668: }

    2670: /*@
    2671:   PCFieldSplitGetISByIndex - Retrieves the elements for a given split as an `IS`

    2673:   Logically Collective

    2675:   Input Parameters:
    2676: + pc    - the preconditioner context
    2677: - index - index of this split

    2679:   Output Parameter:
    2680: . is - the index set that defines the elements in this split

    2682:   Level: intermediate

    2684: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitGetIS()`, `PCFieldSplitSetIS()`
    2685: @*/
    2686: PetscErrorCode PCFieldSplitGetISByIndex(PC pc, PetscInt index, IS *is)
    2687: {
    2688:   PetscFunctionBegin;
    2689:   PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative field %" PetscInt_FMT " requested", index);
    2691:   PetscAssertPointer(is, 3);
    2692:   {
    2693:     PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    2694:     PC_FieldSplitLink ilink = jac->head;
    2695:     PetscInt          i     = 0;
    2696:     PetscCheck(index < jac->nsplits, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " requested but only %" PetscInt_FMT " exist", index, jac->nsplits);

    2698:     while (i < index) {
    2699:       ilink = ilink->next;
    2700:       ++i;
    2701:     }
    2702:     PetscCall(PCFieldSplitGetIS(pc, ilink->splitname, is));
    2703:   }
    2704:   PetscFunctionReturn(PETSC_SUCCESS);
    2705: }

    2707: /*@
    2708:   PCFieldSplitSetBlockSize - Sets the block size for defining where fields start in the
    2709:   fieldsplit preconditioner when calling `PCFieldSplitSetFields()`. If not set the matrix block size is used.

    2711:   Logically Collective

    2713:   Input Parameters:
    2714: + pc - the preconditioner context
    2715: - bs - the block size

    2717:   Level: intermediate

    2719:   Note:
    2720:   If the matrix is a `MATNEST` then the `is_rows[]` passed to `MatCreateNest()` determines the fields.

    2722: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
    2723: @*/
    2724: PetscErrorCode PCFieldSplitSetBlockSize(PC pc, PetscInt bs)
    2725: {
    2726:   PetscFunctionBegin;
    2729:   PetscTryMethod(pc, "PCFieldSplitSetBlockSize_C", (PC, PetscInt), (pc, bs));
    2730:   PetscFunctionReturn(PETSC_SUCCESS);
    2731: }

    2733: /*@C
    2734:   PCFieldSplitGetSubKSP - Gets the `KSP` contexts for all splits

    2736:   Collective

    2738:   Input Parameter:
    2739: . pc - the preconditioner context

    2741:   Output Parameters:
    2742: + n      - the number of splits
    2743: - subksp - the array of `KSP` contexts

    2745:   Level: advanced

    2747:   Notes:
    2748:   After `PCFieldSplitGetSubKSP()` the array of `KSP`s is to be freed by the user with `PetscFree()`
    2749:   (not the `KSP`, just the array that contains them).

    2751:   You must call `PCSetUp()` before calling `PCFieldSplitGetSubKSP()`.

    2753:   If the fieldsplit is of type `PC_COMPOSITE_SCHUR`, it returns the `KSP` object used inside the
    2754:   Schur complement and the `KSP` object used to iterate over the Schur complement.
    2755:   To access all the `KSP` objects used in `PC_COMPOSITE_SCHUR`, use `PCFieldSplitSchurGetSubKSP()`.

    2757:   If the fieldsplit is of type `PC_COMPOSITE_GKB`, it returns the `KSP` object used to solve the
    2758:   inner linear system defined by the matrix H in each loop.

    2760:   Fortran Note:
    2761:   Call `PCFieldSplitRestoreSubKSP()` when the array of `KSP` is no longer needed

    2763:   Developer Notes:
    2764:   There should be a `PCFieldSplitRestoreSubKSP()` instead of requiring the user to call `PetscFree()`

    2766: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitSchurGetSubKSP()`
    2767: @*/
    2768: PetscErrorCode PCFieldSplitGetSubKSP(PC pc, PetscInt *n, KSP *subksp[])
    2769: {
    2770:   PetscFunctionBegin;
    2772:   if (n) PetscAssertPointer(n, 2);
    2773:   PetscUseMethod(pc, "PCFieldSplitGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp));
    2774:   PetscFunctionReturn(PETSC_SUCCESS);
    2775: }

    2777: /*@C
    2778:   PCFieldSplitSchurGetSubKSP - Gets the `KSP` contexts used inside the Schur complement based `PCFIELDSPLIT`

    2780:   Collective

    2782:   Input Parameter:
    2783: . pc - the preconditioner context

    2785:   Output Parameters:
    2786: + n      - the number of splits
    2787: - subksp - the array of `KSP` contexts

    2789:   Level: advanced

    2791:   Notes:
    2792:   After `PCFieldSplitSchurGetSubKSP()` the array of `KSP`s is to be freed by the user with `PetscFree()`
    2793:   (not the `KSP` just the array that contains them).

    2795:   You must call `PCSetUp()` before calling `PCFieldSplitSchurGetSubKSP()`.

    2797:   If the fieldsplit type is of type `PC_COMPOSITE_SCHUR`, it returns (in order)
    2798: +  1  - the `KSP` used for the (1,1) block
    2799: .  2  - the `KSP` used for the Schur complement (not the one used for the interior Schur solver)
    2800: -  3  - the `KSP` used for the (1,1) block in the upper triangular factor (if different from that of the (1,1) block).

    2802:   It returns a null array if the fieldsplit is not of type `PC_COMPOSITE_SCHUR`; in this case, you should use `PCFieldSplitGetSubKSP()`.

    2804:   Fortran Note:
    2805:   Call `PCFieldSplitSchurRestoreSubKSP()` when the array of `KSP` is no longer needed

    2807:   Developer Notes:
    2808:   There should be a `PCFieldSplitRestoreSubKSP()` instead of requiring the user to call `PetscFree()`

    2810:   Should the functionality of `PCFieldSplitSchurGetSubKSP()` and `PCFieldSplitGetSubKSP()` be merged?

    2812: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitGetSubKSP()`
    2813: @*/
    2814: PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc, PetscInt *n, KSP *subksp[])
    2815: {
    2816:   PetscFunctionBegin;
    2818:   if (n) PetscAssertPointer(n, 2);
    2819:   PetscUseMethod(pc, "PCFieldSplitSchurGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp));
    2820:   PetscFunctionReturn(PETSC_SUCCESS);
    2821: }

    2823: /*@
    2824:   PCFieldSplitSetSchurPre -  Indicates from what operator the preconditioner is constructed for the Schur complement.
    2825:   The default is the A11 matrix.

    2827:   Collective

    2829:   Input Parameters:
    2830: + pc    - the preconditioner context
    2831: . ptype - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11` (default),
    2832:               `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_SCHUR_PRE_USER`,
    2833:               `PC_FIELDSPLIT_SCHUR_PRE_SELFP`, and `PC_FIELDSPLIT_SCHUR_PRE_FULL`
    2834: - pre   - matrix to use for preconditioning, or `NULL`

    2836:   Options Database Keys:
    2837: + -pc_fieldsplit_schur_precondition (self|selfp|user|a11|full) - default is `a11`. See notes for meaning of various arguments
    2838: - -fieldsplit_1_pc_type pctype                                 - the preconditioner algorithm that is used to construct the preconditioner from the operator

    2840:   Level: intermediate

    2842:   Notes:
    2843:   If ptype is
    2844: +     a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
    2845:   matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
    2846: .     self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
    2847:   The only preconditioners that currently work with this symbolic representation matrix object are `PCLSC` and `PCHPDDM`
    2848: .     user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
    2849:   to this function).
    2850: .     selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation $ Sp = A11 - A10 inv(diag(A00)) A01 $
    2851:   This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
    2852:   lumped before extracting the diagonal using the additional option `-fieldsplit_1_mat_schur_complement_ainv_type lump`
    2853: -     full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation
    2854:   computed internally by `PCFIELDSPLIT` (this is expensive)
    2855:   useful mostly as a test that the Schur complement approach can work for your problem

    2857:   When solving a saddle point problem, where the A11 block is identically zero, using `a11` as the ptype only makes sense
    2858:   with the additional option `-fieldsplit_1_pc_type none`. Usually for saddle point problems one would use a `ptype` of `self` and
    2859:   `-fieldsplit_1_pc_type lsc` which uses the least squares commutator to compute a preconditioner for the Schur complement.

    2861:   Developer Note:
    2862:   The name of this function and the option `-pc_fieldsplit_schur_precondition` are inconsistent; precondition should be used everywhere.

    2864: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`,
    2865:           `MatSchurComplementSetAinvType()`, `PCLSC`, `PCFieldSplitSetSchurFactType()`
    2866: @*/
    2867: PetscErrorCode PCFieldSplitSetSchurPre(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
    2868: {
    2869:   PetscFunctionBegin;
    2871:   PetscTryMethod(pc, "PCFieldSplitSetSchurPre_C", (PC, PCFieldSplitSchurPreType, Mat), (pc, ptype, pre));
    2872:   PetscFunctionReturn(PETSC_SUCCESS);
    2873: }

    2875: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
    2876: {
    2877:   return PCFieldSplitSetSchurPre(pc, ptype, pre);
    2878: } /* Deprecated name */

    2880: /*@
    2881:   PCFieldSplitGetSchurPre - For Schur complement fieldsplit, determine how the Schur complement will be
    2882:   preconditioned.  See `PCFieldSplitSetSchurPre()` for details.

    2884:   Logically Collective

    2886:   Input Parameter:
    2887: . pc - the preconditioner context

    2889:   Output Parameters:
    2890: + ptype - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11`, `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_SCHUR_PRE_USER`
    2891: - pre   - matrix to use for preconditioning (with `PC_FIELDSPLIT_SCHUR_PRE_USER`), or `NULL`

    2893:   Level: intermediate

    2895: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCLSC`
    2896: @*/
    2897: PetscErrorCode PCFieldSplitGetSchurPre(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre)
    2898: {
    2899:   PetscFunctionBegin;
    2901:   PetscUseMethod(pc, "PCFieldSplitGetSchurPre_C", (PC, PCFieldSplitSchurPreType *, Mat *), (pc, ptype, pre));
    2902:   PetscFunctionReturn(PETSC_SUCCESS);
    2903: }

    2905: /*@
    2906:   PCFieldSplitSchurGetS -  extract the `MATSCHURCOMPLEMENT` object used by this `PCFIELDSPLIT` in case it needs to be configured separately

    2908:   Not Collective

    2910:   Input Parameter:
    2911: . pc - the preconditioner context

    2913:   Output Parameter:
    2914: . S - the Schur complement matrix

    2916:   Level: advanced

    2918:   Note:
    2919:   This matrix should not be destroyed using `MatDestroy()`; rather, use `PCFieldSplitSchurRestoreS()`.

    2921: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MATSCHURCOMPLEMENT`, `PCFieldSplitSchurRestoreS()`,
    2922:           `MatCreateSchurComplement()`, `MatSchurComplementGetKSP()`, `MatSchurComplementComputeExplicitOperator()`, `MatGetSchurComplement()`
    2923: @*/
    2924: PetscErrorCode PCFieldSplitSchurGetS(PC pc, Mat *S)
    2925: {
    2926:   const char    *t;
    2927:   PetscBool      isfs;
    2928:   PC_FieldSplit *jac;

    2930:   PetscFunctionBegin;
    2932:   PetscCall(PetscObjectGetType((PetscObject)pc, &t));
    2933:   PetscCall(PetscStrcmp(t, PCFIELDSPLIT, &isfs));
    2934:   PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PC of type PCFIELDSPLIT, got %s instead", t);
    2935:   jac = (PC_FieldSplit *)pc->data;
    2936:   PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PCFIELDSPLIT of type SCHUR, got %d instead", jac->type);
    2937:   if (S) *S = jac->schur;
    2938:   PetscFunctionReturn(PETSC_SUCCESS);
    2939: }

    2941: /*@
    2942:   PCFieldSplitSchurRestoreS -  returns the `MATSCHURCOMPLEMENT` matrix used by this `PC`

    2944:   Not Collective

    2946:   Input Parameters:
    2947: + pc - the preconditioner context
    2948: - S  - the Schur complement matrix

    2950:   Level: advanced

    2952: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MatSchurComplement`, `PCFieldSplitSchurGetS()`
    2953: @*/
    2954: PetscErrorCode PCFieldSplitSchurRestoreS(PC pc, Mat *S)
    2955: {
    2956:   const char    *t;
    2957:   PetscBool      isfs;
    2958:   PC_FieldSplit *jac;

    2960:   PetscFunctionBegin;
    2962:   PetscCall(PetscObjectGetType((PetscObject)pc, &t));
    2963:   PetscCall(PetscStrcmp(t, PCFIELDSPLIT, &isfs));
    2964:   PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PC of type PCFIELDSPLIT, got %s instead", t);
    2965:   jac = (PC_FieldSplit *)pc->data;
    2966:   PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PCFIELDSPLIT of type SCHUR, got %d instead", jac->type);
    2967:   PetscCheck(S && (*S == jac->schur), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatSchurComplement restored is not the same as gotten");
    2968:   PetscFunctionReturn(PETSC_SUCCESS);
    2969: }

    2971: static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
    2972: {
    2973:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    2975:   PetscFunctionBegin;
    2976:   jac->schurpre = ptype;
    2977:   if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
    2978:     PetscCall(MatDestroy(&jac->schur_user));
    2979:     jac->schur_user = pre;
    2980:     PetscCall(PetscObjectReference((PetscObject)jac->schur_user));
    2981:   }
    2982:   PetscFunctionReturn(PETSC_SUCCESS);
    2983: }

    2985: static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre)
    2986: {
    2987:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    2989:   PetscFunctionBegin;
    2990:   if (ptype) *ptype = jac->schurpre;
    2991:   if (pre) *pre = jac->schur_user;
    2992:   PetscFunctionReturn(PETSC_SUCCESS);
    2993: }

    2995: /*@
    2996:   PCFieldSplitSetSchurFactType -  sets which blocks of the approximate block factorization to retain in the preconditioner {cite}`murphy2000note` and {cite}`ipsen2001note`

    2998:   Collective

    3000:   Input Parameters:
    3001: + pc    - the preconditioner context
    3002: - ftype - which blocks of factorization to retain, `PC_FIELDSPLIT_SCHUR_FACT_FULL` is default

    3004:   Options Database Key:
    3005: . -pc_fieldsplit_schur_fact_type (diag|lower|upper|full) - default is `full`

    3007:   Level: intermediate

    3009:   Notes:
    3010:   The `full` factorization is

    3012:   ```{math}
    3013:   \left(\begin{array}{cc} A & B \\
    3014:   C & E \\
    3015:   \end{array}\right) =
    3016:   \left(\begin{array}{cc} I & 0 \\
    3017:   C A^{-1} & I \\
    3018:   \end{array}\right)
    3019:   \left(\begin{array}{cc} A & 0 \\
    3020:   0 & S \\
    3021:   \end{array}\right)
    3022:   \left(\begin{array}{cc} I & A^{-1}B \\
    3023:   0 & I \\
    3024:   \end{array}\right) = L D U,
    3025:   ```

    3027:   where $ S = E - C A^{-1} B $. In practice, the full factorization is applied via block triangular solves with the grouping $L(DU)$. `upper` uses $DU$, `lower` uses $LD$,
    3028:   and `diag` is the diagonal part with the sign of $S$ flipped (because this makes the preconditioner positive definite for many formulations,
    3029:   thus allowing the use of `KSPMINRES)`. Sign flipping of $S$ can be turned off with `PCFieldSplitSetSchurScale()`.

    3031:   If $A$ and $S$ are solved exactly
    3032: +  1 - `full` factorization is a direct solver.
    3033: .  2 - The preconditioned operator with `lower` or `upper` has all eigenvalues equal to 1 and minimal polynomial of degree 2, so `KSPGMRES` converges in 2 iterations.
    3034: -  3 - With `diag`, the preconditioned operator has three distinct nonzero eigenvalues and minimal polynomial of degree at most 4, so `KSPGMRES` converges in at most 4 iterations.

    3036:   If the iteration count is very low, consider using `KSPFGMRES` or `KSPGCR` which can use one less preconditioner
    3037:   application in this case. Note that the preconditioned operator may be highly non-normal, so such fast convergence may not be observed in practice.

    3039:   For symmetric problems in which $A$ is positive definite and $S$ is negative definite, `diag` can be used with `KSPMINRES`.

    3041:   A flexible method like `KSPFGMRES` or `KSPGCR`, [](sec_flexibleksp), must be used if the fieldsplit preconditioner is nonlinear (e.g., a few iterations of a Krylov method is used to solve with $A$ or $S$).

    3043: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurScale()`,
    3044:           [](sec_flexibleksp), `PCFieldSplitSetSchurPre()`
    3045: @*/
    3046: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc, PCFieldSplitSchurFactType ftype)
    3047: {
    3048:   PetscFunctionBegin;
    3050:   PetscTryMethod(pc, "PCFieldSplitSetSchurFactType_C", (PC, PCFieldSplitSchurFactType), (pc, ftype));
    3051:   PetscFunctionReturn(PETSC_SUCCESS);
    3052: }

    3054: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc, PCFieldSplitSchurFactType ftype)
    3055: {
    3056:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3058:   PetscFunctionBegin;
    3059:   jac->schurfactorization = ftype;
    3060:   PetscFunctionReturn(PETSC_SUCCESS);
    3061: }

    3063: /*@
    3064:   PCFieldSplitSetSchurScale -  Controls the sign flip of S for `PC_FIELDSPLIT_SCHUR_FACT_DIAG`.

    3066:   Collective

    3068:   Input Parameters:
    3069: + pc    - the preconditioner context
    3070: - scale - scaling factor for the Schur complement

    3072:   Options Database Key:
    3073: . -pc_fieldsplit_schur_scale scale - default is -1.0

    3075:   Level: intermediate

    3077: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurFactType`, `PCFieldSplitSetSchurFactType()`
    3078: @*/
    3079: PetscErrorCode PCFieldSplitSetSchurScale(PC pc, PetscScalar scale)
    3080: {
    3081:   PetscFunctionBegin;
    3084:   PetscTryMethod(pc, "PCFieldSplitSetSchurScale_C", (PC, PetscScalar), (pc, scale));
    3085:   PetscFunctionReturn(PETSC_SUCCESS);
    3086: }

    3088: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc, PetscScalar scale)
    3089: {
    3090:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3092:   PetscFunctionBegin;
    3093:   jac->schurscale = scale;
    3094:   PetscFunctionReturn(PETSC_SUCCESS);
    3095: }

    3097: /*@C
    3098:   PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement

    3100:   Collective

    3102:   Input Parameter:
    3103: . pc - the preconditioner context

    3105:   Output Parameters:
    3106: + A00 - the (0,0) block
    3107: . A01 - the (0,1) block
    3108: . A10 - the (1,0) block
    3109: - A11 - the (1,1) block

    3111:   Level: advanced

    3113:   Note:
    3114:   Use `NULL` for any unneeded output arguments

    3116: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `MatSchurComplementGetSubMatrices()`, `MatSchurComplementSetSubMatrices()`
    3117: @*/
    3118: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc, Mat *A00, Mat *A01, Mat *A10, Mat *A11)
    3119: {
    3120:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3122:   PetscFunctionBegin;
    3124:   PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
    3125:   if (A00) *A00 = jac->pmat[0];
    3126:   if (A01) *A01 = jac->B;
    3127:   if (A10) *A10 = jac->C;
    3128:   if (A11) *A11 = jac->pmat[1];
    3129:   PetscFunctionReturn(PETSC_SUCCESS);
    3130: }

    3132: /*@
    3133:   PCFieldSplitSetGKBTol -  Sets the solver tolerance for the generalized Golub-Kahan bidiagonalization preconditioner {cite}`arioli2013` in `PCFIELDSPLIT`

    3135:   Collective

    3137:   Input Parameters:
    3138: + pc        - the preconditioner context
    3139: - tolerance - the solver tolerance

    3141:   Options Database Key:
    3142: . -pc_fieldsplit_gkb_tol tolerance - default is 1e-5

    3144:   Level: intermediate

    3146:   Note:
    3147:   The generalized GKB algorithm {cite}`arioli2013` uses a lower bound estimate of the error in energy norm as stopping criterion.
    3148:   It stops once the lower bound estimate undershoots the required solver tolerance. Although the actual error might be bigger than
    3149:   this estimate, the stopping criterion is satisfactory in practical cases.

    3151: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBMaxit()`
    3152: @*/
    3153: PetscErrorCode PCFieldSplitSetGKBTol(PC pc, PetscReal tolerance)
    3154: {
    3155:   PetscFunctionBegin;
    3158:   PetscTryMethod(pc, "PCFieldSplitSetGKBTol_C", (PC, PetscReal), (pc, tolerance));
    3159:   PetscFunctionReturn(PETSC_SUCCESS);
    3160: }

    3162: static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc, PetscReal tolerance)
    3163: {
    3164:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3166:   PetscFunctionBegin;
    3167:   jac->gkbtol = tolerance;
    3168:   PetscFunctionReturn(PETSC_SUCCESS);
    3169: }

    3171: /*@
    3172:   PCFieldSplitSetGKBMaxit -  Sets the maximum number of iterations for the generalized Golub-Kahan bidiagonalization preconditioner {cite}`arioli2013` in `PCFIELDSPLIT`

    3174:   Collective

    3176:   Input Parameters:
    3177: + pc    - the preconditioner context
    3178: - maxit - the maximum number of iterations

    3180:   Options Database Key:
    3181: . -pc_fieldsplit_gkb_maxit maxit - default is 100

    3183:   Level: intermediate

    3185: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBNu()`
    3186: @*/
    3187: PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc, PetscInt maxit)
    3188: {
    3189:   PetscFunctionBegin;
    3192:   PetscTryMethod(pc, "PCFieldSplitSetGKBMaxit_C", (PC, PetscInt), (pc, maxit));
    3193:   PetscFunctionReturn(PETSC_SUCCESS);
    3194: }

    3196: static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc, PetscInt maxit)
    3197: {
    3198:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3200:   PetscFunctionBegin;
    3201:   jac->gkbmaxit = maxit;
    3202:   PetscFunctionReturn(PETSC_SUCCESS);
    3203: }

    3205: /*@
    3206:   PCFieldSplitSetGKBDelay -  Sets the delay in the lower bound error estimate in the generalized Golub-Kahan bidiagonalization {cite}`arioli2013` in `PCFIELDSPLIT`
    3207:   preconditioner.

    3209:   Collective

    3211:   Input Parameters:
    3212: + pc    - the preconditioner context
    3213: - delay - the delay window in the lower bound estimate

    3215:   Options Database Key:
    3216: . -pc_fieldsplit_gkb_delay delay - default is 5

    3218:   Level: intermediate

    3220:   Notes:
    3221:   The algorithm uses a lower bound estimate of the error in energy norm as stopping criterion. The lower bound of the error $ ||u-u^k||_H $
    3222:   is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + `delay`), and thus the algorithm needs
    3223:   at least (`delay` + 1) iterations to stop.

    3225:   For more details on the generalized Golub-Kahan bidiagonalization method and its lower bound stopping criterion, please refer to {cite}`arioli2013`

    3227: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()`
    3228: @*/
    3229: PetscErrorCode PCFieldSplitSetGKBDelay(PC pc, PetscInt delay)
    3230: {
    3231:   PetscFunctionBegin;
    3234:   PetscTryMethod(pc, "PCFieldSplitSetGKBDelay_C", (PC, PetscInt), (pc, delay));
    3235:   PetscFunctionReturn(PETSC_SUCCESS);
    3236: }

    3238: static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc, PetscInt delay)
    3239: {
    3240:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3242:   PetscFunctionBegin;
    3243:   jac->gkbdelay = delay;
    3244:   PetscFunctionReturn(PETSC_SUCCESS);
    3245: }

    3247: /*@
    3248:   PCFieldSplitSetGKBNu -  Sets the scalar value nu >= 0 in the transformation H = A00 + nu*A01*A01' of the (1,1) block in the
    3249:   Golub-Kahan bidiagonalization preconditioner {cite}`arioli2013` in `PCFIELDSPLIT`

    3251:   Collective

    3253:   Input Parameters:
    3254: + pc - the preconditioner context
    3255: - nu - the shift parameter

    3257:   Options Database Key:
    3258: . -pc_fieldsplit_gkb_nu nu - default is 1

    3260:   Level: intermediate

    3262:   Notes:
    3263:   This shift is in general done to obtain better convergence properties for the outer loop of the algorithm. This is often achieved by choosing `nu` sufficiently large. However,
    3264:   if `nu` is chosen too large, the matrix H might be badly conditioned and the solution of the linear system $Hx = b$ in the inner loop becomes difficult. It is therefore
    3265:   necessary to find a good balance in between the convergence of the inner and outer loop.

    3267:   For `nu` = 0, no shift is done. In this case A00 has to be positive definite. The matrix N in {cite}`arioli2013` is then chosen as identity.

    3269: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()`
    3270: @*/
    3271: PetscErrorCode PCFieldSplitSetGKBNu(PC pc, PetscReal nu)
    3272: {
    3273:   PetscFunctionBegin;
    3276:   PetscTryMethod(pc, "PCFieldSplitSetGKBNu_C", (PC, PetscReal), (pc, nu));
    3277:   PetscFunctionReturn(PETSC_SUCCESS);
    3278: }

    3280: static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc, PetscReal nu)
    3281: {
    3282:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3284:   PetscFunctionBegin;
    3285:   jac->gkbnu = nu;
    3286:   PetscFunctionReturn(PETSC_SUCCESS);
    3287: }

    3289: static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc, PCCompositeType type)
    3290: {
    3291:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3293:   PetscFunctionBegin;
    3294:   jac->type = type;
    3295:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL));
    3296:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL));
    3297:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL));
    3298:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL));
    3299:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL));
    3300:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL));
    3301:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL));
    3302:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL));
    3303:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL));

    3305:   if (type == PC_COMPOSITE_SCHUR) {
    3306:     pc->ops->apply          = PCApply_FieldSplit_Schur;
    3307:     pc->ops->applytranspose = PCApplyTranspose_FieldSplit_Schur;
    3308:     pc->ops->matapply       = PCMatApply_FieldSplit_Schur;
    3309:     pc->ops->view           = PCView_FieldSplit_Schur;
    3310:     pc->ops->setuponblocks  = PCSetUpOnBlocks_FieldSplit_Schur;

    3312:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit_Schur));
    3313:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", PCFieldSplitSetSchurPre_FieldSplit));
    3314:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", PCFieldSplitGetSchurPre_FieldSplit));
    3315:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", PCFieldSplitSetSchurFactType_FieldSplit));
    3316:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", PCFieldSplitSetSchurScale_FieldSplit));
    3317:   } else if (type == PC_COMPOSITE_GKB) {
    3318:     pc->ops->apply          = PCApply_FieldSplit_GKB;
    3319:     pc->ops->applytranspose = NULL;
    3320:     pc->ops->matapply       = NULL;
    3321:     pc->ops->view           = PCView_FieldSplit_GKB;
    3322:     pc->ops->setuponblocks  = PCSetUpOnBlocks_FieldSplit_GKB;

    3324:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit));
    3325:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", PCFieldSplitSetGKBTol_FieldSplit));
    3326:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", PCFieldSplitSetGKBMaxit_FieldSplit));
    3327:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", PCFieldSplitSetGKBNu_FieldSplit));
    3328:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", PCFieldSplitSetGKBDelay_FieldSplit));
    3329:   } else {
    3330:     pc->ops->apply          = PCApply_FieldSplit;
    3331:     pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
    3332:     pc->ops->matapply       = PCMatApply_FieldSplit;
    3333:     pc->ops->view           = PCView_FieldSplit;
    3334:     pc->ops->setuponblocks  = PCSetUpOnBlocks_FieldSplit;

    3336:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit));
    3337:   }
    3338:   PetscFunctionReturn(PETSC_SUCCESS);
    3339: }

    3341: static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc, PetscInt bs)
    3342: {
    3343:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3345:   PetscFunctionBegin;
    3346:   PetscCheck(bs >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Blocksize must be positive, you gave %" PetscInt_FMT, bs);
    3347:   PetscCheck(jac->bs <= 0 || jac->bs == bs, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Cannot change fieldsplit blocksize from %" PetscInt_FMT " to %" PetscInt_FMT " after it has been set", jac->bs, bs);
    3348:   jac->bs = bs;
    3349:   PetscFunctionReturn(PETSC_SUCCESS);
    3350: }

    3352: static PetscErrorCode PCSetCoordinates_FieldSplit(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
    3353: {
    3354:   PC_FieldSplit    *jac           = (PC_FieldSplit *)pc->data;
    3355:   PC_FieldSplitLink ilink_current = jac->head;
    3356:   IS                is_owned;

    3358:   PetscFunctionBegin;
    3359:   jac->coordinates_set = PETSC_TRUE; // Internal flag
    3360:   PetscCall(MatGetOwnershipIS(pc->mat, &is_owned, NULL));

    3362:   while (ilink_current) {
    3363:     // For each IS, embed it to get local coords indces
    3364:     IS              is_coords;
    3365:     PetscInt        ndofs_block;
    3366:     const PetscInt *block_dofs_enumeration; // Numbering of the dofs relevant to the current block

    3368:     // Setting drop to true for safety. It should make no difference.
    3369:     PetscCall(ISEmbed(ilink_current->is, is_owned, PETSC_TRUE, &is_coords));
    3370:     PetscCall(ISGetLocalSize(is_coords, &ndofs_block));
    3371:     PetscCall(ISGetIndices(is_coords, &block_dofs_enumeration));

    3373:     // Allocate coordinates vector and set it directly
    3374:     PetscCall(PetscMalloc1(ndofs_block * dim, &ilink_current->coords));
    3375:     for (PetscInt dof = 0; dof < ndofs_block; ++dof) {
    3376:       for (PetscInt d = 0; d < dim; ++d) (ilink_current->coords)[dim * dof + d] = coords[dim * block_dofs_enumeration[dof] + d];
    3377:     }
    3378:     ilink_current->dim   = dim;
    3379:     ilink_current->ndofs = ndofs_block;
    3380:     PetscCall(ISRestoreIndices(is_coords, &block_dofs_enumeration));
    3381:     PetscCall(ISDestroy(&is_coords));
    3382:     ilink_current = ilink_current->next;
    3383:   }
    3384:   PetscCall(ISDestroy(&is_owned));
    3385:   PetscFunctionReturn(PETSC_SUCCESS);
    3386: }

    3388: /*@
    3389:   PCFieldSplitSetType - Sets the type, `PCCompositeType`, of a `PCFIELDSPLIT`

    3391:   Collective

    3393:   Input Parameters:
    3394: + pc   - the preconditioner context
    3395: - type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`,
    3396:          `PC_COMPOSITE_GKB`

    3398:   Options Database Key:
    3399: . -pc_fieldsplit_type (multiplicative|additive|symmetric_multiplicative|special|schur) - Sets fieldsplit preconditioner type

    3401:   Level: intermediate

    3403: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCCompositeType`, `PCCompositeGetType()`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`,
    3404:           `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`, `PCFieldSplitSetSchurFactType()`
    3405: @*/
    3406: PetscErrorCode PCFieldSplitSetType(PC pc, PCCompositeType type)
    3407: {
    3408:   PetscFunctionBegin;
    3410:   PetscTryMethod(pc, "PCFieldSplitSetType_C", (PC, PCCompositeType), (pc, type));
    3411:   PetscFunctionReturn(PETSC_SUCCESS);
    3412: }

    3414: /*@
    3415:   PCFieldSplitGetType - Gets the type, `PCCompositeType`, of a `PCFIELDSPLIT`

    3417:   Not collective

    3419:   Input Parameter:
    3420: . pc - the preconditioner context

    3422:   Output Parameter:
    3423: . type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`

    3425:   Level: intermediate

    3427: .seealso: [](sec_block_matrices), `PC`, `PCCompositeSetType()`, `PCFIELDSPLIT`, `PCCompositeType`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`,
    3428:           `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`
    3429: @*/
    3430: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
    3431: {
    3432:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3434:   PetscFunctionBegin;
    3436:   PetscAssertPointer(type, 2);
    3437:   *type = jac->type;
    3438:   PetscFunctionReturn(PETSC_SUCCESS);
    3439: }

    3441: /*@
    3442:   PCFieldSplitSetDMSplits - Flags whether `DMCreateFieldDecomposition()` should be used to define the splits in a `PCFIELDSPLIT`, whenever possible.

    3444:   Logically Collective

    3446:   Input Parameters:
    3447: + pc  - the preconditioner context
    3448: - flg - boolean indicating whether to use field splits defined by the `DM`

    3450:   Options Database Key:
    3451: . -pc_fieldsplit_dm_splits (true|false) - use the field splits defined by the `DM`

    3453:   Level: intermediate

    3455:   Developer Note:
    3456:   The name should be `PCFieldSplitSetUseDMSplits()`, similar change to options database

    3458: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDMSplits()`, `DMCreateFieldDecomposition()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
    3459: @*/
    3460: PetscErrorCode PCFieldSplitSetDMSplits(PC pc, PetscBool flg)
    3461: {
    3462:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
    3463:   PetscBool      isfs;

    3465:   PetscFunctionBegin;
    3468:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
    3469:   if (isfs) jac->dm_splits = flg;
    3470:   PetscFunctionReturn(PETSC_SUCCESS);
    3471: }

    3473: /*@
    3474:   PCFieldSplitGetDMSplits - Returns flag indicating whether `DMCreateFieldDecomposition()` should be used to define the splits in a `PCFIELDSPLIT`, whenever possible.

    3476:   Logically Collective

    3478:   Input Parameter:
    3479: . pc - the preconditioner context

    3481:   Output Parameter:
    3482: . flg - boolean indicating whether to use field splits defined by the `DM`

    3484:   Level: intermediate

    3486:   Developer Note:
    3487:   The name should be `PCFieldSplitGetUseDMSplits()`

    3489: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDMSplits()`, `DMCreateFieldDecomposition()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
    3490: @*/
    3491: PetscErrorCode PCFieldSplitGetDMSplits(PC pc, PetscBool *flg)
    3492: {
    3493:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
    3494:   PetscBool      isfs;

    3496:   PetscFunctionBegin;
    3498:   PetscAssertPointer(flg, 2);
    3499:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
    3500:   if (isfs) {
    3501:     if (flg) *flg = jac->dm_splits;
    3502:   }
    3503:   PetscFunctionReturn(PETSC_SUCCESS);
    3504: }

    3506: /*@
    3507:   PCFieldSplitGetDetectSaddlePoint - Returns flag indicating whether `PCFIELDSPLIT` will attempt to automatically determine fields based on zero diagonal entries.

    3509:   Logically Collective

    3511:   Input Parameter:
    3512: . pc - the preconditioner context

    3514:   Output Parameter:
    3515: . flg - boolean indicating whether to detect fields or not

    3517:   Level: intermediate

    3519: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDetectSaddlePoint()`
    3520: @*/
    3521: PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc, PetscBool *flg)
    3522: {
    3523:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3525:   PetscFunctionBegin;
    3526:   *flg = jac->detect;
    3527:   PetscFunctionReturn(PETSC_SUCCESS);
    3528: }

    3530: /*@
    3531:   PCFieldSplitSetDetectSaddlePoint - Sets flag indicating whether `PCFIELDSPLIT` will attempt to automatically determine fields based on zero diagonal entries.

    3533:   Logically Collective

    3535:   Input Parameter:
    3536: . pc - the preconditioner context

    3538:   Output Parameter:
    3539: . flg - boolean indicating whether to detect fields or not

    3541:   Options Database Key:
    3542: . -pc_fieldsplit_detect_saddle_point (true|false) - detect and use the saddle point

    3544:   Level: intermediate

    3546:   Note:
    3547:   Also sets the split type to `PC_COMPOSITE_SCHUR` (see `PCFieldSplitSetType()`) and the Schur preconditioner type to `PC_FIELDSPLIT_SCHUR_PRE_SELF` (see `PCFieldSplitSetSchurPre()`).

    3549: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDetectSaddlePoint()`, `PCFieldSplitSetType()`, `PCFieldSplitSetSchurPre()`, `PC_FIELDSPLIT_SCHUR_PRE_SELF`
    3550: @*/
    3551: PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc, PetscBool flg)
    3552: {
    3553:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3555:   PetscFunctionBegin;
    3556:   jac->detect = flg;
    3557:   if (jac->detect) {
    3558:     PetscCall(PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR));
    3559:     PetscCall(PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_SELF, NULL));
    3560:   }
    3561:   PetscFunctionReturn(PETSC_SUCCESS);
    3562: }

    3564: /*MC
    3565:   PCFIELDSPLIT - Preconditioner created by combining separate preconditioners for individual
    3566:   collections of variables (that may overlap) called fields or splits. Each field often represents a different continuum variable
    3567:   represented on a grid, such as velocity, pressure, or temperature.
    3568:   In the literature these are sometimes called block preconditioners; but should not be confused with `PCBJACOBI`.
    3569:   See [the users manual section on "Solving Block Matrices"](sec_block_matrices) for more details.

    3571:   Options Database Keys:
    3572: +   -pc_fieldsplit_%d_fields a,b,...                                                 - indicates the fields to be used in the `%d`'th split
    3573: .   -pc_fieldsplit_default                                                           - automatically add any fields to additional splits that have not
    3574:                                                                                        been supplied explicitly by `-pc_fieldsplit_%d_fields`
    3575: .   -pc_fieldsplit_block_size bs                                                     - size of block that defines fields (i.e. there are bs fields)
    3576:                                                                                        when the matrix is not of `MatType` `MATNEST`
    3577: .   -pc_fieldsplit_type (additive|multiplicative|symmetric_multiplicative|schur|gkb) - type of relaxation or factorization splitting
    3578: .   -pc_fieldsplit_schur_precondition (self|selfp|user|a11|full)                     - default is `a11`; see `PCFieldSplitSetSchurPre()`
    3579: .   -pc_fieldsplit_schur_fact_type (diag|lower|upper|full)                           - set factorization type when using `-pc_fieldsplit_type schur`;
    3580:                                                                                        see `PCFieldSplitSetSchurFactType()`
    3581: .   -pc_fieldsplit_dm_splits (true|false) (default is true)                          - Whether to use `DMCreateFieldDecomposition()` for splits
    3582: -   -pc_fieldsplit_detect_saddle_point (true|false)                                  - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver

    3584:   Options prefixes for inner solvers when using the Schur complement preconditioner are `-fieldsplit_0_` and `-fieldsplit_1_` .
    3585:   The options prefix for the inner solver when using the Golub-Kahan biadiagonalization preconditioner is `-fieldsplit_0_`
    3586:   For all other solvers they are `-fieldsplit_%d_` for the `%d`'th field; use `-fieldsplit_` for all fields.

    3588:   To set options on the solvers for all blocks, prepend `-fieldsplit_` to all the `PC`
    3589:   options database keys. For example, `-fieldsplit_pc_type ilu` `-fieldsplit_pc_factor_levels 1`.

    3591:   To set the options on the solvers separate for each block call `PCFieldSplitGetSubKSP()`
    3592:   and set the options directly on the resulting `KSP` object

    3594:   Level: intermediate

    3596:   Notes:
    3597:   Use `PCFieldSplitSetFields()` to set splits defined by "strided" entries or with a `MATNEST` and `PCFieldSplitSetIS()`
    3598:   to define a split by an arbitrary collection of entries.

    3600:   If no splits are set, the default is used. If a `DM` is associated with the `PC` and it supports
    3601:   `DMCreateFieldDecomposition()`, then that is used for the default. Otherwise if the matrix is not `MATNEST`, the splits are defined by entries strided by bs,
    3602:   beginning at 0 then 1, etc to bs-1. The block size can be set with `PCFieldSplitSetBlockSize()`,
    3603:   if this is not called the block size defaults to the blocksize of the second matrix passed
    3604:   to `KSPSetOperators()`/`PCSetOperators()`.

    3606:   For the Schur complement preconditioner if
    3607:   ```{math}
    3608:     J = \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{10} & A_{11} \end{array}\right]
    3609:   ```

    3611:   the preconditioner using `full` factorization is logically
    3612:   ```{math}
    3613:     \left[\begin{array}{cc} I & -\text{ksp}(A_{00}) A_{01} \\ 0 & I \end{array}\right] \left[\begin{array}{cc} \text{ksp}(A_{00}) & 0 \\ 0 & \text{ksp}(S) \end{array}\right] \left[\begin{array}{cc} I & 0 \\ -A_{10} \text{ksp}(A_{00}) & I \end{array}\right]
    3614:       ```
    3615:   where the action of $\text{ksp}(A_{00})$ is applied using the `KSP` solver with prefix `-fieldsplit_0_`.  $S$ is the Schur complement
    3616:   ```{math}
    3617:      S = A_{11} - A_{10} \text{ksp}(A_{00}) A_{01}
    3618:   ```
    3619:   which is usually dense and not stored explicitly.  The action of $\text{ksp}(S)$ is computed using the `KSP` solver with prefix `-fieldsplit_splitname_` (where `splitname`
    3620:   was given in providing the SECOND split or 1 if not given). Accordingly, if using `PCFieldSplitGetSubKSP()`, the array of sub-`KSP` contexts will hold two `KSP`s: at its
    3621:   0th index, the `KSP` associated with `-fieldsplit_0_`, and at its 1st index, the `KSP` corresponding to `-fieldsplit_1_`.
    3622:   By default, $A_{11}$ is used to construct a preconditioner for $S$, use `PCFieldSplitSetSchurPre()` for all the possible ways to construct the preconditioner for $S$.

    3624:   The factorization type is set using `-pc_fieldsplit_schur_fact_type <diag, lower, upper, full>`. `full` is shown above,
    3625:   `diag` gives
    3626:   ```{math}
    3627:     \left[\begin{array}{cc} \text{ksp}(A_{00}) & 0 \\  0 & -\text{ksp}(S) \end{array}\right]
    3628:   ```
    3629:   Note that, slightly counter intuitively, there is a negative in front of the $\text{ksp}(S)$  so that the preconditioner is positive definite. For SPD matrices $J$, the sign flip
    3630:   can be turned off with `PCFieldSplitSetSchurScale()` or by command line `-pc_fieldsplit_schur_scale 1.0`. The `lower` factorization is the inverse of
    3631:   ```{math}
    3632:     \left[\begin{array}{cc} A_{00} & 0 \\  A_{10} & S \end{array}\right]
    3633:   ```
    3634:   where the inverses of $A_{00}$ and $S$ are applied using `KSP`s. The upper factorization is the inverse of
    3635:   ```{math}
    3636:     \left[\begin{array}{cc} A_{00} & A_{01} \\  0 & S \end{array}\right]
    3637:   ```
    3638:   where again the inverses of $A_{00}$ and $S$ are applied using `KSP`s.

    3640:   If only one set of indices (one `IS`) is provided with `PCFieldSplitSetIS()` then the complement of that `IS`
    3641:   is used automatically for a second submatrix.

    3643:   The fieldsplit preconditioner cannot currently be used with the `MATBAIJ` or `MATSBAIJ` data formats if the blocksize is larger than 1.
    3644:   Generally it should be used with the `MATAIJ` or `MATNEST` `MatType`

    3646:   The forms of these preconditioners are closely related, if not identical, to forms derived as "Distributive Iterations", see,
    3647:   for example, page 294 in "Principles of Computational Fluid Dynamics" by Pieter Wesseling {cite}`wesseling2009`.
    3648:   One can also use `PCFIELDSPLIT` inside a smoother resulting in "Distributive Smoothers".

    3650:   See "A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible Navier-Stokes equations" {cite}`elman2008tcp`.

    3652:   The Constrained Pressure Preconditioner (CPR) can be implemented using `PCCOMPOSITE` with `PCGALERKIN`. CPR first solves an $R A P$ subsystem, updates the
    3653:   residual on all variables (`PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE)`), and then applies a simple ILU like preconditioner on all the variables.

    3655:   The generalized Golub-Kahan bidiagonalization preconditioner (GKB) can be applied to symmetric $2 \times 2$ block matrices of the shape
    3656:   ```{math}
    3657:     \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{01}' & 0 \end{array}\right]
    3658:   ```
    3659:   with $A_{00}$ positive semi-definite. The implementation follows {cite}`arioli2013`. Therein, we choose $N := 1/\nu * I$ and the $(1,1)$-block of the matrix is modified to $H = _{A00} + \nu*A_{01}*A_{01}'$.
    3660:   A linear system $Hx = b$ has to be solved in each iteration of the GKB algorithm. This solver is chosen with the option prefix `-fieldsplit_0_`.

    3662:   Some `PCFIELDSPLIT` variants are called physics-based preconditioners, since the preconditioner takes into account the underlying physics of the
    3663:   problem. But this nomenclature is not well-defined.

    3665:   Developer Note:
    3666:   The Schur complement functionality of `PCFIELDSPLIT` should likely be factored into its own `PC` thus simplifying the implementation of the preconditioners and their
    3667:   user API.

    3669: .seealso: [](sec_block_matrices), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCLSC`,
    3670:           `PCFieldSplitGetSubKSP()`, `PCFieldSplitSchurGetSubKSP()`, `PCFieldSplitSetFields()`,
    3671:           `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitSetSchurFactType()`,
    3672:           `MatSchurComplementSetAinvType()`, `PCFieldSplitSetSchurScale()`, `PCFieldSplitSetDetectSaddlePoint()`
    3673: M*/

    3675: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
    3676: {
    3677:   PC_FieldSplit *jac;

    3679:   PetscFunctionBegin;
    3680:   PetscCall(PetscNew(&jac));

    3682:   jac->bs                 = -1;
    3683:   jac->type               = PC_COMPOSITE_MULTIPLICATIVE;
    3684:   jac->schurpre           = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
    3685:   jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
    3686:   jac->schurscale         = -1.0;
    3687:   jac->dm_splits          = PETSC_TRUE;
    3688:   jac->gkbtol             = 1e-5;
    3689:   jac->gkbdelay           = 5;
    3690:   jac->gkbnu              = 1;
    3691:   jac->gkbmaxit           = 100;

    3693:   pc->data = (void *)jac;

    3695:   pc->ops->setup           = PCSetUp_FieldSplit;
    3696:   pc->ops->reset           = PCReset_FieldSplit;
    3697:   pc->ops->destroy         = PCDestroy_FieldSplit;
    3698:   pc->ops->setfromoptions  = PCSetFromOptions_FieldSplit;
    3699:   pc->ops->applyrichardson = NULL;

    3701:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", PCFieldSplitSchurGetSubKSP_FieldSplit));
    3702:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", PCFieldSplitSetFields_FieldSplit));
    3703:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", PCFieldSplitSetIS_FieldSplit));
    3704:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", PCFieldSplitSetType_FieldSplit));
    3705:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", PCFieldSplitSetBlockSize_FieldSplit));
    3706:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", PCFieldSplitRestrictIS_FieldSplit));
    3707:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_FieldSplit));

    3709:   /* Initialize function pointers */
    3710:   PetscCall(PCFieldSplitSetType(pc, jac->type));
    3711:   PetscFunctionReturn(PETSC_SUCCESS);
    3712: }