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 KSPHPDDMPrecisionTypes[] = {"HALF", "SINGLE", "DOUBLE", "QUADRUPLE", "KSPHPDDMPrecisionType", "KSP_HPDDM_PRECISION_", NULL};
  6: const char *const HPDDMOrthogonalization[] = {"cgs", "mgs"};
  7: const char *const HPDDMQR[]                = {"cholqr", "cgs", "mgs"};
  8: const char *const HPDDMVariant[]           = {"left", "right", "flexible"};
  9: const char *const HPDDMRecycleTarget[]     = {"SM", "LM", "SR", "LR", "SI", "LI"};
 10: const char *const HPDDMRecycleStrategy[]   = {"A", "B"};

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

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

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

 34:   PetscFunctionBegin;
 35:   PetscOptionsHeadBegin(PetscOptionsObject, "KSPHPDDM options, cf. https://github.com/hpddm/hpddm");
 36:   i = (data->cntl[0] == static_cast<char>(PETSC_DECIDE) ? HPDDM_KRYLOV_METHOD_GMRES : data->cntl[0]);
 37:   PetscCall(PetscOptionsEList("-ksp_hpddm_type", "Type of Krylov method", "KSPHPDDMGetType", KSPHPDDMTypes, PETSC_STATIC_ARRAY_LENGTH(KSPHPDDMTypes), KSPHPDDMTypes[HPDDM_KRYLOV_METHOD_GMRES], &i, NULL));
 38:   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 */
 39:   data->cntl[0] = i;
 40:   PetscCall(PetscOptionsEnum("-ksp_hpddm_precision", "Precision in which Krylov bases are stored", "KSPHPDDM", KSPHPDDMPrecisionTypes, (PetscEnum)data->precision, (PetscEnum *)&data->precision, NULL));
 41:   PetscCheck(data->precision != KSP_HPDDM_PRECISION_QUADRUPLE || PetscDefined(HAVE_REAL___FLOAT128), PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP_SYS, "Unsupported %s precision", KSPHPDDMPrecisionTypes[data->precision]);
 42:   PetscCheck(std::abs(data->precision - PETSC_KSPHPDDM_DEFAULT_PRECISION) <= 1, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Unhandled mixed %s and %s precisions", KSPHPDDMPrecisionTypes[data->precision], KSPHPDDMPrecisionTypes[PETSC_KSPHPDDM_DEFAULT_PRECISION]);
 43:   if (data->cntl[0] != HPDDM_KRYLOV_METHOD_NONE) {
 44:     if (data->cntl[0] != HPDDM_KRYLOV_METHOD_BCG && data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG) {
 45:       i = (data->cntl[1] == static_cast<char>(PETSC_DECIDE) ? HPDDM_VARIANT_LEFT : data->cntl[1]);
 46:       if (ksp->pc_side_set == PC_SIDE_DEFAULT)
 47:         PetscCall(PetscOptionsEList("-ksp_hpddm_variant", "Left, right, or variable preconditioning", "KSPHPDDM", HPDDMVariant, PETSC_STATIC_ARRAY_LENGTH(HPDDMVariant), HPDDMVariant[HPDDM_VARIANT_LEFT], &i, NULL));
 48:       else if (ksp->pc_side_set == PC_RIGHT) i = HPDDM_VARIANT_RIGHT;
 49:       data->cntl[1] = i;
 50:       if (i > 0) PetscCall(KSPSetPCSide(ksp, PC_RIGHT));
 51:     }
 52:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
 53:       data->rcntl[0] = (PetscAbsReal(data->rcntl[0] - static_cast<PetscReal>(PETSC_DECIDE)) < PETSC_SMALL ? -1.0 : data->rcntl[0]);
 54:       PetscCall(PetscOptionsReal("-ksp_hpddm_deflation_tol", "Tolerance when deflating right-hand sides inside block methods", "KSPHPDDM", data->rcntl[0], data->rcntl, NULL));
 55:       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]));
 56:       PetscCall(PetscOptionsRangeInt("-ksp_hpddm_enlarge_krylov_subspace", "Split the initial right-hand side into multiple vectors", "KSPHPDDM", i, &i, NULL, 1, std::numeric_limits<unsigned short>::max() - 1));
 57:       data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG] = i;
 58:     } else data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BCG] = 0;
 59:     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) {
 60:       i = (data->cntl[2] == static_cast<char>(PETSC_DECIDE) ? HPDDM_ORTHOGONALIZATION_CGS : data->cntl[2] & 3);
 61:       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, NULL));
 62:       j = (data->cntl[2] == static_cast<char>(PETSC_DECIDE) ? HPDDM_QR_CHOLQR : ((data->cntl[2] >> 2) & 7));
 63:       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, NULL));
 64:       data->cntl[2] = static_cast<char>(i) + (static_cast<char>(j) << 2);
 65:       i             = (data->scntl[0] == static_cast<unsigned short>(PETSC_DECIDE) ? PetscMin(30, ksp->max_it) : data->scntl[0]);
 66:       PetscCall(PetscOptionsRangeInt("-ksp_gmres_restart", "Maximum number of Arnoldi vectors generated per cycle", "KSPHPDDM", i, &i, NULL, PetscMin(1, ksp->max_it), PetscMin(ksp->max_it, std::numeric_limits<unsigned short>::max() - 1)));
 67:       data->scntl[0] = i;
 68:     }
 69:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BCG || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
 70:       j = (data->cntl[1] == static_cast<char>(PETSC_DECIDE) ? HPDDM_QR_CHOLQR : data->cntl[1]);
 71:       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, NULL));
 72:       data->cntl[1] = j;
 73:     }
 74:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
 75:       i = (data->icntl[0] == static_cast<int>(PETSC_DECIDE) ? PetscMin(20, data->scntl[0] - 1) : data->icntl[0]);
 76:       PetscCall(PetscOptionsRangeInt("-ksp_hpddm_recycle", "Number of harmonic Ritz vectors to compute", "KSPHPDDM", i, &i, NULL, 1, data->scntl[0] - 1));
 77:       data->icntl[0] = i;
 78:       if (!PetscDefined(HAVE_SLEPC) || !PetscDefined(USE_SHARED_LIBRARIES) || data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR) {
 79:         i = (data->cntl[3] == static_cast<char>(PETSC_DECIDE) ? HPDDM_RECYCLE_TARGET_SM : data->cntl[3]);
 80:         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, NULL));
 81:         data->cntl[3] = i;
 82:       } else {
 83:         PetscCheck(data->precision == PETSC_KSPHPDDM_DEFAULT_PRECISION, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "Cannot use SLEPc with a different precision than PETSc for harmonic Ritz eigensolves");
 84:         PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ksp), &size));
 85:         i = (data->cntl[3] == static_cast<char>(PETSC_DECIDE) ? 1 : data->cntl[3]);
 86:         PetscCall(PetscOptionsRangeInt("-ksp_hpddm_recycle_redistribute", "Number of processes used to solve eigenvalue problems when recycling in BGCRODR", "KSPHPDDM", i, &i, NULL, 1, PetscMin(size, 192)));
 87:         data->cntl[3] = i;
 88:       }
 89:       i = (data->cntl[4] == static_cast<char>(PETSC_DECIDE) ? HPDDM_RECYCLE_STRATEGY_A : data->cntl[4]);
 90:       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, NULL));
 91:       data->cntl[4] = i;
 92:     }
 93:   } else {
 94:     data->cntl[0]  = HPDDM_KRYLOV_METHOD_NONE;
 95:     data->scntl[1] = 1;
 96:   }
 97:   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",
 98:              ksp->nmax);
 99:   data->icntl[1] = static_cast<int>(ksp->nmax);
100:   PetscOptionsHeadEnd();
101:   PetscFunctionReturn(PETSC_SUCCESS);
102: }

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

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

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

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

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

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

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

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

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

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

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

327: static PetscErrorCode KSPSolve_HPDDM(KSP ksp)
328: {
329:   KSP_HPDDM         *data = (KSP_HPDDM *)ksp->data;
330:   Mat                A, B;
331:   PetscScalar       *x, *bt = NULL, **ptr;
332:   const PetscScalar *b;
333:   PetscInt           i, j, n;
334:   PetscBool          flg;
335:   PetscMemType       type[2];

337:   PetscFunctionBegin;
338:   PetscCall(PetscCitationsRegister(HPDDMCitation, &HPDDMCite));
339:   PetscCall(KSPGetOperators(ksp, &A, NULL));
340:   PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &flg, MATSEQKAIJ, MATMPIKAIJ, ""));
341:   PetscCall(VecGetArrayWriteAndMemType(ksp->vec_sol, &x, type));
342:   PetscCall(VecGetArrayReadAndMemType(ksp->vec_rhs, &b, type + 1));
343:   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]));
344:   if (!flg) {
345:     if (PetscMemTypeCUDA(type[0])) PetscCall(KSPSolve_HPDDM_Private<PETSC_MEMTYPE_CUDA>(ksp, b, x, 1));
346:     else {
347:       PetscCheck(PetscMemTypeHost(type[0]), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "PetscMemType (%s) is neither PETSC_MEMTYPE_HOST nor PETSC_MEMTYPE_CUDA", PetscMemTypeToString(type[0]));
348:       PetscCall(KSPSolve_HPDDM_Private(ksp, b, x, 1));
349:     }
350:   } else {
351:     PetscCheck(PetscMemTypeHost(type[0]), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "PetscMemType (%s) is not PETSC_MEMTYPE_HOST", PetscMemTypeToString(type[0]));
352:     PetscCall(MatKAIJGetScaledIdentity(A, &flg));
353:     PetscCall(MatKAIJGetAIJ(A, &B));
354:     PetscCall(MatGetBlockSize(A, &n));
355:     PetscCall(MatGetLocalSize(B, &i, NULL));
356:     j = data->op->getDof();
357:     if (!flg) i *= n; /* S and T are not scaled identities, cannot use block methods */
358:     if (i != j) {     /* switching between block and standard methods */
359:       delete data->op;
360:       data->op = new HPDDM::PETScOperator(ksp, i);
361:     }
362:     if (flg && n > 1) {
363:       PetscCall(PetscMalloc1(i * n, &bt));
364:       /* from row- to column-major to be consistent with HPDDM */
365:       HPDDM::Wrapper<PetscScalar>::omatcopy<'T'>(i, n, b, n, bt, i);
366:       ptr = const_cast<PetscScalar **>(&b);
367:       std::swap(*ptr, bt);
368:       HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(i, n, x, n, i);
369:     }
370:     PetscCall(KSPSolve_HPDDM_Private(ksp, b, x, flg ? n : 1));
371:     if (flg && n > 1) {
372:       std::swap(*ptr, bt);
373:       PetscCall(PetscFree(bt));
374:       /* from column- to row-major to be consistent with MatKAIJ format */
375:       HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(n, i, x, i, n);
376:     }
377:   }
378:   PetscCall(VecRestoreArrayReadAndMemType(ksp->vec_rhs, &b));
379:   PetscCall(VecRestoreArrayWriteAndMemType(ksp->vec_sol, &x));
380:   PetscFunctionReturn(PETSC_SUCCESS);
381: }

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

387:    Input Parameters:
388: +     ksp - iterative context
389: -     U - deflation space to be used during KSPSolve()

391:    Level: intermediate

393: .seealso: [](chapter_ksp), `KSPHPDDM`, `KSPCreate()`, `KSPType`, `KSPHPDDMGetDeflationMat()`
394: @*/
395: PetscErrorCode KSPHPDDMSetDeflationMat(KSP ksp, Mat U)
396: {
397:   PetscFunctionBegin;
400:   PetscCheckSameComm(ksp, 1, U, 2);
401:   PetscUseMethod(ksp, "KSPHPDDMSetDeflationMat_C", (KSP, Mat), (ksp, U));
402:   PetscFunctionReturn(PETSC_SUCCESS);
403: }

405: /*@
406:      KSPHPDDMGetDeflationMat - Gets the deflation space computed by Krylov methods in `KSPHPDDM`  with recycling or NULL if `KSPSolve()` has not been called yet.
407:      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`.

409:    Input Parameter:
410: .     ksp - iterative context

412:    Output Parameter:
413: .     U - deflation space generated during `KSPSolve()`

415:    Level: intermediate

417: .seealso: [](chapter_ksp), `KSPHPDDM`, `KSPCreate()`, `KSPType`, `KSPHPDDMSetDeflationMat()`
418: @*/
419: PetscErrorCode KSPHPDDMGetDeflationMat(KSP ksp, Mat *U)
420: {
421:   PetscFunctionBegin;
423:   if (U) {
425:     PetscUseMethod(ksp, "KSPHPDDMGetDeflationMat_C", (KSP, Mat *), (ksp, U));
426:   }
427:   PetscFunctionReturn(PETSC_SUCCESS);
428: }

430: static PetscErrorCode KSPHPDDMSetDeflationMat_HPDDM(KSP ksp, Mat U)
431: {
432:   KSP_HPDDM            *data = (KSP_HPDDM *)ksp->data;
433:   HPDDM::PETScOperator *op   = data->op;
434:   Mat                   A;
435:   const PetscScalar    *array;
436:   PetscScalar          *copy;
437:   PetscInt              m1, M1, m2, M2, n2, N2, ldu;
438:   PetscBool             match;

440:   PetscFunctionBegin;
441:   if (!op) {
442:     PetscCall(KSPSetUp(ksp));
443:     op = data->op;
444:   }
445:   PetscCheck(data->precision == PETSC_KSPHPDDM_DEFAULT_PRECISION, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s != %s", KSPHPDDMPrecisionTypes[data->precision], KSPHPDDMPrecisionTypes[PETSC_KSPHPDDM_DEFAULT_PRECISION]);
446:   PetscCall(KSPGetOperators(ksp, &A, NULL));
447:   PetscCall(MatGetLocalSize(A, &m1, NULL));
448:   PetscCall(MatGetLocalSize(U, &m2, &n2));
449:   PetscCall(MatGetSize(A, &M1, NULL));
450:   PetscCall(MatGetSize(U, &M2, &N2));
451:   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);
452:   PetscCall(PetscObjectTypeCompareAny((PetscObject)U, &match, MATSEQDENSE, MATMPIDENSE, ""));
453:   PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided deflation space not stored in a dense Mat");
454:   PetscCall(MatDenseGetArrayRead(U, &array));
455:   copy = op->allocate(m2, 1, N2);
456:   PetscCheck(copy, PETSC_COMM_SELF, PETSC_ERR_POINTER, "Memory allocation error");
457:   PetscCall(MatDenseGetLDA(U, &ldu));
458:   HPDDM::Wrapper<PetscScalar>::omatcopy<'N'>(N2, m2, array, ldu, copy, m2);
459:   PetscCall(MatDenseRestoreArrayRead(U, &array));
460:   PetscFunctionReturn(PETSC_SUCCESS);
461: }

463: static PetscErrorCode KSPHPDDMGetDeflationMat_HPDDM(KSP ksp, Mat *U)
464: {
465:   KSP_HPDDM            *data = (KSP_HPDDM *)ksp->data;
466:   HPDDM::PETScOperator *op   = data->op;
467:   Mat                   A;
468:   const PetscScalar    *array;
469:   PetscScalar          *copy;
470:   PetscInt              m1, M1, N2;

472:   PetscFunctionBegin;
473:   if (!op) {
474:     PetscCall(KSPSetUp(ksp));
475:     op = data->op;
476:   }
477:   PetscCheck(data->precision == PETSC_KSPHPDDM_DEFAULT_PRECISION, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s != %s", KSPHPDDMPrecisionTypes[data->precision], KSPHPDDMPrecisionTypes[PETSC_KSPHPDDM_DEFAULT_PRECISION]);
478:   array = op->storage();
479:   N2    = op->k().first * op->k().second;
480:   if (!array) *U = NULL;
481:   else {
482:     PetscCall(KSPGetOperators(ksp, &A, NULL));
483:     PetscCall(MatGetLocalSize(A, &m1, NULL));
484:     PetscCall(MatGetSize(A, &M1, NULL));
485:     PetscCall(MatCreateDense(PetscObjectComm((PetscObject)ksp), m1, PETSC_DECIDE, M1, N2, NULL, U));
486:     PetscCall(MatDenseGetArrayWrite(*U, &copy));
487:     PetscCall(PetscArraycpy(copy, array, m1 * N2));
488:     PetscCall(MatDenseRestoreArrayWrite(*U, &copy));
489:   }
490:   PetscFunctionReturn(PETSC_SUCCESS);
491: }

493: static PetscErrorCode KSPMatSolve_HPDDM(KSP ksp, Mat B, Mat X)
494: {
495:   KSP_HPDDM         *data = (KSP_HPDDM *)ksp->data;
496:   Mat                A;
497:   const PetscScalar *b;
498:   PetscScalar       *x;
499:   PetscInt           n, lda;
500:   PetscMemType       type[2];

502:   PetscFunctionBegin;
503:   PetscCall(PetscCitationsRegister(HPDDMCitation, &HPDDMCite));
504:   if (!data->op) PetscCall(KSPSetUp(ksp));
505:   PetscCall(KSPGetOperators(ksp, &A, NULL));
506:   PetscCall(MatGetLocalSize(B, &n, NULL));
507:   PetscCall(MatDenseGetLDA(B, &lda));
508:   PetscCheck(n == lda, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Unhandled leading dimension lda = %" PetscInt_FMT " with n = %" PetscInt_FMT, lda, n);
509:   PetscCall(MatGetLocalSize(A, &n, NULL));
510:   PetscCall(MatDenseGetLDA(X, &lda));
511:   PetscCheck(n == lda, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Unhandled leading dimension lda = %" PetscInt_FMT " with n = %" PetscInt_FMT, lda, n);
512:   PetscCall(MatGetSize(X, NULL, &n));
513:   PetscCall(MatDenseGetArrayWriteAndMemType(X, &x, type));
514:   PetscCall(MatDenseGetArrayReadAndMemType(B, &b, type + 1));
515:   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]));
516:   if (PetscMemTypeCUDA(type[0])) PetscCall(KSPSolve_HPDDM_Private<PETSC_MEMTYPE_CUDA>(ksp, b, x, n));
517:   else {
518:     PetscCheck(PetscMemTypeHost(type[0]), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "PetscMemType (%s) is neither PETSC_MEMTYPE_HOST nor PETSC_MEMTYPE_CUDA", PetscMemTypeToString(type[0]));
519:     PetscCall(KSPSolve_HPDDM_Private(ksp, b, x, n));
520:   }
521:   PetscCall(MatDenseRestoreArrayReadAndMemType(B, &b));
522:   PetscCall(MatDenseRestoreArrayWriteAndMemType(X, &x));
523:   PetscFunctionReturn(PETSC_SUCCESS);
524: }

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

529:    Collective

531:    Input Parameters:
532: +     ksp - iterative context
533: -     type - any of gmres, bgmres, cg, bcg, gcrodr, bgcrodr, bfbcg, or preonly

535:    Level: intermediate

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

540:      As an example, in the following sequence:
541: .vb
542:      KSPHPDDMSetType(ksp, KSPGCRODR);
543:      KSPSolve(ksp, b, x);
544:      KSPHPDDMSetType(ksp, KSPGMRES);
545:      KSPHPDDMSetType(ksp, KSPGCRODR);
546:      KSPSolve(ksp, b, x);
547: .ve
548:     the recycled space is reused in the second `KSPSolve()`.

550: .seealso: [](chapter_ksp), `KSPCreate()`, `KSPType`, `KSPHPDDMType`, `KSPHPDDMGetType()`
551: @*/
552: PetscErrorCode KSPHPDDMSetType(KSP ksp, KSPHPDDMType type)
553: {
554:   PetscFunctionBegin;
557:   PetscUseMethod(ksp, "KSPHPDDMSetType_C", (KSP, KSPHPDDMType), (ksp, type));
558:   PetscFunctionReturn(PETSC_SUCCESS);
559: }

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

564:    Input Parameter:
565: .     ksp - iterative context

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

570:    Level: intermediate

572: .seealso: [](chapter_ksp), `KSPCreate()`, `KSPType`, `KSPHPDDMType`, `KSPHPDDMSetType()`
573: @*/
574: PetscErrorCode KSPHPDDMGetType(KSP ksp, KSPHPDDMType *type)
575: {
576:   PetscFunctionBegin;
578:   if (type) {
580:     PetscUseMethod(ksp, "KSPHPDDMGetType_C", (KSP, KSPHPDDMType *), (ksp, type));
581:   }
582:   PetscFunctionReturn(PETSC_SUCCESS);
583: }

585: static PetscErrorCode KSPHPDDMSetType_HPDDM(KSP ksp, KSPHPDDMType type)
586: {
587:   KSP_HPDDM *data = (KSP_HPDDM *)ksp->data;
588:   PetscInt   i;
589:   PetscBool  flg = PETSC_FALSE;

591:   PetscFunctionBegin;
592:   for (i = 0; i < static_cast<PetscInt>(PETSC_STATIC_ARRAY_LENGTH(KSPHPDDMTypes)); ++i) {
593:     PetscCall(PetscStrcmp(KSPHPDDMTypes[type], KSPHPDDMTypes[i], &flg));
594:     if (flg) break;
595:   }
596:   PetscCheck(i != PETSC_STATIC_ARRAY_LENGTH(KSPHPDDMTypes), PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown KSPHPDDMType %d", type);
597:   if (data->cntl[0] != static_cast<char>(PETSC_DECIDE) && data->cntl[0] != i) PetscCall(KSPReset_HPDDM_Private(ksp));
598:   data->cntl[0] = i;
599:   PetscFunctionReturn(PETSC_SUCCESS);
600: }

602: static PetscErrorCode KSPHPDDMGetType_HPDDM(KSP ksp, KSPHPDDMType *type)
603: {
604:   KSP_HPDDM *data = (KSP_HPDDM *)ksp->data;

606:   PetscFunctionBegin;
607:   PetscCheck(data->cntl[0] != static_cast<char>(PETSC_DECIDE), PETSC_COMM_SELF, PETSC_ERR_ORDER, "KSPHPDDMType not set yet");
608:   /* need to shift by -1 for HPDDM_KRYLOV_METHOD_NONE */
609:   *type = static_cast<KSPHPDDMType>(PetscMin(data->cntl[0], static_cast<char>(PETSC_STATIC_ARRAY_LENGTH(KSPHPDDMTypes) - 1)));
610:   PetscFunctionReturn(PETSC_SUCCESS);
611: }

613: /*MC
614:      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.,
615:      GCRODR [2006], a recycled Krylov method which is similar to `KSPLGMRES`, see [2016] for a comparison. ex75.c shows how to reproduce the results
616:      from the aforementioned paper [2006]. A chronological bibliography of relevant publications linked with `KSP` available in HPDDM through `KSPHPDDM`,
617:      and not available directly in PETSc, may be found below. The interface is explained in details in [2021].

619:    Options Database Keys:
620: +   -ksp_gmres_restart <restart, default=30> - see `KSPGMRES`
621: .   -ksp_hpddm_type <type, default=gmres> - any of gmres, bgmres, cg, bcg, gcrodr, bgcrodr, bfbcg, or preonly, see `KSPHPDDMType`
622: .   -ksp_hpddm_precision <value, default=same as PetscScalar> - any of single or double, see `KSPHPDDMPrecision`
623: .   -ksp_hpddm_deflation_tol <eps, default=\-1.0> - tolerance when deflating right-hand sides inside block methods (no deflation by default, only relevant with block methods)
624: .   -ksp_hpddm_enlarge_krylov_subspace <p, default=1> - split the initial right-hand side into multiple vectors (only relevant with nonblock methods)
625: .   -ksp_hpddm_orthogonalization <type, default=cgs> - any of cgs or mgs, see KSPGMRES
626: .   -ksp_hpddm_qr <type, default=cholqr> - distributed QR factorizations with any of cholqr, cgs, or mgs (only relevant with block methods)
627: .   -ksp_hpddm_variant <type, default=left> - any of left, right, or flexible (this option is superseded by `KSPSetPCSide()`)
628: .   -ksp_hpddm_recycle <n, default=0> - number of harmonic Ritz vectors to compute (only relevant with GCRODR or BGCRODR)
629: .   -ksp_hpddm_recycle_target <type, default=SM> - criterion to select harmonic Ritz vectors using either SM, LM, SR, LR, SI, or LI (only relevant with GCRODR or BGCRODR).
630:      For BGCRODR, if PETSc is compiled with SLEPc, this option is not relevant, since SLEPc is used instead. Options are set with the prefix -ksp_hpddm_recycle_eps_
631: .   -ksp_hpddm_recycle_strategy <type, default=A> - generalized eigenvalue problem A or B to solve for recycling (only relevant with flexible GCRODR or BGCRODR)
632: -   -ksp_hpddm_recycle_symmetric <true, default=false> - symmetric generalized eigenproblems in BGCRODR, useful to switch to distributed solvers like EPSELEMENTAL or EPSSCALAPACK
633:      (only relevant when PETSc is compiled with SLEPc)

635:    Level: intermediate

637:    References:
638: +   1980 - The block conjugate gradient algorithm and related methods. O'Leary. Linear Algebra and its Applications.
639: .   2006 - Recycling Krylov subspaces for sequences of linear systems. Parks, de Sturler, Mackey, Johnson, and Maiti. SIAM Journal on Scientific Computing
640: .   2013 - A modified block flexible GMRES method with deflation at each iteration for the solution of non-Hermitian linear systems with multiple right-hand sides.
641:            Calandra, Gratton, Lago, Vasseur, and Carvalho. SIAM Journal on Scientific Computing.
642: .   2016 - Block iterative methods and recycling for improved scalability of linear solvers. Jolivet and Tournier. SC16.
643: .   2017 - A breakdown-free block conjugate gradient method. Ji and Li. BIT Numerical Mathematics.
644: -   2021 - KSPHPDDM and PCHPDDM: extending PETSc with advanced Krylov methods and robust multilevel overlapping Schwarz preconditioners. Jolivet, Roman, and Zampini.
645:            Computer & Mathematics with Applications.

647: .seealso: [](chapter_ksp), [](sec_flexibleksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPGMRES`, `KSPCG`, `KSPLGMRES`, `KSPDGMRES`
648: M*/

650: PETSC_EXTERN PetscErrorCode KSPCreate_HPDDM(KSP ksp)
651: {
652:   KSP_HPDDM  *data;
653:   PetscInt    i;
654:   const char *common[] = {KSPGMRES, KSPCG, KSPPREONLY};
655:   PetscBool   flg      = PETSC_FALSE;

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