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, ©));
503: PetscCall(PetscArraycpy(copy, array, m1 * N2));
504: PetscCall(MatDenseRestoreArrayWrite(*U, ©));
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, default=30> - see `KSPGMRES`
639: . -ksp_hpddm_type <type, default=gmres> - any of gmres, bgmres, cg, bcg, gcrodr, bgcrodr, bfbcg, or preonly, see `KSPHPDDMType`
640: . -ksp_hpddm_precision <value, default=same as PetscScalar> - any of __fp16, single, double or __float128, see `PetscPrecision`
641: . -ksp_hpddm_deflation_tol <eps, default=\-1.0> - tolerance when deflating right-hand sides inside block methods (no deflation by default,
642: only relevant with block methods)
643: . -ksp_hpddm_enlarge_krylov_subspace <p, default=1> - split the initial right-hand side into multiple vectors (only relevant with nonblock methods)
644: . -ksp_hpddm_orthogonalization <type, default=cgs> - any of cgs or mgs, see KSPGMRES
645: . -ksp_hpddm_qr <type, default=cholqr> - distributed QR factorizations with any of cholqr, cgs, or mgs (only relevant with block methods)
646: . -ksp_hpddm_variant <type, default=left> - any of left, right, or flexible (this option is superseded by `KSPSetPCSide()`)
647: . -ksp_hpddm_recycle <n, default=0> - number of harmonic Ritz vectors to compute (only relevant with GCRODR or BGCRODR)
648: . -ksp_hpddm_recycle_target <type, default=SM> - criterion to select harmonic Ritz vectors using either SM, LM, SR, LR, SI, or LI
649: (only relevant with GCRODR or BGCRODR). For BGCRODR, if PETSc is compiled with SLEPc,
650: this option is not relevant, since SLEPc is used instead. Options are set with the prefix -ksp_hpddm_recycle_eps_
651: . -ksp_hpddm_recycle_strategy <type, default=A> - generalized eigenvalue problem A or B to solve for recycling (only relevant with flexible GCRODR or BGCRODR)
652: - -ksp_hpddm_recycle_symmetric <true, default=false> - symmetric generalized eigenproblems in BGCRODR, useful to switch to distributed solvers like
653: `EPSELEMENTAL` or `EPSSCALAPACK` (only relevant when PETSc is compiled with SLEPc)
655: Level: intermediate
657: .seealso: [](ch_ksp), [](sec_flexibleksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPGMRES`, `KSPCG`, `KSPLGMRES`, `KSPDGMRES`
658: M*/
660: PETSC_EXTERN PetscErrorCode KSPCreate_HPDDM(KSP ksp)
661: {
662: KSP_HPDDM *data;
663: PetscInt i;
664: const char *common[] = {KSPGMRES, KSPCG, KSPPREONLY};
665: PetscBool flg = PETSC_FALSE;
667: PetscFunctionBegin;
668: PetscCall(PetscNew(&data));
669: ksp->data = (void *)data;
670: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
671: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 1));
672: ksp->ops->solve = KSPSolve_HPDDM;
673: ksp->ops->matsolve = KSPMatSolve_HPDDM;
674: ksp->ops->setup = KSPSetUp_HPDDM;
675: ksp->ops->setfromoptions = KSPSetFromOptions_HPDDM;
676: ksp->ops->destroy = KSPDestroy_HPDDM;
677: ksp->ops->view = KSPView_HPDDM;
678: ksp->ops->reset = KSPReset_HPDDM;
679: PetscCall(KSPReset_HPDDM_Private(ksp));
680: for (i = 0; i < static_cast<PetscInt>(PETSC_STATIC_ARRAY_LENGTH(common)); ++i) {
681: PetscCall(PetscStrcmp(((PetscObject)ksp)->type_name, common[i], &flg));
682: if (flg) break;
683: }
684: if (!i) data->cntl[0] = HPDDM_KRYLOV_METHOD_GMRES;
685: else if (i == 1) data->cntl[0] = HPDDM_KRYLOV_METHOD_CG;
686: else if (i == 2) data->cntl[0] = HPDDM_KRYLOV_METHOD_NONE;
687: if (data->cntl[0] != static_cast<char>(PETSC_DECIDE)) PetscCall(PetscInfo(ksp, "Using the previously set KSPType %s\n", common[i]));
688: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetDeflationMat_C", KSPHPDDMSetDeflationMat_HPDDM));
689: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetDeflationMat_C", KSPHPDDMGetDeflationMat_HPDDM));
690: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetType_C", KSPHPDDMSetType_HPDDM));
691: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetType_C", KSPHPDDMGetType_HPDDM));
692: #if PetscDefined(HAVE_SLEPC) && PetscDefined(HAVE_DYNAMIC_LIBRARIES) && PetscDefined(USE_SHARED_LIBRARIES)
693: if (!loadedDL) PetscCall(HPDDMLoadDL_Private(&loadedDL));
694: #endif
695: data->precision = PETSC_SCALAR_PRECISION;
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }