Actual source code: fieldsplit.c

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

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

     15: 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;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    475:       if (jac->detect) {
    476:         IS       zerodiags, rest;
    477:         PetscInt nmin, nmax;

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

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

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

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

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

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

    552: static PetscErrorCode MatGolubKahanComputeExplicitOperator(Mat A, Mat B, Mat C, Mat *H, PetscReal gkbnu)
    553: {
    554:   Mat       BT, T;
    555:   PetscReal nrmT, nrmB;

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

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

    570:   PetscCall(MatDestroy(&BT));
    571:   PetscCall(MatDestroy(&T));
    572:   PetscFunctionReturn(PETSC_SUCCESS);
    573: }

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

    577: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
    578: {
    579:   PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
    580:   PC_FieldSplitLink ilink;
    581:   PetscInt          i, nsplit;
    582:   PetscBool         sorted, sorted_col, matnest = PETSC_FALSE;

    584:   PetscFunctionBegin;
    585:   pc->failedreason = PC_NOERROR;
    586:   PetscCall(PCFieldSplitSetDefaults(pc));
    587:   nsplit = jac->nsplits;
    588:   ilink  = jac->head;
    589:   if (pc->pmat) PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &matnest));

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

    595:     jac->issetup = PETSC_TRUE;

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

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

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

    611:     if (!matnest) { /* use the matrix blocksize and stride IS to determine the index sets that define the submatrices */
    612:       bs = jac->bs;
    613:       PetscCall(MatGetOwnershipRange(pc->pmat, &rstart, &rend));
    614:       nslots = (rend - rstart) / bs;
    615:       for (i = 0; i < nsplit; i++) {
    616:         if (jac->defaultsplit) {
    617:           PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + i, nsplit, &ilink->is));
    618:           PetscCall(ISDuplicate(ilink->is, &ilink->is_col));
    619:         } else if (!ilink->is) {
    620:           if (ilink->nfields > 1) {
    621:             PetscInt *ii, *jj, j, k, nfields = ilink->nfields, *fields = ilink->fields, *fields_col = ilink->fields_col;

    623:             PetscCall(PetscMalloc1(ilink->nfields * nslots, &ii));
    624:             PetscCall(PetscMalloc1(ilink->nfields * nslots, &jj));
    625:             for (j = 0; j < nslots; j++) {
    626:               for (k = 0; k < nfields; k++) {
    627:                 ii[nfields * j + k] = rstart + bs * j + fields[k];
    628:                 jj[nfields * j + k] = rstart + bs * j + fields_col[k];
    629:               }
    630:             }
    631:             PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, ii, PETSC_OWN_POINTER, &ilink->is));
    632:             PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), nslots * nfields, jj, PETSC_OWN_POINTER, &ilink->is_col));
    633:             PetscCall(ISSetBlockSize(ilink->is, nfields));
    634:             PetscCall(ISSetBlockSize(ilink->is_col, nfields));
    635:           } else {
    636:             PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields[0], bs, &ilink->is));
    637:             PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), nslots, rstart + ilink->fields_col[0], bs, &ilink->is_col));
    638:           }
    639:         }
    640:         PetscCall(ISSorted(ilink->is, &sorted));
    641:         if (ilink->is_col) PetscCall(ISSorted(ilink->is_col, &sorted_col));
    642:         PetscCheck(sorted && sorted_col, PETSC_COMM_SELF, PETSC_ERR_USER, "Fields must be sorted when creating split");
    643:         ilink = ilink->next;
    644:       }
    645:     } else { /* use the IS that define the MATNEST to determine the index sets that define the submatrices */
    646:       IS      *rowis, *colis, *ises = NULL;
    647:       PetscInt mis, nis;

    649:       PetscCall(MatNestGetSize(pc->pmat, &mis, &nis));
    650:       PetscCall(PetscMalloc2(mis, &rowis, nis, &colis));
    651:       PetscCall(MatNestGetISs(pc->pmat, rowis, colis));
    652:       if (!jac->defaultsplit) PetscCall(PetscMalloc1(mis, &ises));

    654:       for (i = 0; i < nsplit; i++) {
    655:         if (jac->defaultsplit) {
    656:           PetscCall(ISDuplicate(rowis[i], &ilink->is));
    657:           PetscCall(ISDuplicate(ilink->is, &ilink->is_col));
    658:         } else if (!ilink->is) {
    659:           if (ilink->nfields > 1) {
    660:             for (PetscInt j = 0; j < ilink->nfields; j++) ises[j] = rowis[ilink->fields[j]];
    661:             PetscCall(ISConcatenate(PetscObjectComm((PetscObject)pc), ilink->nfields, ises, &ilink->is));
    662:           } else {
    663:             PetscCall(ISDuplicate(rowis[ilink->fields[0]], &ilink->is));
    664:           }
    665:           PetscCall(ISDuplicate(ilink->is, &ilink->is_col));
    666:         }
    667:         ilink = ilink->next;
    668:       }
    669:       PetscCall(PetscFree2(rowis, colis));
    670:       PetscCall(PetscFree(ises));
    671:     }
    672:   }

    674:   ilink = jac->head;
    675:   if (!jac->pmat) {
    676:     Vec xtmp;

    678:     PetscCall(MatCreateVecs(pc->pmat, &xtmp, NULL));
    679:     PetscCall(PetscMalloc1(nsplit, &jac->pmat));
    680:     PetscCall(PetscMalloc2(nsplit, &jac->x, nsplit, &jac->y));
    681:     for (i = 0; i < nsplit; i++) {
    682:       MatNullSpace sp;

    684:       /* Check for preconditioning matrix attached to IS */
    685:       PetscCall(PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&jac->pmat[i]));
    686:       if (jac->pmat[i]) {
    687:         PetscCall(PetscObjectReference((PetscObject)jac->pmat[i]));
    688:         if (jac->type == PC_COMPOSITE_SCHUR) {
    689:           jac->schur_user = jac->pmat[i];

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

    720:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
    721:       PetscCall(MatGetNullSpaces(nsplit, jac->pmat, &nullsp));
    722:       for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->pmat[i]));
    723:       scall = MAT_INITIAL_MATRIX;
    724:     } else scall = MAT_REUSE_MATRIX;

    726:     for (i = 0; i < nsplit; i++) {
    727:       Mat pmat;

    729:       /* Check for preconditioning matrix attached to IS */
    730:       PetscCall(PetscObjectQuery((PetscObject)ilink->is, "pmat", (PetscObject *)&pmat));
    731:       if (!pmat) PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ilink->is_col, scall, &jac->pmat[i]));
    732:       ilink = ilink->next;
    733:     }
    734:     if (nullsp) PetscCall(MatRestoreNullSpaces(nsplit, jac->pmat, &nullsp));
    735:   }
    736:   if (jac->diag_use_amat) {
    737:     ilink = jac->head;
    738:     if (!jac->mat) {
    739:       PetscCall(PetscMalloc1(nsplit, &jac->mat));
    740:       for (i = 0; i < nsplit; i++) {
    741:         PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, MAT_INITIAL_MATRIX, &jac->mat[i]));
    742:         ilink = ilink->next;
    743:       }
    744:     } else {
    745:       MatReuse      scall;
    746:       MatNullSpace *nullsp = NULL;

    748:       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
    749:         PetscCall(MatGetNullSpaces(nsplit, jac->mat, &nullsp));
    750:         for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->mat[i]));
    751:         scall = MAT_INITIAL_MATRIX;
    752:       } else scall = MAT_REUSE_MATRIX;

    754:       for (i = 0; i < nsplit; i++) {
    755:         PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ilink->is_col, scall, &jac->mat[i]));
    756:         ilink = ilink->next;
    757:       }
    758:       if (nullsp) PetscCall(MatRestoreNullSpaces(nsplit, jac->mat, &nullsp));
    759:     }
    760:   } else {
    761:     jac->mat = jac->pmat;
    762:   }

    764:   /* Check for null space attached to IS */
    765:   ilink = jac->head;
    766:   for (i = 0; i < nsplit; i++) {
    767:     MatNullSpace sp;

    769:     PetscCall(PetscObjectQuery((PetscObject)ilink->is, "nullspace", (PetscObject *)&sp));
    770:     if (sp) PetscCall(MatSetNullSpace(jac->mat[i], sp));
    771:     ilink = ilink->next;
    772:   }

    774:   if (jac->type != PC_COMPOSITE_ADDITIVE && jac->type != PC_COMPOSITE_SCHUR && jac->type != PC_COMPOSITE_GKB) {
    775:     /* extract the rows of the matrix associated with each field: used for efficient computation of residual inside algorithm */
    776:     /* FIXME: Can/should we reuse jac->mat whenever (jac->diag_use_amat) is true? */
    777:     ilink = jac->head;
    778:     if (nsplit == 2 && jac->type == PC_COMPOSITE_MULTIPLICATIVE) {
    779:       /* 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 */
    780:       if (!jac->Afield) {
    781:         PetscCall(PetscCalloc1(nsplit, &jac->Afield));
    782:         if (jac->offdiag_use_amat) {
    783:           PetscCall(MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1]));
    784:         } else {
    785:           PetscCall(MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, MAT_INITIAL_MATRIX, &jac->Afield[1]));
    786:         }
    787:       } else {
    788:         MatReuse scall;

    790:         if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
    791:           PetscCall(MatDestroy(&jac->Afield[1]));
    792:           scall = MAT_INITIAL_MATRIX;
    793:         } else scall = MAT_REUSE_MATRIX;

    795:         if (jac->offdiag_use_amat) {
    796:           PetscCall(MatCreateSubMatrix(pc->mat, ilink->next->is, ilink->is, scall, &jac->Afield[1]));
    797:         } else {
    798:           PetscCall(MatCreateSubMatrix(pc->pmat, ilink->next->is, ilink->is, scall, &jac->Afield[1]));
    799:         }
    800:       }
    801:     } else {
    802:       if (!jac->Afield) {
    803:         PetscCall(PetscMalloc1(nsplit, &jac->Afield));
    804:         for (i = 0; i < nsplit; i++) {
    805:           if (jac->offdiag_use_amat) {
    806:             PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i]));
    807:           } else {
    808:             PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, NULL, MAT_INITIAL_MATRIX, &jac->Afield[i]));
    809:           }
    810:           ilink = ilink->next;
    811:         }
    812:       } else {
    813:         MatReuse scall;
    814:         if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
    815:           for (i = 0; i < nsplit; i++) PetscCall(MatDestroy(&jac->Afield[i]));
    816:           scall = MAT_INITIAL_MATRIX;
    817:         } else scall = MAT_REUSE_MATRIX;

    819:         for (i = 0; i < nsplit; i++) {
    820:           if (jac->offdiag_use_amat) {
    821:             PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, NULL, scall, &jac->Afield[i]));
    822:           } else {
    823:             PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, NULL, scall, &jac->Afield[i]));
    824:           }
    825:           ilink = ilink->next;
    826:         }
    827:       }
    828:     }
    829:   }

    831:   if (jac->type == PC_COMPOSITE_SCHUR) {
    832:     IS          ccis;
    833:     PetscBool   isset, isspd;
    834:     PetscInt    rstart, rend;
    835:     char        lscname[256];
    836:     PetscObject LSC_L;
    837:     PetscBool   set, flg;

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

    841:     /* If pc->mat is SPD, don't scale by -1 the Schur complement */
    842:     if (jac->schurscale == (PetscScalar)-1.0) {
    843:       PetscCall(MatIsSPDKnown(pc->pmat, &isset, &isspd));
    844:       jac->schurscale = (isset && isspd) ? 1.0 : -1.0;
    845:     }

    847:     /* When extracting off-diagonal submatrices, we take complements from this range */
    848:     PetscCall(MatGetOwnershipRangeColumn(pc->mat, &rstart, &rend));
    849:     PetscCall(PetscObjectTypeCompareAny(jac->offdiag_use_amat ? (PetscObject)pc->mat : (PetscObject)pc->pmat, &flg, MATSEQSBAIJ, MATMPISBAIJ, ""));

    851:     if (jac->schur) {
    852:       KSP      kspA = jac->head->ksp, kspInner = NULL, kspUpper = jac->kspupper;
    853:       MatReuse scall;

    855:       if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
    856:         scall = MAT_INITIAL_MATRIX;
    857:         PetscCall(MatDestroy(&jac->B));
    858:         PetscCall(MatDestroy(&jac->C));
    859:       } else scall = MAT_REUSE_MATRIX;

    861:       PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner));
    862:       ilink = jac->head;
    863:       PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis));
    864:       if (jac->offdiag_use_amat) {
    865:         PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, scall, &jac->B));
    866:       } else {
    867:         PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, scall, &jac->B));
    868:       }
    869:       PetscCall(ISDestroy(&ccis));
    870:       if (!flg) {
    871:         ilink = ilink->next;
    872:         PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis));
    873:         if (jac->offdiag_use_amat) {
    874:           PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, scall, &jac->C));
    875:         } else {
    876:           PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, scall, &jac->C));
    877:         }
    878:         PetscCall(ISDestroy(&ccis));
    879:       } else {
    880:         PetscCall(MatIsHermitianKnown(jac->offdiag_use_amat ? pc->mat : pc->pmat, &set, &flg));
    881:         if (set && flg) PetscCall(MatCreateHermitianTranspose(jac->B, &jac->C));
    882:         else PetscCall(MatCreateTranspose(jac->B, &jac->C));
    883:       }
    884:       PetscCall(MatSchurComplementUpdateSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1]));
    885:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) {
    886:         PetscCall(MatDestroy(&jac->schurp));
    887:         PetscCall(MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp));
    888:       } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
    889:         PetscCall(MatDestroy(&jac->schur_user));
    890:         if (jac->kspupper == jac->head->ksp) {
    891:           Mat AinvB;

    893:           PetscCall(MatCreate(PetscObjectComm((PetscObject)jac->schur), &AinvB));
    894:           PetscCall(PetscObjectCompose((PetscObject)jac->schur, "AinvB", (PetscObject)AinvB));
    895:           PetscCall(MatDestroy(&AinvB));
    896:         }
    897:         PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user));
    898:       }
    899:       if (kspA != kspInner) PetscCall(KSPSetOperators(kspA, jac->mat[0], jac->pmat[0]));
    900:       if (kspUpper != kspA) PetscCall(KSPSetOperators(kspUpper, jac->mat[0], jac->pmat[0]));
    901:       PetscCall(KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac)));
    902:     } else {
    903:       const char  *Dprefix;
    904:       char         schurprefix[256], schurmatprefix[256];
    905:       char         schurtestoption[256];
    906:       MatNullSpace sp;
    907:       KSP          kspt;

    909:       /* extract the A01 and A10 matrices */
    910:       ilink = jac->head;
    911:       PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis));
    912:       if (jac->offdiag_use_amat) {
    913:         PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B));
    914:       } else {
    915:         PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B));
    916:       }
    917:       PetscCall(ISDestroy(&ccis));
    918:       ilink = ilink->next;
    919:       if (!flg) {
    920:         PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis));
    921:         if (jac->offdiag_use_amat) {
    922:           PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C));
    923:         } else {
    924:           PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C));
    925:         }
    926:         PetscCall(ISDestroy(&ccis));
    927:       } else {
    928:         PetscCall(MatIsHermitianKnown(jac->offdiag_use_amat ? pc->mat : pc->pmat, &set, &flg));
    929:         if (set && flg) PetscCall(MatCreateHermitianTranspose(jac->B, &jac->C));
    930:         else PetscCall(MatCreateTranspose(jac->B, &jac->C));
    931:       }
    932:       /* Use mat[0] (diagonal block of Amat) preconditioned by pmat[0] to define Schur complement */
    933:       PetscCall(MatCreate(((PetscObject)jac->mat[0])->comm, &jac->schur));
    934:       PetscCall(MatSetType(jac->schur, MATSCHURCOMPLEMENT));
    935:       PetscCall(MatSchurComplementSetSubMatrices(jac->schur, jac->mat[0], jac->pmat[0], jac->B, jac->C, jac->mat[1]));
    936:       PetscCall(PetscSNPrintf(schurmatprefix, sizeof(schurmatprefix), "%sfieldsplit_%s_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
    937:       PetscCall(MatSetOptionsPrefix(jac->schur, schurmatprefix));
    938:       PetscCall(MatSchurComplementGetKSP(jac->schur, &kspt));
    939:       PetscCall(KSPSetOptionsPrefix(kspt, schurmatprefix));

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

    945:       PetscCall(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_inner_", ilink->splitname));
    946:       PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, NULL, &flg));
    947:       if (flg) {
    948:         DM  dmInner;
    949:         KSP kspInner;
    950:         PC  pcInner;

    952:         PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner));
    953:         PetscCall(KSPReset(kspInner));
    954:         PetscCall(KSPSetOperators(kspInner, jac->mat[0], jac->pmat[0]));
    955:         PetscCall(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_inner_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
    956:         /* Indent this deeper to emphasize the "inner" nature of this solver. */
    957:         PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspInner, (PetscObject)pc, 2));
    958:         PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspInner->pc, (PetscObject)pc, 2));
    959:         PetscCall(KSPSetOptionsPrefix(kspInner, schurprefix));

    961:         /* Set DM for new solver */
    962:         PetscCall(KSPGetDM(jac->head->ksp, &dmInner));
    963:         PetscCall(KSPSetDM(kspInner, dmInner));
    964:         PetscCall(KSPSetDMActive(kspInner, PETSC_FALSE));

    966:         /* Defaults to PCKSP as preconditioner */
    967:         PetscCall(KSPGetPC(kspInner, &pcInner));
    968:         PetscCall(PCSetType(pcInner, PCKSP));
    969:         PetscCall(PCKSPSetKSP(pcInner, jac->head->ksp));
    970:       } else {
    971:         /* Use the outer solver for the inner solve, but revert the KSPPREONLY from PCFieldSplitSetFields_FieldSplit or
    972:           * PCFieldSplitSetIS_FieldSplit. We don't want KSPPREONLY because it makes the Schur complement inexact,
    973:           * preventing Schur complement reduction to be an accurate solve. Usually when an iterative solver is used for
    974:           * S = D - C A_inner^{-1} B, we expect S to be defined using an accurate definition of A_inner^{-1}, so we make
    975:           * GMRES the default. Note that it is also common to use PREONLY for S, in which case S may not be used
    976:           * directly, and the user is responsible for setting an inexact method for fieldsplit's A^{-1}. */
    977:         PetscCall(KSPSetType(jac->head->ksp, KSPGMRES));
    978:         PetscCall(MatSchurComplementSetKSP(jac->schur, jac->head->ksp));
    979:       }
    980:       PetscCall(KSPSetOperators(jac->head->ksp, jac->mat[0], jac->pmat[0]));
    981:       PetscCall(KSPSetFromOptions(jac->head->ksp));
    982:       PetscCall(MatSetFromOptions(jac->schur));

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

    989:         PetscCall(MatSchurComplementGetKSP(jac->schur, &kspInner));
    990:         PetscCall(KSPGetPC(kspInner, &pcInner));
    991:         PetscCall(PetscObjectTypeCompare((PetscObject)pcInner, PCKSP, &flg));
    992:         if (flg) {
    993:           KSP ksp;

    995:           PetscCall(PCKSPGetKSP(pcInner, &ksp));
    996:           if (ksp == jac->head->ksp) PetscCall(PCSetUseAmat(pcInner, PETSC_TRUE));
    997:         }
    998:       }
    999:       PetscCall(PetscSNPrintf(schurtestoption, sizeof(schurtestoption), "-fieldsplit_%s_upper_", ilink->splitname));
    1000:       PetscCall(PetscOptionsFindPairPrefix_Private(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, schurtestoption, NULL, NULL, &flg));
    1001:       if (flg) {
    1002:         DM dmInner;

    1004:         PetscCall(PetscSNPrintf(schurprefix, sizeof(schurprefix), "%sfieldsplit_%s_upper_", ((PetscObject)pc)->prefix ? ((PetscObject)pc)->prefix : "", ilink->splitname));
    1005:         PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspupper));
    1006:         PetscCall(KSPSetNestLevel(jac->kspupper, pc->kspnestlevel));
    1007:         PetscCall(KSPSetErrorIfNotConverged(jac->kspupper, pc->erroriffailure));
    1008:         PetscCall(KSPSetOptionsPrefix(jac->kspupper, schurprefix));
    1009:         PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper, (PetscObject)pc, 1));
    1010:         PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspupper->pc, (PetscObject)pc, 1));
    1011:         PetscCall(KSPGetDM(jac->head->ksp, &dmInner));
    1012:         PetscCall(KSPSetDM(jac->kspupper, dmInner));
    1013:         PetscCall(KSPSetDMActive(jac->kspupper, PETSC_FALSE));
    1014:         PetscCall(KSPSetFromOptions(jac->kspupper));
    1015:         PetscCall(KSPSetOperators(jac->kspupper, jac->mat[0], jac->pmat[0]));
    1016:         PetscCall(VecDuplicate(jac->head->x, &jac->head->z));
    1017:       } else {
    1018:         jac->kspupper = jac->head->ksp;
    1019:         PetscCall(PetscObjectReference((PetscObject)jac->head->ksp));
    1020:       }

    1022:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELFP) PetscCall(MatSchurComplementGetPmat(jac->schur, MAT_INITIAL_MATRIX, &jac->schurp));
    1023:       PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->kspschur));
    1024:       PetscCall(KSPSetNestLevel(jac->kspschur, pc->kspnestlevel));
    1025:       PetscCall(KSPSetErrorIfNotConverged(jac->kspschur, pc->erroriffailure));
    1026:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->kspschur, (PetscObject)pc, 1));
    1027:       if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_SELF) {
    1028:         PC pcschur;
    1029:         PetscCall(KSPGetPC(jac->kspschur, &pcschur));
    1030:         PetscCall(PCSetType(pcschur, PCNONE));
    1031:         /* Note: This is bad if there exist preconditioners for MATSCHURCOMPLEMENT */
    1032:       } else if (jac->schurpre == PC_FIELDSPLIT_SCHUR_PRE_FULL) {
    1033:         if (jac->schurfactorization == PC_FIELDSPLIT_SCHUR_FACT_FULL && jac->kspupper == jac->head->ksp) {
    1034:           Mat AinvB;

    1036:           PetscCall(MatCreate(PetscObjectComm((PetscObject)jac->schur), &AinvB));
    1037:           PetscCall(PetscObjectCompose((PetscObject)jac->schur, "AinvB", (PetscObject)AinvB));
    1038:           PetscCall(MatDestroy(&AinvB));
    1039:         }
    1040:         PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user));
    1041:       }
    1042:       PetscCall(KSPSetOperators(jac->kspschur, jac->schur, FieldSplitSchurPre(jac)));
    1043:       PetscCall(KSPGetOptionsPrefix(jac->head->next->ksp, &Dprefix));
    1044:       PetscCall(KSPSetOptionsPrefix(jac->kspschur, Dprefix));
    1045:       /* propagate DM */
    1046:       {
    1047:         DM sdm;
    1048:         PetscCall(KSPGetDM(jac->head->next->ksp, &sdm));
    1049:         if (sdm) {
    1050:           PetscCall(KSPSetDM(jac->kspschur, sdm));
    1051:           PetscCall(KSPSetDMActive(jac->kspschur, PETSC_FALSE));
    1052:         }
    1053:       }
    1054:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
    1055:       /* need to call this every time, since the jac->kspschur is freshly created, otherwise its options never get set */
    1056:       PetscCall(KSPSetFromOptions(jac->kspschur));
    1057:     }
    1058:     PetscCall(MatAssemblyBegin(jac->schur, MAT_FINAL_ASSEMBLY));
    1059:     PetscCall(MatAssemblyEnd(jac->schur, MAT_FINAL_ASSEMBLY));

    1061:     /* HACK: special support to forward L and Lp matrices that might be used by PCLSC */
    1062:     PetscCall(PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_L", ilink->splitname));
    1063:     PetscCall(PetscObjectQuery((PetscObject)pc->mat, lscname, (PetscObject *)&LSC_L));
    1064:     if (!LSC_L) PetscCall(PetscObjectQuery((PetscObject)pc->pmat, lscname, (PetscObject *)&LSC_L));
    1065:     if (LSC_L) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "LSC_L", (PetscObject)LSC_L));
    1066:     PetscCall(PetscSNPrintf(lscname, sizeof(lscname), "%s_LSC_Lp", ilink->splitname));
    1067:     PetscCall(PetscObjectQuery((PetscObject)pc->pmat, lscname, (PetscObject *)&LSC_L));
    1068:     if (!LSC_L) PetscCall(PetscObjectQuery((PetscObject)pc->mat, lscname, (PetscObject *)&LSC_L));
    1069:     if (LSC_L) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "LSC_Lp", (PetscObject)LSC_L));
    1070:   } else if (jac->type == PC_COMPOSITE_GKB) {
    1071:     IS       ccis;
    1072:     PetscInt rstart, rend;

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

    1076:     ilink = jac->head;

    1078:     /* When extracting off-diagonal submatrices, we take complements from this range */
    1079:     PetscCall(MatGetOwnershipRangeColumn(pc->mat, &rstart, &rend));

    1081:     PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis));
    1082:     if (jac->offdiag_use_amat) {
    1083:       PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B));
    1084:     } else {
    1085:       PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->B));
    1086:     }
    1087:     PetscCall(ISDestroy(&ccis));
    1088:     /* Create work vectors for GKB algorithm */
    1089:     PetscCall(VecDuplicate(ilink->x, &jac->u));
    1090:     PetscCall(VecDuplicate(ilink->x, &jac->Hu));
    1091:     PetscCall(VecDuplicate(ilink->x, &jac->w2));
    1092:     ilink = ilink->next;
    1093:     PetscCall(ISComplement(ilink->is_col, rstart, rend, &ccis));
    1094:     if (jac->offdiag_use_amat) {
    1095:       PetscCall(MatCreateSubMatrix(pc->mat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C));
    1096:     } else {
    1097:       PetscCall(MatCreateSubMatrix(pc->pmat, ilink->is, ccis, MAT_INITIAL_MATRIX, &jac->C));
    1098:     }
    1099:     PetscCall(ISDestroy(&ccis));
    1100:     /* Create work vectors for GKB algorithm */
    1101:     PetscCall(VecDuplicate(ilink->x, &jac->v));
    1102:     PetscCall(VecDuplicate(ilink->x, &jac->d));
    1103:     PetscCall(VecDuplicate(ilink->x, &jac->w1));
    1104:     PetscCall(MatGolubKahanComputeExplicitOperator(jac->mat[0], jac->B, jac->C, &jac->H, jac->gkbnu));
    1105:     PetscCall(PetscCalloc1(jac->gkbdelay, &jac->vecz));

    1107:     ilink = jac->head;
    1108:     PetscCall(KSPSetOperators(ilink->ksp, jac->H, jac->H));
    1109:     if (!jac->suboptionsset) PetscCall(KSPSetFromOptions(ilink->ksp));
    1110:     /* Create gkb_monitor context */
    1111:     if (jac->gkbmonitor) {
    1112:       PetscInt tablevel;
    1113:       PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &jac->gkbviewer));
    1114:       PetscCall(PetscViewerSetType(jac->gkbviewer, PETSCVIEWERASCII));
    1115:       PetscCall(PetscObjectGetTabLevel((PetscObject)ilink->ksp, &tablevel));
    1116:       PetscCall(PetscViewerASCIISetTab(jac->gkbviewer, tablevel));
    1117:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)ilink->ksp, 1));
    1118:     }
    1119:   } else {
    1120:     /* set up the individual splits' PCs */
    1121:     i     = 0;
    1122:     ilink = jac->head;
    1123:     while (ilink) {
    1124:       PetscCall(KSPSetOperators(ilink->ksp, jac->mat[i], jac->pmat[i]));
    1125:       /* really want setfromoptions called in PCSetFromOptions_FieldSplit(), but it is not ready yet */
    1126:       if (!jac->suboptionsset) PetscCall(KSPSetFromOptions(ilink->ksp));
    1127:       i++;
    1128:       ilink = ilink->next;
    1129:     }
    1130:   }

    1132:   /* Set coordinates to the sub PC objects whenever these are set */
    1133:   if (jac->coordinates_set) {
    1134:     PC pc_coords;
    1135:     if (jac->type == PC_COMPOSITE_SCHUR) {
    1136:       // Head is first block.
    1137:       PetscCall(KSPGetPC(jac->head->ksp, &pc_coords));
    1138:       PetscCall(PCSetCoordinates(pc_coords, jac->head->dim, jac->head->ndofs, jac->head->coords));
    1139:       // Second one is Schur block, but its KSP object is in kspschur.
    1140:       PetscCall(KSPGetPC(jac->kspschur, &pc_coords));
    1141:       PetscCall(PCSetCoordinates(pc_coords, jac->head->next->dim, jac->head->next->ndofs, jac->head->next->coords));
    1142:     } else if (jac->type == PC_COMPOSITE_GKB) {
    1143:       PetscCall(PetscInfo(pc, "Warning: Setting coordinates does nothing for the GKB Fieldpslit preconditioner\n"));
    1144:     } else {
    1145:       ilink = jac->head;
    1146:       while (ilink) {
    1147:         PetscCall(KSPGetPC(ilink->ksp, &pc_coords));
    1148:         PetscCall(PCSetCoordinates(pc_coords, ilink->dim, ilink->ndofs, ilink->coords));
    1149:         ilink = ilink->next;
    1150:       }
    1151:     }
    1152:   }

    1154:   jac->suboptionsset = PETSC_TRUE;
    1155:   PetscFunctionReturn(PETSC_SUCCESS);
    1156: }

    1158: #define FieldSplitSplitSolveAdd(ilink, xx, yy) \
    1159:   ((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) || \
    1160:                     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) || \
    1161:                     VecScatterEnd(ilink->sctx, ilink->y, yy, ADD_VALUES, SCATTER_REVERSE)))

    1163: static PetscErrorCode PCSetUpOnBlocks_FieldSplit_Schur(PC pc)
    1164: {
    1165:   PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
    1166:   PC_FieldSplitLink ilinkA = jac->head;
    1167:   KSP               kspA = ilinkA->ksp, kspUpper = jac->kspupper;

    1169:   PetscFunctionBegin;
    1170:   if (jac->schurfactorization == PC_FIELDSPLIT_SCHUR_FACT_FULL && kspUpper != kspA) {
    1171:     PetscCall(KSPSetUp(kspUpper));
    1172:     PetscCall(KSPSetUpOnBlocks(kspUpper));
    1173:   }
    1174:   PetscCall(KSPSetUp(kspA));
    1175:   PetscCall(KSPSetUpOnBlocks(kspA));
    1176:   if (jac->schurpre != PC_FIELDSPLIT_SCHUR_PRE_FULL) {
    1177:     PetscCall(KSPSetUp(jac->kspschur));
    1178:     PetscCall(KSPSetUpOnBlocks(jac->kspschur));
    1179:   } else if (kspUpper == kspA) {
    1180:     Mat      AinvB, A;
    1181:     PetscInt m, M, N;

    1183:     PetscCall(PetscObjectQuery((PetscObject)jac->schur, "AinvB", (PetscObject *)&AinvB));
    1184:     if (AinvB) {
    1185:       PetscCall(MatGetSize(AinvB, NULL, &N));
    1186:       if (N == -1) { // first time PCSetUpOnBlocks_FieldSplit_Schur() is called
    1187:         VecType      vtype;
    1188:         PetscMemType mtype;
    1189:         PetscScalar *array;

    1191:         PetscCall(MatGetSize(jac->B, &M, &N));
    1192:         PetscCall(MatGetLocalSize(jac->B, &m, NULL));
    1193:         PetscCall(MatGetVecType(jac->B, &vtype));
    1194:         PetscCall(VecGetArrayAndMemType(ilinkA->x, &array, &mtype));
    1195:         PetscCall(VecRestoreArrayAndMemType(ilinkA->x, &array));
    1196:         if (PetscMemTypeHost(mtype) || (!PetscDefined(HAVE_CUDA) && !PetscDefined(HAVE_HIP))) PetscCall(PetscMalloc1(m * (N + 1), &array));
    1197: #if PetscDefined(HAVE_CUDA)
    1198:         else if (PetscMemTypeCUDA(mtype)) PetscCallCUDA(cudaMalloc((void **)&array, sizeof(PetscScalar) * m * (N + 1)));
    1199: #endif
    1200: #if PetscDefined(HAVE_HIP)
    1201:         else if (PetscMemTypeHIP(mtype)) PetscCallHIP(hipMalloc((void **)&array, sizeof(PetscScalar) * m * (N + 1)));
    1202: #endif
    1203:         PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)jac->schur), vtype, m, PETSC_DECIDE, M, N + 1, -1, array, &A)); // number of columns of the Schur complement plus one
    1204:         PetscCall(MatHeaderReplace(AinvB, &A));
    1205:       }
    1206:     }
    1207:   }
    1208:   PetscFunctionReturn(PETSC_SUCCESS);
    1209: }

    1211: static PetscErrorCode PCSetUpOnBlocks_FieldSplit(PC pc)
    1212: {
    1213:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    1214:   PC_FieldSplitLink ilink = jac->head;

    1216:   PetscFunctionBegin;
    1217:   while (ilink) {
    1218:     PetscCall(KSPSetUp(ilink->ksp));
    1219:     PetscCall(KSPSetUpOnBlocks(ilink->ksp));
    1220:     ilink = ilink->next;
    1221:   }
    1222:   PetscFunctionReturn(PETSC_SUCCESS);
    1223: }

    1225: static PetscErrorCode PCSetUpOnBlocks_FieldSplit_GKB(PC pc)
    1226: {
    1227:   PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
    1228:   PC_FieldSplitLink ilinkA = jac->head;
    1229:   KSP               ksp    = ilinkA->ksp;

    1231:   PetscFunctionBegin;
    1232:   PetscCall(KSPSetUp(ksp));
    1233:   PetscCall(KSPSetUpOnBlocks(ksp));
    1234:   PetscFunctionReturn(PETSC_SUCCESS);
    1235: }

    1237: static PetscErrorCode PCApply_FieldSplit_Schur(PC pc, Vec x, Vec y)
    1238: {
    1239:   PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
    1240:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
    1241:   KSP               kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;
    1242:   Mat               AinvB = NULL;
    1243:   PetscInt          N, P;

    1245:   PetscFunctionBegin;
    1246:   switch (jac->schurfactorization) {
    1247:   case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
    1248:     /* [A00 0; 0 -S], positive definite, suitable for MINRES */
    1249:     PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1250:     PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1251:     PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1252:     PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1253:     PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
    1254:     PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1255:     PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1256:     PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1257:     PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1258:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1259:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1260:     PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
    1261:     PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    1262:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1263:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1264:     PetscCall(VecScale(ilinkD->y, jac->schurscale));
    1265:     PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1266:     PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1267:     PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1268:     break;
    1269:   case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
    1270:     /* [A00 0; A10 S], suitable for left preconditioning */
    1271:     PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1272:     PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1273:     PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1274:     PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
    1275:     PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1276:     PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1277:     PetscCall(MatMult(jac->C, ilinkA->y, ilinkD->x));
    1278:     PetscCall(VecScale(ilinkD->x, -1.));
    1279:     PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    1280:     PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1281:     PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    1282:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1283:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1284:     PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
    1285:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1286:     PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    1287:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1288:     PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1289:     PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1290:     PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1291:     break;
    1292:   case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
    1293:     /* [A00 A01; 0 S], suitable for right preconditioning */
    1294:     PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1295:     PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1296:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1297:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1298:     PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
    1299:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1300:     PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    1301:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1302:     PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->x));
    1303:     PetscCall(VecScale(ilinkA->x, -1.));
    1304:     PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
    1305:     PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1306:     PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
    1307:     PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1308:     PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
    1309:     PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1310:     PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1311:     PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1312:     PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1313:     PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1314:     break;
    1315:   case PC_FIELDSPLIT_SCHUR_FACT_FULL:
    1316:     /* [1 0; A10 A00^{-1} 1] [A00 0; 0 S] [1 A00^{-1}A01; 0 1] */
    1317:     PetscCall(MatGetSize(jac->B, NULL, &P));
    1318:     N = P;
    1319:     PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1320:     PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1321:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL));
    1322:     if (kspUpper == kspA) {
    1323:       PetscCall(PetscObjectQuery((PetscObject)jac->schur, "AinvB", (PetscObject *)&AinvB));
    1324:       if (AinvB) {
    1325:         PetscCall(MatGetSize(AinvB, NULL, &N));
    1326:         if (N > P) { // first time PCApply_FieldSplit_Schur() is called
    1327:           PetscMemType mtype;
    1328:           Vec          c = NULL;
    1329:           PetscScalar *array;
    1330:           PetscInt     m, M;

    1332:           PetscCall(MatGetSize(jac->B, &M, NULL));
    1333:           PetscCall(MatGetLocalSize(jac->B, &m, NULL));
    1334:           PetscCall(MatDenseGetArrayAndMemType(AinvB, &array, &mtype));
    1335:           if (PetscMemTypeHost(mtype) || (!PetscDefined(HAVE_CUDA) && !PetscDefined(HAVE_HIP))) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)jac->schur), 1, m, M, array + m * P, &c));
    1336: #if PetscDefined(HAVE_CUDA)
    1337:           else if (PetscMemTypeCUDA(mtype)) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)jac->schur), 1, m, M, array + m * P, &c));
    1338: #endif
    1339: #if PetscDefined(HAVE_HIP)
    1340:           else if (PetscMemTypeHIP(mtype)) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)jac->schur), 1, m, M, array + m * P, &c));
    1341: #endif
    1342:           PetscCall(MatDenseRestoreArrayAndMemType(AinvB, &array));
    1343:           PetscCall(VecCopy(ilinkA->x, c));
    1344:           PetscCall(MatSchurComplementComputeExplicitOperator(jac->schur, &jac->schur_user));
    1345:           PetscCall(KSPSetOperators(jac->kspschur, jac->schur, jac->schur_user));
    1346:           PetscCall(VecCopy(c, ilinkA->y)); // retrieve the solution as the last column of the composed Mat
    1347:           PetscCall(VecDestroy(&c));
    1348:         }
    1349:       }
    1350:     }
    1351:     if (N == P) PetscCall(KSPSolve(kspLower, ilinkA->x, ilinkA->y));
    1352:     PetscCall(KSPCheckSolve(kspLower, pc, ilinkA->y));
    1353:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->y, NULL));
    1354:     PetscCall(MatMult(jac->C, ilinkA->y, ilinkD->x));
    1355:     PetscCall(VecScale(ilinkD->x, -1.0));
    1356:     PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    1357:     PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));

    1359:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1360:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1361:     PetscCall(KSPSolve(jac->kspschur, ilinkD->x, ilinkD->y));
    1362:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1363:     PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    1364:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1365:     PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));

    1367:     if (kspUpper == kspA) {
    1368:       if (!AinvB) {
    1369:         PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->y));
    1370:         PetscCall(VecAXPY(ilinkA->x, -1.0, ilinkA->y));
    1371:         PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1372:         PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
    1373:         PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1374:         PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1375:       } else PetscCall(MatMultAdd(AinvB, ilinkD->y, ilinkA->y, ilinkA->y));
    1376:     } else {
    1377:       PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1378:       PetscCall(KSPSolve(kspA, ilinkA->x, ilinkA->y));
    1379:       PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1380:       PetscCall(MatMult(jac->B, ilinkD->y, ilinkA->x));
    1381:       PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL));
    1382:       PetscCall(KSPSolve(kspUpper, ilinkA->x, ilinkA->z));
    1383:       PetscCall(KSPCheckSolve(kspUpper, pc, ilinkA->z));
    1384:       PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->z, NULL));
    1385:       PetscCall(VecAXPY(ilinkA->y, -1.0, ilinkA->z));
    1386:     }
    1387:     PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1388:     PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1389:     PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1390:   }
    1391:   PetscFunctionReturn(PETSC_SUCCESS);
    1392: }

    1394: static PetscErrorCode PCApplyTranspose_FieldSplit_Schur(PC pc, Vec x, Vec y)
    1395: {
    1396:   PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
    1397:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
    1398:   KSP               kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;

    1400:   PetscFunctionBegin;
    1401:   switch (jac->schurfactorization) {
    1402:   case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
    1403:     /* [A00 0; 0 -S], positive definite, suitable for MINRES */
    1404:     PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1405:     PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1406:     PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1407:     PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1408:     PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
    1409:     PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1410:     PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1411:     PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1412:     PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1413:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1414:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1415:     PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
    1416:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1417:     PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    1418:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1419:     PetscCall(VecScale(ilinkD->y, jac->schurscale));
    1420:     PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1421:     PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1422:     PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1423:     break;
    1424:   case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
    1425:     PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1426:     PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1427:     PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1428:     PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
    1429:     PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1430:     PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1431:     PetscCall(MatMultTranspose(jac->B, ilinkA->y, ilinkD->x));
    1432:     PetscCall(VecScale(ilinkD->x, -1.));
    1433:     PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    1434:     PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1435:     PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_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(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
    1439:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1440:     PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    1441:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1442:     PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1443:     PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1444:     PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1445:     break;
    1446:   case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
    1447:     PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1448:     PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1449:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1450:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1451:     PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
    1452:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1453:     PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    1454:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1455:     PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->x));
    1456:     PetscCall(VecScale(ilinkA->x, -1.));
    1457:     PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
    1458:     PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1459:     PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
    1460:     PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1461:     PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
    1462:     PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1463:     PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1464:     PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1465:     PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1466:     PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1467:     break;
    1468:   case PC_FIELDSPLIT_SCHUR_FACT_FULL:
    1469:     PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1470:     PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1471:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->y, NULL));
    1472:     PetscCall(KSPSolveTranspose(kspUpper, ilinkA->x, ilinkA->y));
    1473:     PetscCall(KSPCheckSolve(kspUpper, pc, ilinkA->y));
    1474:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->y, NULL));
    1475:     PetscCall(MatMultTranspose(jac->B, ilinkA->y, ilinkD->x));
    1476:     PetscCall(VecScale(ilinkD->x, -1.0));
    1477:     PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
    1478:     PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));

    1480:     PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1481:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, 1));
    1482:     PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
    1483:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)kspA, (PetscObject)kspA, -1));
    1484:     PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
    1485:     PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
    1486:     PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));

    1488:     if (kspLower == kspA) {
    1489:       PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->y));
    1490:       PetscCall(VecAXPY(ilinkA->x, -1.0, ilinkA->y));
    1491:       PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1492:       PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
    1493:       PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1494:       PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1495:     } else {
    1496:       PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
    1497:       PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
    1498:       PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
    1499:       PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->x));
    1500:       PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->z, NULL));
    1501:       PetscCall(KSPSolveTranspose(kspLower, ilinkA->x, ilinkA->z));
    1502:       PetscCall(KSPCheckSolve(kspLower, pc, ilinkA->z));
    1503:       PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->z, NULL));
    1504:       PetscCall(VecAXPY(ilinkA->y, -1.0, ilinkA->z));
    1505:     }
    1506:     PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1507:     PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1508:     PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1509:   }
    1510:   PetscFunctionReturn(PETSC_SUCCESS);
    1511: }

    1513: static PetscErrorCode PCApply_FieldSplit(PC pc, Vec x, Vec y)
    1514: {
    1515:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    1516:   PC_FieldSplitLink ilink = jac->head;
    1517:   PetscInt          cnt, bs;

    1519:   PetscFunctionBegin;
    1520:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
    1521:     PetscBool matnest;

    1523:     PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &matnest));
    1524:     if (jac->defaultsplit && !matnest) {
    1525:       PetscCall(VecGetBlockSize(x, &bs));
    1526:       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);
    1527:       PetscCall(VecGetBlockSize(y, &bs));
    1528:       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);
    1529:       PetscCall(VecStrideGatherAll(x, jac->x, INSERT_VALUES));
    1530:       while (ilink) {
    1531:         PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1532:         PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
    1533:         PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
    1534:         PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1535:         ilink = ilink->next;
    1536:       }
    1537:       PetscCall(VecStrideScatterAll(jac->y, y, INSERT_VALUES));
    1538:     } else {
    1539:       PetscCall(VecSet(y, 0.0));
    1540:       while (ilink) {
    1541:         PetscCall(FieldSplitSplitSolveAdd(ilink, x, y));
    1542:         ilink = ilink->next;
    1543:       }
    1544:     }
    1545:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE && jac->nsplits == 2) {
    1546:     PetscCall(VecSet(y, 0.0));
    1547:     /* solve on first block for first block variables */
    1548:     PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD));
    1549:     PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, INSERT_VALUES, SCATTER_FORWARD));
    1550:     PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1551:     PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
    1552:     PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
    1553:     PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1554:     PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    1555:     PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));

    1557:     /* compute the residual only onto second block variables using first block variables */
    1558:     PetscCall(MatMult(jac->Afield[1], ilink->y, ilink->next->x));
    1559:     ilink = ilink->next;
    1560:     PetscCall(VecScale(ilink->x, -1.0));
    1561:     PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
    1562:     PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));

    1564:     /* solve on second block variables */
    1565:     PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1566:     PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
    1567:     PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
    1568:     PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1569:     PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    1570:     PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    1571:   } else if (jac->type == PC_COMPOSITE_MULTIPLICATIVE || jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
    1572:     if (!jac->w1) {
    1573:       PetscCall(VecDuplicate(x, &jac->w1));
    1574:       PetscCall(VecDuplicate(x, &jac->w2));
    1575:     }
    1576:     PetscCall(VecSet(y, 0.0));
    1577:     PetscCall(FieldSplitSplitSolveAdd(ilink, x, y));
    1578:     cnt = 1;
    1579:     while (ilink->next) {
    1580:       ilink = ilink->next;
    1581:       /* compute the residual only over the part of the vector needed */
    1582:       PetscCall(MatMult(jac->Afield[cnt++], y, ilink->x));
    1583:       PetscCall(VecScale(ilink->x, -1.0));
    1584:       PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
    1585:       PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
    1586:       PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1587:       PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
    1588:       PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
    1589:       PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1590:       PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    1591:       PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    1592:     }
    1593:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
    1594:       cnt -= 2;
    1595:       while (ilink->previous) {
    1596:         ilink = ilink->previous;
    1597:         /* compute the residual only over the part of the vector needed */
    1598:         PetscCall(MatMult(jac->Afield[cnt--], y, ilink->x));
    1599:         PetscCall(VecScale(ilink->x, -1.0));
    1600:         PetscCall(VecScatterBegin(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
    1601:         PetscCall(VecScatterEnd(ilink->sctx, x, ilink->x, ADD_VALUES, SCATTER_FORWARD));
    1602:         PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1603:         PetscCall(KSPSolve(ilink->ksp, ilink->x, ilink->y));
    1604:         PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
    1605:         PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1606:         PetscCall(VecScatterBegin(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    1607:         PetscCall(VecScatterEnd(ilink->sctx, ilink->y, y, ADD_VALUES, SCATTER_REVERSE));
    1608:       }
    1609:     }
    1610:   } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Unsupported or unknown composition %d", (int)jac->type);
    1611:   PetscFunctionReturn(PETSC_SUCCESS);
    1612: }

    1614: static PetscErrorCode PCApply_FieldSplit_GKB(PC pc, Vec x, Vec y)
    1615: {
    1616:   PC_FieldSplit    *jac    = (PC_FieldSplit *)pc->data;
    1617:   PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
    1618:   KSP               ksp = ilinkA->ksp;
    1619:   Vec               u, v, Hu, d, work1, work2;
    1620:   PetscScalar       alpha, z, nrmz2, *vecz;
    1621:   PetscReal         lowbnd, nu, beta;
    1622:   PetscInt          j, iterGKB;

    1624:   PetscFunctionBegin;
    1625:   PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1626:   PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
    1627:   PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
    1628:   PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));

    1630:   u     = jac->u;
    1631:   v     = jac->v;
    1632:   Hu    = jac->Hu;
    1633:   d     = jac->d;
    1634:   work1 = jac->w1;
    1635:   work2 = jac->w2;
    1636:   vecz  = jac->vecz;

    1638:   /* Change RHS to comply with matrix regularization H = A + nu*B*B' */
    1639:   /* Add q = q + nu*B*b */
    1640:   if (jac->gkbnu) {
    1641:     nu = jac->gkbnu;
    1642:     PetscCall(VecScale(ilinkD->x, jac->gkbnu));
    1643:     PetscCall(MatMultAdd(jac->B, ilinkD->x, ilinkA->x, ilinkA->x)); /* q = q + nu*B*b */
    1644:   } else {
    1645:     /* Situation when no augmented Lagrangian is used. Then we set inner  */
    1646:     /* matrix N = I in [Ar13], and thus nu = 1.                           */
    1647:     nu = 1;
    1648:   }

    1650:   /* Transform rhs from [q,tilde{b}] to [0,b] */
    1651:   PetscCall(PetscLogEventBegin(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL));
    1652:   PetscCall(KSPSolve(ksp, ilinkA->x, ilinkA->y));
    1653:   PetscCall(KSPCheckSolve(ksp, pc, ilinkA->y));
    1654:   PetscCall(PetscLogEventEnd(ilinkA->event, ksp, ilinkA->x, ilinkA->y, NULL));
    1655:   PetscCall(MatMultHermitianTranspose(jac->B, ilinkA->y, work1));
    1656:   PetscCall(VecAXPBY(work1, 1.0 / nu, -1.0, ilinkD->x)); /* c = b - B'*x        */

    1658:   /* First step of algorithm */
    1659:   PetscCall(VecNorm(work1, NORM_2, &beta)); /* beta = sqrt(nu*c'*c)*/
    1660:   KSPCheckDot(ksp, beta);
    1661:   beta = PetscSqrtReal(nu) * beta;
    1662:   PetscCall(VecAXPBY(v, nu / beta, 0.0, work1)); /* v = nu/beta *c      */
    1663:   PetscCall(MatMult(jac->B, v, work2));          /* u = H^{-1}*B*v      */
    1664:   PetscCall(PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL));
    1665:   PetscCall(KSPSolve(ksp, work2, u));
    1666:   PetscCall(KSPCheckSolve(ksp, pc, u));
    1667:   PetscCall(PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL));
    1668:   PetscCall(MatMult(jac->H, u, Hu)); /* alpha = u'*H*u      */
    1669:   PetscCall(VecDot(Hu, u, &alpha));
    1670:   KSPCheckDot(ksp, alpha);
    1671:   PetscCheck(PetscRealPart(alpha) > 0.0, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "GKB preconditioner diverged, H is not positive definite");
    1672:   alpha = PetscSqrtReal(PetscAbsScalar(alpha));
    1673:   PetscCall(VecScale(u, 1.0 / alpha));
    1674:   PetscCall(VecAXPBY(d, 1.0 / alpha, 0.0, v)); /* v = nu/beta *c      */

    1676:   z       = beta / alpha;
    1677:   vecz[1] = z;

    1679:   /* Computation of first iterate x(1) and p(1) */
    1680:   PetscCall(VecAXPY(ilinkA->y, z, u));
    1681:   PetscCall(VecCopy(d, ilinkD->y));
    1682:   PetscCall(VecScale(ilinkD->y, -z));

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

    1688:   while (iterGKB < jac->gkbmaxit && lowbnd > jac->gkbtol) {
    1689:     iterGKB += 1;
    1690:     PetscCall(MatMultHermitianTranspose(jac->B, u, work1)); /* v <- nu*(B'*u-alpha/nu*v) */
    1691:     PetscCall(VecAXPBY(v, nu, -alpha, work1));
    1692:     PetscCall(VecNorm(v, NORM_2, &beta)); /* beta = sqrt(nu)*v'*v      */
    1693:     beta = beta / PetscSqrtReal(nu);
    1694:     PetscCall(VecScale(v, 1.0 / beta));
    1695:     PetscCall(MatMult(jac->B, v, work2)); /* u <- H^{-1}*(B*v-beta*H*u) */
    1696:     PetscCall(MatMult(jac->H, u, Hu));
    1697:     PetscCall(VecAXPY(work2, -beta, Hu));
    1698:     PetscCall(PetscLogEventBegin(ilinkA->event, ksp, work2, u, NULL));
    1699:     PetscCall(KSPSolve(ksp, work2, u));
    1700:     PetscCall(KSPCheckSolve(ksp, pc, u));
    1701:     PetscCall(PetscLogEventEnd(ilinkA->event, ksp, work2, u, NULL));
    1702:     PetscCall(MatMult(jac->H, u, Hu)); /* alpha = u'*H*u            */
    1703:     PetscCall(VecDot(Hu, u, &alpha));
    1704:     KSPCheckDot(ksp, alpha);
    1705:     PetscCheck(PetscRealPart(alpha) > 0.0, PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "GKB preconditioner diverged, H is not positive definite");
    1706:     alpha = PetscSqrtReal(PetscAbsScalar(alpha));
    1707:     PetscCall(VecScale(u, 1.0 / alpha));

    1709:     z       = -beta / alpha * z; /* z <- beta/alpha*z     */
    1710:     vecz[0] = z;

    1712:     /* Computation of new iterate x(i+1) and p(i+1) */
    1713:     PetscCall(VecAXPBY(d, 1.0 / alpha, -beta / alpha, v)); /* d = (v-beta*d)/alpha */
    1714:     PetscCall(VecAXPY(ilinkA->y, z, u));                   /* r = r + z*u          */
    1715:     PetscCall(VecAXPY(ilinkD->y, -z, d));                  /* p = p - z*d          */
    1716:     PetscCall(MatMult(jac->H, ilinkA->y, Hu));             /* ||u||_H = u'*H*u     */
    1717:     PetscCall(VecDot(Hu, ilinkA->y, &nrmz2));

    1719:     /* Compute Lower Bound estimate */
    1720:     if (iterGKB > jac->gkbdelay) {
    1721:       lowbnd = 0.0;
    1722:       for (j = 0; j < jac->gkbdelay; j++) lowbnd += PetscAbsScalar(vecz[j] * vecz[j]);
    1723:       lowbnd = PetscSqrtReal(lowbnd / PetscAbsScalar(nrmz2));
    1724:     }

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

    1730:   PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1731:   PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1732:   PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1733:   PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
    1734:   PetscFunctionReturn(PETSC_SUCCESS);
    1735: }

    1737: #define FieldSplitSplitSolveAddTranspose(ilink, xx, yy) \
    1738:   ((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) || \
    1739:                     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) || \
    1740:                     VecScatterEnd(ilink->sctx, ilink->x, yy, ADD_VALUES, SCATTER_REVERSE)))

    1742: static PetscErrorCode PCApplyTranspose_FieldSplit(PC pc, Vec x, Vec y)
    1743: {
    1744:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    1745:   PC_FieldSplitLink ilink = jac->head;
    1746:   PetscInt          bs;

    1748:   PetscFunctionBegin;
    1749:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
    1750:     PetscBool matnest;

    1752:     PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATNEST, &matnest));
    1753:     if (jac->defaultsplit && !matnest) {
    1754:       PetscCall(VecGetBlockSize(x, &bs));
    1755:       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);
    1756:       PetscCall(VecGetBlockSize(y, &bs));
    1757:       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);
    1758:       PetscCall(VecStrideGatherAll(x, jac->x, INSERT_VALUES));
    1759:       while (ilink) {
    1760:         PetscCall(PetscLogEventBegin(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1761:         PetscCall(KSPSolveTranspose(ilink->ksp, ilink->x, ilink->y));
    1762:         PetscCall(KSPCheckSolve(ilink->ksp, pc, ilink->y));
    1763:         PetscCall(PetscLogEventEnd(ilink->event, ilink->ksp, ilink->x, ilink->y, NULL));
    1764:         ilink = ilink->next;
    1765:       }
    1766:       PetscCall(VecStrideScatterAll(jac->y, y, INSERT_VALUES));
    1767:     } else {
    1768:       PetscCall(VecSet(y, 0.0));
    1769:       while (ilink) {
    1770:         PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y));
    1771:         ilink = ilink->next;
    1772:       }
    1773:     }
    1774:   } else {
    1775:     if (!jac->w1) {
    1776:       PetscCall(VecDuplicate(x, &jac->w1));
    1777:       PetscCall(VecDuplicate(x, &jac->w2));
    1778:     }
    1779:     PetscCall(VecSet(y, 0.0));
    1780:     if (jac->type == PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE) {
    1781:       PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y));
    1782:       while (ilink->next) {
    1783:         ilink = ilink->next;
    1784:         PetscCall(MatMultTranspose(pc->mat, y, jac->w1));
    1785:         PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x));
    1786:         PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y));
    1787:       }
    1788:       while (ilink->previous) {
    1789:         ilink = ilink->previous;
    1790:         PetscCall(MatMultTranspose(pc->mat, y, jac->w1));
    1791:         PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x));
    1792:         PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y));
    1793:       }
    1794:     } else {
    1795:       while (ilink->next) { /* get to last entry in linked list */
    1796:         ilink = ilink->next;
    1797:       }
    1798:       PetscCall(FieldSplitSplitSolveAddTranspose(ilink, x, y));
    1799:       while (ilink->previous) {
    1800:         ilink = ilink->previous;
    1801:         PetscCall(MatMultTranspose(pc->mat, y, jac->w1));
    1802:         PetscCall(VecWAXPY(jac->w2, -1.0, jac->w1, x));
    1803:         PetscCall(FieldSplitSplitSolveAddTranspose(ilink, jac->w2, y));
    1804:       }
    1805:     }
    1806:   }
    1807:   PetscFunctionReturn(PETSC_SUCCESS);
    1808: }

    1810: static PetscErrorCode PCReset_FieldSplit(PC pc)
    1811: {
    1812:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    1813:   PC_FieldSplitLink ilink = jac->head, next;

    1815:   PetscFunctionBegin;
    1816:   while (ilink) {
    1817:     PetscCall(KSPDestroy(&ilink->ksp));
    1818:     PetscCall(VecDestroy(&ilink->x));
    1819:     PetscCall(VecDestroy(&ilink->y));
    1820:     PetscCall(VecDestroy(&ilink->z));
    1821:     PetscCall(VecScatterDestroy(&ilink->sctx));
    1822:     PetscCall(ISDestroy(&ilink->is));
    1823:     PetscCall(ISDestroy(&ilink->is_col));
    1824:     PetscCall(PetscFree(ilink->splitname));
    1825:     PetscCall(PetscFree(ilink->fields));
    1826:     PetscCall(PetscFree(ilink->fields_col));
    1827:     next = ilink->next;
    1828:     PetscCall(PetscFree(ilink));
    1829:     ilink = next;
    1830:   }
    1831:   jac->head = NULL;
    1832:   PetscCall(PetscFree2(jac->x, jac->y));
    1833:   if (jac->mat && jac->mat != jac->pmat) {
    1834:     PetscCall(MatDestroyMatrices(jac->nsplits, &jac->mat));
    1835:   } else if (jac->mat) {
    1836:     jac->mat = NULL;
    1837:   }
    1838:   if (jac->pmat) PetscCall(MatDestroyMatrices(jac->nsplits, &jac->pmat));
    1839:   if (jac->Afield) PetscCall(MatDestroyMatrices(jac->nsplits, &jac->Afield));
    1840:   jac->nsplits = 0;
    1841:   PetscCall(VecDestroy(&jac->w1));
    1842:   PetscCall(VecDestroy(&jac->w2));
    1843:   if (jac->schur) PetscCall(PetscObjectCompose((PetscObject)jac->schur, "AinvB", NULL));
    1844:   PetscCall(MatDestroy(&jac->schur));
    1845:   PetscCall(MatDestroy(&jac->schurp));
    1846:   PetscCall(MatDestroy(&jac->schur_user));
    1847:   PetscCall(KSPDestroy(&jac->kspschur));
    1848:   PetscCall(KSPDestroy(&jac->kspupper));
    1849:   PetscCall(MatDestroy(&jac->B));
    1850:   PetscCall(MatDestroy(&jac->C));
    1851:   PetscCall(MatDestroy(&jac->H));
    1852:   PetscCall(VecDestroy(&jac->u));
    1853:   PetscCall(VecDestroy(&jac->v));
    1854:   PetscCall(VecDestroy(&jac->Hu));
    1855:   PetscCall(VecDestroy(&jac->d));
    1856:   PetscCall(PetscFree(jac->vecz));
    1857:   PetscCall(PetscViewerDestroy(&jac->gkbviewer));
    1858:   jac->isrestrict = PETSC_FALSE;
    1859:   PetscFunctionReturn(PETSC_SUCCESS);
    1860: }

    1862: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
    1863: {
    1864:   PetscFunctionBegin;
    1865:   PetscCall(PCReset_FieldSplit(pc));
    1866:   PetscCall(PetscFree(pc->data));
    1867:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
    1868:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", NULL));
    1869:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", NULL));
    1870:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", NULL));
    1871:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", NULL));
    1872:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", NULL));
    1873:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", NULL));
    1874:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL));
    1875:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL));
    1876:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL));
    1877:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL));
    1878:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL));
    1879:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL));
    1880:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL));
    1881:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL));
    1882:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL));
    1883:   PetscFunctionReturn(PETSC_SUCCESS);
    1884: }

    1886: static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc, PetscOptionItems *PetscOptionsObject)
    1887: {
    1888:   PetscInt        bs;
    1889:   PetscBool       flg;
    1890:   PC_FieldSplit  *jac = (PC_FieldSplit *)pc->data;
    1891:   PCCompositeType ctype;

    1893:   PetscFunctionBegin;
    1894:   PetscOptionsHeadBegin(PetscOptionsObject, "FieldSplit options");
    1895:   PetscCall(PetscOptionsBool("-pc_fieldsplit_dm_splits", "Whether to use DMCreateFieldDecomposition() for splits", "PCFieldSplitSetDMSplits", jac->dm_splits, &jac->dm_splits, NULL));
    1896:   PetscCall(PetscOptionsInt("-pc_fieldsplit_block_size", "Blocksize that defines number of fields", "PCFieldSplitSetBlockSize", jac->bs, &bs, &flg));
    1897:   if (flg) PetscCall(PCFieldSplitSetBlockSize(pc, bs));
    1898:   jac->diag_use_amat = pc->useAmat;
    1899:   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));
    1900:   jac->offdiag_use_amat = pc->useAmat;
    1901:   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));
    1902:   PetscCall(PetscOptionsBool("-pc_fieldsplit_detect_saddle_point", "Form 2-way split by detecting zero diagonal entries", "PCFieldSplitSetDetectSaddlePoint", jac->detect, &jac->detect, NULL));
    1903:   PetscCall(PCFieldSplitSetDetectSaddlePoint(pc, jac->detect)); /* Sets split type and Schur PC type */
    1904:   PetscCall(PetscOptionsEnum("-pc_fieldsplit_type", "Type of composition", "PCFieldSplitSetType", PCCompositeTypes, (PetscEnum)jac->type, (PetscEnum *)&ctype, &flg));
    1905:   if (flg) PetscCall(PCFieldSplitSetType(pc, ctype));
    1906:   /* Only setup fields once */
    1907:   if ((jac->bs > 0) && (jac->nsplits == 0)) {
    1908:     /* only allow user to set fields from command line.
    1909:        otherwise user can set them in PCFieldSplitSetDefaults() */
    1910:     PetscCall(PCFieldSplitSetRuntimeSplits_Private(pc));
    1911:     if (jac->splitdefined) PetscCall(PetscInfo(pc, "Splits defined using the options database\n"));
    1912:   }
    1913:   if (jac->type == PC_COMPOSITE_SCHUR) {
    1914:     PetscCall(PetscOptionsGetEnum(((PetscObject)pc)->options, ((PetscObject)pc)->prefix, "-pc_fieldsplit_schur_factorization_type", PCFieldSplitSchurFactTypes, (PetscEnum *)&jac->schurfactorization, &flg));
    1915:     if (flg) PetscCall(PetscInfo(pc, "Deprecated use of -pc_fieldsplit_schur_factorization_type\n"));
    1916:     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));
    1917:     PetscCall(PetscOptionsEnum("-pc_fieldsplit_schur_precondition", "How to build preconditioner for Schur complement", "PCFieldSplitSetSchurPre", PCFieldSplitSchurPreTypes, (PetscEnum)jac->schurpre, (PetscEnum *)&jac->schurpre, NULL));
    1918:     PetscCall(PetscOptionsScalar("-pc_fieldsplit_schur_scale", "Scale Schur complement", "PCFieldSplitSetSchurScale", jac->schurscale, &jac->schurscale, NULL));
    1919:   } else if (jac->type == PC_COMPOSITE_GKB) {
    1920:     PetscCall(PetscOptionsReal("-pc_fieldsplit_gkb_tol", "The tolerance for the lower bound stopping criterion", "PCFieldSplitSetGKBTol", jac->gkbtol, &jac->gkbtol, NULL));
    1921:     PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_delay", "The delay value for lower bound criterion", "PCFieldSplitSetGKBDelay", jac->gkbdelay, &jac->gkbdelay, NULL));
    1922:     PetscCall(PetscOptionsBoundedReal("-pc_fieldsplit_gkb_nu", "Parameter in augmented Lagrangian approach", "PCFieldSplitSetGKBNu", jac->gkbnu, &jac->gkbnu, NULL, 0.0));
    1923:     PetscCall(PetscOptionsInt("-pc_fieldsplit_gkb_maxit", "Maximum allowed number of iterations", "PCFieldSplitSetGKBMaxit", jac->gkbmaxit, &jac->gkbmaxit, NULL));
    1924:     PetscCall(PetscOptionsBool("-pc_fieldsplit_gkb_monitor", "Prints number of GKB iterations and error", "PCFieldSplitGKB", jac->gkbmonitor, &jac->gkbmonitor, NULL));
    1925:   }
    1926:   /*
    1927:     In the initial call to this routine the sub-solver data structures do not exist so we cannot call KSPSetFromOptions() on them yet.
    1928:     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
    1929:     is called on the outer solver in case changes were made in the options database

    1931:     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()
    1932:     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.
    1933:     Without this extra check test p2p1fetidp_olof_full and others fail with incorrect matrix types.

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

    1937:     If one captured the PetscObjectState of the options database one could skip these calls if the database has not changed from the previous call
    1938:   */
    1939:   if (jac->issetup) {
    1940:     PC_FieldSplitLink ilink = jac->head;
    1941:     if (jac->type == PC_COMPOSITE_SCHUR) {
    1942:       if (jac->kspupper && jac->kspupper->totalits > 0) PetscCall(KSPSetFromOptions(jac->kspupper));
    1943:       if (jac->kspschur && jac->kspschur->totalits > 0) PetscCall(KSPSetFromOptions(jac->kspschur));
    1944:     }
    1945:     while (ilink) {
    1946:       if (ilink->ksp->totalits > 0) PetscCall(KSPSetFromOptions(ilink->ksp));
    1947:       ilink = ilink->next;
    1948:     }
    1949:   }
    1950:   PetscOptionsHeadEnd();
    1951:   PetscFunctionReturn(PETSC_SUCCESS);
    1952: }

    1954: static PetscErrorCode PCFieldSplitSetFields_FieldSplit(PC pc, const char splitname[], PetscInt n, const PetscInt *fields, const PetscInt *fields_col)
    1955: {
    1956:   PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
    1957:   PC_FieldSplitLink ilink, next = jac->head;
    1958:   char              prefix[128];
    1959:   PetscInt          i;

    1961:   PetscFunctionBegin;
    1962:   if (jac->splitdefined) {
    1963:     PetscCall(PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname));
    1964:     PetscFunctionReturn(PETSC_SUCCESS);
    1965:   }
    1966:   for (i = 0; i < n; i++) { PetscCheck(fields[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative field %" PetscInt_FMT " requested", fields[i]); }
    1967:   PetscCall(PetscNew(&ilink));
    1968:   if (splitname) {
    1969:     PetscCall(PetscStrallocpy(splitname, &ilink->splitname));
    1970:   } else {
    1971:     PetscCall(PetscMalloc1(3, &ilink->splitname));
    1972:     PetscCall(PetscSNPrintf(ilink->splitname, 2, "%" PetscInt_FMT, jac->nsplits));
    1973:   }
    1974:   ilink->event = jac->nsplits < 5 ? (PetscLogEvent)(KSP_Solve_FS_0 + jac->nsplits) : (PetscLogEvent)(KSP_Solve_FS_0 + 4); /* Splits greater than 4 logged in 4th split */
    1975:   PetscCall(PetscMalloc1(n, &ilink->fields));
    1976:   PetscCall(PetscArraycpy(ilink->fields, fields, n));
    1977:   PetscCall(PetscMalloc1(n, &ilink->fields_col));
    1978:   PetscCall(PetscArraycpy(ilink->fields_col, fields_col, n));

    1980:   ilink->nfields = n;
    1981:   ilink->next    = NULL;
    1982:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp));
    1983:   PetscCall(KSPSetNestLevel(ilink->ksp, pc->kspnestlevel));
    1984:   PetscCall(KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure));
    1985:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1));
    1986:   PetscCall(KSPSetType(ilink->ksp, KSPPREONLY));

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

    1991:   if (!next) {
    1992:     jac->head       = ilink;
    1993:     ilink->previous = NULL;
    1994:   } else {
    1995:     while (next->next) next = next->next;
    1996:     next->next      = ilink;
    1997:     ilink->previous = next;
    1998:   }
    1999:   jac->nsplits++;
    2000:   PetscFunctionReturn(PETSC_SUCCESS);
    2001: }

    2003: static PetscErrorCode PCFieldSplitSchurGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp)
    2004: {
    2005:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    2007:   PetscFunctionBegin;
    2008:   *subksp = NULL;
    2009:   if (n) *n = 0;
    2010:   if (jac->type == PC_COMPOSITE_SCHUR) {
    2011:     PetscInt nn;

    2013:     PetscCheck(jac->schur, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitSchurGetSubKSP()");
    2014:     PetscCheck(jac->nsplits == 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Unexpected number of splits %" PetscInt_FMT " != 2", jac->nsplits);
    2015:     nn = jac->nsplits + (jac->kspupper != jac->head->ksp ? 1 : 0);
    2016:     PetscCall(PetscMalloc1(nn, subksp));
    2017:     (*subksp)[0] = jac->head->ksp;
    2018:     (*subksp)[1] = jac->kspschur;
    2019:     if (jac->kspupper != jac->head->ksp) (*subksp)[2] = jac->kspupper;
    2020:     if (n) *n = nn;
    2021:   }
    2022:   PetscFunctionReturn(PETSC_SUCCESS);
    2023: }

    2025: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit_Schur(PC pc, PetscInt *n, KSP **subksp)
    2026: {
    2027:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    2029:   PetscFunctionBegin;
    2030:   PetscCheck(jac->schur, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() before calling PCFieldSplitGetSubKSP()");
    2031:   PetscCall(PetscMalloc1(jac->nsplits, subksp));
    2032:   PetscCall(MatSchurComplementGetKSP(jac->schur, *subksp));

    2034:   (*subksp)[1] = jac->kspschur;
    2035:   if (n) *n = jac->nsplits;
    2036:   PetscFunctionReturn(PETSC_SUCCESS);
    2037: }

    2039: static PetscErrorCode PCFieldSplitGetSubKSP_FieldSplit(PC pc, PetscInt *n, KSP **subksp)
    2040: {
    2041:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    2042:   PetscInt          cnt   = 0;
    2043:   PC_FieldSplitLink ilink = jac->head;

    2045:   PetscFunctionBegin;
    2046:   PetscCall(PetscMalloc1(jac->nsplits, subksp));
    2047:   while (ilink) {
    2048:     (*subksp)[cnt++] = ilink->ksp;
    2049:     ilink            = ilink->next;
    2050:   }
    2051:   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);
    2052:   if (n) *n = jac->nsplits;
    2053:   PetscFunctionReturn(PETSC_SUCCESS);
    2054: }

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

    2059:   Input Parameters:
    2060: + pc  - the preconditioner context
    2061: - isy - the index set that defines the indices to which the fieldsplit is to be restricted

    2063:   Level: advanced

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

    2069: .seealso: [](sec_block_matrices), `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
    2070: @*/
    2071: PetscErrorCode PCFieldSplitRestrictIS(PC pc, IS isy)
    2072: {
    2073:   PetscFunctionBegin;
    2076:   PetscTryMethod(pc, "PCFieldSplitRestrictIS_C", (PC, IS), (pc, isy));
    2077:   PetscFunctionReturn(PETSC_SUCCESS);
    2078: }

    2080: static PetscErrorCode PCFieldSplitRestrictIS_FieldSplit(PC pc, IS isy)
    2081: {
    2082:   PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    2083:   PC_FieldSplitLink ilink = jac->head, next;
    2084:   PetscInt          localsize, size, sizez, i;
    2085:   const PetscInt   *ind, *indz;
    2086:   PetscInt         *indc, *indcz;
    2087:   PetscBool         flg;

    2089:   PetscFunctionBegin;
    2090:   PetscCall(ISGetLocalSize(isy, &localsize));
    2091:   PetscCallMPI(MPI_Scan(&localsize, &size, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)isy)));
    2092:   size -= localsize;
    2093:   while (ilink) {
    2094:     IS isrl, isr;
    2095:     PC subpc;
    2096:     PetscCall(ISEmbed(ilink->is, isy, PETSC_TRUE, &isrl));
    2097:     PetscCall(ISGetLocalSize(isrl, &localsize));
    2098:     PetscCall(PetscMalloc1(localsize, &indc));
    2099:     PetscCall(ISGetIndices(isrl, &ind));
    2100:     PetscCall(PetscArraycpy(indc, ind, localsize));
    2101:     PetscCall(ISRestoreIndices(isrl, &ind));
    2102:     PetscCall(ISDestroy(&isrl));
    2103:     for (i = 0; i < localsize; i++) *(indc + i) += size;
    2104:     PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)isy), localsize, indc, PETSC_OWN_POINTER, &isr));
    2105:     PetscCall(PetscObjectReference((PetscObject)isr));
    2106:     PetscCall(ISDestroy(&ilink->is));
    2107:     ilink->is = isr;
    2108:     PetscCall(PetscObjectReference((PetscObject)isr));
    2109:     PetscCall(ISDestroy(&ilink->is_col));
    2110:     ilink->is_col = isr;
    2111:     PetscCall(ISDestroy(&isr));
    2112:     PetscCall(KSPGetPC(ilink->ksp, &subpc));
    2113:     PetscCall(PetscObjectTypeCompare((PetscObject)subpc, PCFIELDSPLIT, &flg));
    2114:     if (flg) {
    2115:       IS       iszl, isz;
    2116:       MPI_Comm comm;
    2117:       PetscCall(ISGetLocalSize(ilink->is, &localsize));
    2118:       comm = PetscObjectComm((PetscObject)ilink->is);
    2119:       PetscCall(ISEmbed(isy, ilink->is, PETSC_TRUE, &iszl));
    2120:       PetscCallMPI(MPI_Scan(&localsize, &sizez, 1, MPIU_INT, MPI_SUM, comm));
    2121:       sizez -= localsize;
    2122:       PetscCall(ISGetLocalSize(iszl, &localsize));
    2123:       PetscCall(PetscMalloc1(localsize, &indcz));
    2124:       PetscCall(ISGetIndices(iszl, &indz));
    2125:       PetscCall(PetscArraycpy(indcz, indz, localsize));
    2126:       PetscCall(ISRestoreIndices(iszl, &indz));
    2127:       PetscCall(ISDestroy(&iszl));
    2128:       for (i = 0; i < localsize; i++) *(indcz + i) += sizez;
    2129:       PetscCall(ISCreateGeneral(comm, localsize, indcz, PETSC_OWN_POINTER, &isz));
    2130:       PetscCall(PCFieldSplitRestrictIS(subpc, isz));
    2131:       PetscCall(ISDestroy(&isz));
    2132:     }
    2133:     next  = ilink->next;
    2134:     ilink = next;
    2135:   }
    2136:   jac->isrestrict = PETSC_TRUE;
    2137:   PetscFunctionReturn(PETSC_SUCCESS);
    2138: }

    2140: static PetscErrorCode PCFieldSplitSetIS_FieldSplit(PC pc, const char splitname[], IS is)
    2141: {
    2142:   PC_FieldSplit    *jac = (PC_FieldSplit *)pc->data;
    2143:   PC_FieldSplitLink ilink, next = jac->head;
    2144:   char              prefix[128];

    2146:   PetscFunctionBegin;
    2147:   if (jac->splitdefined) {
    2148:     PetscCall(PetscInfo(pc, "Ignoring new split \"%s\" because the splits have already been defined\n", splitname));
    2149:     PetscFunctionReturn(PETSC_SUCCESS);
    2150:   }
    2151:   PetscCall(PetscNew(&ilink));
    2152:   if (splitname) {
    2153:     PetscCall(PetscStrallocpy(splitname, &ilink->splitname));
    2154:   } else {
    2155:     PetscCall(PetscMalloc1(8, &ilink->splitname));
    2156:     PetscCall(PetscSNPrintf(ilink->splitname, 7, "%" PetscInt_FMT, jac->nsplits));
    2157:   }
    2158:   ilink->event = jac->nsplits < 5 ? (PetscLogEvent)(KSP_Solve_FS_0 + jac->nsplits) : (PetscLogEvent)(KSP_Solve_FS_0 + 4); /* Splits greater than 4 logged in 4th split */
    2159:   PetscCall(PetscObjectReference((PetscObject)is));
    2160:   PetscCall(ISDestroy(&ilink->is));
    2161:   ilink->is = is;
    2162:   PetscCall(PetscObjectReference((PetscObject)is));
    2163:   PetscCall(ISDestroy(&ilink->is_col));
    2164:   ilink->is_col = is;
    2165:   ilink->next   = NULL;
    2166:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &ilink->ksp));
    2167:   PetscCall(KSPSetNestLevel(ilink->ksp, pc->kspnestlevel));
    2168:   PetscCall(KSPSetErrorIfNotConverged(ilink->ksp, pc->erroriffailure));
    2169:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)ilink->ksp, (PetscObject)pc, 1));
    2170:   PetscCall(KSPSetType(ilink->ksp, KSPPREONLY));

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

    2175:   if (!next) {
    2176:     jac->head       = ilink;
    2177:     ilink->previous = NULL;
    2178:   } else {
    2179:     while (next->next) next = next->next;
    2180:     next->next      = ilink;
    2181:     ilink->previous = next;
    2182:   }
    2183:   jac->nsplits++;
    2184:   PetscFunctionReturn(PETSC_SUCCESS);
    2185: }

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

    2190:   Logically Collective

    2192:   Input Parameters:
    2193: + pc         - the preconditioner context
    2194: . splitname  - name of this split, if `NULL` the number of the split is used
    2195: . n          - the number of fields in this split
    2196: . fields     - the fields in this split
    2197: - 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
    2198:                of the matrix and `fields_col` provides the column indices for that block

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

    2203:   Level: intermediate

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

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

    2210:   If the matrix used to construct the preconditioner is not `MATNEST` then
    2211:   `PCFieldSplitSetFields()` is for defining fields as strided blocks (based on the block size provided to the matrix with `MatSetBlockSize()` or
    2212:   to the `PC` with `PCFieldSplitSetBlockSize()`). For example, if the block
    2213:   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
    2214:   0xx3xx6xx9xx12 ... x1xx4xx7xx ... xx2xx5xx8xx.. 01x34x67x... 0x23x56x8.. x12x45x78x....
    2215:   where the numbered entries indicate what is in the split.

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

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

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

    2227: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetIS()`, `PCFieldSplitRestrictIS()`,
    2228:           `MatSetBlockSize()`, `MatCreateNest()`
    2229: @*/
    2230: PetscErrorCode PCFieldSplitSetFields(PC pc, const char splitname[], PetscInt n, const PetscInt fields[], const PetscInt fields_col[])
    2231: {
    2232:   PetscFunctionBegin;
    2234:   PetscAssertPointer(splitname, 2);
    2235:   PetscCheck(n >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Provided number of fields %" PetscInt_FMT " in split \"%s\" not positive", n, splitname);
    2236:   PetscAssertPointer(fields, 4);
    2237:   PetscTryMethod(pc, "PCFieldSplitSetFields_C", (PC, const char[], PetscInt, const PetscInt *, const PetscInt *), (pc, splitname, n, fields, fields_col));
    2238:   PetscFunctionReturn(PETSC_SUCCESS);
    2239: }

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

    2245:   Logically Collective

    2247:   Input Parameters:
    2248: + pc  - the preconditioner object
    2249: - flg - boolean flag indicating whether or not to use Amat to extract the diagonal blocks from

    2251:   Options Database Key:
    2252: . -pc_fieldsplit_diag_use_amat - use the Amat to provide the diagonal blocks

    2254:   Level: intermediate

    2256: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetDiagUseAmat()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFIELDSPLIT`
    2257: @*/
    2258: PetscErrorCode PCFieldSplitSetDiagUseAmat(PC pc, PetscBool flg)
    2259: {
    2260:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
    2261:   PetscBool      isfs;

    2263:   PetscFunctionBegin;
    2265:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
    2266:   PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
    2267:   jac->diag_use_amat = flg;
    2268:   PetscFunctionReturn(PETSC_SUCCESS);
    2269: }

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

    2275:   Logically Collective

    2277:   Input Parameter:
    2278: . pc - the preconditioner object

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

    2283:   Level: intermediate

    2285: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetDiagUseAmat()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFIELDSPLIT`
    2286: @*/
    2287: PetscErrorCode PCFieldSplitGetDiagUseAmat(PC pc, PetscBool *flg)
    2288: {
    2289:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
    2290:   PetscBool      isfs;

    2292:   PetscFunctionBegin;
    2294:   PetscAssertPointer(flg, 2);
    2295:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
    2296:   PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
    2297:   *flg = jac->diag_use_amat;
    2298:   PetscFunctionReturn(PETSC_SUCCESS);
    2299: }

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

    2305:   Logically Collective

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

    2311:   Options Database Key:
    2312: . -pc_fieldsplit_off_diag_use_amat <bool> - use the Amat to extract the off-diagonal blocks

    2314:   Level: intermediate

    2316: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitGetOffDiagUseAmat()`, `PCFieldSplitSetDiagUseAmat()`, `PCFIELDSPLIT`
    2317: @*/
    2318: PetscErrorCode PCFieldSplitSetOffDiagUseAmat(PC pc, PetscBool flg)
    2319: {
    2320:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
    2321:   PetscBool      isfs;

    2323:   PetscFunctionBegin;
    2325:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
    2326:   PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
    2327:   jac->offdiag_use_amat = flg;
    2328:   PetscFunctionReturn(PETSC_SUCCESS);
    2329: }

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

    2335:   Logically Collective

    2337:   Input Parameter:
    2338: . pc - the preconditioner object

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

    2343:   Level: intermediate

    2345: .seealso: [](sec_block_matrices), `PC`, `PCSetOperators()`, `KSPSetOperators()`, `PCFieldSplitSetOffDiagUseAmat()`, `PCFieldSplitGetDiagUseAmat()`, `PCFIELDSPLIT`
    2346: @*/
    2347: PetscErrorCode PCFieldSplitGetOffDiagUseAmat(PC pc, PetscBool *flg)
    2348: {
    2349:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
    2350:   PetscBool      isfs;

    2352:   PetscFunctionBegin;
    2354:   PetscAssertPointer(flg, 2);
    2355:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
    2356:   PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PC not of type %s", PCFIELDSPLIT);
    2357:   *flg = jac->offdiag_use_amat;
    2358:   PetscFunctionReturn(PETSC_SUCCESS);
    2359: }

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

    2364:   Logically Collective

    2366:   Input Parameters:
    2367: + pc        - the preconditioner context
    2368: . splitname - name of this split, if `NULL` the number of the split is used
    2369: - is        - the index set that defines the elements in this split

    2371:   Level: intermediate

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

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

    2379: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetBlockSize()`, `PCFieldSplitSetFields()`
    2380: @*/
    2381: PetscErrorCode PCFieldSplitSetIS(PC pc, const char splitname[], IS is)
    2382: {
    2383:   PetscFunctionBegin;
    2385:   if (splitname) PetscAssertPointer(splitname, 2);
    2387:   PetscTryMethod(pc, "PCFieldSplitSetIS_C", (PC, const char[], IS), (pc, splitname, is));
    2388:   PetscFunctionReturn(PETSC_SUCCESS);
    2389: }

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

    2394:   Logically Collective

    2396:   Input Parameters:
    2397: + pc        - the preconditioner context
    2398: - splitname - name of this split

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

    2403:   Level: intermediate

    2405: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetIS()`, `PCFieldSplitGetISByIndex()`
    2406: @*/
    2407: PetscErrorCode PCFieldSplitGetIS(PC pc, const char splitname[], IS *is)
    2408: {
    2409:   PetscFunctionBegin;
    2411:   PetscAssertPointer(splitname, 2);
    2412:   PetscAssertPointer(is, 3);
    2413:   {
    2414:     PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    2415:     PC_FieldSplitLink ilink = jac->head;
    2416:     PetscBool         found;

    2418:     *is = NULL;
    2419:     while (ilink) {
    2420:       PetscCall(PetscStrcmp(ilink->splitname, splitname, &found));
    2421:       if (found) {
    2422:         *is = ilink->is;
    2423:         break;
    2424:       }
    2425:       ilink = ilink->next;
    2426:     }
    2427:   }
    2428:   PetscFunctionReturn(PETSC_SUCCESS);
    2429: }

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

    2434:   Logically Collective

    2436:   Input Parameters:
    2437: + pc    - the preconditioner context
    2438: - index - index of this split

    2440:   Output Parameter:
    2441: . is - the index set that defines the elements in this split

    2443:   Level: intermediate

    2445: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitGetIS()`, `PCFieldSplitSetIS()`,

    2447: @*/
    2448: PetscErrorCode PCFieldSplitGetISByIndex(PC pc, PetscInt index, IS *is)
    2449: {
    2450:   PetscFunctionBegin;
    2451:   PetscCheck(index >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative field %" PetscInt_FMT " requested", index);
    2453:   PetscAssertPointer(is, 3);
    2454:   {
    2455:     PC_FieldSplit    *jac   = (PC_FieldSplit *)pc->data;
    2456:     PC_FieldSplitLink ilink = jac->head;
    2457:     PetscInt          i     = 0;
    2458:     PetscCheck(index < jac->nsplits, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT " requested but only %" PetscInt_FMT " exist", index, jac->nsplits);

    2460:     while (i < index) {
    2461:       ilink = ilink->next;
    2462:       ++i;
    2463:     }
    2464:     PetscCall(PCFieldSplitGetIS(pc, ilink->splitname, is));
    2465:   }
    2466:   PetscFunctionReturn(PETSC_SUCCESS);
    2467: }

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

    2473:   Logically Collective

    2475:   Input Parameters:
    2476: + pc - the preconditioner context
    2477: - bs - the block size

    2479:   Level: intermediate

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

    2484: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
    2485: @*/
    2486: PetscErrorCode PCFieldSplitSetBlockSize(PC pc, PetscInt bs)
    2487: {
    2488:   PetscFunctionBegin;
    2491:   PetscTryMethod(pc, "PCFieldSplitSetBlockSize_C", (PC, PetscInt), (pc, bs));
    2492:   PetscFunctionReturn(PETSC_SUCCESS);
    2493: }

    2495: /*@C
    2496:   PCFieldSplitGetSubKSP - Gets the `KSP` contexts for all splits

    2498:   Collective

    2500:   Input Parameter:
    2501: . pc - the preconditioner context

    2503:   Output Parameters:
    2504: + n      - the number of splits
    2505: - subksp - the array of `KSP` contexts

    2507:   Level: advanced

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

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

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

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

    2522:   Fortran Notes:
    2523:   You must pass in a `KSP` array that is large enough to contain all the `KSP`s.
    2524:   You can call `PCFieldSplitGetSubKSP`(pc,n,`PETSC_NULL_KSP`,ierr) to determine how large the
    2525:   `KSP` array must be.

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

    2530:   The Fortran interface could be modernized to return directly the array of values.

    2532: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitSchurGetSubKSP()`
    2533: @*/
    2534: PetscErrorCode PCFieldSplitGetSubKSP(PC pc, PetscInt *n, KSP *subksp[])
    2535: {
    2536:   PetscFunctionBegin;
    2538:   if (n) PetscAssertPointer(n, 2);
    2539:   PetscUseMethod(pc, "PCFieldSplitGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp));
    2540:   PetscFunctionReturn(PETSC_SUCCESS);
    2541: }

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

    2546:   Collective

    2548:   Input Parameter:
    2549: . pc - the preconditioner context

    2551:   Output Parameters:
    2552: + n      - the number of splits
    2553: - subksp - the array of `KSP` contexts

    2555:   Level: advanced

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

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

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

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

    2570:   Fortran Notes:
    2571:   You must pass in a `KSP` array that is large enough to contain all the local `KSP`s.
    2572:   You can call `PCFieldSplitSchurGetSubKSP`(pc,n,`PETSC_NULL_KSP`,ierr) to determine how large the
    2573:   `KSP` array must be.

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

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

    2580:   The Fortran interface should be modernized to return directly the array of values.

    2582: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`, `PCFieldSplitGetSubKSP()`
    2583: @*/
    2584: PetscErrorCode PCFieldSplitSchurGetSubKSP(PC pc, PetscInt *n, KSP *subksp[])
    2585: {
    2586:   PetscFunctionBegin;
    2588:   if (n) PetscAssertPointer(n, 2);
    2589:   PetscUseMethod(pc, "PCFieldSplitSchurGetSubKSP_C", (PC, PetscInt *, KSP **), (pc, n, subksp));
    2590:   PetscFunctionReturn(PETSC_SUCCESS);
    2591: }

    2593: /*@
    2594:   PCFieldSplitSetSchurPre -  Indicates from what operator the preconditioner is constructed for the Schur complement.
    2595:   The default is the A11 matrix.

    2597:   Collective

    2599:   Input Parameters:
    2600: + pc    - the preconditioner context
    2601: . ptype - which matrix to use for preconditioning the Schur complement: `PC_FIELDSPLIT_SCHUR_PRE_A11` (default),
    2602:               `PC_FIELDSPLIT_SCHUR_PRE_SELF`, `PC_FIELDSPLIT_SCHUR_PRE_USER`,
    2603:               `PC_FIELDSPLIT_SCHUR_PRE_SELFP`, and `PC_FIELDSPLIT_SCHUR_PRE_FULL`
    2604: - pre   - matrix to use for preconditioning, or `NULL`

    2606:   Options Database Keys:
    2607: + -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full> - default is `a11`. See notes for meaning of various arguments
    2608: - -fieldsplit_1_pc_type <pctype>                               - the preconditioner algorithm that is used to construct the preconditioner from the operator

    2610:   Level: intermediate

    2612:   Notes:
    2613:   If ptype is
    2614: +     a11 - the preconditioner for the Schur complement is generated from the block diagonal part of the preconditioner
    2615:   matrix associated with the Schur complement (i.e. A11), not the Schur complement matrix
    2616: .     self - the preconditioner for the Schur complement is generated from the symbolic representation of the Schur complement matrix:
    2617:   The only preconditioners that currently work with this symbolic representation matrix object are `PCLSC` and `PCHPDDM`
    2618: .     user - the preconditioner for the Schur complement is generated from the user provided matrix (pre argument
    2619:   to this function).
    2620: .     selfp - the preconditioning for the Schur complement is generated from an explicitly-assembled approximation $ Sp = A11 - A10 inv(diag(A00)) A01 $
    2621:   This is only a good preconditioner when diag(A00) is a good preconditioner for A00. Optionally, A00 can be
    2622:   lumped before extracting the diagonal using the additional option `-fieldsplit_1_mat_schur_complement_ainv_type lump`
    2623: -     full - the preconditioner for the Schur complement is generated from the exact Schur complement matrix representation
    2624:   computed internally by `PCFIELDSPLIT` (this is expensive)
    2625:   useful mostly as a test that the Schur complement approach can work for your problem

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

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

    2634: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`,
    2635:           `MatSchurComplementSetAinvType()`, `PCLSC`, `PCFieldSplitSetSchurFactType()`
    2636: @*/
    2637: PetscErrorCode PCFieldSplitSetSchurPre(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
    2638: {
    2639:   PetscFunctionBegin;
    2641:   PetscTryMethod(pc, "PCFieldSplitSetSchurPre_C", (PC, PCFieldSplitSchurPreType, Mat), (pc, ptype, pre));
    2642:   PetscFunctionReturn(PETSC_SUCCESS);
    2643: }

    2645: PetscErrorCode PCFieldSplitSchurPrecondition(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
    2646: {
    2647:   return PCFieldSplitSetSchurPre(pc, ptype, pre);
    2648: } /* Deprecated name */

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

    2654:   Logically Collective

    2656:   Input Parameter:
    2657: . pc - the preconditioner context

    2659:   Output Parameters:
    2660: + 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`
    2661: - pre   - matrix to use for preconditioning (with `PC_FIELDSPLIT_SCHUR_PRE_USER`), or `NULL`

    2663:   Level: intermediate

    2665: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCLSC`
    2666: @*/
    2667: PetscErrorCode PCFieldSplitGetSchurPre(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre)
    2668: {
    2669:   PetscFunctionBegin;
    2671:   PetscUseMethod(pc, "PCFieldSplitGetSchurPre_C", (PC, PCFieldSplitSchurPreType *, Mat *), (pc, ptype, pre));
    2672:   PetscFunctionReturn(PETSC_SUCCESS);
    2673: }

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

    2678:   Not Collective

    2680:   Input Parameter:
    2681: . pc - the preconditioner context

    2683:   Output Parameter:
    2684: . S - the Schur complement matrix

    2686:   Level: advanced

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

    2691: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MATSCHURCOMPLEMENT`, `PCFieldSplitSchurRestoreS()`,
    2692:           `MatCreateSchurComplement()`, `MatSchurComplementGetKSP()`, `MatSchurComplementComputeExplicitOperator()`, `MatGetSchurComplement()`
    2693: @*/
    2694: PetscErrorCode PCFieldSplitSchurGetS(PC pc, Mat *S)
    2695: {
    2696:   const char    *t;
    2697:   PetscBool      isfs;
    2698:   PC_FieldSplit *jac;

    2700:   PetscFunctionBegin;
    2702:   PetscCall(PetscObjectGetType((PetscObject)pc, &t));
    2703:   PetscCall(PetscStrcmp(t, PCFIELDSPLIT, &isfs));
    2704:   PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PC of type PCFIELDSPLIT, got %s instead", t);
    2705:   jac = (PC_FieldSplit *)pc->data;
    2706:   PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PCFIELDSPLIT of type SCHUR, got %d instead", jac->type);
    2707:   if (S) *S = jac->schur;
    2708:   PetscFunctionReturn(PETSC_SUCCESS);
    2709: }

    2711: /*@
    2712:   PCFieldSplitSchurRestoreS -  returns the `MATSCHURCOMPLEMENT` matrix used by this `PC`

    2714:   Not Collective

    2716:   Input Parameters:
    2717: + pc - the preconditioner context
    2718: - S  - the Schur complement matrix

    2720:   Level: advanced

    2722: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurPre()`, `MatSchurComplement`, `PCFieldSplitSchurGetS()`
    2723: @*/
    2724: PetscErrorCode PCFieldSplitSchurRestoreS(PC pc, Mat *S)
    2725: {
    2726:   const char    *t;
    2727:   PetscBool      isfs;
    2728:   PC_FieldSplit *jac;

    2730:   PetscFunctionBegin;
    2732:   PetscCall(PetscObjectGetType((PetscObject)pc, &t));
    2733:   PetscCall(PetscStrcmp(t, PCFIELDSPLIT, &isfs));
    2734:   PetscCheck(isfs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PC of type PCFIELDSPLIT, got %s instead", t);
    2735:   jac = (PC_FieldSplit *)pc->data;
    2736:   PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Expected PCFIELDSPLIT of type SCHUR, got %d instead", jac->type);
    2737:   PetscCheck(S && (*S == jac->schur), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "MatSchurComplement restored is not the same as gotten");
    2738:   PetscFunctionReturn(PETSC_SUCCESS);
    2739: }

    2741: static PetscErrorCode PCFieldSplitSetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType ptype, Mat pre)
    2742: {
    2743:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    2745:   PetscFunctionBegin;
    2746:   jac->schurpre = ptype;
    2747:   if (ptype == PC_FIELDSPLIT_SCHUR_PRE_USER && pre) {
    2748:     PetscCall(MatDestroy(&jac->schur_user));
    2749:     jac->schur_user = pre;
    2750:     PetscCall(PetscObjectReference((PetscObject)jac->schur_user));
    2751:   }
    2752:   PetscFunctionReturn(PETSC_SUCCESS);
    2753: }

    2755: static PetscErrorCode PCFieldSplitGetSchurPre_FieldSplit(PC pc, PCFieldSplitSchurPreType *ptype, Mat *pre)
    2756: {
    2757:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    2759:   PetscFunctionBegin;
    2760:   if (ptype) *ptype = jac->schurpre;
    2761:   if (pre) *pre = jac->schur_user;
    2762:   PetscFunctionReturn(PETSC_SUCCESS);
    2763: }

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

    2768:   Collective

    2770:   Input Parameters:
    2771: + pc    - the preconditioner context
    2772: - ftype - which blocks of factorization to retain, `PC_FIELDSPLIT_SCHUR_FACT_FULL` is default

    2774:   Options Database Key:
    2775: . -pc_fieldsplit_schur_fact_type <diag,lower,upper,full> - default is `full`

    2777:   Level: intermediate

    2779:   Notes:
    2780:   The FULL factorization is

    2782:   ```{math}
    2783:   \left(\begin{array}{cc} A & B \\
    2784:   C & E \\
    2785:   \end{array}\right) =
    2786:   \left(\begin{array}{cc} 1 & 0 \\
    2787:   C*A^{-1} & I \\
    2788:   \end{array}\right)
    2789:   \left(\begin{array}{cc} A & 0 \\
    2790:   0 & S \\
    2791:   \end{array}\right)
    2792:   \left(\begin{array}{cc} I & A^{-1}B \\
    2793:   0 & I \\
    2794:   \end{array}\right) = L D U.
    2795:   ```

    2797:   where $ S = E - C*A^{-1}*B $. In practice, the full factorization is applied via block triangular solves with the grouping $L*(D*U)$. UPPER uses $D*U$, LOWER uses $L*D$,
    2798:   and DIAG is the diagonal part with the sign of $ S $ flipped (because this makes the preconditioner positive definite for many formulations,
    2799:   thus allowing the use of `KSPMINRES)`. Sign flipping of $ S $ can be turned off with `PCFieldSplitSetSchurScale()`.

    2801:   If $A$ and $S$ are solved exactly
    2802: +  1 - FULL factorization is a direct solver.
    2803: .  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.
    2804: -  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.

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

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

    2811:   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).

    2813: .seealso: [](sec_block_matrices), `PC`, `PCFieldSplitGetSubKSP()`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurPreType`, `PCFieldSplitSetSchurScale()`,
    2814:           [](sec_flexibleksp), `PCFieldSplitSetSchurPre()`
    2815: @*/
    2816: PetscErrorCode PCFieldSplitSetSchurFactType(PC pc, PCFieldSplitSchurFactType ftype)
    2817: {
    2818:   PetscFunctionBegin;
    2820:   PetscTryMethod(pc, "PCFieldSplitSetSchurFactType_C", (PC, PCFieldSplitSchurFactType), (pc, ftype));
    2821:   PetscFunctionReturn(PETSC_SUCCESS);
    2822: }

    2824: static PetscErrorCode PCFieldSplitSetSchurFactType_FieldSplit(PC pc, PCFieldSplitSchurFactType ftype)
    2825: {
    2826:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    2828:   PetscFunctionBegin;
    2829:   jac->schurfactorization = ftype;
    2830:   PetscFunctionReturn(PETSC_SUCCESS);
    2831: }

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

    2836:   Collective

    2838:   Input Parameters:
    2839: + pc    - the preconditioner context
    2840: - scale - scaling factor for the Schur complement

    2842:   Options Database Key:
    2843: . -pc_fieldsplit_schur_scale <scale> - default is -1.0

    2845:   Level: intermediate

    2847: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetFields()`, `PCFieldSplitSchurFactType`, `PCFieldSplitSetSchurFactType()`
    2848: @*/
    2849: PetscErrorCode PCFieldSplitSetSchurScale(PC pc, PetscScalar scale)
    2850: {
    2851:   PetscFunctionBegin;
    2854:   PetscTryMethod(pc, "PCFieldSplitSetSchurScale_C", (PC, PetscScalar), (pc, scale));
    2855:   PetscFunctionReturn(PETSC_SUCCESS);
    2856: }

    2858: static PetscErrorCode PCFieldSplitSetSchurScale_FieldSplit(PC pc, PetscScalar scale)
    2859: {
    2860:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    2862:   PetscFunctionBegin;
    2863:   jac->schurscale = scale;
    2864:   PetscFunctionReturn(PETSC_SUCCESS);
    2865: }

    2867: /*@C
    2868:   PCFieldSplitGetSchurBlocks - Gets all matrix blocks for the Schur complement

    2870:   Collective

    2872:   Input Parameter:
    2873: . pc - the preconditioner context

    2875:   Output Parameters:
    2876: + A00 - the (0,0) block
    2877: . A01 - the (0,1) block
    2878: . A10 - the (1,0) block
    2879: - A11 - the (1,1) block

    2881:   Level: advanced

    2883:   Note:
    2884:   Use `NULL` for any unneeded output arguments

    2886: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `MatSchurComplementGetSubMatrices()`, `MatSchurComplementSetSubMatrices()`
    2887: @*/
    2888: PetscErrorCode PCFieldSplitGetSchurBlocks(PC pc, Mat *A00, Mat *A01, Mat *A10, Mat *A11)
    2889: {
    2890:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    2892:   PetscFunctionBegin;
    2894:   PetscCheck(jac->type == PC_COMPOSITE_SCHUR, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "FieldSplit is not using a Schur complement approach.");
    2895:   if (A00) *A00 = jac->pmat[0];
    2896:   if (A01) *A01 = jac->B;
    2897:   if (A10) *A10 = jac->C;
    2898:   if (A11) *A11 = jac->pmat[1];
    2899:   PetscFunctionReturn(PETSC_SUCCESS);
    2900: }

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

    2905:   Collective

    2907:   Input Parameters:
    2908: + pc        - the preconditioner context
    2909: - tolerance - the solver tolerance

    2911:   Options Database Key:
    2912: . -pc_fieldsplit_gkb_tol <tolerance> - default is 1e-5

    2914:   Level: intermediate

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

    2921: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBMaxit()`
    2922: @*/
    2923: PetscErrorCode PCFieldSplitSetGKBTol(PC pc, PetscReal tolerance)
    2924: {
    2925:   PetscFunctionBegin;
    2928:   PetscTryMethod(pc, "PCFieldSplitSetGKBTol_C", (PC, PetscReal), (pc, tolerance));
    2929:   PetscFunctionReturn(PETSC_SUCCESS);
    2930: }

    2932: static PetscErrorCode PCFieldSplitSetGKBTol_FieldSplit(PC pc, PetscReal tolerance)
    2933: {
    2934:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    2936:   PetscFunctionBegin;
    2937:   jac->gkbtol = tolerance;
    2938:   PetscFunctionReturn(PETSC_SUCCESS);
    2939: }

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

    2944:   Collective

    2946:   Input Parameters:
    2947: + pc    - the preconditioner context
    2948: - maxit - the maximum number of iterations

    2950:   Options Database Key:
    2951: . -pc_fieldsplit_gkb_maxit <maxit> - default is 100

    2953:   Level: intermediate

    2955: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBNu()`
    2956: @*/
    2957: PetscErrorCode PCFieldSplitSetGKBMaxit(PC pc, PetscInt maxit)
    2958: {
    2959:   PetscFunctionBegin;
    2962:   PetscTryMethod(pc, "PCFieldSplitSetGKBMaxit_C", (PC, PetscInt), (pc, maxit));
    2963:   PetscFunctionReturn(PETSC_SUCCESS);
    2964: }

    2966: static PetscErrorCode PCFieldSplitSetGKBMaxit_FieldSplit(PC pc, PetscInt maxit)
    2967: {
    2968:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    2970:   PetscFunctionBegin;
    2971:   jac->gkbmaxit = maxit;
    2972:   PetscFunctionReturn(PETSC_SUCCESS);
    2973: }

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

    2979:   Collective

    2981:   Input Parameters:
    2982: + pc    - the preconditioner context
    2983: - delay - the delay window in the lower bound estimate

    2985:   Options Database Key:
    2986: . -pc_fieldsplit_gkb_delay <delay> - default is 5

    2988:   Level: intermediate

    2990:   Notes:
    2991:   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 $
    2992:   is expressed as a truncated sum. The error at iteration k can only be measured at iteration (k + `delay`), and thus the algorithm needs
    2993:   at least (`delay` + 1) iterations to stop.

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

    2997: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBNu()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()`
    2998: @*/
    2999: PetscErrorCode PCFieldSplitSetGKBDelay(PC pc, PetscInt delay)
    3000: {
    3001:   PetscFunctionBegin;
    3004:   PetscTryMethod(pc, "PCFieldSplitSetGKBDelay_C", (PC, PetscInt), (pc, delay));
    3005:   PetscFunctionReturn(PETSC_SUCCESS);
    3006: }

    3008: static PetscErrorCode PCFieldSplitSetGKBDelay_FieldSplit(PC pc, PetscInt delay)
    3009: {
    3010:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3012:   PetscFunctionBegin;
    3013:   jac->gkbdelay = delay;
    3014:   PetscFunctionReturn(PETSC_SUCCESS);
    3015: }

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

    3021:   Collective

    3023:   Input Parameters:
    3024: + pc - the preconditioner context
    3025: - nu - the shift parameter

    3027:   Options Database Key:
    3028: . -pc_fieldsplit_gkb_nu <nu> - default is 1

    3030:   Level: intermediate

    3032:   Notes:
    3033:   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,
    3034:   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
    3035:   necessary to find a good balance in between the convergence of the inner and outer loop.

    3037:   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.

    3039: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetGKBDelay()`, `PCFieldSplitSetGKBTol()`, `PCFieldSplitSetGKBMaxit()`
    3040: @*/
    3041: PetscErrorCode PCFieldSplitSetGKBNu(PC pc, PetscReal nu)
    3042: {
    3043:   PetscFunctionBegin;
    3046:   PetscTryMethod(pc, "PCFieldSplitSetGKBNu_C", (PC, PetscReal), (pc, nu));
    3047:   PetscFunctionReturn(PETSC_SUCCESS);
    3048: }

    3050: static PetscErrorCode PCFieldSplitSetGKBNu_FieldSplit(PC pc, PetscReal nu)
    3051: {
    3052:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3054:   PetscFunctionBegin;
    3055:   jac->gkbnu = nu;
    3056:   PetscFunctionReturn(PETSC_SUCCESS);
    3057: }

    3059: static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc, PCCompositeType type)
    3060: {
    3061:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3063:   PetscFunctionBegin;
    3064:   jac->type = type;
    3065:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", NULL));
    3066:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", NULL));
    3067:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", NULL));
    3068:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", NULL));
    3069:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", NULL));
    3070:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", NULL));
    3071:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", NULL));
    3072:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", NULL));
    3073:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL));

    3075:   if (type == PC_COMPOSITE_SCHUR) {
    3076:     pc->ops->apply          = PCApply_FieldSplit_Schur;
    3077:     pc->ops->applytranspose = PCApplyTranspose_FieldSplit_Schur;
    3078:     pc->ops->view           = PCView_FieldSplit_Schur;
    3079:     pc->ops->setuponblocks  = PCSetUpOnBlocks_FieldSplit_Schur;

    3081:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit_Schur));
    3082:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", PCFieldSplitSetSchurPre_FieldSplit));
    3083:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSchurPre_C", PCFieldSplitGetSchurPre_FieldSplit));
    3084:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurFactType_C", PCFieldSplitSetSchurFactType_FieldSplit));
    3085:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurScale_C", PCFieldSplitSetSchurScale_FieldSplit));
    3086:   } else if (type == PC_COMPOSITE_GKB) {
    3087:     pc->ops->apply          = PCApply_FieldSplit_GKB;
    3088:     pc->ops->applytranspose = NULL;
    3089:     pc->ops->view           = PCView_FieldSplit_GKB;
    3090:     pc->ops->setuponblocks  = PCSetUpOnBlocks_FieldSplit_GKB;

    3092:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit));
    3093:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBTol_C", PCFieldSplitSetGKBTol_FieldSplit));
    3094:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBMaxit_C", PCFieldSplitSetGKBMaxit_FieldSplit));
    3095:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBNu_C", PCFieldSplitSetGKBNu_FieldSplit));
    3096:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", PCFieldSplitSetGKBDelay_FieldSplit));
    3097:   } else {
    3098:     pc->ops->apply          = PCApply_FieldSplit;
    3099:     pc->ops->applytranspose = PCApplyTranspose_FieldSplit;
    3100:     pc->ops->view           = PCView_FieldSplit;
    3101:     pc->ops->setuponblocks  = PCSetUpOnBlocks_FieldSplit;

    3103:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit));
    3104:   }
    3105:   PetscFunctionReturn(PETSC_SUCCESS);
    3106: }

    3108: static PetscErrorCode PCFieldSplitSetBlockSize_FieldSplit(PC pc, PetscInt bs)
    3109: {
    3110:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3112:   PetscFunctionBegin;
    3113:   PetscCheck(bs >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Blocksize must be positive, you gave %" PetscInt_FMT, bs);
    3114:   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);
    3115:   jac->bs = bs;
    3116:   PetscFunctionReturn(PETSC_SUCCESS);
    3117: }

    3119: static PetscErrorCode PCSetCoordinates_FieldSplit(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
    3120: {
    3121:   PC_FieldSplit    *jac           = (PC_FieldSplit *)pc->data;
    3122:   PC_FieldSplitLink ilink_current = jac->head;
    3123:   IS                is_owned;

    3125:   PetscFunctionBegin;
    3126:   jac->coordinates_set = PETSC_TRUE; // Internal flag
    3127:   PetscCall(MatGetOwnershipIS(pc->mat, &is_owned, NULL));

    3129:   while (ilink_current) {
    3130:     // For each IS, embed it to get local coords indces
    3131:     IS              is_coords;
    3132:     PetscInt        ndofs_block;
    3133:     const PetscInt *block_dofs_enumeration; // Numbering of the dofs relevant to the current block

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

    3140:     // Allocate coordinates vector and set it directly
    3141:     PetscCall(PetscMalloc1(ndofs_block * dim, &ilink_current->coords));
    3142:     for (PetscInt dof = 0; dof < ndofs_block; ++dof) {
    3143:       for (PetscInt d = 0; d < dim; ++d) (ilink_current->coords)[dim * dof + d] = coords[dim * block_dofs_enumeration[dof] + d];
    3144:     }
    3145:     ilink_current->dim   = dim;
    3146:     ilink_current->ndofs = ndofs_block;
    3147:     PetscCall(ISRestoreIndices(is_coords, &block_dofs_enumeration));
    3148:     PetscCall(ISDestroy(&is_coords));
    3149:     ilink_current = ilink_current->next;
    3150:   }
    3151:   PetscCall(ISDestroy(&is_owned));
    3152:   PetscFunctionReturn(PETSC_SUCCESS);
    3153: }

    3155: /*@
    3156:   PCFieldSplitSetType - Sets the type, `PCCompositeType`, of a `PCFIELDSPLIT`

    3158:   Collective

    3160:   Input Parameters:
    3161: + pc   - the preconditioner context
    3162: - type - `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE` (default), `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`,
    3163:          `PC_COMPOSITE_GKB`

    3165:   Options Database Key:
    3166: . -pc_fieldsplit_type <one of multiplicative, additive, symmetric_multiplicative, special, schur> - Sets fieldsplit preconditioner type

    3168:   Level: intermediate

    3170: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCCompositeType`, `PCCompositeGetType()`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`,
    3171:           `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`, `PCFieldSplitSetSchurFactType()`
    3172: @*/
    3173: PetscErrorCode PCFieldSplitSetType(PC pc, PCCompositeType type)
    3174: {
    3175:   PetscFunctionBegin;
    3177:   PetscTryMethod(pc, "PCFieldSplitSetType_C", (PC, PCCompositeType), (pc, type));
    3178:   PetscFunctionReturn(PETSC_SUCCESS);
    3179: }

    3181: /*@
    3182:   PCFieldSplitGetType - Gets the type, `PCCompositeType`, of a `PCFIELDSPLIT`

    3184:   Not collective

    3186:   Input Parameter:
    3187: . pc - the preconditioner context

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

    3192:   Level: intermediate

    3194: .seealso: [](sec_block_matrices), `PC`, `PCCompositeSetType()`, `PCFIELDSPLIT`, `PCCompositeType`, `PC_COMPOSITE_ADDITIVE`, `PC_COMPOSITE_MULTIPLICATIVE`,
    3195:           `PC_COMPOSITE_SYMMETRIC_MULTIPLICATIVE`, `PC_COMPOSITE_SPECIAL`, `PC_COMPOSITE_SCHUR`
    3196: @*/
    3197: PetscErrorCode PCFieldSplitGetType(PC pc, PCCompositeType *type)
    3198: {
    3199:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3201:   PetscFunctionBegin;
    3203:   PetscAssertPointer(type, 2);
    3204:   *type = jac->type;
    3205:   PetscFunctionReturn(PETSC_SUCCESS);
    3206: }

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

    3211:   Logically Collective

    3213:   Input Parameters:
    3214: + pc  - the preconditioner context
    3215: - flg - boolean indicating whether to use field splits defined by the `DM`

    3217:   Options Database Key:
    3218: . -pc_fieldsplit_dm_splits <bool> - use the field splits defined by the `DM`

    3220:   Level: intermediate

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

    3225: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDMSplits()`, `DMCreateFieldDecomposition()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
    3226: @*/
    3227: PetscErrorCode PCFieldSplitSetDMSplits(PC pc, PetscBool flg)
    3228: {
    3229:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
    3230:   PetscBool      isfs;

    3232:   PetscFunctionBegin;
    3235:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
    3236:   if (isfs) jac->dm_splits = flg;
    3237:   PetscFunctionReturn(PETSC_SUCCESS);
    3238: }

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

    3243:   Logically Collective

    3245:   Input Parameter:
    3246: . pc - the preconditioner context

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

    3251:   Level: intermediate

    3253:   Developer Note:
    3254:   The name should be `PCFieldSplitGetUseDMSplits()`

    3256: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDMSplits()`, `DMCreateFieldDecomposition()`, `PCFieldSplitSetFields()`, `PCFieldSplitSetIS()`
    3257: @*/
    3258: PetscErrorCode PCFieldSplitGetDMSplits(PC pc, PetscBool *flg)
    3259: {
    3260:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
    3261:   PetscBool      isfs;

    3263:   PetscFunctionBegin;
    3265:   PetscAssertPointer(flg, 2);
    3266:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &isfs));
    3267:   if (isfs) {
    3268:     if (flg) *flg = jac->dm_splits;
    3269:   }
    3270:   PetscFunctionReturn(PETSC_SUCCESS);
    3271: }

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

    3276:   Logically Collective

    3278:   Input Parameter:
    3279: . pc - the preconditioner context

    3281:   Output Parameter:
    3282: . flg - boolean indicating whether to detect fields or not

    3284:   Level: intermediate

    3286: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitSetDetectSaddlePoint()`
    3287: @*/
    3288: PetscErrorCode PCFieldSplitGetDetectSaddlePoint(PC pc, PetscBool *flg)
    3289: {
    3290:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3292:   PetscFunctionBegin;
    3293:   *flg = jac->detect;
    3294:   PetscFunctionReturn(PETSC_SUCCESS);
    3295: }

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

    3300:   Logically Collective

    3302:   Input Parameter:
    3303: . pc - the preconditioner context

    3305:   Output Parameter:
    3306: . flg - boolean indicating whether to detect fields or not

    3308:   Options Database Key:
    3309: . -pc_fieldsplit_detect_saddle_point <bool> - detect and use the saddle point

    3311:   Level: intermediate

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

    3316: .seealso: [](sec_block_matrices), `PC`, `PCFIELDSPLIT`, `PCFieldSplitGetDetectSaddlePoint()`, `PCFieldSplitSetType()`, `PCFieldSplitSetSchurPre()`, `PC_FIELDSPLIT_SCHUR_PRE_SELF`
    3317: @*/
    3318: PetscErrorCode PCFieldSplitSetDetectSaddlePoint(PC pc, PetscBool flg)
    3319: {
    3320:   PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;

    3322:   PetscFunctionBegin;
    3323:   jac->detect = flg;
    3324:   if (jac->detect) {
    3325:     PetscCall(PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR));
    3326:     PetscCall(PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_SELF, NULL));
    3327:   }
    3328:   PetscFunctionReturn(PETSC_SUCCESS);
    3329: }

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

    3338:   Options Database Keys:
    3339: +   -pc_fieldsplit_%d_fields <a,b,..>                                                - indicates the fields to be used in the `%d`'th split
    3340: .   -pc_fieldsplit_default                                                           - automatically add any fields to additional splits that have not
    3341:                                                                                        been supplied explicitly by `-pc_fieldsplit_%d_fields`
    3342: .   -pc_fieldsplit_block_size <bs>                                                   - size of block that defines fields (i.e. there are bs fields)
    3343:                                                                                        when the matrix is not of `MatType` `MATNEST`
    3344: .   -pc_fieldsplit_type <additive,multiplicative,symmetric_multiplicative,schur,gkb> - type of relaxation or factorization splitting
    3345: .   -pc_fieldsplit_schur_precondition <self,selfp,user,a11,full>                     - default is `a11`; see `PCFieldSplitSetSchurPre()`
    3346: .   -pc_fieldsplit_schur_fact_type <diag,lower,upper,full>                           - set factorization type when using `-pc_fieldsplit_type schur`;
    3347:                                                                                        see `PCFieldSplitSetSchurFactType()`
    3348: .   -pc_fieldsplit_dm_splits <true,false> (default is true)                          - Whether to use `DMCreateFieldDecomposition()` for splits
    3349: -   -pc_fieldsplit_detect_saddle_point                                               - automatically finds rows with zero diagonal and uses Schur complement with no preconditioner as the solver

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

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

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

    3361:   Level: intermediate

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

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

    3373:   For the Schur complement preconditioner if
    3374:   ```{math}
    3375:     J = \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{10} & A_{11} \end{array}\right]
    3376:   ```

    3378:   the preconditioner using `full` factorization is logically
    3379:   ```{math}
    3380:     \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]
    3381:       ```
    3382:   where the action of $\text{ksp}(A_{00})$ is applied using the `KSP` solver with prefix `-fieldsplit_0_`.  $S$ is the Schur complement
    3383:   ```{math}
    3384:      S = A_{11} - A_{10} \text{ksp}(A_{00}) A_{01}
    3385:   ```
    3386:   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`
    3387:   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
    3388:   0th index, the `KSP` associated with `-fieldsplit_0_`, and at its 1st index, the `KSP` corresponding to `-fieldsplit_1_`.
    3389:   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$.

    3391:   The factorization type is set using `-pc_fieldsplit_schur_fact_type <diag, lower, upper, full>`. `full` is shown above,
    3392:   `diag` gives
    3393:   ```{math}
    3394:     \left[\begin{array}{cc} \text{ksp}(A_{00}) & 0 \\  0 & -\text{ksp}(S) \end{array}\right]
    3395:   ```
    3396:   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
    3397:   can be turned off with `PCFieldSplitSetSchurScale()` or by command line `-pc_fieldsplit_schur_scale 1.0`. The `lower` factorization is the inverse of
    3398:   ```{math}
    3399:     \left[\begin{array}{cc} A_{00} & 0 \\  A_{10} & S \end{array}\right]
    3400:   ```
    3401:   where the inverses of $A_{00}$ and $S$ are applied using `KSP`s. The upper factorization is the inverse of
    3402:   ```{math}
    3403:     \left[\begin{array}{cc} A_{00} & A_{01} \\  0 & S \end{array}\right]
    3404:   ```
    3405:   where again the inverses of $A_{00}$ and $S$ are applied using `KSP`s.

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

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

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

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

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

    3422:   The generalized Golub-Kahan bidiagonalization preconditioner (GKB) can be applied to symmetric $2 \times 2$ block matrices of the shape
    3423:   ```{math}
    3424:     \left[\begin{array}{cc} A_{00} & A_{01} \\ A_{01}' & 0 \end{array}\right]
    3425:   ```
    3426:   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}'$.
    3427:   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_`.

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

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

    3436: .seealso: [](sec_block_matrices), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCLSC`,
    3437:           `PCFieldSplitGetSubKSP()`, `PCFieldSplitSchurGetSubKSP()`, `PCFieldSplitSetFields()`,
    3438:           `PCFieldSplitSetType()`, `PCFieldSplitSetIS()`, `PCFieldSplitSetSchurPre()`, `PCFieldSplitSetSchurFactType()`,
    3439:           `MatSchurComplementSetAinvType()`, `PCFieldSplitSetSchurScale()`, `PCFieldSplitSetDetectSaddlePoint()`
    3440: M*/

    3442: PETSC_EXTERN PetscErrorCode PCCreate_FieldSplit(PC pc)
    3443: {
    3444:   PC_FieldSplit *jac;

    3446:   PetscFunctionBegin;
    3447:   PetscCall(PetscNew(&jac));

    3449:   jac->bs                 = -1;
    3450:   jac->type               = PC_COMPOSITE_MULTIPLICATIVE;
    3451:   jac->schurpre           = PC_FIELDSPLIT_SCHUR_PRE_USER; /* Try user preconditioner first, fall back on diagonal */
    3452:   jac->schurfactorization = PC_FIELDSPLIT_SCHUR_FACT_FULL;
    3453:   jac->schurscale         = -1.0;
    3454:   jac->dm_splits          = PETSC_TRUE;
    3455:   jac->gkbtol             = 1e-5;
    3456:   jac->gkbdelay           = 5;
    3457:   jac->gkbnu              = 1;
    3458:   jac->gkbmaxit           = 100;

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

    3462:   pc->ops->setup           = PCSetUp_FieldSplit;
    3463:   pc->ops->reset           = PCReset_FieldSplit;
    3464:   pc->ops->destroy         = PCDestroy_FieldSplit;
    3465:   pc->ops->setfromoptions  = PCSetFromOptions_FieldSplit;
    3466:   pc->ops->applyrichardson = NULL;

    3468:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSchurGetSubKSP_C", PCFieldSplitSchurGetSubKSP_FieldSplit));
    3469:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetFields_C", PCFieldSplitSetFields_FieldSplit));
    3470:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetIS_C", PCFieldSplitSetIS_FieldSplit));
    3471:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetType_C", PCFieldSplitSetType_FieldSplit));
    3472:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetBlockSize_C", PCFieldSplitSetBlockSize_FieldSplit));
    3473:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitRestrictIS_C", PCFieldSplitRestrictIS_FieldSplit));
    3474:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_FieldSplit));

    3476:   /* Initialize function pointers */
    3477:   PetscCall(PCFieldSplitSetType(pc, jac->type));
    3478:   PetscFunctionReturn(PETSC_SUCCESS);
    3479: }