Actual source code: hpddm.cxx

      1: #define HPDDM_MIXED_PRECISION 1
      2: #include <petsc/private/petschpddm.h>

      4: const char *const KSPHPDDMTypes[]          = {KSPGMRES, "bgmres", KSPCG, "bcg", "gcrodr", "bgcrodr", "bfbcg", KSPPREONLY};
      5: const char *const HPDDMOrthogonalization[] = {"cgs", "mgs"};
      6: const char *const HPDDMQR[]                = {"cholqr", "cgs", "mgs"};
      7: const char *const HPDDMVariant[]           = {"left", "right", "flexible"};
      8: const char *const HPDDMRecycleTarget[]     = {"SM", "LM", "SR", "LR", "SI", "LI"};
      9: const char *const HPDDMRecycleStrategy[]   = {"A", "B"};

     11: PetscBool  HPDDMCite       = PETSC_FALSE;
     12: const char HPDDMCitation[] = "@article{jolivet2020petsc,\n"
     13:                              "  Author = {Jolivet, Pierre and Roman, Jose E. and Zampini, Stefano},\n"
     14:                              "  Title = {{KSPHPDDM} and {PCHPDDM}: Extending {PETSc} with Robust Overlapping {Schwarz} Preconditioners and Advanced {Krylov} Methods},\n"
     15:                              "  Year = {2021},\n"
     16:                              "  Publisher = {Elsevier},\n"
     17:                              "  Journal = {Computer \\& Mathematics with Applications},\n"
     18:                              "  Volume = {84},\n"
     19:                              "  Pages = {277--295},\n"
     20:                              "  Url = {https://github.com/prj-/jolivet2020petsc}\n"
     21:                              "}\n";

     23: #if PetscDefined(HAVE_SLEPC) && PetscDefined(HAVE_DYNAMIC_LIBRARIES) && PetscDefined(USE_SHARED_LIBRARIES)
     24: static PetscBool loadedDL = PETSC_FALSE;
     25: #endif

     27: static PetscErrorCode KSPSetFromOptions_HPDDM(KSP ksp, PetscOptionItems PetscOptionsObject)
     28: {
     29:   KSP_HPDDM  *data = (KSP_HPDDM *)ksp->data;
     30:   PetscInt    i, j;
     31:   PetscMPIInt size;

     33:   PetscFunctionBegin;
     34:   PetscOptionsHeadBegin(PetscOptionsObject, "KSPHPDDM options, cf. https://github.com/hpddm/hpddm");
     35:   i = (data->cntl[0] == static_cast<char>(PETSC_DECIDE) ? HPDDM_KRYLOV_METHOD_GMRES : data->cntl[0]);
     36:   PetscCall(PetscOptionsEList("-ksp_hpddm_type", "Type of Krylov method", "KSPHPDDMGetType", KSPHPDDMTypes, PETSC_STATIC_ARRAY_LENGTH(KSPHPDDMTypes), KSPHPDDMTypes[HPDDM_KRYLOV_METHOD_GMRES], &i, nullptr));
     37:   if (i == PETSC_STATIC_ARRAY_LENGTH(KSPHPDDMTypes) - 1) i = HPDDM_KRYLOV_METHOD_NONE; /* need to shift the value since HPDDM_KRYLOV_METHOD_RICHARDSON is not registered in PETSc */
     38:   data->cntl[0] = i;
     39:   PetscCall(PetscOptionsEnum("-ksp_hpddm_precision", "Precision in which Krylov bases are stored", "KSPHPDDM", PetscPrecisionTypes, (PetscEnum)data->precision, (PetscEnum *)&data->precision, nullptr));
     40:   PetscCheck(data->precision != PETSC_PRECISION___FLOAT128 || PetscDefined(HAVE_REAL___FLOAT128), PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP_SYS, "Unsupported %s precision", PetscPrecisionTypes[data->precision]);
     41:   PetscCheck(std::abs(data->precision - PETSC_SCALAR_PRECISION) <= 1, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Unsupported mixed %s and %s precisions", PetscPrecisionTypes[data->precision], PetscPrecisionTypes[PETSC_SCALAR_PRECISION]);
     42:   PetscCheck(data->precision != PETSC_PRECISION_INVALID && data->precision != PETSC_PRECISION_BFLOAT16, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Unsupported PetscPrecision %s", PetscPrecisionTypes[data->precision]);
     43:   static_assert(PETSC_PRECISION___FP16 == PETSC_PRECISION_SINGLE - 1 && PETSC_PRECISION_SINGLE == PETSC_PRECISION_DOUBLE - 1 && PETSC_PRECISION_DOUBLE == PETSC_PRECISION___FLOAT128 - 1, "");
     44:   if (data->cntl[0] != HPDDM_KRYLOV_METHOD_NONE) {
     45:     if (data->cntl[0] != HPDDM_KRYLOV_METHOD_BCG && data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG) {
     46:       i = (data->cntl[1] == static_cast<char>(PETSC_DECIDE) ? HPDDM_VARIANT_LEFT : data->cntl[1]);
     47:       if (ksp->pc_side_set == PC_SIDE_DEFAULT)
     48:         PetscCall(PetscOptionsEList("-ksp_hpddm_variant", "Left, right, or variable preconditioning", "KSPHPDDM", HPDDMVariant, PETSC_STATIC_ARRAY_LENGTH(HPDDMVariant), HPDDMVariant[HPDDM_VARIANT_LEFT], &i, nullptr));
     49:       else if (ksp->pc_side_set == PC_RIGHT) i = HPDDM_VARIANT_RIGHT;
     50:       data->cntl[1] = i;
     51:       if (i > 0) PetscCall(KSPSetPCSide(ksp, PC_RIGHT));
     52:     }
     53:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
     54:       data->rcntl[0] = (PetscAbsReal(data->rcntl[0] - static_cast<PetscReal>(PETSC_DECIDE)) < PETSC_SMALL ? -1.0 : data->rcntl[0]);
     55:       PetscCall(PetscOptionsReal("-ksp_hpddm_deflation_tol", "Tolerance when deflating right-hand sides inside block methods", "KSPHPDDM", data->rcntl[0], data->rcntl, nullptr));
     56:       i = (data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG] == static_cast<unsigned short>(PETSC_DECIDE) ? 1 : PetscMax(1, data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG]));
     57:       PetscCall(PetscOptionsRangeInt("-ksp_hpddm_enlarge_krylov_subspace", "Split the initial right-hand side into multiple vectors", "KSPHPDDM", i, &i, nullptr, 1, std::numeric_limits<unsigned short>::max() - 1));
     58:       data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG] = i;
     59:     } else data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BCG] = 0;
     60:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
     61:       i = (data->cntl[2] == static_cast<char>(PETSC_DECIDE) ? HPDDM_ORTHOGONALIZATION_CGS : data->cntl[2] & 3);
     62:       PetscCall(PetscOptionsEList("-ksp_hpddm_orthogonalization", "Classical (faster) or Modified (more robust) Gram--Schmidt process", "KSPHPDDM", HPDDMOrthogonalization, PETSC_STATIC_ARRAY_LENGTH(HPDDMOrthogonalization), HPDDMOrthogonalization[HPDDM_ORTHOGONALIZATION_CGS], &i, nullptr));
     63:       j = (data->cntl[2] == static_cast<char>(PETSC_DECIDE) ? HPDDM_QR_CHOLQR : ((data->cntl[2] >> 2) & 7));
     64:       PetscCall(PetscOptionsEList("-ksp_hpddm_qr", "Distributed QR factorizations computed with Cholesky QR, Classical or Modified Gram--Schmidt process", "KSPHPDDM", HPDDMQR, PETSC_STATIC_ARRAY_LENGTH(HPDDMQR), HPDDMQR[HPDDM_QR_CHOLQR], &j, nullptr));
     65:       data->cntl[2] = static_cast<char>(i) + (static_cast<char>(j) << 2);
     66:       i             = (data->scntl[0] == static_cast<unsigned short>(PETSC_DECIDE) ? PetscMin(30, ksp->max_it) : data->scntl[0]);
     67:       PetscCall(PetscOptionsRangeInt("-ksp_gmres_restart", "Maximum number of Arnoldi vectors generated per cycle", "KSPHPDDM", i, &i, nullptr, PetscMin(1, ksp->max_it), PetscMin(ksp->max_it, std::numeric_limits<unsigned short>::max() - 1)));
     68:       data->scntl[0] = i;
     69:     }
     70:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BCG || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
     71:       j = (data->cntl[1] == static_cast<char>(PETSC_DECIDE) ? HPDDM_QR_CHOLQR : data->cntl[1]);
     72:       PetscCall(PetscOptionsEList("-ksp_hpddm_qr", "Distributed QR factorizations computed with Cholesky QR, Classical or Modified Gram--Schmidt process", "KSPHPDDM", HPDDMQR, PETSC_STATIC_ARRAY_LENGTH(HPDDMQR), HPDDMQR[HPDDM_QR_CHOLQR], &j, nullptr));
     73:       data->cntl[1] = j;
     74:     }
     75:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
     76:       i = (data->icntl[0] == static_cast<int>(PETSC_DECIDE) ? PetscMin(20, data->scntl[0] - 1) : data->icntl[0]);
     77:       PetscCall(PetscOptionsRangeInt("-ksp_hpddm_recycle", "Number of harmonic Ritz vectors to compute", "KSPHPDDM", i, &i, nullptr, 1, data->scntl[0] - 1));
     78:       data->icntl[0] = i;
     79:       if (!PetscDefined(HAVE_SLEPC) || !PetscDefined(USE_SHARED_LIBRARIES) || data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR) {
     80:         i = (data->cntl[3] == static_cast<char>(PETSC_DECIDE) ? HPDDM_RECYCLE_TARGET_SM : data->cntl[3]);
     81:         PetscCall(PetscOptionsEList("-ksp_hpddm_recycle_target", "Criterion to select harmonic Ritz vectors", "KSPHPDDM", HPDDMRecycleTarget, PETSC_STATIC_ARRAY_LENGTH(HPDDMRecycleTarget), HPDDMRecycleTarget[HPDDM_RECYCLE_TARGET_SM], &i, nullptr));
     82:         data->cntl[3] = i;
     83:       } else {
     84:         PetscCheck(data->precision == PETSC_SCALAR_PRECISION, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "Cannot use SLEPc with a different precision than PETSc for harmonic Ritz eigensolves");
     85:         PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ksp), &size));
     86:         i = (data->cntl[3] == static_cast<char>(PETSC_DECIDE) ? 1 : data->cntl[3]);
     87:         PetscCall(PetscOptionsRangeInt("-ksp_hpddm_recycle_redistribute", "Number of processes used to solve eigenvalue problems when recycling in BGCRODR", "KSPHPDDM", i, &i, nullptr, 1, PetscMin(size, 192)));
     88:         data->cntl[3] = i;
     89:       }
     90:       i = (data->cntl[4] == static_cast<char>(PETSC_DECIDE) ? HPDDM_RECYCLE_STRATEGY_A : data->cntl[4]);
     91:       PetscCall(PetscOptionsEList("-ksp_hpddm_recycle_strategy", "Generalized eigenvalue problem to solve for recycling", "KSPHPDDM", HPDDMRecycleStrategy, PETSC_STATIC_ARRAY_LENGTH(HPDDMRecycleStrategy), HPDDMRecycleStrategy[HPDDM_RECYCLE_STRATEGY_A], &i, nullptr));
     92:       data->cntl[4] = i;
     93:     }
     94:   } else {
     95:     data->cntl[0]  = HPDDM_KRYLOV_METHOD_NONE;
     96:     data->scntl[1] = 1;
     97:   }
     98:   PetscCheck(ksp->nmax >= std::numeric_limits<int>::min() && ksp->nmax <= std::numeric_limits<int>::max(), PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "KSPMatSolve() block size %" PetscInt_FMT " not representable by an integer, which is not handled by KSPHPDDM",
     99:              ksp->nmax);
    100:   data->icntl[1] = static_cast<int>(ksp->nmax);
    101:   PetscOptionsHeadEnd();
    102:   PetscFunctionReturn(PETSC_SUCCESS);
    103: }

    105: static PetscErrorCode KSPView_HPDDM(KSP ksp, PetscViewer viewer)
    106: {
    107:   KSP_HPDDM            *data  = (KSP_HPDDM *)ksp->data;
    108:   HPDDM::PETScOperator *op    = data->op;
    109:   const PetscScalar    *array = op ? op->storage() : nullptr;
    110:   PetscBool             ascii;

    112:   PetscFunctionBegin;
    113:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii));
    114:   if (op && ascii) {
    115:     PetscCall(PetscViewerASCIIPrintf(viewer, "HPDDM type: %s%s\n", KSPHPDDMTypes[std::min(static_cast<PetscInt>(data->cntl[0]), static_cast<PetscInt>(PETSC_STATIC_ARRAY_LENGTH(KSPHPDDMTypes) - 1))], data->cntl[1] == HPDDM_VARIANT_FLEXIBLE ? " (with support for variable preconditioning)" : ""));
    116:     PetscCall(PetscViewerASCIIPrintf(viewer, "precision: %s\n", PetscPrecisionTypes[data->precision]));
    117:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
    118:       if (PetscAbsReal(data->rcntl[0] - static_cast<PetscReal>(PETSC_DECIDE)) < PETSC_SMALL) PetscCall(PetscViewerASCIIPrintf(viewer, "no deflation at restarts\n"));
    119:       else PetscCall(PetscViewerASCIIPrintf(viewer, "deflation tolerance: %g\n", static_cast<double>(data->rcntl[0])));
    120:     }
    121:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
    122:       PetscCall(PetscViewerASCIIPrintf(viewer, "deflation subspace attached? %s\n", PetscBools[array ? PETSC_TRUE : PETSC_FALSE]));
    123:       if (!PetscDefined(HAVE_SLEPC) || !PetscDefined(USE_SHARED_LIBRARIES) || data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR) PetscCall(PetscViewerASCIIPrintf(viewer, "deflation target: %s\n", HPDDMRecycleTarget[static_cast<PetscInt>(data->cntl[3])]));
    124:       else PetscCall(PetscViewerASCIIPrintf(viewer, "redistribution size: %d\n", static_cast<PetscMPIInt>(data->cntl[3])));
    125:     }
    126:     if (data->icntl[1] != static_cast<int>(PETSC_DECIDE)) PetscCall(PetscViewerASCIIPrintf(viewer, "  block size is %d\n", data->icntl[1]));
    127:   }
    128:   PetscFunctionReturn(PETSC_SUCCESS);
    129: }

    131: static PetscErrorCode KSPSetUp_HPDDM(KSP ksp)
    132: {
    133:   KSP_HPDDM *data = (KSP_HPDDM *)ksp->data;
    134:   Mat        A;
    135:   PetscInt   n, bs;
    136:   PetscBool  match;

    138:   PetscFunctionBegin;
    139:   PetscCall(KSPGetOperators(ksp, &A, nullptr));
    140:   PetscCall(MatGetLocalSize(A, &n, nullptr));
    141:   PetscCall(MatGetBlockSize(A, &bs));
    142:   PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &match, MATSEQKAIJ, MATMPIKAIJ, ""));
    143:   if (match) n /= bs;
    144:   data->op = new HPDDM::PETScOperator(ksp, n);
    145:   if (PetscUnlikely(!ksp->setfromoptionscalled || data->cntl[0] == static_cast<char>(PETSC_DECIDE))) { /* what follows is basically a copy/paste of KSPSetFromOptions_HPDDM, with no call to PetscOptions() */
    146:     PetscCall(PetscInfo(ksp, "KSPSetFromOptions() not called or uninitialized internal structure, hardwiring default KSPHPDDM options\n"));
    147:     if (data->cntl[0] == static_cast<char>(PETSC_DECIDE)) data->cntl[0] = 0; /* GMRES by default */
    148:     if (data->cntl[0] != HPDDM_KRYLOV_METHOD_NONE) {                         /* following options do not matter with PREONLY */
    149:       if (data->cntl[0] != HPDDM_KRYLOV_METHOD_BCG && data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG) {
    150:         data->cntl[1] = HPDDM_VARIANT_LEFT; /* left preconditioning by default */
    151:         if (ksp->pc_side_set == PC_RIGHT) data->cntl[1] = HPDDM_VARIANT_RIGHT;
    152:         if (data->cntl[1] > 0) PetscCall(KSPSetPCSide(ksp, PC_RIGHT));
    153:       }
    154:       if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
    155:         data->rcntl[0]                                          = -1.0; /* no deflation by default */
    156:         data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG] = 1;    /* Krylov subspace not enlarged by default */
    157:       } else data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BCG] = 0;
    158:       if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
    159:         data->cntl[2]  = static_cast<char>(HPDDM_ORTHOGONALIZATION_CGS) + (static_cast<char>(HPDDM_QR_CHOLQR) << 2); /* CGS and CholQR by default */
    160:         data->scntl[0] = PetscMin(30, ksp->max_it);                                                                  /* restart parameter of 30 by default */
    161:       }
    162:       if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BCG || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) data->cntl[1] = HPDDM_QR_CHOLQR; /* CholQR by default */
    163:       if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
    164:         data->icntl[0] = PetscMin(20, data->scntl[0] - 1); /* recycled subspace of size 20 by default */
    165:         if (!PetscDefined(HAVE_SLEPC) || !PetscDefined(USE_SHARED_LIBRARIES) || data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR) {
    166:           data->cntl[3] = HPDDM_RECYCLE_TARGET_SM; /* default recycling target */
    167:         } else {
    168:           data->cntl[3] = 1; /* redistribution parameter of 1 by default */
    169:         }
    170:         data->cntl[4] = HPDDM_RECYCLE_STRATEGY_A; /* default recycling strategy */
    171:       }
    172:     } else data->scntl[1] = 1;
    173:   }
    174:   PetscCheck(ksp->nmax >= std::numeric_limits<int>::min() && ksp->nmax <= std::numeric_limits<int>::max(), PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "KSPMatSolve() block size %" PetscInt_FMT " not representable by an integer, which is not handled by KSPHPDDM",
    175:              ksp->nmax);
    176:   data->icntl[1] = static_cast<int>(ksp->nmax);
    177:   PetscFunctionReturn(PETSC_SUCCESS);
    178: }

    180: static inline PetscErrorCode KSPReset_HPDDM_Private(KSP ksp)
    181: {
    182:   KSP_HPDDM *data = (KSP_HPDDM *)ksp->data;

    184:   PetscFunctionBegin;
    185:   /* cast PETSC_DECIDE into the appropriate types to avoid compiler warnings */
    186:   std::fill_n(data->rcntl, PETSC_STATIC_ARRAY_LENGTH(data->rcntl), static_cast<PetscReal>(PETSC_DECIDE));
    187:   std::fill_n(data->icntl, PETSC_STATIC_ARRAY_LENGTH(data->icntl), static_cast<int>(PETSC_DECIDE));
    188:   std::fill_n(data->scntl, PETSC_STATIC_ARRAY_LENGTH(data->scntl), static_cast<unsigned short>(PETSC_DECIDE));
    189:   std::fill_n(data->cntl, PETSC_STATIC_ARRAY_LENGTH(data->cntl), static_cast<char>(PETSC_DECIDE));
    190:   data->precision = PETSC_SCALAR_PRECISION;
    191:   PetscFunctionReturn(PETSC_SUCCESS);
    192: }

    194: static PetscErrorCode KSPReset_HPDDM(KSP ksp)
    195: {
    196:   KSP_HPDDM *data = (KSP_HPDDM *)ksp->data;

    198:   PetscFunctionBegin;
    199:   delete data->op;
    200:   data->op = nullptr;
    201:   PetscCall(KSPReset_HPDDM_Private(ksp));
    202:   PetscFunctionReturn(PETSC_SUCCESS);
    203: }

    205: static PetscErrorCode KSPDestroy_HPDDM(KSP ksp)
    206: {
    207:   PetscFunctionBegin;
    208:   PetscCall(KSPReset_HPDDM(ksp));
    209:   PetscCall(KSPDestroyDefault(ksp));
    210:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetDeflationMat_C", nullptr));
    211:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetDeflationMat_C", nullptr));
    212:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetType_C", nullptr));
    213:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetType_C", nullptr));
    214:   PetscFunctionReturn(PETSC_SUCCESS);
    215: }

    217: template <PetscMemType type = PETSC_MEMTYPE_HOST>
    218: static inline PetscErrorCode KSPSolve_HPDDM_Private(KSP ksp, const PetscScalar *b, PetscScalar *x, PetscInt n)
    219: {
    220:   KSP_HPDDM              *data = (KSP_HPDDM *)ksp->data;
    221:   KSPConvergedDefaultCtx *ctx  = (KSPConvergedDefaultCtx *)ksp->cnvP;
    222:   const PetscInt          N    = data->op->getDof() * n;
    223:   PetscBool               flg;
    224: #if !PetscDefined(USE_REAL_DOUBLE) || PetscDefined(HAVE_F2CBLASLAPACK___FLOAT128_BINDINGS)
    225:   HPDDM::upscaled_type<PetscScalar> *high[2];
    226: #endif
    227: #if !PetscDefined(USE_REAL_SINGLE) || PetscDefined(HAVE_F2CBLASLAPACK___FP16_BINDINGS)
    228:   typedef HPDDM::downscaled_type<PetscReal> PetscDownscaledReal PETSC_ATTRIBUTE_MAY_ALIAS;
    229:   #if !PetscDefined(USE_COMPLEX)
    230:   PetscDownscaledReal *low[2];
    231:   #else
    232:   typedef PetscReal PetscAliasedReal   PETSC_ATTRIBUTE_MAY_ALIAS;
    233:   HPDDM::downscaled_type<PetscScalar> *low[2];
    234:   PetscAliasedReal                    *x_r;
    235:   PetscDownscaledReal                 *low_r;
    236:   #endif
    237: #endif
    238: #if PetscDefined(HAVE_CUDA)
    239:   Mat     A;
    240:   VecType vtype;
    241: #endif

    243:   PetscFunctionBegin;
    244: #if PetscDefined(HAVE_CUDA)
    245:   PetscCall(KSPGetOperators(ksp, &A, nullptr));
    246:   PetscCall(MatGetVecType(A, &vtype));
    247:   std::initializer_list<std::string>                 list = {VECCUDA, VECSEQCUDA, VECMPICUDA};
    248:   std::initializer_list<std::string>::const_iterator it   = std::find(list.begin(), list.end(), std::string(vtype));
    249:   PetscCheck(type != PETSC_MEMTYPE_HOST || it == list.end(), PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "MatGetVecType() must return a Vec with the same PetscMemType as the right-hand side and solution, PetscMemType(%s) != %s", vtype, PetscMemTypeToString(type));
    250: #endif
    251:   PetscCall(PCGetDiagonalScale(ksp->pc, &flg));
    252:   PetscCheck(!flg, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
    253:   if (n > 1) {
    254:     if (ksp->converged == KSPConvergedDefault) {
    255:       PetscCheck(!ctx->mininitialrtol, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support KSPConvergedDefaultSetUMIRNorm()", ((PetscObject)ksp)->type_name);
    256:       if (!ctx->initialrtol) {
    257:         PetscCall(PetscInfo(ksp, "Forcing KSPConvergedDefaultSetUIRNorm() since KSPConvergedDefault() cannot handle multiple norms\n"));
    258:         ctx->initialrtol = PETSC_TRUE;
    259:       }
    260:     } else PetscCall(PetscInfo(ksp, "Using a special \"converged\" callback, be careful, it is used in KSPHPDDM to track blocks of residuals\n"));
    261:   }
    262:   /* initial guess is always nonzero with recycling methods if there is a deflation subspace available */
    263:   if ((data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) && data->op->storage()) ksp->guess_zero = PETSC_FALSE;
    264:   ksp->its    = 0;
    265:   ksp->reason = KSP_CONVERGED_ITERATING;
    266:   if (data->precision > PETSC_SCALAR_PRECISION) { /* Krylov basis stored in higher precision than PetscScalar */
    267: #if !PetscDefined(USE_REAL_DOUBLE) || PetscDefined(HAVE_F2CBLASLAPACK___FLOAT128_BINDINGS)
    268:     if (type == PETSC_MEMTYPE_HOST) {
    269:       PetscCall(PetscMalloc2(N, high, N, high + 1));
    270:       HPDDM::copy_n(b, N, high[0]);
    271:       HPDDM::copy_n(x, N, high[1]);
    272:       PetscCall(HPDDM::IterativeMethod::solve(*data->op, high[0], high[1], n, PetscObjectComm((PetscObject)ksp)));
    273:       HPDDM::copy_n(high[1], N, x);
    274:       PetscCall(PetscFree2(high[0], high[1]));
    275:     } else {
    276:       PetscCheck(PetscDefined(HAVE_CUDA) && PetscDefined(USE_REAL_SINGLE), PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "CUDA in PETSc has no support for precisions other than single or double");
    277:   #if PetscDefined(HAVE_CUDA)
    278:     #if PetscDefined(HAVE_HPDDM)
    279:       PetscCall(KSPSolve_HPDDM_CUDA_Private(data, b, x, n, PetscObjectComm((PetscObject)ksp)));
    280:     #else
    281:       SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "No CUDA support with --download-hpddm from SLEPc");
    282:     #endif
    283:   #endif
    284:     }
    285: #else
    286:     PetscCheck(data->precision != PETSC_PRECISION___FLOAT128, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Reconfigure with --download-f2cblaslapack --with-f2cblaslapack-float128-bindings");
    287: #endif
    288:   } else if (data->precision < PETSC_SCALAR_PRECISION) { /* Krylov basis stored in lower precision than PetscScalar */
    289: #if !PetscDefined(USE_REAL_SINGLE) || PetscDefined(HAVE_F2CBLASLAPACK___FP16_BINDINGS)
    290:     if (type == PETSC_MEMTYPE_HOST) {
    291:       PetscCall(PetscMalloc1(N, low));
    292:   #if !PetscDefined(USE_COMPLEX)
    293:       low[1] = reinterpret_cast<PetscDownscaledReal *>(x);
    294:   #else
    295:       low[1] = reinterpret_cast<HPDDM::downscaled_type<PetscScalar> *>(x);
    296:   #endif
    297:       std::copy_n(b, N, low[0]);
    298:       for (PetscInt i = 0; i < N; ++i) low[1][i] = x[i];
    299:       PetscCall(HPDDM::IterativeMethod::solve(*data->op, low[0], low[1], n, PetscObjectComm((PetscObject)ksp)));
    300:   #if !PetscDefined(USE_COMPLEX)
    301:       for (PetscInt i = N; i-- > 0;) x[i] = static_cast<PetscScalar>(low[1][i]);
    302:   #else
    303:       x_r = reinterpret_cast<PetscAliasedReal *>(x), low_r = reinterpret_cast<PetscDownscaledReal *>(x_r);
    304:       for (PetscInt i = 2 * N; i-- > 0;) x_r[i] = static_cast<PetscReal>(low_r[i]);
    305:   #endif
    306:       PetscCall(PetscFree(low[0]));
    307:     } else {
    308:       PetscCheck(PetscDefined(HAVE_CUDA) && PetscDefined(USE_REAL_DOUBLE), PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "CUDA in PETSc has no support for precisions other than single or double");
    309:   #if PetscDefined(HAVE_CUDA)
    310:     #if PetscDefined(HAVE_HPDDM)
    311:       PetscCall(KSPSolve_HPDDM_CUDA_Private(data, b, x, n, PetscObjectComm((PetscObject)ksp)));
    312:     #else
    313:       SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "No CUDA support with --download-hpddm from SLEPc");
    314:     #endif
    315:   #endif
    316:     }
    317: #else
    318:     PetscCheck(data->precision != PETSC_PRECISION___FP16, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Reconfigure with --download-f2cblaslapack --with-f2cblaslapack-fp16-bindings");
    319: #endif
    320:   } else { /* Krylov basis stored in the same precision as PetscScalar */
    321:     if (type == PETSC_MEMTYPE_HOST) PetscCall(HPDDM::IterativeMethod::solve(*data->op, b, x, n, PetscObjectComm((PetscObject)ksp)));
    322:     else {
    323:       PetscCheck(PetscDefined(USE_REAL_SINGLE) || PetscDefined(USE_REAL_DOUBLE), PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "CUDA in PETSc has no support for precisions other than single or double");
    324: #if PetscDefined(HAVE_CUDA)
    325:   #if PetscDefined(HAVE_HPDDM)
    326:       PetscCall(KSPSolve_HPDDM_CUDA_Private(data, b, x, n, PetscObjectComm((PetscObject)ksp)));
    327:   #else
    328:       SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "No CUDA support with --download-hpddm from SLEPc");
    329:   #endif
    330: #endif
    331:     }
    332:   }
    333:   if (!ksp->reason) { /* KSPConvergedDefault() is still returning 0 (= KSP_CONVERGED_ITERATING) */
    334:     if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
    335:     else ksp->reason = KSP_CONVERGED_RTOL; /* early exit by HPDDM, which only happens on breakdowns or convergence */
    336:   }
    337:   ksp->its = PetscMin(ksp->its, ksp->max_it);
    338:   PetscFunctionReturn(PETSC_SUCCESS);
    339: }

    341: static PetscErrorCode KSPSolve_HPDDM(KSP ksp)
    342: {
    343:   KSP_HPDDM         *data = (KSP_HPDDM *)ksp->data;
    344:   Mat                A, B;
    345:   PetscScalar       *x, *bt = nullptr, **ptr;
    346:   const PetscScalar *b;
    347:   PetscInt           i, j, n;
    348:   PetscBool          flg;
    349:   PetscMemType       type[2];

    351:   PetscFunctionBegin;
    352:   PetscCall(PetscCitationsRegister(HPDDMCitation, &HPDDMCite));
    353:   PetscCall(KSPGetOperators(ksp, &A, nullptr));
    354:   PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &flg, MATSEQKAIJ, MATMPIKAIJ, ""));
    355:   PetscCall(VecGetArrayWriteAndMemType(ksp->vec_sol, &x, type));
    356:   PetscCall(VecGetArrayReadAndMemType(ksp->vec_rhs, &b, type + 1));
    357:   PetscCheck(type[0] == type[1], PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Right-hand side and solution vectors must have the same PetscMemType, %s != %s", PetscMemTypeToString(type[0]), PetscMemTypeToString(type[1]));
    358:   if (!flg) {
    359:     if (PetscMemTypeCUDA(type[0])) PetscCall(KSPSolve_HPDDM_Private<PETSC_MEMTYPE_CUDA>(ksp, b, x, 1));
    360:     else {
    361:       PetscCheck(PetscMemTypeHost(type[0]), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "PetscMemType (%s) is neither PETSC_MEMTYPE_HOST nor PETSC_MEMTYPE_CUDA", PetscMemTypeToString(type[0]));
    362:       PetscCall(KSPSolve_HPDDM_Private(ksp, b, x, 1));
    363:     }
    364:   } else {
    365:     PetscCheck(PetscMemTypeHost(type[0]), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "PetscMemType (%s) is not PETSC_MEMTYPE_HOST", PetscMemTypeToString(type[0]));
    366:     PetscCall(MatKAIJGetScaledIdentity(A, &flg));
    367:     PetscCall(MatKAIJGetAIJ(A, &B));
    368:     PetscCall(MatGetBlockSize(A, &n));
    369:     PetscCall(MatGetLocalSize(B, &i, nullptr));
    370:     j = data->op->getDof();
    371:     if (!flg) i *= n; /* S and T are not scaled identities, cannot use block methods */
    372:     if (i != j) {     /* switching between block and standard methods */
    373:       delete data->op;
    374:       data->op = new HPDDM::PETScOperator(ksp, i);
    375:     }
    376:     if (flg && n > 1) {
    377:       PetscCall(PetscMalloc1(i * n, &bt));
    378:       /* from row- to column-major to be consistent with HPDDM */
    379:       HPDDM::Wrapper<PetscScalar>::omatcopy<'T'>(i, n, b, n, bt, i);
    380:       ptr = const_cast<PetscScalar **>(&b);
    381:       std::swap(*ptr, bt);
    382:       HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(i, n, x, n, i);
    383:     }
    384:     PetscCall(KSPSolve_HPDDM_Private(ksp, b, x, flg ? n : 1));
    385:     if (flg && n > 1) {
    386:       std::swap(*ptr, bt);
    387:       PetscCall(PetscFree(bt));
    388:       /* from column- to row-major to be consistent with MatKAIJ format */
    389:       HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(n, i, x, i, n);
    390:     }
    391:   }
    392:   PetscCall(VecRestoreArrayReadAndMemType(ksp->vec_rhs, &b));
    393:   PetscCall(VecRestoreArrayWriteAndMemType(ksp->vec_sol, &x));
    394:   PetscFunctionReturn(PETSC_SUCCESS);
    395: }

    397: /*@
    398:   KSPHPDDMSetDeflationMat - Sets the deflation space used by Krylov methods in `KSPHPDDM` with recycling. This space is viewed as a set of vectors stored in
    399:   a `MATDENSE` (column major).

    401:   Input Parameters:
    402: + ksp - iterative context
    403: - U   - deflation space to be used during `KSPSolve()`

    405:   Level: intermediate

    407: .seealso: [](ch_ksp), `KSPHPDDM`, `KSPCreate()`, `KSPType`, `KSPHPDDMGetDeflationMat()`
    408: @*/
    409: PetscErrorCode KSPHPDDMSetDeflationMat(KSP ksp, Mat U)
    410: {
    411:   PetscFunctionBegin;
    414:   PetscCheckSameComm(ksp, 1, U, 2);
    415:   PetscUseMethod(ksp, "KSPHPDDMSetDeflationMat_C", (KSP, Mat), (ksp, U));
    416:   PetscFunctionReturn(PETSC_SUCCESS);
    417: }

    419: /*@
    420:   KSPHPDDMGetDeflationMat - Gets the deflation space computed by Krylov methods in `KSPHPDDM`  with recycling or `NULL` if `KSPSolve()` has not been called yet.

    422:   Input Parameter:
    423: . ksp - iterative context

    425:   Output Parameter:
    426: . U - deflation space generated during `KSPSolve()`

    428:   Level: intermediate

    430:   Note:
    431:   This space is viewed as a set of vectors stored in a `MATDENSE` (column major). It is the responsibility of the user to free the returned `Mat`.

    433: .seealso: [](ch_ksp), `KSPHPDDM`, `KSPCreate()`, `KSPType`, `KSPHPDDMSetDeflationMat()`
    434: @*/
    435: PetscErrorCode KSPHPDDMGetDeflationMat(KSP ksp, Mat *U)
    436: {
    437:   PetscFunctionBegin;
    439:   if (U) {
    440:     PetscAssertPointer(U, 2);
    441:     PetscUseMethod(ksp, "KSPHPDDMGetDeflationMat_C", (KSP, Mat *), (ksp, U));
    442:   }
    443:   PetscFunctionReturn(PETSC_SUCCESS);
    444: }

    446: static PetscErrorCode KSPHPDDMSetDeflationMat_HPDDM(KSP ksp, Mat U)
    447: {
    448:   KSP_HPDDM            *data = (KSP_HPDDM *)ksp->data;
    449:   HPDDM::PETScOperator *op   = data->op;
    450:   Mat                   A;
    451:   const PetscScalar    *array;
    452:   PetscScalar          *copy;
    453:   PetscInt              m1, M1, m2, M2, n2, N2, ldu;
    454:   PetscBool             match;

    456:   PetscFunctionBegin;
    457:   if (!op) {
    458:     PetscCall(KSPSetUp(ksp));
    459:     op = data->op;
    460:   }
    461:   PetscCheck(data->precision == PETSC_SCALAR_PRECISION, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s != %s", PetscPrecisionTypes[data->precision], PetscPrecisionTypes[PETSC_SCALAR_PRECISION]);
    462:   PetscCall(KSPGetOperators(ksp, &A, nullptr));
    463:   PetscCall(MatGetLocalSize(A, &m1, nullptr));
    464:   PetscCall(MatGetLocalSize(U, &m2, &n2));
    465:   PetscCall(MatGetSize(A, &M1, nullptr));
    466:   PetscCall(MatGetSize(U, &M2, &N2));
    467:   PetscCheck(m1 == m2 && M1 == M2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot use a deflation space with (m2,M2) = (%" PetscInt_FMT ",%" PetscInt_FMT ") for a linear system with (m1,M1) = (%" PetscInt_FMT ",%" PetscInt_FMT ")", m2, M2, m1, M1);
    468:   PetscCall(PetscObjectTypeCompareAny((PetscObject)U, &match, MATSEQDENSE, MATMPIDENSE, ""));
    469:   PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided deflation space not stored in a dense Mat");
    470:   PetscCall(MatDenseGetArrayRead(U, &array));
    471:   copy = op->allocate(m2, 1, N2);
    472:   PetscCheck(copy, PETSC_COMM_SELF, PETSC_ERR_POINTER, "Memory allocation error");
    473:   PetscCall(MatDenseGetLDA(U, &ldu));
    474:   HPDDM::Wrapper<PetscScalar>::omatcopy<'N'>(N2, m2, array, ldu, copy, m2);
    475:   PetscCall(MatDenseRestoreArrayRead(U, &array));
    476:   PetscFunctionReturn(PETSC_SUCCESS);
    477: }

    479: static PetscErrorCode KSPHPDDMGetDeflationMat_HPDDM(KSP ksp, Mat *U)
    480: {
    481:   KSP_HPDDM            *data = (KSP_HPDDM *)ksp->data;
    482:   HPDDM::PETScOperator *op   = data->op;
    483:   Mat                   A;
    484:   const PetscScalar    *array;
    485:   PetscScalar          *copy;
    486:   PetscInt              m1, M1, N2;

    488:   PetscFunctionBegin;
    489:   if (!op) {
    490:     PetscCall(KSPSetUp(ksp));
    491:     op = data->op;
    492:   }
    493:   PetscCheck(data->precision == PETSC_SCALAR_PRECISION, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s != %s", PetscPrecisionTypes[data->precision], PetscPrecisionTypes[PETSC_SCALAR_PRECISION]);
    494:   array = op->storage();
    495:   N2    = op->k().first * op->k().second;
    496:   if (!array) *U = nullptr;
    497:   else {
    498:     PetscCall(KSPGetOperators(ksp, &A, nullptr));
    499:     PetscCall(MatGetLocalSize(A, &m1, nullptr));
    500:     PetscCall(MatGetSize(A, &M1, nullptr));
    501:     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)ksp), m1, PETSC_DECIDE, M1, N2, nullptr, U));
    502:     PetscCall(MatDenseGetArrayWrite(*U, &copy));
    503:     PetscCall(PetscArraycpy(copy, array, m1 * N2));
    504:     PetscCall(MatDenseRestoreArrayWrite(*U, &copy));
    505:   }
    506:   PetscFunctionReturn(PETSC_SUCCESS);
    507: }

    509: static PetscErrorCode KSPMatSolve_HPDDM(KSP ksp, Mat B, Mat X)
    510: {
    511:   KSP_HPDDM         *data = (KSP_HPDDM *)ksp->data;
    512:   Mat                A;
    513:   const PetscScalar *b;
    514:   PetscScalar       *x;
    515:   PetscInt           n, lda;
    516:   PetscMemType       type[2];

    518:   PetscFunctionBegin;
    519:   PetscCall(PetscCitationsRegister(HPDDMCitation, &HPDDMCite));
    520:   if (!data->op) PetscCall(KSPSetUp(ksp));
    521:   PetscCall(KSPGetOperators(ksp, &A, nullptr));
    522:   PetscCall(MatGetLocalSize(B, &n, nullptr));
    523:   PetscCall(MatDenseGetLDA(B, &lda));
    524:   PetscCheck(n == lda, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Unsupported leading dimension lda = %" PetscInt_FMT " with n = %" PetscInt_FMT, lda, n);
    525:   PetscCall(MatGetLocalSize(A, &n, nullptr));
    526:   PetscCall(MatDenseGetLDA(X, &lda));
    527:   PetscCheck(n == lda, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Unsupported leading dimension lda = %" PetscInt_FMT " with n = %" PetscInt_FMT, lda, n);
    528:   PetscCall(MatGetSize(X, nullptr, &n));
    529:   PetscCall(MatDenseGetArrayWriteAndMemType(X, &x, type));
    530:   PetscCall(MatDenseGetArrayReadAndMemType(B, &b, type + 1));
    531:   PetscCheck(type[0] == type[1], PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Right-hand side and solution matrices must have the same PetscMemType, %s != %s", PetscMemTypeToString(type[0]), PetscMemTypeToString(type[1]));
    532:   if (PetscMemTypeCUDA(type[0])) PetscCall(KSPSolve_HPDDM_Private<PETSC_MEMTYPE_CUDA>(ksp, b, x, n));
    533:   else {
    534:     PetscCheck(PetscMemTypeHost(type[0]), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "PetscMemType (%s) is neither PETSC_MEMTYPE_HOST nor PETSC_MEMTYPE_CUDA", PetscMemTypeToString(type[0]));
    535:     PetscCall(KSPSolve_HPDDM_Private(ksp, b, x, n));
    536:   }
    537:   PetscCall(MatDenseRestoreArrayReadAndMemType(B, &b));
    538:   PetscCall(MatDenseRestoreArrayWriteAndMemType(X, &x));
    539:   PetscFunctionReturn(PETSC_SUCCESS);
    540: }

    542: /*@
    543:   KSPHPDDMSetType - Sets the type of Krylov method used in `KSPHPDDM`.

    545:   Collective

    547:   Input Parameters:
    548: + ksp  - iterative context
    549: - type - any of gmres, bgmres, cg, bcg, gcrodr, bgcrodr, bfbcg, or preonly

    551:   Level: intermediate

    553:   Notes:
    554:   Unlike `KSPReset()`, this function does not destroy any deflation space attached to the `KSP`.

    556:   As an example, in the following sequence\:
    557: .vb
    558:      KSPHPDDMSetType(ksp, KSPGCRODR);
    559:      KSPSolve(ksp, b, x);
    560:      KSPHPDDMSetType(ksp, KSPGMRES);
    561:      KSPHPDDMSetType(ksp, KSPGCRODR);
    562:      KSPSolve(ksp, b, x);
    563: .ve
    564:   the recycled space is reused in the second `KSPSolve()`.

    566: .seealso: [](ch_ksp), `KSPCreate()`, `KSPType`, `KSPHPDDMType`, `KSPHPDDMGetType()`
    567: @*/
    568: PetscErrorCode KSPHPDDMSetType(KSP ksp, KSPHPDDMType type)
    569: {
    570:   PetscFunctionBegin;
    573:   PetscUseMethod(ksp, "KSPHPDDMSetType_C", (KSP, KSPHPDDMType), (ksp, type));
    574:   PetscFunctionReturn(PETSC_SUCCESS);
    575: }

    577: /*@
    578:   KSPHPDDMGetType - Gets the type of Krylov method used in `KSPHPDDM`.

    580:   Input Parameter:
    581: . ksp - iterative context

    583:   Output Parameter:
    584: . type - any of gmres, bgmres, cg, bcg, gcrodr, bgcrodr, bfbcg, or preonly

    586:   Level: intermediate

    588: .seealso: [](ch_ksp), `KSPCreate()`, `KSPType`, `KSPHPDDMType`, `KSPHPDDMSetType()`
    589: @*/
    590: PetscErrorCode KSPHPDDMGetType(KSP ksp, KSPHPDDMType *type)
    591: {
    592:   PetscFunctionBegin;
    594:   if (type) {
    595:     PetscAssertPointer(type, 2);
    596:     PetscUseMethod(ksp, "KSPHPDDMGetType_C", (KSP, KSPHPDDMType *), (ksp, type));
    597:   }
    598:   PetscFunctionReturn(PETSC_SUCCESS);
    599: }

    601: static PetscErrorCode KSPHPDDMSetType_HPDDM(KSP ksp, KSPHPDDMType type)
    602: {
    603:   KSP_HPDDM *data = (KSP_HPDDM *)ksp->data;
    604:   PetscInt   i;
    605:   PetscBool  flg = PETSC_FALSE;

    607:   PetscFunctionBegin;
    608:   for (i = 0; i < static_cast<PetscInt>(PETSC_STATIC_ARRAY_LENGTH(KSPHPDDMTypes)); ++i) {
    609:     PetscCall(PetscStrcmp(KSPHPDDMTypes[type], KSPHPDDMTypes[i], &flg));
    610:     if (flg) break;
    611:   }
    612:   PetscCheck(i != PETSC_STATIC_ARRAY_LENGTH(KSPHPDDMTypes), PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown KSPHPDDMType %d", type);
    613:   if (data->cntl[0] != static_cast<char>(PETSC_DECIDE) && data->cntl[0] != i) PetscCall(KSPReset_HPDDM_Private(ksp));
    614:   data->cntl[0] = i;
    615:   PetscFunctionReturn(PETSC_SUCCESS);
    616: }

    618: static PetscErrorCode KSPHPDDMGetType_HPDDM(KSP ksp, KSPHPDDMType *type)
    619: {
    620:   KSP_HPDDM *data = (KSP_HPDDM *)ksp->data;

    622:   PetscFunctionBegin;
    623:   PetscCheck(data->cntl[0] != static_cast<char>(PETSC_DECIDE), PETSC_COMM_SELF, PETSC_ERR_ORDER, "KSPHPDDMType not set yet");
    624:   /* need to shift by -1 for HPDDM_KRYLOV_METHOD_NONE */
    625:   *type = static_cast<KSPHPDDMType>(PetscMin(data->cntl[0], static_cast<char>(PETSC_STATIC_ARRAY_LENGTH(KSPHPDDMTypes) - 1)));
    626:   PetscFunctionReturn(PETSC_SUCCESS);
    627: }

    629: /*MC
    630:    KSPHPDDM - Interface with the HPDDM library. This `KSP` may be used to further select methods that are currently not implemented natively in PETSc, e.g.,
    631:    GCRODR {cite}`parks2006recycling`, a recycled Krylov method which is similar to `KSPLGMRES`, see {cite}`jolivet2016block` for a comparison.
    632:    ex75.c shows how to reproduce the results
    633:    from the aforementioned paper {cite}`parks2006recycling`. A chronological bibliography of relevant publications linked with `KSP` available in HPDDM through `KSPHPDDM`,
    634:    and not available directly in PETSc, may be found below. The interface is explained in details in {cite}`jolivetromanzampini2020`.
    635:    See also {cite}`o1980block`, {cite}`ji2017breakdown` and {cite}`calandra2013modified`

    637:    Options Database Keys:
    638: +   -ksp_gmres_restart restart                             - see `KSPGMRES`, default is 30
    639: .   -ksp_hpddm_type type                                   - see `KSPHPDDMType`, default is `gmres`
    640: .   -ksp_hpddm_precision (__fp16|single|double|__float128) - see `PetscPrecision`, default is the same as `PetscScalar`
    641: .   -ksp_hpddm_deflation_tol eps                           - tolerance when deflating right-hand sides inside block methods (only relevant with block methods), default is -1.0
    642: .   -ksp_hpddm_enlarge_krylov_subspace p                   - split the initial right-hand side into multiple vectors (only relevant with nonblock methods), default is 1
    643: .   -ksp_hpddm_orthogonalization (cgs|mgs)                 - see `KSPGMRES`, default is `cgs`
    644: .   -ksp_hpddm_qr (cholqr|cgs|mgs)                         - distributed QR factorizations, only relevant with block methods, default is `cholqr`
    645: .   -ksp_hpddm_variant (left|right|flexible)               - this option is superseded by `KSPSetPCSide()`, default is `left`
    646: .   -ksp_hpddm_recycle n                                   - number of harmonic Ritz vectors to compute (only relevant with GCRODR or BGCRODR), default is 0
    647: .   -ksp_hpddm_recycle_target (SM|LM|SR|LR|SI|LI)          - criterion to select harmonic Ritz vectors, default is `SM`
    648:                                                              (only relevant with GCRODR or BGCRODR). For BGCRODR, if PETSc is compiled with SLEPc,
    649:                                                              this option is not relevant, since SLEPc is used instead. Options are set with the prefix `-ksp_hpddm_recycle_eps_`
    650: .   -ksp_hpddm_recycle_strategy (A|B)                      - generalized eigenvalue problem to solve for recycling (only relevant with flexible GCRODR or BGCRODR), default is `A`
    651: -   -ksp_hpddm_recycle_symmetric (true|false)              - symmetric generalized eigenproblems in BGCRODR, useful to switch to distributed solvers like
    652:                                                              `EPSELEMENTAL` or `EPSSCALAPACK` (only relevant when PETSc is compiled with SLEPc), default is `false`

    654:    Level: intermediate

    656: .seealso: [](ch_ksp), [](sec_flexibleksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPGMRES`, `KSPCG`, `KSPLGMRES`, `KSPDGMRES`
    657: M*/

    659: PETSC_EXTERN PetscErrorCode KSPCreate_HPDDM(KSP ksp)
    660: {
    661:   KSP_HPDDM  *data;
    662:   PetscInt    i;
    663:   const char *common[] = {KSPGMRES, KSPCG, KSPPREONLY};
    664:   PetscBool   flg      = PETSC_FALSE;

    666:   PetscFunctionBegin;
    667:   PetscCall(PetscNew(&data));
    668:   ksp->data = (void *)data;
    669:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
    670:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 1));
    671:   ksp->ops->solve          = KSPSolve_HPDDM;
    672:   ksp->ops->matsolve       = KSPMatSolve_HPDDM;
    673:   ksp->ops->setup          = KSPSetUp_HPDDM;
    674:   ksp->ops->setfromoptions = KSPSetFromOptions_HPDDM;
    675:   ksp->ops->destroy        = KSPDestroy_HPDDM;
    676:   ksp->ops->view           = KSPView_HPDDM;
    677:   ksp->ops->reset          = KSPReset_HPDDM;
    678:   PetscCall(KSPReset_HPDDM_Private(ksp));
    679:   for (i = 0; i < static_cast<PetscInt>(PETSC_STATIC_ARRAY_LENGTH(common)); ++i) {
    680:     PetscCall(PetscStrcmp(((PetscObject)ksp)->type_name, common[i], &flg));
    681:     if (flg) break;
    682:   }
    683:   if (!i) data->cntl[0] = HPDDM_KRYLOV_METHOD_GMRES;
    684:   else if (i == 1) data->cntl[0] = HPDDM_KRYLOV_METHOD_CG;
    685:   else if (i == 2) data->cntl[0] = HPDDM_KRYLOV_METHOD_NONE;
    686:   if (data->cntl[0] != static_cast<char>(PETSC_DECIDE)) PetscCall(PetscInfo(ksp, "Using the previously set KSPType %s\n", common[i]));
    687:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetDeflationMat_C", KSPHPDDMSetDeflationMat_HPDDM));
    688:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetDeflationMat_C", KSPHPDDMGetDeflationMat_HPDDM));
    689:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetType_C", KSPHPDDMSetType_HPDDM));
    690:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetType_C", KSPHPDDMGetType_HPDDM));
    691: #if PetscDefined(HAVE_SLEPC) && PetscDefined(HAVE_DYNAMIC_LIBRARIES) && PetscDefined(USE_SHARED_LIBRARIES)
    692:   if (!loadedDL) PetscCall(HPDDMLoadDL_Private(&loadedDL));
    693: #endif
    694:   data->precision = PETSC_SCALAR_PRECISION;
    695:   PetscFunctionReturn(PETSC_SUCCESS);
    696: }