Actual source code: hpddm.cxx

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

  4: /* static array length */
  5: #define ALEN(a) (sizeof(a)/sizeof((a)[0]))

  7: const char *const KSPHPDDMTypes[]          = { KSPGMRES, "bgmres", KSPCG, "bcg", "gcrodr", "bgcrodr", "bfbcg", KSPPREONLY };
  8: const char *const KSPHPDDMPrecisionTypes[] = { "HALF", "SINGLE", "DOUBLE", "QUADRUPLE", "KSPHPDDMPrecisionType", "KSP_HPDDM_PRECISION_", NULL };
  9: const char *const HPDDMOrthogonalization[] = { "cgs", "mgs" };
 10: const char *const HPDDMQR[]                = { "cholqr", "cgs", "mgs" };
 11: const char *const HPDDMVariant[]           = { "left", "right", "flexible" };
 12: const char *const HPDDMRecycleTarget[]     = { "SM", "LM", "SR", "LR", "SI", "LI" };
 13: const char *const HPDDMRecycleStrategy[]   = { "A", "B" };

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

 27: #if defined(PETSC_HAVE_SLEPC) && defined(PETSC_USE_SHARED_LIBRARIES)
 28: static PetscBool loadedDL = PETSC_FALSE;
 29: #endif

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

 37:   PetscOptionsHead(PetscOptionsObject, "KSPHPDDM options, cf. https://github.com/hpddm/hpddm");
 38:   i = (data->cntl[0] == static_cast<char>(PETSC_DECIDE) ? HPDDM_KRYLOV_METHOD_GMRES : data->cntl[0]);
 39:   PetscOptionsEList("-ksp_hpddm_type", "Type of Krylov method", "KSPHPDDMGetType", KSPHPDDMTypes, ALEN(KSPHPDDMTypes), KSPHPDDMTypes[HPDDM_KRYLOV_METHOD_GMRES], &i, NULL);
 40:   if (i == ALEN(KSPHPDDMTypes) - 1)
 41:     i = HPDDM_KRYLOV_METHOD_NONE; /* need to shift the value since HPDDM_KRYLOV_METHOD_RICHARDSON is not registered in PETSc */
 42:   data->cntl[0] = i;
 43:   PetscOptionsEnum("-ksp_hpddm_precision", "Precision in which Krylov bases are stored", "KSPHPDDM", KSPHPDDMPrecisionTypes, (PetscEnum)data->precision, (PetscEnum*)&data->precision, NULL);
 45:   if (data->cntl[0] != HPDDM_KRYLOV_METHOD_NONE) {
 46:     if (data->cntl[0] != HPDDM_KRYLOV_METHOD_BCG && data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG) {
 47:       i = (data->cntl[1] == static_cast<char>(PETSC_DECIDE) ? HPDDM_VARIANT_LEFT : data->cntl[1]);
 48:       if (ksp->pc_side_set == PC_SIDE_DEFAULT) PetscOptionsEList("-ksp_hpddm_variant", "Left, right, or variable preconditioning", "KSPHPDDM", HPDDMVariant, ALEN(HPDDMVariant), HPDDMVariant[HPDDM_VARIANT_LEFT], &i, NULL);
 49:       else if (ksp->pc_side_set == PC_RIGHT) i = HPDDM_VARIANT_RIGHT;
 51:       data->cntl[1] = i;
 52:       if (i > 0) KSPSetPCSide(ksp, PC_RIGHT);
 53:     }
 54:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
 55:       data->rcntl[0] = (std::abs(data->rcntl[0] - static_cast<PetscReal>(PETSC_DECIDE)) < PETSC_SMALL ? -1.0 : data->rcntl[0]);
 56:       PetscOptionsReal("-ksp_hpddm_deflation_tol", "Tolerance when deflating right-hand sides inside block methods", "KSPHPDDM", data->rcntl[0], data->rcntl, NULL);
 57:       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]));
 58:       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);
 59:       data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG] = i;
 60:     } else data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BCG] = 0;
 61:     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) {
 62:       i = (data->cntl[2] == static_cast<char>(PETSC_DECIDE) ? HPDDM_ORTHOGONALIZATION_CGS : data->cntl[2] & 3);
 63:       PetscOptionsEList("-ksp_hpddm_orthogonalization", "Classical (faster) or Modified (more robust) Gram--Schmidt process", "KSPHPDDM", HPDDMOrthogonalization, ALEN(HPDDMOrthogonalization), HPDDMOrthogonalization[HPDDM_ORTHOGONALIZATION_CGS], &i, NULL);
 64:       j = (data->cntl[2] == static_cast<char>(PETSC_DECIDE) ? HPDDM_QR_CHOLQR : ((data->cntl[2] >> 2) & 7));
 65:       PetscOptionsEList("-ksp_hpddm_qr", "Distributed QR factorizations computed with Cholesky QR, Classical or Modified Gram--Schmidt process", "KSPHPDDM", HPDDMQR, ALEN(HPDDMQR), HPDDMQR[HPDDM_QR_CHOLQR], &j, NULL);
 66:       data->cntl[2] = static_cast<char>(i) + (static_cast<char>(j) << 2);
 67:       i = (data->scntl[0] == static_cast<unsigned short>(PETSC_DECIDE) ? PetscMin(30, ksp->max_it) : data->scntl[0]);
 68:       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));
 69:       data->scntl[0] = i;
 70:     }
 71:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BCG || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
 72:       j = (data->cntl[1] == static_cast<char>(PETSC_DECIDE) ? HPDDM_QR_CHOLQR : data->cntl[1]);
 73:       PetscOptionsEList("-ksp_hpddm_qr", "Distributed QR factorizations computed with Cholesky QR, Classical or Modified Gram--Schmidt process", "KSPHPDDM", HPDDMQR, ALEN(HPDDMQR), HPDDMQR[HPDDM_QR_CHOLQR], &j, NULL);
 74:       data->cntl[1] = j;
 75:     }
 76:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
 77:       i = (data->icntl[0] == static_cast<int>(PETSC_DECIDE) ? PetscMin(20, data->scntl[0] - 1) : data->icntl[0]);
 78:       PetscOptionsRangeInt("-ksp_hpddm_recycle", "Number of harmonic Ritz vectors to compute", "KSPHPDDM", i, &i, NULL, 1, data->scntl[0] - 1);
 79:       data->icntl[0] = i;
 80:       if (!PetscDefined(HAVE_SLEPC) || !PetscDefined(USE_SHARED_LIBRARIES) || data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR) {
 81:         i = (data->cntl[3] == static_cast<char>(PETSC_DECIDE) ? HPDDM_RECYCLE_TARGET_SM : data->cntl[3]);
 82:         PetscOptionsEList("-ksp_hpddm_recycle_target", "Criterion to select harmonic Ritz vectors", "KSPHPDDM", HPDDMRecycleTarget, ALEN(HPDDMRecycleTarget), HPDDMRecycleTarget[HPDDM_RECYCLE_TARGET_SM], &i, NULL);
 83:         data->cntl[3] = i;
 84:       } else {
 86:         MPI_Comm_size(PetscObjectComm((PetscObject)ksp), &size);
 87:         i = (data->cntl[3] == static_cast<char>(PETSC_DECIDE) ? 1 : data->cntl[3]);
 88:         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));
 89:         data->cntl[3] = i;
 90:       }
 91:       i = (data->cntl[4] == static_cast<char>(PETSC_DECIDE) ? HPDDM_RECYCLE_STRATEGY_A : data->cntl[4]);
 92:       PetscOptionsEList("-ksp_hpddm_recycle_strategy", "Generalized eigenvalue problem to solve for recycling", "KSPHPDDM", HPDDMRecycleStrategy, ALEN(HPDDMRecycleStrategy), HPDDMRecycleStrategy[HPDDM_RECYCLE_STRATEGY_A], &i, NULL);
 93:       data->cntl[4] = i;
 94:     }
 95:   } else {
 96:     data->cntl[0] = HPDDM_KRYLOV_METHOD_NONE;
 97:     data->scntl[1] = 1;
 98:   }
100:   else data->icntl[1] = static_cast<int>(ksp->nmax);
101:   PetscOptionsTail();
102:   return 0;
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() : NULL;
110:   PetscBool            ascii;

112:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii);
113:   if (op && ascii) {
114:     PetscViewerASCIIPrintf(viewer, "HPDDM type: %s\n", KSPHPDDMTypes[std::min(static_cast<PetscInt>(data->cntl[0]), static_cast<PetscInt>(ALEN(KSPHPDDMTypes) - 1))]);
115:     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 (std::abs(data->rcntl[0] - static_cast<PetscReal>(PETSC_DECIDE)) < PETSC_SMALL) PetscViewerASCIIPrintf(viewer, "no deflation at restarts\n", PetscBools[array ? PETSC_TRUE : PETSC_FALSE]);
118:       else PetscViewerASCIIPrintf(viewer, "deflation tolerance: %g\n", data->rcntl[0]);
119:     }
120:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
121:       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) PetscViewerASCIIPrintf(viewer, "deflation target: %s\n", HPDDMRecycleTarget[static_cast<PetscInt>(data->cntl[3])]);
123:       else PetscViewerASCIIPrintf(viewer, "redistribution size: %d\n", static_cast<PetscMPIInt>(data->cntl[3]));
124:     }
125:     if (data->icntl[1] != static_cast<int>(PETSC_DECIDE)) PetscViewerASCIIPrintf(viewer, "  block size is %d\n", data->icntl[1]);
126:   }
127:   return 0;
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:   KSPGetOperators(ksp, &A, NULL);
138:   MatGetLocalSize(A, &n, NULL);
139:   MatGetBlockSize(A, &bs);
140:   PetscObjectTypeCompareAny((PetscObject)A, &match, MATSEQKAIJ, MATMPIKAIJ, "");
141:   if (match) n /= bs;
142:   data->op = new HPDDM::PETScOperator(ksp, n);
143:   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() */
144:     PetscInfo(ksp, "KSPSetFromOptions() not called or uninitialized internal structure, hardwiring default KSPHPDDM options\n");
145:     if (data->cntl[0] == static_cast<char>(PETSC_DECIDE))
146:       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;
152:         if (data->cntl[1] > 0) 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) {
163:         data->cntl[1] = HPDDM_QR_CHOLQR; /* CholQR by default */
164:       }
165:       if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
166:         data->icntl[0] = PetscMin(20, data->scntl[0] - 1); /* recycled subspace of size 20 by default */
167:         if (!PetscDefined(HAVE_SLEPC) || !PetscDefined(USE_SHARED_LIBRARIES) || data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR) {
168:           data->cntl[3] = HPDDM_RECYCLE_TARGET_SM; /* default recycling target */
169:         } else {
170:           data->cntl[3] = 1; /* redistribution parameter of 1 by default */
171:         }
172:         data->cntl[4] = HPDDM_RECYCLE_STRATEGY_A; /* default recycling strategy */
173:       }
174:     } else data->scntl[1] = 1;
175:   }
177:   else data->icntl[1] = static_cast<int>(ksp->nmax);
178:   return 0;
179: }

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

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

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

198:   if (data->op) {
199:     delete data->op;
200:     data->op = NULL;
201:   }
202:   KSPHPDDMReset_Private(ksp);
203:   return 0;
204: }

206: static PetscErrorCode KSPDestroy_HPDDM(KSP ksp)
207: {
208:   KSPReset_HPDDM(ksp);
209:   KSPDestroyDefault(ksp);
210:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetDeflationSpace_C", NULL);
211:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetDeflationSpace_C", NULL);
212:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetType_C", NULL);
213:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetType_C", NULL);
214:   return 0;
215: }

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:   HPDDM::upscaled_type<PetscScalar>   *dbl[2];
222:   HPDDM::downscaled_type<PetscScalar> *sgl[2];
223:   const PetscInt                      N = data->op->getDof() * n;
224:   PetscBool                           scale;

226:   PCGetDiagonalScale(ksp->pc, &scale);
228:   if (n > 1) {
229:     if (ksp->converged == KSPConvergedDefault) {
231:       if (!ctx->initialrtol) {
232:         PetscInfo(ksp, "Forcing KSPConvergedDefaultSetUIRNorm() since KSPConvergedDefault() cannot handle multiple norms\n");
233:         ctx->initialrtol = PETSC_TRUE;
234:       }
235:     } else PetscInfo(ksp, "Using a special \"converged\" callback, be careful, it is used in KSPHPDDM to track blocks of residuals\n");
236:   }
237:   /* initial guess is always nonzero with recycling methods if there is a deflation subspace available */
238:   if ((data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) && data->op->storage()) ksp->guess_zero = PETSC_FALSE;
239:   ksp->its = 0;
240:   ksp->reason = KSP_CONVERGED_ITERATING;
241:   if (data->precision == KSP_HPDDM_PRECISION_DOUBLE && PetscDefined(USE_REAL_SINGLE)) {
242:     PetscMalloc2(N, dbl, N, dbl + 1);
243:     std::copy_n(b, N, dbl[0]);
244:     std::copy_n(x, N, dbl[1]);
245:     static_cast<PetscErrorCode>(HPDDM::IterativeMethod::solve(*data->op, dbl[0], dbl[1], n, PetscObjectComm((PetscObject)ksp)));
246:     std::copy_n(dbl[1], N, x);
247:     PetscFree2(dbl[0], dbl[1]);
248:   } else if (data->precision == KSP_HPDDM_PRECISION_SINGLE && PetscDefined(USE_REAL_DOUBLE)) {
249:     PetscMalloc1(N, sgl);
250:     sgl[1] = reinterpret_cast<HPDDM::downscaled_type<PetscScalar>*>(x);
251:     std::copy_n(b, N, sgl[0]);
252:     for (PetscInt i = 0; i < N; ++i) sgl[1][i] = x[i];
253:     static_cast<PetscErrorCode>(HPDDM::IterativeMethod::solve(*data->op, sgl[0], sgl[1], n, PetscObjectComm((PetscObject)ksp)));
254:     if (N) {
255:       sgl[0][0] = sgl[1][0];
256:       std::copy_backward(sgl[1] + 1, sgl[1] + N, x + N);
257:       x[0] = sgl[0][0];
258:     }
259:     PetscFree(sgl[0]);
260:   } else static_cast<PetscErrorCode>(HPDDM::IterativeMethod::solve(*data->op, b, x, n, PetscObjectComm((PetscObject)ksp)));
261:   if (!ksp->reason) { /* KSPConvergedDefault() is still returning 0 (= KSP_CONVERGED_ITERATING) */
262:     if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
263:     else ksp->reason = KSP_CONVERGED_RTOL; /* early exit by HPDDM, which only happens on breakdowns or convergence */
264:   }
265:   ksp->its = PetscMin(ksp->its, ksp->max_it);
266:   return 0;
267: }

269: static PetscErrorCode KSPSolve_HPDDM(KSP ksp)
270: {
271:   KSP_HPDDM         *data = (KSP_HPDDM*)ksp->data;
272:   Mat               A, B;
273:   PetscScalar       *x, *bt = NULL, **ptr;
274:   const PetscScalar *b;
275:   PetscInt          i, j, n;
276:   PetscBool         flg;

278:   PetscCitationsRegister(HPDDMCitation, &HPDDMCite);
279:   KSPGetOperators(ksp, &A, NULL);
280:   PetscObjectTypeCompareAny((PetscObject)A, &flg, MATSEQKAIJ, MATMPIKAIJ, "");
281:   VecGetArrayWrite(ksp->vec_sol, &x);
282:   VecGetArrayRead(ksp->vec_rhs, &b);
283:   if (!flg) KSPSolve_HPDDM_Private(ksp, b, x, 1);
284:   else {
285:     MatKAIJGetScaledIdentity(A, &flg);
286:     MatKAIJGetAIJ(A, &B);
287:     MatGetBlockSize(A, &n);
288:     MatGetLocalSize(B, &i, NULL);
289:     j = data->op->getDof();
290:     if (!flg) i *= n; /* S and T are not scaled identities, cannot use block methods */
291:     if (i != j) { /* switching between block and standard methods */
292:       delete data->op;
293:       data->op = new HPDDM::PETScOperator(ksp, i);
294:     }
295:     if (flg && n > 1) {
296:       PetscMalloc1(i * n, &bt);
297:       /* from row- to column-major to be consistent with HPDDM */
298:       HPDDM::Wrapper<PetscScalar>::omatcopy<'T'>(i, n, b, n, bt, i);
299:       ptr = const_cast<PetscScalar**>(&b);
300:       std::swap(*ptr, bt);
301:       HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(i, n, x, n, i);
302:     }
303:     KSPSolve_HPDDM_Private(ksp, b, x, flg ? n : 1);
304:     if (flg && n > 1) {
305:       std::swap(*ptr, bt);
306:       PetscFree(bt);
307:       /* from column- to row-major to be consistent with MatKAIJ format */
308:       HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(n, i, x, i, n);
309:     }
310:   }
311:   VecRestoreArrayRead(ksp->vec_rhs, &b);
312:   VecRestoreArrayWrite(ksp->vec_sol, &x);
313:   return 0;
314: }

316: /*@
317:      KSPHPDDMSetDeflationSpace - Sets the deflation space used by Krylov methods with recycling. This space is viewed as a set of vectors stored in a MATDENSE (column major).

319:    Input Parameters:
320: +     ksp - iterative context
321: -     U - deflation space to be used during KSPSolve()

323:    Level: intermediate

325: .seealso:  KSPCreate(), KSPType (for list of available types), KSPHPDDMGetDeflationSpace()
326: @*/
327: PetscErrorCode KSPHPDDMSetDeflationSpace(KSP ksp, Mat U)
328: {
332:   PetscUseMethod(ksp, "KSPHPDDMSetDeflationSpace_C", (KSP, Mat), (ksp, U));
333:   return 0;
334: }

336: /*@
337:      KSPHPDDMGetDeflationSpace - Gets the deflation space computed by Krylov methods with recycling or NULL if KSPSolve() has not been called yet. 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.

339:    Input Parameter:
340: .     ksp - iterative context

342:    Output Parameter:
343: .     U - deflation space generated during KSPSolve()

345:    Level: intermediate

347: .seealso:  KSPCreate(), KSPType (for list of available types), KSPHPDDMSetDeflationSpace()
348: @*/
349: PetscErrorCode KSPHPDDMGetDeflationSpace(KSP ksp, Mat *U)
350: {
352:   if (U) {
354:     PetscUseMethod(ksp, "KSPHPDDMGetDeflationSpace_C", (KSP, Mat*), (ksp, U));
355:   }
356:   return 0;
357: }

359: static PetscErrorCode KSPHPDDMSetDeflationSpace_HPDDM(KSP ksp, Mat U)
360: {
361:   KSP_HPDDM            *data = (KSP_HPDDM*)ksp->data;
362:   HPDDM::PETScOperator *op = data->op;
363:   Mat                  A;
364:   const PetscScalar    *array;
365:   PetscScalar          *copy;
366:   PetscInt             m1, M1, m2, M2, n2, N2, ldu;
367:   PetscBool            match;

369:   if (!op) {
370:     KSPSetUp(ksp);
371:     op = data->op;
372:   }
374:   KSPGetOperators(ksp, &A, NULL);
375:   MatGetLocalSize(A, &m1, NULL);
376:   MatGetLocalSize(U, &m2, &n2);
377:   MatGetSize(A, &M1, NULL);
378:   MatGetSize(U, &M2, &N2);
380:   PetscObjectTypeCompareAny((PetscObject)U, &match, MATSEQDENSE, MATMPIDENSE, "");
382:   MatDenseGetArrayRead(U, &array);
383:   copy = op->allocate(m2, 1, N2);
385:   MatDenseGetLDA(U, &ldu);
386:   HPDDM::Wrapper<PetscScalar>::omatcopy<'N'>(N2, m2, array, ldu, copy, m2);
387:   MatDenseRestoreArrayRead(U, &array);
388:   return 0;
389: }

391: static PetscErrorCode KSPHPDDMGetDeflationSpace_HPDDM(KSP ksp, Mat *U)
392: {
393:   KSP_HPDDM            *data = (KSP_HPDDM*)ksp->data;
394:   HPDDM::PETScOperator *op = data->op;
395:   Mat                  A;
396:   const PetscScalar    *array;
397:   PetscScalar          *copy;
398:   PetscInt             m1, M1, N2;

400:   if (!op) {
401:     KSPSetUp(ksp);
402:     op = data->op;
403:   }
405:   array = op->storage();
406:   N2 = op->k().first * op->k().second;
407:   if (!array) *U = NULL;
408:   else {
409:     KSPGetOperators(ksp, &A, NULL);
410:     MatGetLocalSize(A, &m1, NULL);
411:     MatGetSize(A, &M1, NULL);
412:     MatCreateDense(PetscObjectComm((PetscObject)ksp), m1, PETSC_DECIDE, M1, N2, NULL, U);
413:     MatDenseGetArrayWrite(*U, &copy);
414:     PetscArraycpy(copy, array, m1 * N2);
415:     MatDenseRestoreArrayWrite(*U, &copy);
416:   }
417:   return 0;
418: }

420: static PetscErrorCode KSPMatSolve_HPDDM(KSP ksp, Mat B, Mat X)
421: {
422:   KSP_HPDDM            *data = (KSP_HPDDM*)ksp->data;
423:   HPDDM::PETScOperator *op = data->op;
424:   Mat                  A;
425:   const PetscScalar    *b;
426:   PetscScalar          *x;
427:   PetscInt             n, lda;

429:   PetscCitationsRegister(HPDDMCitation, &HPDDMCite);
430:   if (!op) {
431:     KSPSetUp(ksp);
432:     op = data->op;
433:   }
434:   KSPGetOperators(ksp, &A, NULL);
435:   MatGetLocalSize(B, &n, NULL);
436:   MatDenseGetLDA(B, &lda);
438:   MatGetLocalSize(A, &n, NULL);
439:   MatDenseGetLDA(X, &lda);
441:   MatDenseGetArrayRead(B, &b);
442:   MatDenseGetArrayWrite(X, &x);
443:   MatGetSize(X, NULL, &n);
444:   KSPSolve_HPDDM_Private(ksp, b, x, n);
445:   MatDenseRestoreArrayWrite(X, &x);
446:   MatDenseRestoreArrayRead(B, &b);
447:   return 0;
448: }

450: /*@
451:      KSPHPDDMSetType - Sets the type of Krylov method used in KSPHPDDM.

453:    Input Parameters:
454: +     ksp - iterative context
455: -     type - any of gmres, bgmres, cg, bcg, gcrodr, bgcrodr, bfbcg, or preonly

457:    Level: intermediate

459:    Notes:
460:      Unlike KSPReset(), this function does not destroy any deflation space attached to the KSP.
461:      As an example, in the following sequence: KSPHPDDMSetType(ksp, KSPGCRODR); KSPSolve(ksp, b, x); KSPHPDDMSetType(ksp, KSPGMRES); KSPHPDDMSetType(ksp, KSPGCRODR); KSPSolve(ksp, b, x); the recycled space is reused in the second KSPSolve().

463: .seealso:  KSPCreate(), KSPType (for list of available types), KSPHPDDMType, KSPHPDDMGetType()
464: @*/
465: PetscErrorCode KSPHPDDMSetType(KSP ksp, KSPHPDDMType type)
466: {
468:   PetscUseMethod(ksp, "KSPHPDDMSetType_C", (KSP, KSPHPDDMType), (ksp, type));
469:   return 0;
470: }

472: /*@
473:      KSPHPDDMGetType - Gets the type of Krylov method used in KSPHPDDM.

475:    Input Parameter:
476: .     ksp - iterative context

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

481:    Level: intermediate

483: .seealso:  KSPCreate(), KSPType (for list of available types), KSPHPDDMType, KSPHPDDMSetType()
484: @*/
485: PetscErrorCode KSPHPDDMGetType(KSP ksp, KSPHPDDMType *type)
486: {
488:   if (type) {
490:     PetscUseMethod(ksp, "KSPHPDDMGetType_C", (KSP, KSPHPDDMType*), (ksp, type));
491:   }
492:   return 0;
493: }

495: static PetscErrorCode KSPHPDDMSetType_HPDDM(KSP ksp, KSPHPDDMType type)
496: {
497:   KSP_HPDDM *data = (KSP_HPDDM*)ksp->data;
498:   PetscInt  i;
499:   PetscBool flg = PETSC_FALSE;

501:   for (i = 0; i < static_cast<PetscInt>(ALEN(KSPHPDDMTypes)); ++i) {
502:     PetscStrcmp(KSPHPDDMTypes[type], KSPHPDDMTypes[i], &flg);
503:     if (flg) break;
504:   }
506:   if (data->cntl[0] != static_cast<char>(PETSC_DECIDE) && data->cntl[0] != i) KSPHPDDMReset_Private(ksp);
507:   data->cntl[0] = i;
508:   return 0;
509: }

511: static PetscErrorCode KSPHPDDMGetType_HPDDM(KSP ksp, KSPHPDDMType *type)
512: {
513:   KSP_HPDDM *data = (KSP_HPDDM*)ksp->data;

516:   /* need to shift by -1 for HPDDM_KRYLOV_METHOD_NONE */
517:   *type = static_cast<KSPHPDDMType>(PetscMin(data->cntl[0], static_cast<char>(ALEN(KSPHPDDMTypes) - 1)));
518:   return 0;
519: }

521: /*MC
522:      KSPHPDDM - Interface with the HPDDM library.

524:    This KSP may be used to further select methods that are currently not implemented natively in PETSc, e.g., GCRODR [2006], a recycled Krylov method which is similar to KSPLGMRES, see [2016] for a comparison. ex75.c shows how to reproduce the results from the aforementioned paper [2006]. A chronological bibliography of relevant publications linked with KSP available in HPDDM through KSPHPDDM, and not available directly in PETSc, may be found below. The interface is explained in details in [2021].

526:    Options Database Keys:
527: +   -ksp_gmres_restart <restart, default=30> - see KSPGMRES
528: .   -ksp_hpddm_type <type, default=gmres> - any of gmres, bgmres, cg, bcg, gcrodr, bgcrodr, bfbcg, or preonly, see KSPHPDDMType
529: .   -ksp_hpddm_precision <value, default=same as PetscScalar> - any of single or double, see KSPHPDDMPrecision
530: .   -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)
531: .   -ksp_hpddm_enlarge_krylov_subspace <p, default=1> - split the initial right-hand side into multiple vectors (only relevant with nonblock methods)
532: .   -ksp_hpddm_orthogonalization <type, default=cgs> - any of cgs or mgs, see KSPGMRES
533: .   -ksp_hpddm_qr <type, default=cholqr> - distributed QR factorizations with any of cholqr, cgs, or mgs (only relevant with block methods)
534: .   -ksp_hpddm_variant <type, default=left> - any of left, right, or flexible (this option is superseded by KSPSetPCSide())
535: .   -ksp_hpddm_recycle <n, default=0> - number of harmonic Ritz vectors to compute (only relevant with GCRODR or BGCRODR)
536: .   -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). 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_
537: .   -ksp_hpddm_recycle_strategy <type, default=A> - generalized eigenvalue problem A or B to solve for recycling (only relevant with flexible GCRODR or BGCRODR)
538: -   -ksp_hpddm_recycle_symmetric <true, default=false> - symmetric generalized eigenproblems in BGCRODR, useful to switch to distributed solvers like EPSELEMENTAL or EPSSCALAPACK (only relevant when PETSc is compiled with SLEPc)

540:    References:
541: +   1980 - The block conjugate gradient algorithm and related methods. O'Leary. Linear Algebra and its Applications.
542: .   2006 - Recycling Krylov subspaces for sequences of linear systems. Parks, de Sturler, Mackey, Johnson, and Maiti. SIAM Journal on Scientific Computing
543: .   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. Calandra, Gratton, Lago, Vasseur, and Carvalho. SIAM Journal on Scientific Computing.
544: .   2016 - Block iterative methods and recycling for improved scalability of linear solvers. Jolivet and Tournier. SC16.
545: .   2017 - A breakdown-free block conjugate gradient method. Ji and Li. BIT Numerical Mathematics.
546: -   2021 - KSPHPDDM and PCHPDDM: extending PETSc with advanced Krylov methods and robust multilevel overlapping Schwarz preconditioners. Jolivet, Roman, and Zampini. Computer & Mathematics with Applications.

548:    Level: intermediate

550: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPGMRES, KSPCG, KSPLGMRES, KSPDGMRES
551: M*/
552: PETSC_EXTERN PetscErrorCode KSPCreate_HPDDM(KSP ksp)
553: {
554:   KSP_HPDDM  *data;
555:   PetscInt   i;
556:   const char *common[] = { KSPGMRES, KSPCG, KSPPREONLY };
557:   PetscBool  flg = PETSC_FALSE;

559:   PetscNewLog(ksp, &data);
560:   ksp->data = (void*)data;
561:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2);
562:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 1);
563:   ksp->ops->solve          = KSPSolve_HPDDM;
564:   ksp->ops->matsolve       = KSPMatSolve_HPDDM;
565:   ksp->ops->setup          = KSPSetUp_HPDDM;
566:   ksp->ops->setfromoptions = KSPSetFromOptions_HPDDM;
567:   ksp->ops->destroy        = KSPDestroy_HPDDM;
568:   ksp->ops->view           = KSPView_HPDDM;
569:   ksp->ops->reset          = KSPReset_HPDDM;
570:   KSPHPDDMReset_Private(ksp);
571:   for (i = 0; i < static_cast<PetscInt>(ALEN(common)); ++i) {
572:     PetscStrcmp(((PetscObject)ksp)->type_name, common[i], &flg);
573:     if (flg) break;
574:   }
575:   if (!i) data->cntl[0] = HPDDM_KRYLOV_METHOD_GMRES;
576:   else if (i == 1) data->cntl[0] = HPDDM_KRYLOV_METHOD_CG;
577:   else if (i == 2) data->cntl[0] = HPDDM_KRYLOV_METHOD_NONE;
578:   if (data->cntl[0] != static_cast<char>(PETSC_DECIDE)) PetscInfo(ksp, "Using the previously set KSPType %s\n", common[i]);
579:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetDeflationSpace_C", KSPHPDDMSetDeflationSpace_HPDDM);
580:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetDeflationSpace_C", KSPHPDDMGetDeflationSpace_HPDDM);
581:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetType_C", KSPHPDDMSetType_HPDDM);
582:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetType_C", KSPHPDDMGetType_HPDDM);
583: #if defined(PETSC_HAVE_SLEPC) && defined(PETSC_USE_SHARED_LIBRARIES)
584:   if (!loadedDL) HPDDMLoadDL_Private(&loadedDL);
585: #endif
586:   data->precision = PETSC_KSPHPDDM_DEFAULT_PRECISION;
587:   return 0;
588: }