Actual source code: fetidp.c

      1: #include <petsc/private/kspimpl.h>
      2: #include <petsc/private/pcbddcimpl.h>
      3: #include <petsc/private/pcbddcprivateimpl.h>
      4: #include <petscdm.h>

      6: static PetscBool  cited       = PETSC_FALSE;
      7: static PetscBool  cited2      = PETSC_FALSE;
      8: static const char citation[]  = "@article{ZampiniPCBDDC,\n"
      9:                                 "author = {Stefano Zampini},\n"
     10:                                 "title = {{PCBDDC}: A Class of Robust Dual-Primal Methods in {PETS}c},\n"
     11:                                 "journal = {SIAM Journal on Scientific Computing},\n"
     12:                                 "volume = {38},\n"
     13:                                 "number = {5},\n"
     14:                                 "pages = {S282-S306},\n"
     15:                                 "year = {2016},\n"
     16:                                 "doi = {10.1137/15M1025785},\n"
     17:                                 "URL = {http://dx.doi.org/10.1137/15M1025785},\n"
     18:                                 "eprint = {http://dx.doi.org/10.1137/15M1025785}\n"
     19:                                 "}\n"
     20:                                 "@article{ZampiniDualPrimal,\n"
     21:                                 "author = {Stefano Zampini},\n"
     22:                                 "title = {{D}ual-{P}rimal methods for the cardiac {B}idomain model},\n"
     23:                                 "volume = {24},\n"
     24:                                 "number = {04},\n"
     25:                                 "pages = {667-696},\n"
     26:                                 "year = {2014},\n"
     27:                                 "doi = {10.1142/S0218202513500632},\n"
     28:                                 "URL = {https://www.worldscientific.com/doi/abs/10.1142/S0218202513500632},\n"
     29:                                 "eprint = {https://www.worldscientific.com/doi/pdf/10.1142/S0218202513500632}\n"
     30:                                 "}\n";
     31: static const char citation2[] = "@article{li2013nonoverlapping,\n"
     32:                                 "title={A nonoverlapping domain decomposition method for incompressible Stokes equations with continuous pressures},\n"
     33:                                 "author={Li, Jing and Tu, Xuemin},\n"
     34:                                 "journal={SIAM Journal on Numerical Analysis},\n"
     35:                                 "volume={51},\n"
     36:                                 "number={2},\n"
     37:                                 "pages={1235--1253},\n"
     38:                                 "year={2013},\n"
     39:                                 "publisher={Society for Industrial and Applied Mathematics}\n"
     40:                                 "}\n";

     42: /*
     43:     This file implements the FETI-DP method in PETSc as part of KSP.
     44: */
     45: typedef struct {
     46:   KSP parentksp;
     47: } KSP_FETIDPMon;

     49: typedef struct {
     50:   KSP              innerksp;        /* the KSP for the Lagrange multipliers */
     51:   PC               innerbddc;       /* the inner BDDC object */
     52:   PetscBool        fully_redundant; /* true for using a fully redundant set of multipliers */
     53:   PetscBool        userbddc;        /* true if the user provided the PCBDDC object */
     54:   PetscBool        saddlepoint;     /* support for saddle point problems */
     55:   IS               pP;              /* index set for pressure variables */
     56:   Vec              rhs_flip;        /* see KSPFETIDPSetUpOperators */
     57:   KSP_FETIDPMon   *monctx;          /* monitor context, used to pass user defined monitors
     58:                                         in the physical space */
     59:   PetscObjectState matstate;        /* these are needed just in the saddle point case */
     60:   PetscObjectState matnnzstate;     /* where we are going to use MatZeroRows on pmat */
     61:   PetscBool        statechanged;
     62:   PetscBool        check;
     63: } KSP_FETIDP;

     65: static PetscErrorCode KSPFETIDPSetPressureOperator_FETIDP(KSP ksp, Mat P)
     66: {
     67:   KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;

     69:   PetscFunctionBegin;
     70:   if (P) fetidp->saddlepoint = PETSC_TRUE;
     71:   PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_PPmat", (PetscObject)P));
     72:   PetscFunctionReturn(PETSC_SUCCESS);
     73: }

     75: /*@
     76:   KSPFETIDPSetPressureOperator - Sets the operator used to set up the pressure preconditioner for the saddle point `KSPFETIDP` solver,

     78:   Collective

     80:   Input Parameters:
     81: + ksp - the `KSPFETIDP` solver
     82: - P   - the linear operator to be preconditioned, usually the mass matrix.

     84:   Level: advanced

     86:   Notes:
     87:   The operator can be either passed in
     88: .vb
     89:   a) monolithic global ordering,
     90:   b) pressure-only global ordering, or
     91:   c) interface pressure ordering (if `-ksp_fetidp_pressure_all false`).
     92: .ve
     93:   In cases b) and c), the pressure ordering of dofs needs to satisfy
     94:   pid_1 < pid_2  iff  gid_1 < gid_2
     95:   where pid_1 and pid_2 are two different pressure dof numbers and gid_1 and gid_2 the corresponding
     96:   id in the monolithic global ordering.

     98: .seealso: [](ch_ksp), `KSPFETIDP`, `MATIS`, `PCBDDC`, `KSPFETIDPGetInnerBDDC()`, `KSPFETIDPGetInnerKSP()`, `KSPSetOperators()`
     99: @*/
    100: PetscErrorCode KSPFETIDPSetPressureOperator(KSP ksp, Mat P)
    101: {
    102:   PetscFunctionBegin;
    105:   PetscTryMethod(ksp, "KSPFETIDPSetPressureOperator_C", (KSP, Mat), (ksp, P));
    106:   PetscFunctionReturn(PETSC_SUCCESS);
    107: }

    109: static PetscErrorCode KSPFETIDPGetInnerKSP_FETIDP(KSP ksp, KSP *innerksp)
    110: {
    111:   KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;

    113:   PetscFunctionBegin;
    114:   *innerksp = fetidp->innerksp;
    115:   PetscFunctionReturn(PETSC_SUCCESS);
    116: }

    118: /*@
    119:   KSPFETIDPGetInnerKSP - Gets the `KSP` object for the Lagrange multipliers from inside a `KSPFETIDP`

    121:   Input Parameter:
    122: . ksp - the `KSPFETIDP`

    124:   Output Parameter:
    125: . innerksp - the `KSP` for the multipliers

    127:   Level: advanced

    129: .seealso: [](ch_ksp), `KSPFETIDP`, `MATIS`, `PCBDDC`, `KSPFETIDPSetInnerBDDC()`, `KSPFETIDPGetInnerBDDC()`
    130: @*/
    131: PetscErrorCode KSPFETIDPGetInnerKSP(KSP ksp, KSP *innerksp)
    132: {
    133:   PetscFunctionBegin;
    135:   PetscAssertPointer(innerksp, 2);
    136:   PetscUseMethod(ksp, "KSPFETIDPGetInnerKSP_C", (KSP, KSP *), (ksp, innerksp));
    137:   PetscFunctionReturn(PETSC_SUCCESS);
    138: }

    140: static PetscErrorCode KSPFETIDPGetInnerBDDC_FETIDP(KSP ksp, PC *pc)
    141: {
    142:   KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;

    144:   PetscFunctionBegin;
    145:   *pc = fetidp->innerbddc;
    146:   PetscFunctionReturn(PETSC_SUCCESS);
    147: }

    149: /*@
    150:   KSPFETIDPGetInnerBDDC - Gets the `PCBDDC` preconditioner used to set up the `KSPFETIDP` matrix for the Lagrange multipliers

    152:   Input Parameter:
    153: . ksp - the `KSPFETIDP` Krylov solver

    155:   Output Parameter:
    156: . pc - the `PCBDDC` preconditioner

    158:   Level: advanced

    160: .seealso: [](ch_ksp), `MATIS`, `PCBDDC`, `KSPFETIDP`, `KSPFETIDPSetInnerBDDC()`, `KSPFETIDPGetInnerKSP()`
    161: @*/
    162: PetscErrorCode KSPFETIDPGetInnerBDDC(KSP ksp, PC *pc)
    163: {
    164:   PetscFunctionBegin;
    166:   PetscAssertPointer(pc, 2);
    167:   PetscUseMethod(ksp, "KSPFETIDPGetInnerBDDC_C", (KSP, PC *), (ksp, pc));
    168:   PetscFunctionReturn(PETSC_SUCCESS);
    169: }

    171: static PetscErrorCode KSPFETIDPSetInnerBDDC_FETIDP(KSP ksp, PC pc)
    172: {
    173:   KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;

    175:   PetscFunctionBegin;
    176:   PetscCall(PetscObjectReference((PetscObject)pc));
    177:   PetscCall(PCDestroy(&fetidp->innerbddc));
    178:   fetidp->innerbddc = pc;
    179:   fetidp->userbddc  = PETSC_TRUE;
    180:   PetscFunctionReturn(PETSC_SUCCESS);
    181: }

    183: /*@
    184:   KSPFETIDPSetInnerBDDC - Provides the `PCBDDC` preconditioner used to set up the `KSPFETIDP` matrix for the Lagrange multipliers

    186:   Collective

    188:   Input Parameters:
    189: + ksp - the `KSPFETIDP` Krylov solver
    190: - pc  - the `PCBDDC` preconditioner

    192:   Level: advanced

    194:   Note:
    195:   A `PC` is automatically created for the `KSPFETIDP` and can be accessed to change options with `KSPFETIDPGetInnerBDDC()` hence this routine is rarely needed

    197: .seealso: [](ch_ksp), `MATIS`, `PCBDDC`, `KSPFETIDPGetInnerBDDC()`, `KSPFETIDPGetInnerKSP()`
    198: @*/
    199: PetscErrorCode KSPFETIDPSetInnerBDDC(KSP ksp, PC pc)
    200: {
    201:   PetscBool isbddc;

    203:   PetscFunctionBegin;
    206:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &isbddc));
    207:   PetscCheck(isbddc, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_WRONG, "KSPFETIDPSetInnerBDDC need a PCBDDC preconditioner");
    208:   PetscTryMethod(ksp, "KSPFETIDPSetInnerBDDC_C", (KSP, PC), (ksp, pc));
    209:   PetscFunctionReturn(PETSC_SUCCESS);
    210: }

    212: static PetscErrorCode KSPBuildSolution_FETIDP(KSP ksp, Vec v, Vec *V)
    213: {
    214:   KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
    215:   Mat         F;
    216:   Vec         Xl;

    218:   PetscFunctionBegin;
    219:   PetscCall(KSPGetOperators(fetidp->innerksp, &F, NULL));
    220:   PetscCall(KSPBuildSolution(fetidp->innerksp, NULL, &Xl));
    221:   if (v) {
    222:     PetscCall(PCBDDCMatFETIDPGetSolution(F, Xl, v));
    223:     *V = v;
    224:   } else {
    225:     PetscCall(PCBDDCMatFETIDPGetSolution(F, Xl, *V));
    226:   }
    227:   PetscFunctionReturn(PETSC_SUCCESS);
    228: }

    230: static PetscErrorCode KSPMonitor_FETIDP(KSP ksp, PetscInt it, PetscReal rnorm, PetscCtx ctx)
    231: {
    232:   KSP_FETIDPMon *monctx = (KSP_FETIDPMon *)ctx;

    234:   PetscFunctionBegin;
    235:   PetscCall(KSPMonitor(monctx->parentksp, it, rnorm));
    236:   PetscFunctionReturn(PETSC_SUCCESS);
    237: }

    239: static PetscErrorCode KSPComputeEigenvalues_FETIDP(KSP ksp, PetscInt nmax, PetscReal *r, PetscReal *c, PetscInt *neig)
    240: {
    241:   KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;

    243:   PetscFunctionBegin;
    244:   PetscCall(KSPComputeEigenvalues(fetidp->innerksp, nmax, r, c, neig));
    245:   PetscFunctionReturn(PETSC_SUCCESS);
    246: }

    248: static PetscErrorCode KSPComputeExtremeSingularValues_FETIDP(KSP ksp, PetscReal *emax, PetscReal *emin)
    249: {
    250:   KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;

    252:   PetscFunctionBegin;
    253:   PetscCall(KSPComputeExtremeSingularValues(fetidp->innerksp, emax, emin));
    254:   PetscFunctionReturn(PETSC_SUCCESS);
    255: }

    257: static PetscErrorCode KSPFETIDPCheckOperators(KSP ksp, PetscViewer viewer)
    258: {
    259:   KSP_FETIDP     *fetidp = (KSP_FETIDP *)ksp->data;
    260:   PC_BDDC        *pcbddc = (PC_BDDC *)fetidp->innerbddc->data;
    261:   PC_IS          *pcis   = (PC_IS *)fetidp->innerbddc->data;
    262:   Mat_IS         *matis  = (Mat_IS *)fetidp->innerbddc->pmat->data;
    263:   Mat             F;
    264:   FETIDPMat_ctx   fetidpmat_ctx;
    265:   Vec             test_vec, test_vec_p = NULL, fetidp_global;
    266:   IS              dirdofs, isvert;
    267:   MPI_Comm        comm = PetscObjectComm((PetscObject)ksp);
    268:   PetscScalar     sval, *array;
    269:   PetscReal       val, rval;
    270:   const PetscInt *vertex_indices;
    271:   PetscInt        i, n_vertices;
    272:   PetscBool       isascii;

    274:   PetscFunctionBegin;
    275:   PetscCheckSameComm(ksp, 1, viewer, 2);
    276:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    277:   PetscCheck(isascii, comm, PETSC_ERR_SUP, "Unsupported viewer");
    278:   PetscCall(PetscViewerASCIIPrintf(viewer, "----------FETI-DP MAT  --------------\n"));
    279:   PetscCall(PetscViewerASCIIAddTab(viewer, 2));
    280:   PetscCall(KSPGetOperators(fetidp->innerksp, &F, NULL));
    281:   PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
    282:   PetscCall(MatView(F, viewer));
    283:   PetscCall(PetscViewerPopFormat(viewer));
    284:   PetscCall(PetscViewerASCIISubtractTab(viewer, 2));
    285:   PetscCall(MatShellGetContext(F, &fetidpmat_ctx));
    286:   PetscCall(PetscViewerASCIIPrintf(viewer, "----------FETI-DP TESTS--------------\n"));
    287:   PetscCall(PetscViewerASCIIPrintf(viewer, "All tests should return zero!\n"));
    288:   PetscCall(PetscViewerASCIIPrintf(viewer, "FETIDP MAT context in the "));
    289:   if (fetidp->fully_redundant) {
    290:     PetscCall(PetscViewerASCIIPrintf(viewer, "fully redundant case for lagrange multipliers.\n"));
    291:   } else {
    292:     PetscCall(PetscViewerASCIIPrintf(viewer, "Non-fully redundant case for lagrange multiplier.\n"));
    293:   }
    294:   PetscCall(PetscViewerFlush(viewer));

    296:   /* Get Vertices used to define the BDDC */
    297:   PetscCall(PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph, NULL, NULL, NULL, NULL, &isvert));
    298:   PetscCall(ISGetLocalSize(isvert, &n_vertices));
    299:   PetscCall(ISGetIndices(isvert, &vertex_indices));

    301:   /******************************************************************/
    302:   /* TEST A/B: Test numbering of global fetidp dofs                 */
    303:   /******************************************************************/
    304:   PetscCall(MatCreateVecs(F, &fetidp_global, NULL));
    305:   PetscCall(VecDuplicate(fetidpmat_ctx->lambda_local, &test_vec));
    306:   PetscCall(VecSet(fetidp_global, 1.0));
    307:   PetscCall(VecSet(test_vec, 1.));
    308:   PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidp_global, fetidpmat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
    309:   PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidp_global, fetidpmat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
    310:   if (fetidpmat_ctx->l2g_p) {
    311:     PetscCall(VecDuplicate(fetidpmat_ctx->vP, &test_vec_p));
    312:     PetscCall(VecSet(test_vec_p, 1.));
    313:     PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_p, fetidp_global, fetidpmat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
    314:     PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_p, fetidp_global, fetidpmat_ctx->vP, INSERT_VALUES, SCATTER_REVERSE));
    315:   }
    316:   PetscCall(VecAXPY(test_vec, -1.0, fetidpmat_ctx->lambda_local));
    317:   PetscCall(VecNorm(test_vec, NORM_INFINITY, &val));
    318:   PetscCall(VecDestroy(&test_vec));
    319:   PetscCallMPI(MPI_Reduce(&val, &rval, 1, MPIU_REAL, MPIU_MAX, 0, comm));
    320:   PetscCall(PetscViewerASCIIPrintf(viewer, "A: CHECK glob to loc: % 1.14e\n", (double)rval));

    322:   if (fetidpmat_ctx->l2g_p) {
    323:     PetscCall(VecAXPY(test_vec_p, -1.0, fetidpmat_ctx->vP));
    324:     PetscCall(VecNorm(test_vec_p, NORM_INFINITY, &val));
    325:     PetscCallMPI(MPI_Reduce(&val, &rval, 1, MPIU_REAL, MPIU_MAX, 0, comm));
    326:     PetscCall(PetscViewerASCIIPrintf(viewer, "A: CHECK glob to loc (p): % 1.14e\n", (double)rval));
    327:   }

    329:   if (fetidp->fully_redundant) {
    330:     PetscCall(VecSet(fetidp_global, 0.0));
    331:     PetscCall(VecSet(fetidpmat_ctx->lambda_local, 0.5));
    332:     PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
    333:     PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
    334:     PetscCall(VecSum(fetidp_global, &sval));
    335:     val = PetscRealPart(sval) - fetidpmat_ctx->n_lambda;
    336:     PetscCallMPI(MPI_Reduce(&val, &rval, 1, MPIU_REAL, MPIU_MAX, 0, comm));
    337:     PetscCall(PetscViewerASCIIPrintf(viewer, "B: CHECK loc to glob: % 1.14e\n", (double)rval));
    338:   }

    340:   if (fetidpmat_ctx->l2g_p) {
    341:     PetscCall(VecSet(pcis->vec1_N, 1.0));
    342:     PetscCall(VecSet(pcis->vec1_global, 0.0));
    343:     PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_N, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
    344:     PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_N, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));

    346:     PetscCall(VecSet(fetidp_global, 0.0));
    347:     PetscCall(VecSet(fetidpmat_ctx->vP, -1.0));
    348:     PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_p, fetidpmat_ctx->vP, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
    349:     PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_p, fetidpmat_ctx->vP, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
    350:     PetscCall(VecScatterBegin(fetidpmat_ctx->g2g_p, fetidp_global, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
    351:     PetscCall(VecScatterEnd(fetidpmat_ctx->g2g_p, fetidp_global, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
    352:     PetscCall(VecScatterBegin(fetidpmat_ctx->g2g_p, pcis->vec1_global, fetidp_global, INSERT_VALUES, SCATTER_FORWARD));
    353:     PetscCall(VecScatterEnd(fetidpmat_ctx->g2g_p, pcis->vec1_global, fetidp_global, INSERT_VALUES, SCATTER_FORWARD));
    354:     PetscCall(VecSum(fetidp_global, &sval));
    355:     val = PetscRealPart(sval);
    356:     PetscCallMPI(MPI_Reduce(&val, &rval, 1, MPIU_REAL, MPIU_MAX, 0, comm));
    357:     PetscCall(PetscViewerASCIIPrintf(viewer, "B: CHECK loc to glob (p): % 1.14e\n", (double)rval));
    358:   }

    360:   /******************************************************************/
    361:   /* TEST C: It should hold B_delta*w=0, w\in\widehat{W}            */
    362:   /* This is the meaning of the B matrix                            */
    363:   /******************************************************************/

    365:   PetscCall(VecSetRandom(pcis->vec1_N, NULL));
    366:   PetscCall(VecSet(pcis->vec1_global, 0.0));
    367:   PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_N, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
    368:   PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_N, pcis->vec1_global, ADD_VALUES, SCATTER_REVERSE));
    369:   PetscCall(VecScatterBegin(matis->rctx, pcis->vec1_global, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD));
    370:   PetscCall(VecScatterEnd(matis->rctx, pcis->vec1_global, pcis->vec1_N, INSERT_VALUES, SCATTER_FORWARD));
    371:   PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
    372:   PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
    373:   /* Action of B_delta */
    374:   PetscCall(MatMult(fetidpmat_ctx->B_delta, pcis->vec1_B, fetidpmat_ctx->lambda_local));
    375:   PetscCall(VecSet(fetidp_global, 0.0));
    376:   PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
    377:   PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
    378:   PetscCall(VecNorm(fetidp_global, NORM_INFINITY, &val));
    379:   PetscCall(PetscViewerASCIIPrintf(viewer, "C: CHECK infty norm of B_delta*w (w continuous): % 1.14e\n", (double)val));

    381:   /******************************************************************/
    382:   /* TEST D: It should hold E_Dw = w - P_Dw w\in\widetilde{W}       */
    383:   /* E_D = R_D^TR                                                   */
    384:   /* P_D = B_{D,delta}^T B_{delta}                                  */
    385:   /* eq.44 Mandel Tezaur and Dohrmann 2005                          */
    386:   /******************************************************************/

    388:   /* compute a random vector in \widetilde{W} */
    389:   PetscCall(VecSetRandom(pcis->vec1_N, NULL));
    390:   /* set zero at vertices and essential dofs */
    391:   PetscCall(VecGetArray(pcis->vec1_N, &array));
    392:   for (i = 0; i < n_vertices; i++) array[vertex_indices[i]] = 0.0;
    393:   PetscCall(PCBDDCGraphGetDirichletDofs(pcbddc->mat_graph, &dirdofs));
    394:   if (dirdofs) {
    395:     const PetscInt *idxs;
    396:     PetscInt        ndir;

    398:     PetscCall(ISGetLocalSize(dirdofs, &ndir));
    399:     PetscCall(ISGetIndices(dirdofs, &idxs));
    400:     for (i = 0; i < ndir; i++) array[idxs[i]] = 0.0;
    401:     PetscCall(ISRestoreIndices(dirdofs, &idxs));
    402:   }
    403:   PetscCall(VecRestoreArray(pcis->vec1_N, &array));
    404:   /* store w for final comparison */
    405:   PetscCall(VecDuplicate(pcis->vec1_B, &test_vec));
    406:   PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_N, test_vec, INSERT_VALUES, SCATTER_FORWARD));
    407:   PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_N, test_vec, INSERT_VALUES, SCATTER_FORWARD));

    409:   /* Jump operator P_D : results stored in pcis->vec1_B */
    410:   /* Action of B_delta */
    411:   PetscCall(MatMult(fetidpmat_ctx->B_delta, test_vec, fetidpmat_ctx->lambda_local));
    412:   PetscCall(VecSet(fetidp_global, 0.0));
    413:   PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
    414:   PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
    415:   /* Action of B_Ddelta^T */
    416:   PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidp_global, fetidpmat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
    417:   PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidp_global, fetidpmat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
    418:   PetscCall(MatMultTranspose(fetidpmat_ctx->B_Ddelta, fetidpmat_ctx->lambda_local, pcis->vec1_B));

    420:   /* Average operator E_D : results stored in pcis->vec2_B */
    421:   PetscCall(PCBDDCScalingExtension(fetidpmat_ctx->pc, test_vec, pcis->vec1_global));
    422:   PetscCall(VecScatterBegin(pcis->global_to_B, pcis->vec1_global, pcis->vec2_B, INSERT_VALUES, SCATTER_FORWARD));
    423:   PetscCall(VecScatterEnd(pcis->global_to_B, pcis->vec1_global, pcis->vec2_B, INSERT_VALUES, SCATTER_FORWARD));

    425:   /* test E_D=I-P_D */
    426:   PetscCall(VecAXPY(pcis->vec1_B, 1.0, pcis->vec2_B));
    427:   PetscCall(VecAXPY(pcis->vec1_B, -1.0, test_vec));
    428:   PetscCall(VecNorm(pcis->vec1_B, NORM_INFINITY, &val));
    429:   PetscCall(VecDestroy(&test_vec));
    430:   PetscCallMPI(MPI_Reduce(&val, &rval, 1, MPIU_REAL, MPIU_MAX, 0, comm));
    431:   PetscCall(PetscViewerASCIIPrintf(viewer, "%d: CHECK infty norm of E_D + P_D - I: %1.14e\n", PetscGlobalRank, (double)val));

    433:   /******************************************************************/
    434:   /* TEST E: It should hold R_D^TP_Dw=0 w\in\widetilde{W}           */
    435:   /* eq.48 Mandel Tezaur and Dohrmann 2005                          */
    436:   /******************************************************************/

    438:   PetscCall(VecSetRandom(pcis->vec1_N, NULL));
    439:   /* set zero at vertices and essential dofs */
    440:   PetscCall(VecGetArray(pcis->vec1_N, &array));
    441:   for (i = 0; i < n_vertices; i++) array[vertex_indices[i]] = 0.0;
    442:   if (dirdofs) {
    443:     const PetscInt *idxs;
    444:     PetscInt        ndir;

    446:     PetscCall(ISGetLocalSize(dirdofs, &ndir));
    447:     PetscCall(ISGetIndices(dirdofs, &idxs));
    448:     for (i = 0; i < ndir; i++) array[idxs[i]] = 0.0;
    449:     PetscCall(ISRestoreIndices(dirdofs, &idxs));
    450:   }
    451:   PetscCall(VecRestoreArray(pcis->vec1_N, &array));

    453:   /* Jump operator P_D : results stored in pcis->vec1_B */

    455:   PetscCall(VecScatterBegin(pcis->N_to_B, pcis->vec1_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
    456:   PetscCall(VecScatterEnd(pcis->N_to_B, pcis->vec1_N, pcis->vec1_B, INSERT_VALUES, SCATTER_FORWARD));
    457:   /* Action of B_delta */
    458:   PetscCall(MatMult(fetidpmat_ctx->B_delta, pcis->vec1_B, fetidpmat_ctx->lambda_local));
    459:   PetscCall(VecSet(fetidp_global, 0.0));
    460:   PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
    461:   PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, fetidp_global, ADD_VALUES, SCATTER_FORWARD));
    462:   /* Action of B_Ddelta^T */
    463:   PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidp_global, fetidpmat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
    464:   PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidp_global, fetidpmat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
    465:   PetscCall(MatMultTranspose(fetidpmat_ctx->B_Ddelta, fetidpmat_ctx->lambda_local, pcis->vec1_B));
    466:   /* scaling */
    467:   PetscCall(PCBDDCScalingExtension(fetidpmat_ctx->pc, pcis->vec1_B, pcis->vec1_global));
    468:   PetscCall(VecNorm(pcis->vec1_global, NORM_INFINITY, &val));
    469:   PetscCall(PetscViewerASCIIPrintf(viewer, "E: CHECK infty norm of R^T_D P_D: % 1.14e\n", (double)val));

    471:   if (!fetidp->fully_redundant) {
    472:     /******************************************************************/
    473:     /* TEST F: It should holds B_{delta}B^T_{D,delta}=I               */
    474:     /* Corollary thm 14 Mandel Tezaur and Dohrmann 2005               */
    475:     /******************************************************************/
    476:     PetscCall(VecDuplicate(fetidp_global, &test_vec));
    477:     PetscCall(VecSetRandom(fetidp_global, NULL));
    478:     if (fetidpmat_ctx->l2g_p) {
    479:       PetscCall(VecSet(fetidpmat_ctx->vP, 0.));
    480:       PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_p, fetidpmat_ctx->vP, fetidp_global, INSERT_VALUES, SCATTER_FORWARD));
    481:       PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_p, fetidpmat_ctx->vP, fetidp_global, INSERT_VALUES, SCATTER_FORWARD));
    482:     }
    483:     /* Action of B_Ddelta^T */
    484:     PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidp_global, fetidpmat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
    485:     PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidp_global, fetidpmat_ctx->lambda_local, INSERT_VALUES, SCATTER_REVERSE));
    486:     PetscCall(MatMultTranspose(fetidpmat_ctx->B_Ddelta, fetidpmat_ctx->lambda_local, pcis->vec1_B));
    487:     /* Action of B_delta */
    488:     PetscCall(MatMult(fetidpmat_ctx->B_delta, pcis->vec1_B, fetidpmat_ctx->lambda_local));
    489:     PetscCall(VecSet(test_vec, 0.0));
    490:     PetscCall(VecScatterBegin(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, test_vec, ADD_VALUES, SCATTER_FORWARD));
    491:     PetscCall(VecScatterEnd(fetidpmat_ctx->l2g_lambda, fetidpmat_ctx->lambda_local, test_vec, ADD_VALUES, SCATTER_FORWARD));
    492:     PetscCall(VecAXPY(fetidp_global, -1., test_vec));
    493:     PetscCall(VecNorm(fetidp_global, NORM_INFINITY, &val));
    494:     PetscCall(PetscViewerASCIIPrintf(viewer, "E: CHECK infty norm of P^T_D - I: % 1.14e\n", (double)val));
    495:     PetscCall(VecDestroy(&test_vec));
    496:   }
    497:   PetscCall(PetscViewerASCIIPrintf(viewer, "-------------------------------------\n"));
    498:   PetscCall(PetscViewerFlush(viewer));
    499:   PetscCall(VecDestroy(&test_vec_p));
    500:   PetscCall(ISDestroy(&dirdofs));
    501:   PetscCall(VecDestroy(&fetidp_global));
    502:   PetscCall(ISRestoreIndices(isvert, &vertex_indices));
    503:   PetscCall(PCBDDCGraphRestoreCandidatesIS(pcbddc->mat_graph, NULL, NULL, NULL, NULL, &isvert));
    504:   PetscFunctionReturn(PETSC_SUCCESS);
    505: }

    507: static PetscErrorCode KSPFETIDPSetUpOperators(KSP ksp)
    508: {
    509:   KSP_FETIDP      *fetidp = (KSP_FETIDP *)ksp->data;
    510:   PC_BDDC         *pcbddc = (PC_BDDC *)fetidp->innerbddc->data;
    511:   Mat              A, Ap;
    512:   PetscInt         fidp[8] = {-1}, nfp = 8;
    513:   PetscMPIInt      size;
    514:   PetscBool        ismatis, pisz = PETSC_FALSE, allp = PETSC_FALSE, schp = PETSC_FALSE;
    515:   PetscBool        flip = PETSC_FALSE; /* Usually, Stokes is written (B = -\int_\Omega \nabla \cdot u q)
    516:                            | A B'| | v | = | f |
    517:                            | B 0 | | p | = | g |
    518:                             If -ksp_fetidp_saddlepoint_flip is true, the code assumes it is written as
    519:                            | A B'| | v | = | f |
    520:                            |-B 0 | | p | = |-g |
    521:                          */
    522:   PetscObjectState matstate, matnnzstate;

    524:   PetscFunctionBegin;
    525:   PetscOptionsBegin(PetscObjectComm((PetscObject)ksp), ((PetscObject)ksp)->prefix, "FETI-DP options", "PC");
    526:   PetscCall(PetscOptionsIntArray("-ksp_fetidp_pressure_field", "Field id for pressures for saddle-point problems", NULL, fidp, &nfp, NULL));
    527:   PetscCall(PetscOptionsBool("-ksp_fetidp_pressure_all", "Use the whole pressure set instead of just that at the interface", NULL, allp, &allp, NULL));
    528:   PetscCall(PetscOptionsBool("-ksp_fetidp_saddlepoint_flip", "Flip the sign of the pressure-velocity (lower-left) block", NULL, flip, &flip, NULL));
    529:   PetscCall(PetscOptionsBool("-ksp_fetidp_pressure_schur", "Use a BDDC solver for pressure", NULL, schp, &schp, NULL));
    530:   PetscOptionsEnd();

    532:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ksp), &size));
    533:   fetidp->saddlepoint = (nfp > 0 ? PETSC_TRUE : fetidp->saddlepoint);
    534:   if (size == 1) fetidp->saddlepoint = PETSC_FALSE;

    536:   PetscCall(KSPGetOperators(ksp, &A, &Ap));
    537:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATIS, &ismatis));
    538:   PetscCheck(ismatis, PetscObjectComm((PetscObject)ksp), PETSC_ERR_USER, "Amat should be of type MATIS");

    540:   /* Quiet return if the matrix states are unchanged.
    541:      Needed only for the saddle point case since it uses MatZeroRows
    542:      on a matrix that may not have changed */
    543:   PetscCall(PetscObjectStateGet((PetscObject)A, &matstate));
    544:   PetscCall(MatGetNonzeroState(A, &matnnzstate));
    545:   if (matstate == fetidp->matstate && matnnzstate == fetidp->matnnzstate) PetscFunctionReturn(PETSC_SUCCESS);
    546:   fetidp->matstate     = matstate;
    547:   fetidp->matnnzstate  = matnnzstate;
    548:   fetidp->statechanged = fetidp->saddlepoint;

    550:   /* see if we have some fields attached */
    551:   if (!pcbddc->n_ISForDofsLocal && !pcbddc->n_ISForDofs) {
    552:     DM             dm;
    553:     PetscContainer c;

    555:     PetscCall(KSPGetDM(ksp, &dm));
    556:     PetscCall(PetscObjectQuery((PetscObject)A, "_convert_nest_lfields", (PetscObject *)&c));
    557:     if (dm) {
    558:       IS      *fields;
    559:       PetscInt nf, i;

    561:       PetscCall(DMCreateFieldDecomposition(dm, &nf, NULL, &fields, NULL));
    562:       PetscCall(PCBDDCSetDofsSplitting(fetidp->innerbddc, nf, fields));
    563:       for (i = 0; i < nf; i++) PetscCall(ISDestroy(&fields[i]));
    564:       PetscCall(PetscFree(fields));
    565:     } else if (c) {
    566:       MatISLocalFields lf;

    568:       PetscCall(PetscContainerGetPointer(c, &lf));
    569:       PetscCall(PCBDDCSetDofsSplittingLocal(fetidp->innerbddc, lf->nr, lf->rf));
    570:     }
    571:   }

    573:   if (!fetidp->saddlepoint) {
    574:     PetscCall(PCSetOperators(fetidp->innerbddc, A, A));
    575:   } else {
    576:     Mat          nA, lA, PPmat;
    577:     MatNullSpace nnsp;
    578:     IS           pP;
    579:     PetscInt     totP;

    581:     PetscCall(MatISGetLocalMat(A, &lA));
    582:     PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_lA", (PetscObject)lA));

    584:     pP = fetidp->pP;
    585:     if (!pP) { /* first time, need to compute pressure dofs */
    586:       PC_IS                 *pcis  = (PC_IS *)fetidp->innerbddc->data;
    587:       Mat_IS                *matis = (Mat_IS *)A->data;
    588:       ISLocalToGlobalMapping l2g;
    589:       IS                     lP = NULL, II, pII, lPall, Pall, is1, is2;
    590:       const PetscInt        *idxs;
    591:       PetscInt               nl, ni, *widxs;
    592:       PetscInt               i, j, n_neigh, *neigh, *n_shared, **shared, *count;
    593:       PetscInt               rst, ren, n;
    594:       PetscBool              ploc;

    596:       PetscCall(MatGetLocalSize(A, &nl, NULL));
    597:       PetscCall(MatGetOwnershipRange(A, &rst, &ren));
    598:       PetscCall(MatGetLocalSize(lA, &n, NULL));
    599:       PetscCall(MatISGetLocalToGlobalMapping(A, &l2g, NULL));

    601:       if (!pcis->is_I_local) { /* need to compute interior dofs */
    602:         PetscCall(PetscCalloc1(n, &count));
    603:         PetscCall(ISLocalToGlobalMappingGetInfo(l2g, &n_neigh, &neigh, &n_shared, &shared));
    604:         for (i = 1; i < n_neigh; i++)
    605:           for (j = 0; j < n_shared[i]; j++) count[shared[i][j]] += 1;
    606:         for (i = 0, j = 0; i < n; i++)
    607:           if (!count[i]) count[j++] = i;
    608:         PetscCall(ISLocalToGlobalMappingRestoreInfo(l2g, &n_neigh, &neigh, &n_shared, &shared));
    609:         PetscCall(ISCreateGeneral(PETSC_COMM_SELF, j, count, PETSC_OWN_POINTER, &II));
    610:       } else {
    611:         PetscCall(PetscObjectReference((PetscObject)pcis->is_I_local));
    612:         II = pcis->is_I_local;
    613:       }

    615:       /* interior dofs in layout */
    616:       PetscCall(PetscArrayzero(matis->sf_leafdata, n));
    617:       PetscCall(PetscArrayzero(matis->sf_rootdata, nl));
    618:       PetscCall(ISGetLocalSize(II, &ni));
    619:       PetscCall(ISGetIndices(II, &idxs));
    620:       for (i = 0; i < ni; i++) matis->sf_leafdata[idxs[i]] = 1;
    621:       PetscCall(ISRestoreIndices(II, &idxs));
    622:       PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, MPI_REPLACE));
    623:       PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, MPI_REPLACE));
    624:       PetscCall(PetscMalloc1(PetscMax(nl, n), &widxs));
    625:       for (i = 0, ni = 0; i < nl; i++)
    626:         if (matis->sf_rootdata[i]) widxs[ni++] = i + rst;
    627:       PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)ksp), ni, widxs, PETSC_COPY_VALUES, &pII));

    629:       /* pressure dofs */
    630:       Pall  = NULL;
    631:       lPall = NULL;
    632:       ploc  = PETSC_FALSE;
    633:       if (nfp == 0) { /* zero pressure block */
    634:         PetscInt np;

    636:         PetscCall(MatFindZeroDiagonals(A, &Pall));
    637:         PetscCall(ISGetSize(Pall, &np));
    638:         if (!np) { /* zero-block not found, defaults to last field (if set) */
    639:           nfp     = 1;
    640:           fidp[0] = pcbddc->n_ISForDofsLocal ? pcbddc->n_ISForDofsLocal - 1 : pcbddc->n_ISForDofs - 1;
    641:           PetscCall(ISDestroy(&Pall));
    642:         } else if (!pcbddc->n_ISForDofsLocal && !pcbddc->n_ISForDofs) {
    643:           PetscCall(PCBDDCSetDofsSplitting(fetidp->innerbddc, 1, &Pall));
    644:         }
    645:       }
    646:       if (!Pall) { /* look for registered fields when no zero block has been found */
    647:         IS *tis;

    649:         PetscCall(PetscMalloc1(nfp, &tis));
    650:         if (pcbddc->n_ISForDofsLocal) {
    651:           for (PetscInt i = 0; i < nfp; i++) {
    652:             PetscInt fid = fidp[i];

    654:             PetscCheck(fid >= 0 && fid < pcbddc->n_ISForDofsLocal, PetscObjectComm((PetscObject)ksp), PETSC_ERR_USER, "Invalid field id for pressure %" PetscInt_FMT ", max %" PetscInt_FMT, fid, pcbddc->n_ISForDofsLocal);
    655:             /* need a sequential IS */
    656:             PetscCall(ISOnComm(pcbddc->ISForDofsLocal[fid], PETSC_COMM_SELF, PETSC_COPY_VALUES, &tis[i]));
    657:           }
    658:           PetscCall(ISConcatenate(PETSC_COMM_SELF, nfp, tis, &lPall));
    659:           ploc = PETSC_TRUE;
    660:         } else if (pcbddc->n_ISForDofs) {
    661:           for (PetscInt i = 0; i < nfp; i++) {
    662:             PetscInt fid = fidp[i];

    664:             PetscCheck(fid >= 0 && fid < pcbddc->n_ISForDofs, PetscObjectComm((PetscObject)ksp), PETSC_ERR_USER, "Invalid field id for pressure %" PetscInt_FMT ", max %" PetscInt_FMT, fid, pcbddc->n_ISForDofs);
    665:             PetscCall(PetscObjectReference((PetscObject)pcbddc->ISForDofs[fid]));
    666:             tis[i] = pcbddc->ISForDofs[fid];
    667:           }
    668:           PetscCall(ISConcatenate(PetscObjectComm((PetscObject)ksp), nfp, tis, &Pall));
    669:           PetscCall(ISSort(Pall));
    670:         } else SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_USER, "Cannot detect pressure field! Use KSPFETIDPGetInnerBDDC() + PCBDDCSetDofsSplitting or PCBDDCSetDofsSplittingLocal");
    671:         for (PetscInt i = 0; i < nfp; i++) PetscCall(ISDestroy(&tis[i]));
    672:         PetscCall(PetscFree(tis));
    673:       }

    675:       /* if the user requested the entire pressure,
    676:          remove the interior pressure dofs from II (or pII) */
    677:       if (allp) {
    678:         if (ploc) {
    679:           IS nII;
    680:           PetscCall(ISDifference(II, lPall, &nII));
    681:           PetscCall(ISDestroy(&II));
    682:           II = nII;
    683:         } else {
    684:           IS nII;
    685:           PetscCall(ISDifference(pII, Pall, &nII));
    686:           PetscCall(ISDestroy(&pII));
    687:           pII = nII;
    688:         }
    689:       }
    690:       if (ploc) {
    691:         PetscCall(ISDifference(lPall, II, &lP));
    692:         PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_lP", (PetscObject)lP));
    693:       } else {
    694:         PetscCall(ISDifference(Pall, pII, &pP));
    695:         PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_pP", (PetscObject)pP));
    696:         /* need all local pressure dofs */
    697:         PetscCall(PetscArrayzero(matis->sf_leafdata, n));
    698:         PetscCall(PetscArrayzero(matis->sf_rootdata, nl));
    699:         PetscCall(ISGetLocalSize(Pall, &ni));
    700:         PetscCall(ISGetIndices(Pall, &idxs));
    701:         for (i = 0; i < ni; i++) matis->sf_rootdata[idxs[i] - rst] = 1;
    702:         PetscCall(ISRestoreIndices(Pall, &idxs));
    703:         PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
    704:         PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
    705:         for (i = 0, ni = 0; i < n; i++)
    706:           if (matis->sf_leafdata[i]) widxs[ni++] = i;
    707:         PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ni, widxs, PETSC_COPY_VALUES, &lPall));
    708:       }

    710:       if (!Pall) {
    711:         PetscCall(PetscArrayzero(matis->sf_leafdata, n));
    712:         PetscCall(PetscArrayzero(matis->sf_rootdata, nl));
    713:         PetscCall(ISGetLocalSize(lPall, &ni));
    714:         PetscCall(ISGetIndices(lPall, &idxs));
    715:         for (i = 0; i < ni; i++) matis->sf_leafdata[idxs[i]] = 1;
    716:         PetscCall(ISRestoreIndices(lPall, &idxs));
    717:         PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, MPI_REPLACE));
    718:         PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, MPI_REPLACE));
    719:         for (i = 0, ni = 0; i < nl; i++)
    720:           if (matis->sf_rootdata[i]) widxs[ni++] = i + rst;
    721:         PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)ksp), ni, widxs, PETSC_COPY_VALUES, &Pall));
    722:       }
    723:       PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_aP", (PetscObject)Pall));

    725:       if (flip) {
    726:         PetscInt npl;
    727:         PetscCall(ISGetLocalSize(Pall, &npl));
    728:         PetscCall(ISGetIndices(Pall, &idxs));
    729:         PetscCall(MatCreateVecs(A, NULL, &fetidp->rhs_flip));
    730:         PetscCall(VecSet(fetidp->rhs_flip, 1.));
    731:         PetscCall(VecSetOption(fetidp->rhs_flip, VEC_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE));
    732:         for (i = 0; i < npl; i++) PetscCall(VecSetValue(fetidp->rhs_flip, idxs[i], -1., INSERT_VALUES));
    733:         PetscCall(VecAssemblyBegin(fetidp->rhs_flip));
    734:         PetscCall(VecAssemblyEnd(fetidp->rhs_flip));
    735:         PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_flip", (PetscObject)fetidp->rhs_flip));
    736:         PetscCall(ISRestoreIndices(Pall, &idxs));
    737:       }
    738:       PetscCall(ISDestroy(&Pall));
    739:       PetscCall(ISDestroy(&pII));

    741:       /* local selected pressures in subdomain-wise and global ordering */
    742:       PetscCall(PetscArrayzero(matis->sf_leafdata, n));
    743:       PetscCall(PetscArrayzero(matis->sf_rootdata, nl));
    744:       if (!ploc) {
    745:         PetscInt *widxs2;

    747:         PetscCheck(pP, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "Missing parallel pressure IS");
    748:         PetscCall(ISGetLocalSize(pP, &ni));
    749:         PetscCall(ISGetIndices(pP, &idxs));
    750:         for (i = 0; i < ni; i++) matis->sf_rootdata[idxs[i] - rst] = 1;
    751:         PetscCall(ISRestoreIndices(pP, &idxs));
    752:         PetscCall(PetscSFBcastBegin(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
    753:         PetscCall(PetscSFBcastEnd(matis->sf, MPIU_INT, matis->sf_rootdata, matis->sf_leafdata, MPI_REPLACE));
    754:         for (i = 0, ni = 0; i < n; i++)
    755:           if (matis->sf_leafdata[i]) widxs[ni++] = i;
    756:         PetscCall(PetscMalloc1(ni, &widxs2));
    757:         PetscCall(ISLocalToGlobalMappingApply(l2g, ni, widxs, widxs2));
    758:         PetscCall(ISCreateGeneral(PETSC_COMM_SELF, ni, widxs, PETSC_COPY_VALUES, &lP));
    759:         PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_lP", (PetscObject)lP));
    760:         PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)ksp), ni, widxs2, PETSC_OWN_POINTER, &is1));
    761:         PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_gP", (PetscObject)is1));
    762:         PetscCall(ISDestroy(&is1));
    763:       } else {
    764:         PetscCheck(lP, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Missing sequential pressure IS");
    765:         PetscCall(ISGetLocalSize(lP, &ni));
    766:         PetscCall(ISGetIndices(lP, &idxs));
    767:         for (i = 0; i < ni; i++)
    768:           if (idxs[i] >= 0 && idxs[i] < n) matis->sf_leafdata[idxs[i]] = 1;
    769:         PetscCall(ISRestoreIndices(lP, &idxs));
    770:         PetscCall(PetscSFReduceBegin(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, MPI_REPLACE));
    771:         PetscCall(ISLocalToGlobalMappingApply(l2g, ni, idxs, widxs));
    772:         PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)ksp), ni, widxs, PETSC_COPY_VALUES, &is1));
    773:         PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_gP", (PetscObject)is1));
    774:         PetscCall(ISDestroy(&is1));
    775:         PetscCall(PetscSFReduceEnd(matis->sf, MPIU_INT, matis->sf_leafdata, matis->sf_rootdata, MPI_REPLACE));
    776:         for (i = 0, ni = 0; i < nl; i++)
    777:           if (matis->sf_rootdata[i]) widxs[ni++] = i + rst;
    778:         PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)ksp), ni, widxs, PETSC_COPY_VALUES, &pP));
    779:         PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_pP", (PetscObject)pP));
    780:       }
    781:       PetscCall(PetscFree(widxs));

    783:       /* If there's any "interior pressure",
    784:          we may want to use a discrete harmonic solver instead
    785:          of a Stokes harmonic for the Dirichlet preconditioner
    786:          Need to extract the interior velocity dofs in interior dofs ordering (iV)
    787:          and interior pressure dofs in local ordering (iP) */
    788:       if (!allp) {
    789:         ISLocalToGlobalMapping l2g_t;

    791:         PetscCall(ISDifference(lPall, lP, &is1));
    792:         PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_iP", (PetscObject)is1));
    793:         PetscCall(ISDifference(II, is1, &is2));
    794:         PetscCall(ISDestroy(&is1));
    795:         PetscCall(ISLocalToGlobalMappingCreateIS(II, &l2g_t));
    796:         PetscCall(ISGlobalToLocalMappingApplyIS(l2g_t, IS_GTOLM_DROP, is2, &is1));
    797:         PetscCall(ISGetLocalSize(is1, &i));
    798:         PetscCall(ISGetLocalSize(is2, &j));
    799:         PetscCheck(i == j, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent local sizes %" PetscInt_FMT " and %" PetscInt_FMT " for iV", i, j);
    800:         PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_iV", (PetscObject)is1));
    801:         PetscCall(ISLocalToGlobalMappingDestroy(&l2g_t));
    802:         PetscCall(ISDestroy(&is1));
    803:         PetscCall(ISDestroy(&is2));
    804:       }

    806:       /* exclude selected pressures from the inner BDDC */
    807:       if (pcbddc->DirichletBoundariesLocal) {
    808:         IS list[2], plP, isout;

    810:         /* need a parallel IS */
    811:         PetscCall(ISOnComm(lP, PetscObjectComm((PetscObject)ksp), PETSC_COPY_VALUES, &plP));
    812:         list[0] = plP;
    813:         list[1] = pcbddc->DirichletBoundariesLocal;
    814:         PetscCall(ISConcatenate(PetscObjectComm((PetscObject)ksp), 2, list, &isout));
    815:         PetscCall(ISSortRemoveDups(isout));
    816:         PetscCall(ISDestroy(&plP));
    817:         PetscCall(PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc, isout));
    818:         PetscCall(ISDestroy(&isout));
    819:       } else if (pcbddc->DirichletBoundaries) {
    820:         IS list[2], isout;

    822:         list[0] = pP;
    823:         list[1] = pcbddc->DirichletBoundaries;
    824:         PetscCall(ISConcatenate(PetscObjectComm((PetscObject)ksp), 2, list, &isout));
    825:         PetscCall(ISSortRemoveDups(isout));
    826:         PetscCall(PCBDDCSetDirichletBoundaries(fetidp->innerbddc, isout));
    827:         PetscCall(ISDestroy(&isout));
    828:       } else {
    829:         IS plP;

    831:         /* need a parallel IS */
    832:         PetscCall(ISOnComm(lP, PetscObjectComm((PetscObject)ksp), PETSC_COPY_VALUES, &plP));
    833:         PetscCall(PCBDDCSetDirichletBoundariesLocal(fetidp->innerbddc, plP));
    834:         PetscCall(ISDestroy(&plP));
    835:       }

    837:       /* Need to zero output of interface BDDC for lP */
    838:       {
    839:         IS BB, lP_I, lP_B;

    841:         PetscCall(ISComplement(II, 0, n, &BB));
    842:         PetscCall(ISEmbed(lP, II, PETSC_TRUE, &lP_I));
    843:         PetscCall(ISEmbed(lP, BB, PETSC_TRUE, &lP_B));
    844:         PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_lP_I", (PetscObject)lP_I));
    845:         PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_lP_B", (PetscObject)lP_B));
    846:         PetscCall(ISDestroy(&BB));
    847:         PetscCall(ISDestroy(&lP_I));
    848:         PetscCall(ISDestroy(&lP_B));
    849:       }
    850:       PetscCall(ISDestroy(&II));

    852:       /* save CSR information for the pressure BDDC solver (if any) */
    853:       if (schp) {
    854:         PetscInt np, nt;

    856:         PetscCall(MatGetSize(matis->A, &nt, NULL));
    857:         PetscCall(ISGetLocalSize(lP, &np));
    858:         if (np) {
    859:           PetscInt *xadj = pcbddc->mat_graph->xadj;
    860:           PetscInt *adjn = pcbddc->mat_graph->adjncy;
    861:           PetscInt  nv   = pcbddc->mat_graph->nvtxs_csr;

    863:           if (nv && nv == nt) {
    864:             ISLocalToGlobalMapping pmap;
    865:             PetscInt              *schp_csr, *schp_xadj, *schp_adjn, p;
    866:             PetscContainer         c;

    868:             PetscCall(ISLocalToGlobalMappingCreateIS(lPall, &pmap));
    869:             PetscCall(ISGetIndices(lPall, &idxs));
    870:             for (p = 0, nv = 0; p < np; p++) {
    871:               PetscInt x, n = idxs[p];

    873:               PetscCall(ISGlobalToLocalMappingApply(pmap, IS_GTOLM_DROP, xadj[n + 1] - xadj[n], adjn + xadj[n], &x, NULL));
    874:               nv += x;
    875:             }
    876:             PetscCall(PetscMalloc1(np + 1 + nv, &schp_csr));
    877:             schp_xadj = schp_csr;
    878:             schp_adjn = schp_csr + np + 1;
    879:             for (p = 0, schp_xadj[0] = 0; p < np; p++) {
    880:               PetscInt x, n = idxs[p];

    882:               PetscCall(ISGlobalToLocalMappingApply(pmap, IS_GTOLM_DROP, xadj[n + 1] - xadj[n], adjn + xadj[n], &x, schp_adjn + schp_xadj[p]));
    883:               schp_xadj[p + 1] = schp_xadj[p] + x;
    884:             }
    885:             PetscCall(ISRestoreIndices(lPall, &idxs));
    886:             PetscCall(ISLocalToGlobalMappingDestroy(&pmap));
    887:             PetscCall(PetscContainerCreate(PETSC_COMM_SELF, &c));
    888:             PetscCall(PetscContainerSetPointer(c, schp_csr));
    889:             PetscCall(PetscContainerSetCtxDestroy(c, PetscCtxDestroyDefault));
    890:             PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_pCSR", (PetscObject)c));
    891:             PetscCall(PetscContainerDestroy(&c));
    892:           }
    893:         }
    894:       }
    895:       PetscCall(ISDestroy(&lPall));
    896:       PetscCall(ISDestroy(&lP));
    897:       fetidp->pP = pP;
    898:     }

    900:     /* total number of selected pressure dofs */
    901:     PetscCall(ISGetSize(fetidp->pP, &totP));

    903:     /* Set operator for inner BDDC */
    904:     if (totP || fetidp->rhs_flip) {
    905:       PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &nA));
    906:     } else {
    907:       PetscCall(PetscObjectReference((PetscObject)A));
    908:       nA = A;
    909:     }
    910:     if (fetidp->rhs_flip) {
    911:       PetscCall(MatDiagonalScale(nA, fetidp->rhs_flip, NULL));
    912:       if (totP) {
    913:         Mat lA2;

    915:         PetscCall(MatISGetLocalMat(nA, &lA));
    916:         PetscCall(MatDuplicate(lA, MAT_COPY_VALUES, &lA2));
    917:         PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_lA", (PetscObject)lA2));
    918:         PetscCall(MatDestroy(&lA2));
    919:       }
    920:     }

    922:     /* assign operator to compute inner bddc */
    923:     if (totP) {
    924:       /* in this case, we remove all the used pressure couplings */
    925:       PetscCall(MatSetOption(nA, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_FALSE));
    926:       PetscCall(MatZeroRowsColumnsIS(nA, fetidp->pP, 1., NULL, NULL));
    927:     } else {
    928:       PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_lA", NULL));
    929:     }
    930:     PetscCall(MatGetNearNullSpace(Ap, &nnsp));
    931:     if (!nnsp) PetscCall(MatGetNullSpace(Ap, &nnsp));
    932:     if (!nnsp) PetscCall(MatGetNearNullSpace(A, &nnsp));
    933:     if (!nnsp) PetscCall(MatGetNullSpace(A, &nnsp));
    934:     PetscCall(MatSetNearNullSpace(nA, nnsp));
    935:     PetscCall(PCSetOperators(fetidp->innerbddc, nA, nA));
    936:     PetscCall(MatDestroy(&nA));

    938:     /* non-zero rhs on interior dofs when applying the preconditioner */
    939:     if (totP) pcbddc->switch_static = PETSC_TRUE;

    941:     /* if there are no interface pressures, set inner bddc flag for benign saddle point */
    942:     if (!totP) {
    943:       pcbddc->benign_saddle_point = PETSC_TRUE;
    944:       pcbddc->compute_nonetflux   = PETSC_TRUE;
    945:     }

    947:     /* Operators for pressure preconditioner */
    948:     if (totP) {
    949:       /* Extract pressure block if needed */
    950:       if (!pisz) {
    951:         Mat C;
    952:         IS  nzrows = NULL;

    954:         PetscCall(MatCreateSubMatrix(A, fetidp->pP, fetidp->pP, MAT_INITIAL_MATRIX, &C));
    955:         PetscCall(MatFindNonzeroRows(C, &nzrows));
    956:         if (nzrows) {
    957:           PetscInt i;

    959:           PetscCall(ISGetSize(nzrows, &i));
    960:           PetscCall(ISDestroy(&nzrows));
    961:           if (!i) pisz = PETSC_TRUE;
    962:         }
    963:         if (!pisz) {
    964:           PetscCall(MatScale(C, -1.)); /* i.e. Almost Incompressible Elasticity, Stokes discretized with Q1xQ1_stabilized, Biot... */
    965:           PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_C", (PetscObject)C));
    966:         }
    967:         PetscCall(MatDestroy(&C));
    968:       }
    969:       /* Divergence mat */
    970:       if (!pcbddc->divudotp) {
    971:         Mat       B;
    972:         IS        P;
    973:         IS        l2l = NULL;
    974:         PetscBool save;

    976:         PetscCall(PetscObjectQuery((PetscObject)fetidp->innerbddc, "__KSPFETIDP_aP", (PetscObject *)&P));
    977:         if (!pisz) {
    978:           IS       F, V, Ps;
    979:           PetscInt m, M;

    981:           PetscCall(MatGetOwnershipRange(A, &m, &M));
    982:           PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A), M - m, m, 1, &F));
    983:           PetscCall(ISDuplicate(P, &Ps));
    984:           PetscCall(ISSort(Ps));
    985:           PetscCall(ISComplement(Ps, m, M, &V));
    986:           PetscCall(ISDestroy(&Ps));
    987:           PetscCall(MatCreateSubMatrix(A, P, V, MAT_INITIAL_MATRIX, &B));
    988:           {
    989:             Mat_IS *Bmatis = (Mat_IS *)B->data;
    990:             PetscCall(PetscObjectReference((PetscObject)Bmatis->getsub_cis));
    991:             l2l = Bmatis->getsub_cis;
    992:           }
    993:           PetscCall(ISDestroy(&V));
    994:           PetscCall(ISDestroy(&F));
    995:         } else {
    996:           PetscCall(MatCreateSubMatrix(A, P, NULL, MAT_INITIAL_MATRIX, &B));
    997:         }
    998:         save = pcbddc->compute_nonetflux; /* SetDivergenceMat activates nonetflux computation */
    999:         PetscCall(PCBDDCSetDivergenceMat(fetidp->innerbddc, B, PETSC_FALSE, l2l));
    1000:         pcbddc->compute_nonetflux = save;
    1001:         PetscCall(MatDestroy(&B));
    1002:         PetscCall(ISDestroy(&l2l));
    1003:       }
    1004:       if (A != Ap) { /* user has provided a different Pmat, this always supersedes the setter (TODO: is it OK?) */
    1005:         /* use monolithic operator, we restrict later */
    1006:         PetscCall(KSPFETIDPSetPressureOperator(ksp, Ap));
    1007:       }
    1008:       PetscCall(PetscObjectQuery((PetscObject)fetidp->innerbddc, "__KSPFETIDP_PPmat", (PetscObject *)&PPmat));

    1010:       /* PPmat not present, use some default choice */
    1011:       if (!PPmat) {
    1012:         Mat C;

    1014:         PetscCall(PetscObjectQuery((PetscObject)fetidp->innerbddc, "__KSPFETIDP_C", (PetscObject *)&C));
    1015:         if (!schp && C) { /* non-zero pressure block, most likely Almost Incompressible Elasticity */
    1016:           PetscCall(KSPFETIDPSetPressureOperator(ksp, C));
    1017:         } else if (!pisz && schp) { /* we need the whole pressure mass matrix to define the interface BDDC */
    1018:           IS P;

    1020:           PetscCall(PetscObjectQuery((PetscObject)fetidp->innerbddc, "__KSPFETIDP_aP", (PetscObject *)&P));
    1021:           PetscCall(MatCreateSubMatrix(A, P, P, MAT_INITIAL_MATRIX, &C));
    1022:           PetscCall(MatScale(C, -1.));
    1023:           PetscCall(KSPFETIDPSetPressureOperator(ksp, C));
    1024:           PetscCall(MatDestroy(&C));
    1025:         } else { /* identity (need to be scaled properly by the user using e.g. a Richardson method */
    1026:           PetscInt nl;

    1028:           PetscCall(ISGetLocalSize(fetidp->pP, &nl));
    1029:           PetscCall(MatCreate(PetscObjectComm((PetscObject)ksp), &C));
    1030:           PetscCall(MatSetSizes(C, nl, nl, totP, totP));
    1031:           PetscCall(MatSetType(C, MATAIJ));
    1032:           PetscCall(MatMPIAIJSetPreallocation(C, 1, NULL, 0, NULL));
    1033:           PetscCall(MatSeqAIJSetPreallocation(C, 1, NULL));
    1034:           PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
    1035:           PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
    1036:           PetscCall(MatShift(C, 1.));
    1037:           PetscCall(KSPFETIDPSetPressureOperator(ksp, C));
    1038:           PetscCall(MatDestroy(&C));
    1039:         }
    1040:       }

    1042:       /* Preconditioned operator for the pressure block */
    1043:       PetscCall(PetscObjectQuery((PetscObject)fetidp->innerbddc, "__KSPFETIDP_PPmat", (PetscObject *)&PPmat));
    1044:       if (PPmat) {
    1045:         Mat      C;
    1046:         IS       Pall;
    1047:         PetscInt AM, PAM, PAN, pam, pan, am, an, pl, pIl, pAg, pIg;

    1049:         PetscCall(PetscObjectQuery((PetscObject)fetidp->innerbddc, "__KSPFETIDP_aP", (PetscObject *)&Pall));
    1050:         PetscCall(MatGetSize(A, &AM, NULL));
    1051:         PetscCall(MatGetSize(PPmat, &PAM, &PAN));
    1052:         PetscCall(ISGetSize(Pall, &pAg));
    1053:         PetscCall(ISGetSize(fetidp->pP, &pIg));
    1054:         PetscCall(MatGetLocalSize(PPmat, &pam, &pan));
    1055:         PetscCall(MatGetLocalSize(A, &am, &an));
    1056:         PetscCall(ISGetLocalSize(Pall, &pIl));
    1057:         PetscCall(ISGetLocalSize(fetidp->pP, &pl));
    1058:         PetscCheck(PAM == PAN, PetscObjectComm((PetscObject)ksp), PETSC_ERR_USER, "Pressure matrix must be square, unsupported %" PetscInt_FMT " x %" PetscInt_FMT, PAM, PAN);
    1059:         PetscCheck(pam == pan, PetscObjectComm((PetscObject)ksp), PETSC_ERR_USER, "Local sizes of pressure matrix must be equal, unsupported %" PetscInt_FMT " x %" PetscInt_FMT, pam, pan);
    1060:         PetscCheck(pam == am || pam == pl || pam == pIl, PETSC_COMM_SELF, PETSC_ERR_USER, "Invalid number of local rows %" PetscInt_FMT " for pressure matrix! Supported are %" PetscInt_FMT ", %" PetscInt_FMT " or %" PetscInt_FMT, pam, am, pl, pIl);
    1061:         PetscCheck(pan == an || pan == pl || pan == pIl, PETSC_COMM_SELF, PETSC_ERR_USER, "Invalid number of local columns %" PetscInt_FMT " for pressure matrix! Supported are %" PetscInt_FMT ", %" PetscInt_FMT " or %" PetscInt_FMT, pan, an, pl, pIl);
    1062:         if (PAM == AM) { /* monolithic ordering, restrict to pressure */
    1063:           if (schp) {
    1064:             PetscCall(MatCreateSubMatrix(PPmat, Pall, Pall, MAT_INITIAL_MATRIX, &C));
    1065:             PetscCall(MatNullSpacePropagateAny_Private(PPmat, Pall, C));
    1066:           } else {
    1067:             PetscCall(MatCreateSubMatrix(PPmat, fetidp->pP, fetidp->pP, MAT_INITIAL_MATRIX, &C));
    1068:           }
    1069:         } else if (pAg == PAM) { /* global ordering for pressure only */
    1070:           if (!allp && !schp) {  /* solving for interface pressure only */
    1071:             IS restr;

    1073:             PetscCall(ISRenumber(fetidp->pP, NULL, NULL, &restr));
    1074:             PetscCall(MatCreateSubMatrix(PPmat, restr, restr, MAT_INITIAL_MATRIX, &C));
    1075:             PetscCall(ISDestroy(&restr));
    1076:           } else {
    1077:             PetscCall(PetscObjectReference((PetscObject)PPmat));
    1078:             C = PPmat;
    1079:           }
    1080:         } else if (pIg == PAM) { /* global ordering for selected pressure only */
    1081:           PetscCheck(!schp, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "Need the entire matrix");
    1082:           PetscCall(PetscObjectReference((PetscObject)PPmat));
    1083:           C = PPmat;
    1084:         } else SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_USER, "Unable to use the pressure matrix");

    1086:         PetscCall(KSPFETIDPSetPressureOperator(ksp, C));
    1087:         PetscCall(MatDestroy(&C));
    1088:       } else SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "Missing Pmat for pressure block");
    1089:     } else { /* totP == 0 */
    1090:       PetscCall(PetscObjectCompose((PetscObject)fetidp->innerbddc, "__KSPFETIDP_pP", NULL));
    1091:     }
    1092:   }
    1093:   PetscFunctionReturn(PETSC_SUCCESS);
    1094: }

    1096: static PetscErrorCode KSPSetUp_FETIDP(KSP ksp)
    1097: {
    1098:   KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
    1099:   PC_BDDC    *pcbddc = (PC_BDDC *)fetidp->innerbddc->data;
    1100:   PetscBool   flg;

    1102:   PetscFunctionBegin;
    1103:   PetscCall(KSPFETIDPSetUpOperators(ksp));
    1104:   /* set up BDDC */
    1105:   PetscCall(PCSetErrorIfFailure(fetidp->innerbddc, ksp->errorifnotconverged));
    1106:   PetscCall(PCSetUp(fetidp->innerbddc));
    1107:   /* FETI-DP as it is implemented needs an exact coarse solver */
    1108:   if (pcbddc->coarse_ksp) {
    1109:     PetscCall(KSPSetTolerances(pcbddc->coarse_ksp, PETSC_SMALL, PETSC_SMALL, PETSC_CURRENT, 1000));
    1110:     PetscCall(KSPSetNormType(pcbddc->coarse_ksp, KSP_NORM_DEFAULT));
    1111:   }
    1112:   /* FETI-DP as it is implemented needs exact local Neumann solvers */
    1113:   PetscCall(KSPSetTolerances(pcbddc->ksp_R, PETSC_SMALL, PETSC_SMALL, PETSC_CURRENT, 1000));
    1114:   PetscCall(KSPSetNormType(pcbddc->ksp_R, KSP_NORM_DEFAULT));

    1116:   /* setup FETI-DP operators
    1117:      If fetidp->statechanged is true, we need to update the operators
    1118:      needed in the saddle-point case. This should be replaced
    1119:      by a better logic when the FETI-DP matrix and preconditioner will
    1120:      have their own classes */
    1121:   if (pcbddc->new_primal_space || fetidp->statechanged) {
    1122:     Mat F; /* the FETI-DP matrix */
    1123:     PC  D; /* the FETI-DP preconditioner */
    1124:     PetscCall(KSPReset(fetidp->innerksp));
    1125:     PetscCall(PCBDDCCreateFETIDPOperators(fetidp->innerbddc, fetidp->fully_redundant, ((PetscObject)ksp)->prefix, &F, &D));
    1126:     PetscCall(KSPSetOperators(fetidp->innerksp, F, F));
    1127:     PetscCall(KSPSetTolerances(fetidp->innerksp, ksp->rtol, ksp->abstol, ksp->divtol, ksp->max_it));
    1128:     PetscCall(KSPSetPC(fetidp->innerksp, D));
    1129:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)D, (PetscObject)fetidp->innerksp, 0));
    1130:     PetscCall(KSPSetFromOptions(fetidp->innerksp));
    1131:     PetscCall(MatCreateVecs(F, &fetidp->innerksp->vec_rhs, &fetidp->innerksp->vec_sol));
    1132:     PetscCall(MatDestroy(&F));
    1133:     PetscCall(PCDestroy(&D));
    1134:     if (fetidp->check) {
    1135:       PetscViewer viewer;

    1137:       if (!pcbddc->dbg_viewer) {
    1138:         viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp));
    1139:       } else {
    1140:         viewer = pcbddc->dbg_viewer;
    1141:       }
    1142:       PetscCall(KSPFETIDPCheckOperators(ksp, viewer));
    1143:     }
    1144:   }
    1145:   fetidp->statechanged     = PETSC_FALSE;
    1146:   pcbddc->new_primal_space = PETSC_FALSE;

    1148:   /* propagate settings to the inner solve */
    1149:   PetscCall(KSPGetComputeSingularValues(ksp, &flg));
    1150:   PetscCall(KSPSetComputeSingularValues(fetidp->innerksp, flg));
    1151:   if (ksp->res_hist) PetscCall(KSPSetResidualHistory(fetidp->innerksp, ksp->res_hist, ksp->res_hist_max, ksp->res_hist_reset));
    1152:   PetscCall(KSPSetErrorIfNotConverged(fetidp->innerksp, ksp->errorifnotconverged));
    1153:   PetscCall(KSPSetUp(fetidp->innerksp));
    1154:   PetscFunctionReturn(PETSC_SUCCESS);
    1155: }

    1157: static PetscErrorCode KSPSolve_FETIDP(KSP ksp)
    1158: {
    1159:   Mat                F, A;
    1160:   MatNullSpace       nsp;
    1161:   Vec                X, B, Xl, Bl;
    1162:   KSP_FETIDP        *fetidp = (KSP_FETIDP *)ksp->data;
    1163:   PC_BDDC           *pcbddc = (PC_BDDC *)fetidp->innerbddc->data;
    1164:   KSPConvergedReason reason;
    1165:   PC                 pc;
    1166:   PCFailedReason     pcreason;
    1167:   PetscInt           hist_len;
    1168:   int                flg;

    1170:   PetscFunctionBegin;
    1171:   PetscCall(PetscCitationsRegister(citation, &cited));
    1172:   if (fetidp->saddlepoint) PetscCall(PetscCitationsRegister(citation2, &cited2));
    1173:   PetscCall(KSPGetOperators(ksp, &A, NULL));
    1174:   PetscCall(KSPGetRhs(ksp, &B));
    1175:   PetscCall(KSPGetSolution(ksp, &X));
    1176:   PetscCall(KSPGetOperators(fetidp->innerksp, &F, NULL));
    1177:   PetscCall(KSPGetRhs(fetidp->innerksp, &Bl));
    1178:   PetscCall(KSPGetSolution(fetidp->innerksp, &Xl));
    1179:   PetscCall(PCBDDCMatFETIDPGetRHS(F, B, Bl));
    1180:   if (ksp->transpose_solve) {
    1181:     PetscCall(KSPSolveTranspose(fetidp->innerksp, Bl, Xl));
    1182:   } else {
    1183:     PetscCall(KSPSolve(fetidp->innerksp, Bl, Xl));
    1184:   }
    1185:   PetscCall(KSPGetConvergedReason(fetidp->innerksp, &reason));
    1186:   PetscCall(KSPGetPC(fetidp->innerksp, &pc));
    1187:   PetscCall(PCGetFailedReason(pc, &pcreason));
    1188:   flg = (reason < 0 && reason != KSP_DIVERGED_ITS) || pcreason;
    1189:   PetscCall(VecFlag(Xl, flg));
    1190:   if (flg) {
    1191:     PetscInt its;

    1193:     PetscCall(KSPGetIterationNumber(fetidp->innerksp, &its));
    1194:     ksp->reason = KSP_DIVERGED_PC_FAILED;
    1195:     PetscCall(PetscInfo(ksp, "Inner KSP solve failed: %s %s at iteration %" PetscInt_FMT "\n", KSPConvergedReasons[reason], PCFailedReasons[pcreason], its));
    1196:   }
    1197:   PetscCall(PCBDDCMatFETIDPGetSolution(F, Xl, X));
    1198:   PetscCall(MatGetNullSpace(A, &nsp));
    1199:   if (nsp) PetscCall(MatNullSpaceRemove(nsp, X));
    1200:   /* update ksp with stats from inner ksp */
    1201:   PetscCall(KSPGetConvergedReason(fetidp->innerksp, &ksp->reason));
    1202:   PetscCall(KSPGetIterationNumber(fetidp->innerksp, &ksp->its));
    1203:   ksp->totalits += ksp->its;
    1204:   PetscCall(KSPGetResidualHistory(fetidp->innerksp, NULL, &hist_len));
    1205:   ksp->res_hist_len = (size_t)hist_len;
    1206:   /* restore defaults for inner BDDC (Pre/PostSolve flags) */
    1207:   pcbddc->temp_solution_used        = PETSC_FALSE;
    1208:   pcbddc->rhs_change                = PETSC_FALSE;
    1209:   pcbddc->exact_dirichlet_trick_app = PETSC_FALSE;
    1210:   PetscFunctionReturn(PETSC_SUCCESS);
    1211: }

    1213: static PetscErrorCode KSPReset_FETIDP(KSP ksp)
    1214: {
    1215:   KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
    1216:   PC_BDDC    *pcbddc;

    1218:   PetscFunctionBegin;
    1219:   PetscCall(ISDestroy(&fetidp->pP));
    1220:   PetscCall(VecDestroy(&fetidp->rhs_flip));
    1221:   /* avoid PCReset that does not take into account ref counting */
    1222:   PetscCall(PCDestroy(&fetidp->innerbddc));
    1223:   PetscCall(PCCreate(PetscObjectComm((PetscObject)ksp), &fetidp->innerbddc));
    1224:   PetscCall(PCSetType(fetidp->innerbddc, PCBDDC));
    1225:   pcbddc                   = (PC_BDDC *)fetidp->innerbddc->data;
    1226:   pcbddc->symmetric_primal = PETSC_FALSE;
    1227:   PetscCall(KSPDestroy(&fetidp->innerksp));
    1228:   fetidp->saddlepoint  = PETSC_FALSE;
    1229:   fetidp->matstate     = -1;
    1230:   fetidp->matnnzstate  = -1;
    1231:   fetidp->statechanged = PETSC_TRUE;
    1232:   PetscFunctionReturn(PETSC_SUCCESS);
    1233: }

    1235: static PetscErrorCode KSPDestroy_FETIDP(KSP ksp)
    1236: {
    1237:   KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;

    1239:   PetscFunctionBegin;
    1240:   PetscCall(KSPReset_FETIDP(ksp));
    1241:   PetscCall(PCDestroy(&fetidp->innerbddc));
    1242:   PetscCall(KSPDestroy(&fetidp->innerksp));
    1243:   PetscCall(PetscFree(fetidp->monctx));
    1244:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFETIDPSetInnerBDDC_C", NULL));
    1245:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFETIDPGetInnerBDDC_C", NULL));
    1246:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFETIDPGetInnerKSP_C", NULL));
    1247:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFETIDPSetPressureOperator_C", NULL));
    1248:   PetscCall(PetscFree(ksp->data));
    1249:   PetscFunctionReturn(PETSC_SUCCESS);
    1250: }

    1252: static PetscErrorCode KSPView_FETIDP(KSP ksp, PetscViewer viewer)
    1253: {
    1254:   KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;
    1255:   PetscBool   isascii;

    1257:   PetscFunctionBegin;
    1258:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    1259:   if (isascii) {
    1260:     PetscCall(PetscViewerASCIIPrintf(viewer, "  fully redundant: %d\n", fetidp->fully_redundant));
    1261:     PetscCall(PetscViewerASCIIPrintf(viewer, "  saddle point:    %d\n", fetidp->saddlepoint));
    1262:     PetscCall(PetscViewerASCIIPrintf(viewer, "Inner KSP solver details\n"));
    1263:   }
    1264:   PetscCall(PetscViewerASCIIPushTab(viewer));
    1265:   PetscCall(KSPView(fetidp->innerksp, viewer));
    1266:   PetscCall(PetscViewerASCIIPopTab(viewer));
    1267:   if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "Inner BDDC solver details\n"));
    1268:   PetscCall(PetscViewerASCIIPushTab(viewer));
    1269:   PetscCall(PCView(fetidp->innerbddc, viewer));
    1270:   PetscCall(PetscViewerASCIIPopTab(viewer));
    1271:   PetscFunctionReturn(PETSC_SUCCESS);
    1272: }

    1274: static PetscErrorCode KSPSetFromOptions_FETIDP(KSP ksp, PetscOptionItems PetscOptionsObject)
    1275: {
    1276:   KSP_FETIDP *fetidp = (KSP_FETIDP *)ksp->data;

    1278:   PetscFunctionBegin;
    1279:   /* set options prefixes for the inner objects, since the parent prefix will be valid at this point */
    1280:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerksp, ((PetscObject)ksp)->prefix));
    1281:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerksp, "fetidp_"));
    1282:   if (!fetidp->userbddc) {
    1283:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)fetidp->innerbddc, ((PetscObject)ksp)->prefix));
    1284:     PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)fetidp->innerbddc, "fetidp_bddc_"));
    1285:   }
    1286:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP FETIDP options");
    1287:   PetscCall(PetscOptionsBool("-ksp_fetidp_fullyredundant", "Use fully redundant multipliers", "none", fetidp->fully_redundant, &fetidp->fully_redundant, NULL));
    1288:   PetscCall(PetscOptionsBool("-ksp_fetidp_saddlepoint", "Activates support for saddle-point problems", NULL, fetidp->saddlepoint, &fetidp->saddlepoint, NULL));
    1289:   PetscCall(PetscOptionsBool("-ksp_fetidp_check", "Activates verbose debugging output FETI-DP operators", NULL, fetidp->check, &fetidp->check, NULL));
    1290:   PetscOptionsHeadEnd();
    1291:   PetscCall(PCSetFromOptions(fetidp->innerbddc));
    1292:   PetscFunctionReturn(PETSC_SUCCESS);
    1293: }

    1295: /*MC
    1296:      KSPFETIDP - The FETI-DP method {cite}`farhat2001feti`

    1298:    Options Database Keys:
    1299: +   -ksp_fetidp_fullyredundant (true|false)   - use a fully redundant set of Lagrange multipliers
    1300: .   -ksp_fetidp_saddlepoint (true|false)      - activates support for saddle point problems, see {cite}`tu2015feti`
    1301: .   -ksp_fetidp_saddlepoint_flip (true|false) - see note below
    1302: .   -ksp_fetidp_pressure_field -1             - activates support for saddle point problems, and identifies the pressure field id.
    1303:                                                 If this information is not provided, the pressure field is detected by using `MatFindZeroDiagonals()`.
    1304: -   -ksp_fetidp_pressure_all (true|false)     - if false, uses the interface pressures, as described in [2]. If true, uses the entire pressure field.

    1306:    Level: Advanced

    1308:    Notes:
    1309:    The matrix for the `KSP` must be of type `MATIS`.

    1311:    Usually, an incompressible Stokes problem is written as
    1312: .vb
    1313:    | A B^T | | v | = | f |
    1314:    | B 0   | | p | = | g |
    1315: .ve
    1316:    with B representing $ -\int_\Omega \nabla \cdot u q $. If -ksp_fetidp_saddlepoint_flip is true, the code assumes that the user provides it as
    1317: .vb
    1318:    | A B^T | | v | = | f |
    1319:    |-B 0   | | p | = |-g |
    1320: .ve

    1322:    The FETI-DP linear system (automatically generated constructing an internal `PCBDDC` object) is solved using an internal `KSP` object.

    1324:    Options for the inner `KSP` and for the customization of the `PCBDDC` object can be specified at command line by using the prefixes `-fetidp_` and `-fetidp_bddc_`. E.g.,
    1325: .vb
    1326:    -fetidp_ksp_type gmres -fetidp_bddc_pc_bddc_symmetric false
    1327: .ve
    1328:    will use `KSPGMRES` for the solution of the linear system on the Lagrange multipliers, generated using a non-symmetric `PCBDDC`.

    1330:    For saddle point problems with continuous pressures, the preconditioned operator for the pressure solver can be specified with `KSPFETIDPSetPressureOperator()`.
    1331:    Alternatively, the pressure operator is extracted from the precondioned matrix (if it is different from the linear solver matrix).
    1332:    If none of the above, an identity matrix will be created; the user then needs to scale it through a Richardson solver.
    1333:    Options for the pressure solver can be prefixed with `-fetidp_fielsplit_p_`, E.g.
    1334: .vb
    1335:    -fetidp_fielsplit_p_ksp_type preonly -fetidp_fielsplit_p_pc_type lu -fetidp_fielsplit_p_pc_factor_mat_solver_type mumps
    1336: .ve
    1337:    In order to use the deluxe version of FETI-DP, you must customize the inner `PCBDDC` operator with -fetidp_bddc_pc_bddc_use_deluxe_scaling -fetidp_bddc_pc_bddc_deluxe_singlemat and use
    1338:    non-redundant multipliers, i.e. `-ksp_fetidp_fullyredundant false`. Options for the scaling solver are prefixed by `-fetidp_bddelta_`, E.g.
    1339: .vb
    1340:    -fetidp_bddelta_pc_factor_mat_solver_type mumps -fetidp_bddelta_pc_type lu
    1341: .ve

    1343:    Some of the basic options such as the maximum number of iterations and tolerances are automatically passed from this `KSP` to the inner `KSP` that actually performs the iterations.

    1345:    The converged reason and number of iterations computed are passed from the inner `KSP` to this `KSP` at the end of the solution.

    1347:    Developer Note:
    1348:    Even though this method does not directly use any norms, the user is allowed to set the `KSPNormType` to any value.
    1349:    This is so users do not have to change `KSPNormType` options when they switch from other `KSP` methods to this one.

    1351: .seealso: [](ch_ksp), `MATIS`, `PCBDDC`, `KSPFETIDPSetInnerBDDC()`, `KSPFETIDPGetInnerBDDC()`, `KSPFETIDPGetInnerKSP()`
    1352: M*/
    1353: PETSC_EXTERN PetscErrorCode KSPCreate_FETIDP(KSP ksp)
    1354: {
    1355:   KSP_FETIDP    *fetidp;
    1356:   KSP_FETIDPMon *monctx;
    1357:   PC_BDDC       *pcbddc;
    1358:   PC             pc;

    1360:   PetscFunctionBegin;
    1361:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 3));
    1362:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 2));
    1363:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
    1364:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_RIGHT, 2));
    1365:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
    1366:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
    1367:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 2));

    1369:   PetscCall(PetscNew(&fetidp));
    1370:   fetidp->matstate     = -1;
    1371:   fetidp->matnnzstate  = -1;
    1372:   fetidp->statechanged = PETSC_TRUE;

    1374:   ksp->data                              = (void *)fetidp;
    1375:   ksp->ops->setup                        = KSPSetUp_FETIDP;
    1376:   ksp->ops->solve                        = KSPSolve_FETIDP;
    1377:   ksp->ops->destroy                      = KSPDestroy_FETIDP;
    1378:   ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_FETIDP;
    1379:   ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_FETIDP;
    1380:   ksp->ops->view                         = KSPView_FETIDP;
    1381:   ksp->ops->setfromoptions               = KSPSetFromOptions_FETIDP;
    1382:   ksp->ops->buildsolution                = KSPBuildSolution_FETIDP;
    1383:   ksp->ops->buildresidual                = KSPBuildResidualDefault;
    1384:   PetscCall(KSPGetPC(ksp, &pc));
    1385:   PetscCall(PCSetType(pc, PCNONE));
    1386:   /* create the inner KSP for the Lagrange multipliers */
    1387:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)ksp), &fetidp->innerksp));
    1388:   PetscCall(KSPGetPC(fetidp->innerksp, &pc));
    1389:   PetscCall(PCSetType(pc, PCNONE));
    1390:   /* monitor */
    1391:   PetscCall(PetscNew(&monctx));
    1392:   monctx->parentksp = ksp;
    1393:   fetidp->monctx    = monctx;
    1394:   PetscCall(KSPMonitorSet(fetidp->innerksp, KSPMonitor_FETIDP, fetidp->monctx, NULL));
    1395:   /* create the inner BDDC */
    1396:   PetscCall(PCCreate(PetscObjectComm((PetscObject)ksp), &fetidp->innerbddc));
    1397:   PetscCall(PCSetType(fetidp->innerbddc, PCBDDC));
    1398:   /* make sure we always obtain a consistent FETI-DP matrix
    1399:      for symmetric problems, the user can always customize it through the command line */
    1400:   pcbddc                   = (PC_BDDC *)fetidp->innerbddc->data;
    1401:   pcbddc->symmetric_primal = PETSC_FALSE;
    1402:   /* composed functions */
    1403:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFETIDPSetInnerBDDC_C", KSPFETIDPSetInnerBDDC_FETIDP));
    1404:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFETIDPGetInnerBDDC_C", KSPFETIDPGetInnerBDDC_FETIDP));
    1405:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFETIDPGetInnerKSP_C", KSPFETIDPGetInnerKSP_FETIDP));
    1406:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPFETIDPSetPressureOperator_C", KSPFETIDPSetPressureOperator_FETIDP));
    1407:   /* need to call KSPSetUp_FETIDP even with KSP_SETUP_NEWMATRIX */
    1408:   ksp->setupnewmatrix = PETSC_TRUE;
    1409:   PetscFunctionReturn(PETSC_SUCCESS);
    1410: }