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 defined(PETSC_HAVE_SLEPC) && defined(PETSC_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:   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:   PetscOptionsEList("-ksp_hpddm_type", "Type of Krylov method", "KSPHPDDMGetType", KSPHPDDMTypes, PETSC_STATIC_ARRAY_LENGTH(KSPHPDDMTypes), KSPHPDDMTypes[HPDDM_KRYLOV_METHOD_GMRES], &i, NULL);
 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:   PetscOptionsEnum("-ksp_hpddm_precision", "Precision in which Krylov bases are stored", "KSPHPDDM", KSPHPDDMPrecisionTypes, (PetscEnum)data->precision, (PetscEnum *)&data->precision, NULL);
 42:   if (data->cntl[0] != HPDDM_KRYLOV_METHOD_NONE) {
 43:     if (data->cntl[0] != HPDDM_KRYLOV_METHOD_BCG && data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG) {
 44:       i = (data->cntl[1] == static_cast<char>(PETSC_DECIDE) ? HPDDM_VARIANT_LEFT : data->cntl[1]);
 45:       if (ksp->pc_side_set == PC_SIDE_DEFAULT)
 46:         PetscOptionsEList("-ksp_hpddm_variant", "Left, right, or variable preconditioning", "KSPHPDDM", HPDDMVariant, PETSC_STATIC_ARRAY_LENGTH(HPDDMVariant), HPDDMVariant[HPDDM_VARIANT_LEFT], &i, NULL);
 47:       else if (ksp->pc_side_set == PC_RIGHT) i = HPDDM_VARIANT_RIGHT;
 48:       data->cntl[1] = i;
 49:       if (i > 0) KSPSetPCSide(ksp, PC_RIGHT);
 50:     }
 51:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
 52:       data->rcntl[0] = (PetscAbsReal(data->rcntl[0] - static_cast<PetscReal>(PETSC_DECIDE)) < PETSC_SMALL ? -1.0 : data->rcntl[0]);
 53:       PetscOptionsReal("-ksp_hpddm_deflation_tol", "Tolerance when deflating right-hand sides inside block methods", "KSPHPDDM", data->rcntl[0], data->rcntl, NULL);
 54:       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]));
 55:       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);
 56:       data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG] = i;
 57:     } else data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BCG] = 0;
 58:     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) {
 59:       i = (data->cntl[2] == static_cast<char>(PETSC_DECIDE) ? HPDDM_ORTHOGONALIZATION_CGS : data->cntl[2] & 3);
 60:       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);
 61:       j = (data->cntl[2] == static_cast<char>(PETSC_DECIDE) ? HPDDM_QR_CHOLQR : ((data->cntl[2] >> 2) & 7));
 62:       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);
 63:       data->cntl[2] = static_cast<char>(i) + (static_cast<char>(j) << 2);
 64:       i             = (data->scntl[0] == static_cast<unsigned short>(PETSC_DECIDE) ? PetscMin(30, ksp->max_it) : data->scntl[0]);
 65:       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));
 66:       data->scntl[0] = i;
 67:     }
 68:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BCG || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
 69:       j = (data->cntl[1] == static_cast<char>(PETSC_DECIDE) ? HPDDM_QR_CHOLQR : data->cntl[1]);
 70:       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);
 71:       data->cntl[1] = j;
 72:     }
 73:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
 74:       i = (data->icntl[0] == static_cast<int>(PETSC_DECIDE) ? PetscMin(20, data->scntl[0] - 1) : data->icntl[0]);
 75:       PetscOptionsRangeInt("-ksp_hpddm_recycle", "Number of harmonic Ritz vectors to compute", "KSPHPDDM", i, &i, NULL, 1, data->scntl[0] - 1);
 76:       data->icntl[0] = i;
 77:       if (!PetscDefined(HAVE_SLEPC) || !PetscDefined(USE_SHARED_LIBRARIES) || data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR) {
 78:         i = (data->cntl[3] == static_cast<char>(PETSC_DECIDE) ? HPDDM_RECYCLE_TARGET_SM : data->cntl[3]);
 79:         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);
 80:         data->cntl[3] = i;
 81:       } else {
 83:         MPI_Comm_size(PetscObjectComm((PetscObject)ksp), &size);
 84:         i = (data->cntl[3] == static_cast<char>(PETSC_DECIDE) ? 1 : data->cntl[3]);
 85:         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));
 86:         data->cntl[3] = i;
 87:       }
 88:       i = (data->cntl[4] == static_cast<char>(PETSC_DECIDE) ? HPDDM_RECYCLE_STRATEGY_A : data->cntl[4]);
 89:       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);
 90:       data->cntl[4] = i;
 91:     }
 92:   } else {
 93:     data->cntl[0]  = HPDDM_KRYLOV_METHOD_NONE;
 94:     data->scntl[1] = 1;
 95:   }
 97:              ksp->nmax);
 98:   data->icntl[1] = static_cast<int>(ksp->nmax);
 99:   PetscOptionsHeadEnd();
100:   return 0;
101: }

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

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

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

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

176: static inline PetscErrorCode KSPHPDDMReset_Private(KSP ksp)
177: {
178:   KSP_HPDDM *data = (KSP_HPDDM *)ksp->data;

180:   /* cast PETSC_DECIDE into the appropriate types to avoid compiler warnings */
181:   std::fill_n(data->rcntl, PETSC_STATIC_ARRAY_LENGTH(data->rcntl), static_cast<PetscReal>(PETSC_DECIDE));
182:   std::fill_n(data->icntl, PETSC_STATIC_ARRAY_LENGTH(data->icntl), static_cast<int>(PETSC_DECIDE));
183:   std::fill_n(data->scntl, PETSC_STATIC_ARRAY_LENGTH(data->scntl), static_cast<unsigned short>(PETSC_DECIDE));
184:   std::fill_n(data->cntl, PETSC_STATIC_ARRAY_LENGTH(data->cntl), static_cast<char>(PETSC_DECIDE));
185:   data->precision = PETSC_KSPHPDDM_DEFAULT_PRECISION;
186:   return 0;
187: }

189: static PetscErrorCode KSPReset_HPDDM(KSP ksp)
190: {
191:   KSP_HPDDM *data = (KSP_HPDDM *)ksp->data;

193:   if (data->op) {
194:     delete data->op;
195:     data->op = NULL;
196:   }
197:   KSPHPDDMReset_Private(ksp);
198:   return 0;
199: }

201: static PetscErrorCode KSPDestroy_HPDDM(KSP ksp)
202: {
203:   KSPReset_HPDDM(ksp);
204:   KSPDestroyDefault(ksp);
205:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetDeflationMat_C", NULL);
206:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetDeflationMat_C", NULL);
207:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetType_C", NULL);
208:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetType_C", NULL);
209:   return 0;
210: }

212: static inline PetscErrorCode KSPSolve_HPDDM_Private(KSP ksp, const PetscScalar *b, PetscScalar *x, PetscInt n)
213: {
214:   KSP_HPDDM              *data = (KSP_HPDDM *)ksp->data;
215:   KSPConvergedDefaultCtx *ctx  = (KSPConvergedDefaultCtx *)ksp->cnvP;
216:   const PetscInt          N    = data->op->getDof() * n;
217:   PetscBool               scale;
218: #if !PetscDefined(USE_REAL_DOUBLE) || PetscDefined(HAVE_F2CBLASLAPACK___FLOAT128_BINDINGS)
219:   HPDDM::upscaled_type<PetscScalar> *high[2];
220: #endif
221: #if !PetscDefined(USE_REAL_SINGLE) || PetscDefined(HAVE_F2CBLASLAPACK___FP16_BINDINGS)
222:   HPDDM::downscaled_type<PetscScalar> *low[2];
223: #endif

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

276: static PetscErrorCode KSPSolve_HPDDM(KSP ksp)
277: {
278:   KSP_HPDDM         *data = (KSP_HPDDM *)ksp->data;
279:   Mat                A, B;
280:   PetscScalar       *x, *bt = NULL, **ptr;
281:   const PetscScalar *b;
282:   PetscInt           i, j, n;
283:   PetscBool          flg;

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

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

327:    Input Parameters:
328: +     ksp - iterative context
329: -     U - deflation space to be used during KSPSolve()

331:    Level: intermediate

333: .seealso: [](chapter_ksp), `KSPHPDDM`, `KSPCreate()`, `KSPType`, `KSPHPDDMGetDeflationMat()`
334: @*/
335: PetscErrorCode KSPHPDDMSetDeflationMat(KSP ksp, Mat U)
336: {
340:   PetscUseMethod(ksp, "KSPHPDDMSetDeflationMat_C", (KSP, Mat), (ksp, U));
341:   return 0;
342: }

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

348:    Input Parameter:
349: .     ksp - iterative context

351:    Output Parameter:
352: .     U - deflation space generated during `KSPSolve()`

354:    Level: intermediate

356: .seealso: [](chapter_ksp), `KSPHPDDM`, `KSPCreate()`, `KSPType`, `KSPHPDDMSetDeflationMat()`
357: @*/
358: PetscErrorCode KSPHPDDMGetDeflationMat(KSP ksp, Mat *U)
359: {
361:   if (U) {
363:     PetscUseMethod(ksp, "KSPHPDDMGetDeflationMat_C", (KSP, Mat *), (ksp, U));
364:   }
365:   return 0;
366: }

368: static PetscErrorCode KSPHPDDMSetDeflationMat_HPDDM(KSP ksp, Mat U)
369: {
370:   KSP_HPDDM            *data = (KSP_HPDDM *)ksp->data;
371:   HPDDM::PETScOperator *op   = data->op;
372:   Mat                   A;
373:   const PetscScalar    *array;
374:   PetscScalar          *copy;
375:   PetscInt              m1, M1, m2, M2, n2, N2, ldu;
376:   PetscBool             match;

378:   if (!op) {
379:     KSPSetUp(ksp);
380:     op = data->op;
381:   }
383:   KSPGetOperators(ksp, &A, NULL);
384:   MatGetLocalSize(A, &m1, NULL);
385:   MatGetLocalSize(U, &m2, &n2);
386:   MatGetSize(A, &M1, NULL);
387:   MatGetSize(U, &M2, &N2);
389:   PetscObjectTypeCompareAny((PetscObject)U, &match, MATSEQDENSE, MATMPIDENSE, "");
391:   MatDenseGetArrayRead(U, &array);
392:   copy = op->allocate(m2, 1, N2);
394:   MatDenseGetLDA(U, &ldu);
395:   HPDDM::Wrapper<PetscScalar>::omatcopy<'N'>(N2, m2, array, ldu, copy, m2);
396:   MatDenseRestoreArrayRead(U, &array);
397:   return 0;
398: }

400: static PetscErrorCode KSPHPDDMGetDeflationMat_HPDDM(KSP ksp, Mat *U)
401: {
402:   KSP_HPDDM            *data = (KSP_HPDDM *)ksp->data;
403:   HPDDM::PETScOperator *op   = data->op;
404:   Mat                   A;
405:   const PetscScalar    *array;
406:   PetscScalar          *copy;
407:   PetscInt              m1, M1, N2;

409:   if (!op) {
410:     KSPSetUp(ksp);
411:     op = data->op;
412:   }
414:   array = op->storage();
415:   N2    = op->k().first * op->k().second;
416:   if (!array) *U = NULL;
417:   else {
418:     KSPGetOperators(ksp, &A, NULL);
419:     MatGetLocalSize(A, &m1, NULL);
420:     MatGetSize(A, &M1, NULL);
421:     MatCreateDense(PetscObjectComm((PetscObject)ksp), m1, PETSC_DECIDE, M1, N2, NULL, U);
422:     MatDenseGetArrayWrite(*U, &copy);
423:     PetscArraycpy(copy, array, m1 * N2);
424:     MatDenseRestoreArrayWrite(*U, &copy);
425:   }
426:   return 0;
427: }

429: static PetscErrorCode KSPMatSolve_HPDDM(KSP ksp, Mat B, Mat X)
430: {
431:   KSP_HPDDM            *data = (KSP_HPDDM *)ksp->data;
432:   HPDDM::PETScOperator *op   = data->op;
433:   Mat                   A;
434:   const PetscScalar    *b;
435:   PetscScalar          *x;
436:   PetscInt              n, lda;

438:   PetscCitationsRegister(HPDDMCitation, &HPDDMCite);
439:   if (!op) {
440:     KSPSetUp(ksp);
441:     op = data->op;
442:   }
443:   KSPGetOperators(ksp, &A, NULL);
444:   MatGetLocalSize(B, &n, NULL);
445:   MatDenseGetLDA(B, &lda);
447:   MatGetLocalSize(A, &n, NULL);
448:   MatDenseGetLDA(X, &lda);
450:   MatDenseGetArrayRead(B, &b);
451:   MatDenseGetArrayWrite(X, &x);
452:   MatGetSize(X, NULL, &n);
453:   KSPSolve_HPDDM_Private(ksp, b, x, n);
454:   MatDenseRestoreArrayWrite(X, &x);
455:   MatDenseRestoreArrayRead(B, &b);
456:   return 0;
457: }

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

462:    Collective

464:    Input Parameters:
465: +     ksp - iterative context
466: -     type - any of gmres, bgmres, cg, bcg, gcrodr, bgcrodr, bfbcg, or preonly

468:    Level: intermediate

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

473:      As an example, in the following sequence:
474: .vb
475:      KSPHPDDMSetType(ksp, KSPGCRODR);
476:      KSPSolve(ksp, b, x);
477:      KSPHPDDMSetType(ksp, KSPGMRES);
478:      KSPHPDDMSetType(ksp, KSPGCRODR);
479:      KSPSolve(ksp, b, x);
480: .ve
481:     the recycled space is reused in the second `KSPSolve()`.

483: .seealso: [](chapter_ksp), `KSPCreate()`, `KSPType`, `KSPHPDDMType`, `KSPHPDDMGetType()`
484: @*/
485: PetscErrorCode KSPHPDDMSetType(KSP ksp, KSPHPDDMType type)
486: {
489:   PetscUseMethod(ksp, "KSPHPDDMSetType_C", (KSP, KSPHPDDMType), (ksp, type));
490:   return 0;
491: }

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

496:    Input Parameter:
497: .     ksp - iterative context

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

502:    Level: intermediate

504: .seealso: [](chapter_ksp), `KSPCreate()`, `KSPType`, `KSPHPDDMType`, `KSPHPDDMSetType()`
505: @*/
506: PetscErrorCode KSPHPDDMGetType(KSP ksp, KSPHPDDMType *type)
507: {
509:   if (type) {
511:     PetscUseMethod(ksp, "KSPHPDDMGetType_C", (KSP, KSPHPDDMType *), (ksp, type));
512:   }
513:   return 0;
514: }

516: static PetscErrorCode KSPHPDDMSetType_HPDDM(KSP ksp, KSPHPDDMType type)
517: {
518:   KSP_HPDDM *data = (KSP_HPDDM *)ksp->data;
519:   PetscInt   i;
520:   PetscBool  flg = PETSC_FALSE;

522:   for (i = 0; i < static_cast<PetscInt>(PETSC_STATIC_ARRAY_LENGTH(KSPHPDDMTypes)); ++i) {
523:     PetscStrcmp(KSPHPDDMTypes[type], KSPHPDDMTypes[i], &flg);
524:     if (flg) break;
525:   }
527:   if (data->cntl[0] != static_cast<char>(PETSC_DECIDE) && data->cntl[0] != i) KSPHPDDMReset_Private(ksp);
528:   data->cntl[0] = i;
529:   return 0;
530: }

532: static PetscErrorCode KSPHPDDMGetType_HPDDM(KSP ksp, KSPHPDDMType *type)
533: {
534:   KSP_HPDDM *data = (KSP_HPDDM *)ksp->data;

537:   /* need to shift by -1 for HPDDM_KRYLOV_METHOD_NONE */
538:   *type = static_cast<KSPHPDDMType>(PetscMin(data->cntl[0], static_cast<char>(PETSC_STATIC_ARRAY_LENGTH(KSPHPDDMTypes) - 1)));
539:   return 0;
540: }

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

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

564:    Level: intermediate

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

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

579: PETSC_EXTERN PetscErrorCode KSPCreate_HPDDM(KSP ksp)
580: {
581:   KSP_HPDDM  *data;
582:   PetscInt    i;
583:   const char *common[] = {KSPGMRES, KSPCG, KSPPREONLY};
584:   PetscBool   flg      = PETSC_FALSE;

586:   PetscNew(&data);
587:   ksp->data = (void *)data;
588:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2);
589:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 1);
590:   ksp->ops->solve          = KSPSolve_HPDDM;
591:   ksp->ops->matsolve       = KSPMatSolve_HPDDM;
592:   ksp->ops->setup          = KSPSetUp_HPDDM;
593:   ksp->ops->setfromoptions = KSPSetFromOptions_HPDDM;
594:   ksp->ops->destroy        = KSPDestroy_HPDDM;
595:   ksp->ops->view           = KSPView_HPDDM;
596:   ksp->ops->reset          = KSPReset_HPDDM;
597:   KSPHPDDMReset_Private(ksp);
598:   for (i = 0; i < static_cast<PetscInt>(PETSC_STATIC_ARRAY_LENGTH(common)); ++i) {
599:     PetscStrcmp(((PetscObject)ksp)->type_name, common[i], &flg);
600:     if (flg) break;
601:   }
602:   if (!i) data->cntl[0] = HPDDM_KRYLOV_METHOD_GMRES;
603:   else if (i == 1) data->cntl[0] = HPDDM_KRYLOV_METHOD_CG;
604:   else if (i == 2) data->cntl[0] = HPDDM_KRYLOV_METHOD_NONE;
605:   if (data->cntl[0] != static_cast<char>(PETSC_DECIDE)) PetscInfo(ksp, "Using the previously set KSPType %s\n", common[i]);
606:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetDeflationMat_C", KSPHPDDMSetDeflationMat_HPDDM);
607:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetDeflationMat_C", KSPHPDDMGetDeflationMat_HPDDM);
608:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetType_C", KSPHPDDMSetType_HPDDM);
609:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetType_C", KSPHPDDMGetType_HPDDM);
610: #if defined(PETSC_HAVE_SLEPC) && PetscDefined(HAVE_DYNAMIC_LIBRARIES) && defined(PETSC_USE_SHARED_LIBRARIES)
611:   if (!loadedDL) HPDDMLoadDL_Private(&loadedDL);
612: #endif
613:   data->precision = PETSC_KSPHPDDM_DEFAULT_PRECISION;
614:   return 0;
615: }