Actual source code: hpddm.cxx

petsc-master 2020-07-09
Report Typos and Errors
  1:  #include <petsc/private/petschpddm.h>
  2: /* access to same_local_solves */
  3:  #include <../src/ksp/pc/impls/asm/asm.h>

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

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

 15: static PetscBool citeKSP = PETSC_FALSE;
 16: static const char hpddmCitationKSP[] = "@inproceedings{jolivet2016block,\n\tTitle = {{Block Iterative Methods and Recycling for Improved Scalability of Linear Solvers}},\n\tAuthor = {Jolivet, Pierre and Tournier, Pierre-Henri},\n\tOrganization = {IEEE},\n\tYear = {2016},\n\tSeries = {SC16},\n\tBooktitle = {Proceedings of the 2016 International Conference for High Performance Computing, Networking, Storage and Analysis}\n}\n";

 18: #if defined(PETSC_HAVE_SLEPC) && defined(PETSC_USE_SHARED_LIBRARIES)
 19: static PetscBool loadedDL = PETSC_FALSE;
 20: #endif

 22: static PetscErrorCode KSPSetFromOptions_HPDDM(PetscOptionItems *PetscOptionsObject, KSP ksp)
 23: {
 24:   KSP_HPDDM      *data = (KSP_HPDDM*)ksp->data;
 25:   PetscInt       i, j;
 26:   PetscMPIInt    size;

 30:   PetscOptionsHead(PetscOptionsObject, "KSPHPDDM options, cf. https://github.com/hpddm/hpddm");
 31:   i = (data->cntl[0] == static_cast<char>(PETSC_DECIDE) ? HPDDM_KRYLOV_METHOD_GMRES : data->cntl[0]);
 32:   PetscOptionsEList("-ksp_hpddm_type", "Type of Krylov method", "KSPHPDDM", KSPHPDDMTypes, ALEN(KSPHPDDMTypes), KSPHPDDMTypes[HPDDM_KRYLOV_METHOD_GMRES], &i, NULL);
 33:   data->cntl[0] = i;
 34:   /* cannot use HPDDM_KRYLOV_METHOD_NONE because HPDDM_KRYLOV_METHOD_RICHARDSON is not registered in PETSc */
 35:   if (data->cntl[0] != 7) {
 36:     if (data->cntl[0] != HPDDM_KRYLOV_METHOD_BCG && data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG) {
 37:       i = (data->cntl[1] == static_cast<char>(PETSC_DECIDE) ? HPDDM_VARIANT_LEFT : data->cntl[1]);
 38:       if (ksp->pc_side_set == PC_SIDE_DEFAULT) {
 39:         PetscOptionsEList("-ksp_hpddm_variant", "Left, right, or variable preconditioning", "KSPHPDDM", HPDDMVariant, ALEN(HPDDMVariant), HPDDMVariant[HPDDM_VARIANT_LEFT], &i, NULL);
 40:       } else if (ksp->pc_side_set == PC_RIGHT) i = HPDDM_VARIANT_RIGHT;
 41:       else if (ksp->pc_side_set == PC_SYMMETRIC) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Symmetric preconditioning not implemented");
 42:       data->cntl[1] = i;
 43:       if (i > 0) {
 44:         KSPSetPCSide(ksp, PC_RIGHT);
 45:       }
 46:     }
 47:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
 48:       data->rcntl[0] = (std::abs(data->rcntl[0] - static_cast<PetscReal>(PETSC_DECIDE)) < PETSC_SMALL ? -1.0 : data->rcntl[0]);
 49:       PetscOptionsReal("-ksp_hpddm_deflation_tol", "Tolerance when deflating right-hand sides inside block methods", "KSPHPDDM", data->rcntl[0], data->rcntl, NULL);
 50:       i = (data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG] == static_cast<unsigned short>(PETSC_DECIDE) ? 1 : data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG]);
 51:       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);
 52:       data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BFBCG] = i;
 53:     } else data->scntl[data->cntl[0] != HPDDM_KRYLOV_METHOD_BCG] = 0;
 54:     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) {
 55:       i = (data->cntl[2] == static_cast<char>(PETSC_DECIDE) ? HPDDM_ORTHOGONALIZATION_CGS : data->cntl[2] & 3);
 56:       PetscOptionsEList("-ksp_hpddm_orthogonalization", "Classical (faster) or Modified (more robust) Gram--Schmidt process", "KSPHPDDM", HPDDMOrthogonalization, ALEN(HPDDMOrthogonalization), HPDDMOrthogonalization[HPDDM_ORTHOGONALIZATION_CGS], &i, NULL);
 57:       j = (data->cntl[2] == static_cast<char>(PETSC_DECIDE) ? HPDDM_QR_CHOLQR : ((data->cntl[2] >> 2) & 7));
 58:       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);
 59:       data->cntl[2] = static_cast<char>(i) + (static_cast<char>(j) << 2);
 60:       i = (data->scntl[0] == static_cast<unsigned short>(PETSC_DECIDE) ? PetscMin(30, ksp->max_it) : data->scntl[0]);
 61:       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));
 62:       data->scntl[0] = i;
 63:     }
 64:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_BCG || data->cntl[0] == HPDDM_KRYLOV_METHOD_BFBCG) {
 65:       j = (data->cntl[1] == static_cast<char>(PETSC_DECIDE) ? HPDDM_QR_CHOLQR : data->cntl[1]);
 66:       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);
 67:       data->cntl[1] = j;
 68:     }
 69:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
 70:       i = (data->icntl[0] == static_cast<int>(PETSC_DECIDE) ? PetscMin(20, data->scntl[0] - 1) : data->icntl[0]);
 71:       PetscOptionsRangeInt("-ksp_hpddm_recycle", "Number of harmonic Ritz vectors to compute", "KSPHPDDM", i, &i, NULL, 1, data->scntl[0] - 1);
 72:       data->icntl[0] = i;
 73:       if (!PetscDefined(HAVE_SLEPC) || !PetscDefined(USE_SHARED_LIBRARIES) || data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR) {
 74:         i = (data->cntl[3] == static_cast<char>(PETSC_DECIDE) ? HPDDM_RECYCLE_TARGET_SM : data->cntl[3]);
 75:         PetscOptionsEList("-ksp_hpddm_recycle_target", "Criterion to select harmonic Ritz vectors", "KSPHPDDM", HPDDMRecycleTarget, ALEN(HPDDMRecycleTarget), HPDDMRecycleTarget[HPDDM_RECYCLE_TARGET_SM], &i, NULL);
 76:         data->cntl[3] = i;
 77:       } else {
 78:         MPI_Comm_size(PetscObjectComm((PetscObject)ksp), &size);
 79:         i = (data->cntl[3] == static_cast<char>(PETSC_DECIDE) ? 1 : data->cntl[3]);
 80:         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));
 81:         data->cntl[3] = i;
 82:       }
 83:       i = (data->cntl[4] == static_cast<char>(PETSC_DECIDE) ? HPDDM_RECYCLE_STRATEGY_A : data->cntl[4]);
 84:       PetscOptionsEList("-ksp_hpddm_recycle_strategy", "Generalized eigenvalue problem to solve for recycling", "KSPHPDDM", HPDDMRecycleStrategy, ALEN(HPDDMRecycleStrategy), HPDDMRecycleStrategy[HPDDM_RECYCLE_STRATEGY_A], &i, NULL);
 85:       data->cntl[4] = i;
 86:     }
 87:   } else data->cntl[0] = HPDDM_KRYLOV_METHOD_NONE;
 88:   PetscOptionsTail();
 89:   return(0);
 90: }

 92: static PetscErrorCode KSPView_HPDDM(KSP ksp, PetscViewer viewer)
 93: {
 94:   KSP_HPDDM            *data = (KSP_HPDDM*)ksp->data;
 95:   HPDDM::PETScOperator *op = data->op;
 96:   const PetscScalar    *array = op ? op->storage() : NULL;
 97:   PetscBool            ascii;
 98:   PetscErrorCode       ierr;

101:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii);
102:   if (op && ascii) {
103:     PetscViewerASCIIPrintf(viewer, "HPDDM type: %s\n", KSPHPDDMTypes[std::min(static_cast<PetscInt>(data->cntl[0]), static_cast<PetscInt>(ALEN(KSPHPDDMTypes) - 1))]);
104:     if (data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) {
105:       PetscViewerASCIIPrintf(viewer, "deflation subspace attached? %s\n", PetscBools[array ? PETSC_TRUE : PETSC_FALSE]);
106:       if (!PetscDefined(HAVE_SLEPC) || !PetscDefined(USE_SHARED_LIBRARIES) || data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR) {
107:         PetscViewerASCIIPrintf(viewer, "deflation target: %s\n", HPDDMRecycleTarget[static_cast<PetscInt>(data->cntl[3])]);
108:       } else {
109:         PetscViewerASCIIPrintf(viewer, "redistribution size: %D\n", static_cast<PetscInt>(data->cntl[3]));
110:       }
111:     }
112:     if (data->icntl[1] != static_cast<int>(PETSC_DECIDE)) {
113:       PetscViewerASCIIPrintf(viewer, "  block size is %d\n", data->icntl[1]);
114:     }
115:   }
116:   return(0);
117: }

119: static PetscErrorCode KSPSetUp_HPDDM(KSP ksp)
120: {
121:   KSP_HPDDM      *data = (KSP_HPDDM*)ksp->data;
122:   Mat            A;
123:   PetscInt       n, bs;
124:   PetscBool      match;

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

180: PETSC_STATIC_INLINE PetscErrorCode KSPHPDDMReset_Private(KSP ksp)
181: {
182:   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:   return(0);
191: }

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

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

207: static PetscErrorCode KSPDestroy_HPDDM(KSP ksp)
208: {

212:   KSPReset_HPDDM(ksp);
213:   KSPDestroyDefault(ksp);
214:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetDeflationSpace_C", NULL);
215:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetDeflationSpace_C", NULL);
216:   PetscObjectComposeFunction((PetscObject)ksp, "KSPSetMatSolveBlockSize_C", NULL);
217:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGetMatSolveBlockSize_C", NULL);
218:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetType_C", NULL);
219:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetType_C", NULL);
220:   return(0);
221: }

223: PETSC_STATIC_INLINE PetscErrorCode KSPSolve_HPDDM_Private(KSP ksp, const PetscScalar *b, PetscScalar *x, PetscInt n)
224: {
225:   KSP_HPDDM              *data = (KSP_HPDDM*)ksp->data;
226:   KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*)ksp->cnvP;
227:   PetscBool              scale;
228:   PetscErrorCode         ierr;

231:   PCGetDiagonalScale(ksp->pc, &scale);
232:   if (scale) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
233:   if (n > 1) {
234:     if (ksp->converged == KSPConvergedDefault) {
235:       if (ctx->mininitialrtol) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support KSPConvergedDefaultSetUMIRNorm()", ((PetscObject)ksp)->type_name);
236:       if (!ctx->initialrtol) {
237:         PetscInfo(ksp, "Forcing KSPConvergedDefaultSetUIRNorm() since KSPConvergedDefault() cannot handle multiple norms\n");
238:         ctx->initialrtol = PETSC_TRUE;
239:       }
240:     } else {
241:       PetscInfo(ksp, "Using a special \"converged\" callback, be careful, it is used in KSPHPDDM to track blocks of residuals\n");
242:     }
243:   }
244:   /* initial guess is always nonzero with recycling methods if there is a deflation subspace available */
245:   if ((data->cntl[0] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[0] == HPDDM_KRYLOV_METHOD_BGCRODR) && data->op->storage()) ksp->guess_zero = PETSC_FALSE;
246:   ksp->its = 0;
247:   ksp->reason = KSP_CONVERGED_ITERATING;
248:   static_cast<PetscErrorCode>(HPDDM::IterativeMethod::solve(*data->op, b, x, n, PetscObjectComm((PetscObject)ksp)));
249:   if (!ksp->reason) { /* KSPConvergedDefault() is still returning 0 (= KSP_CONVERGED_ITERATING) */
250:     if (ksp->its >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
251:     else ksp->reason = KSP_CONVERGED_RTOL; /* early exit by HPDDM, which only happens on breakdowns or convergence */
252:   }
253:   ksp->its = PetscMin(ksp->its, ksp->max_it);
254:   return(0);
255: }

257: static PetscErrorCode KSPSolve_HPDDM(KSP ksp)
258: {
259:   KSP_HPDDM         *data = (KSP_HPDDM*)ksp->data;
260:   Mat               A, B;
261:   PetscScalar       *x, *bt = NULL, **ptr;
262:   const PetscScalar *b;
263:   PetscInt          i, j, n;
264:   PetscBool         flg;
265:   PetscErrorCode    ierr;

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

315: /*@
316:      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).

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

322:    Level: intermediate

324: .seealso:  KSPCreate(), KSPType (for list of available types), KSPHPDDMGetDeflationSpace()
325: @*/
326: PetscErrorCode KSPHPDDMSetDeflationSpace(KSP ksp, Mat U)
327: {

334:   PetscUseMethod(ksp, "KSPHPDDMSetDeflationSpace_C", (KSP, Mat), (ksp, U));
335:   return(0);
336: }

338: /*@
339:      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.

341:    Input Parameter:
342: .     ksp - iterative context

344:    Output Parameter:
345: .     U - deflation space generated during KSPSolve()

347:    Level: intermediate

349: .seealso:  KSPCreate(), KSPType (for list of available types), KSPHPDDMSetDeflationSpace()
350: @*/
351: PetscErrorCode KSPHPDDMGetDeflationSpace(KSP ksp, Mat *U)
352: {

357:   PetscUseMethod(ksp, "KSPHPDDMGetDeflationSpace_C", (KSP, Mat*), (ksp, U));
358:   return(0);
359: }

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

373:   if (!op) {
374:     KSPSetUp(ksp);
375:     op = data->op;
376:   }
377:   KSPGetOperators(ksp, &A, NULL);
378:   MatGetLocalSize(A, &m1, NULL);
379:   MatGetLocalSize(U, &m2, &n2);
380:   MatGetSize(A, &M1, NULL);
381:   MatGetSize(U, &M2, &N2);
382:   if (m1 != m2 || M1 != M2) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot use a deflation space with (m2,M2) = (%D,%D) for a linear system with (m1,M1) = (%D,%D)", m2, M2, m1, M1);
383:   PetscObjectTypeCompareAny((PetscObject)U, &match, MATSEQDENSE, MATMPIDENSE, "");
384:   if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided deflation space not stored in a dense Mat");
385:   MatDenseGetArrayRead(U, &array);
386:   copy = op->allocate(m2, 1, N2);
387:   if (!copy) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_POINTER, "Memory allocation error");
388:   MatDenseGetLDA(U, &ldu);
389:   HPDDM::Wrapper<PetscScalar>::omatcopy<'N'>(N2, m2, array, ldu, copy, m2);
390:   MatDenseRestoreArrayRead(U, &array);
391:   return(0);
392: }

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

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

424: static PetscErrorCode KSPMatSolve_HPDDM(KSP ksp, Mat B, Mat X)
425: {
426:   KSP_HPDDM            *data = (KSP_HPDDM*)ksp->data;
427:   PC                   pc;
428:   PC_ASM               *osm = NULL;
429:   HPDDM::PETScOperator *op = data->op;
430:   Mat                  A;
431:   const PetscScalar    *b;
432:   PetscScalar          *x;
433:   PetscInt             n, lda;
434:   PetscBool            same_local_solves = PETSC_FALSE;
435:   PetscErrorCode       ierr;

438:   if (!op) {
439:     KSPSetUp(ksp);
440:     op = data->op;
441:   }
442:   KSPGetOperators(ksp, &A, NULL);
443:   MatGetLocalSize(B, &n, NULL);
444:   MatDenseGetLDA(B, &lda);
445:   if (n != lda) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Unhandled leading dimension lda = %D with n = %D", lda, n);
446:   MatGetLocalSize(A, &n, NULL);
447:   MatDenseGetLDA(X, &lda);
448:   if (n != lda) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Unhandled leading dimension lda = %D with n = %D", lda, n);
449:   MatDenseGetArrayRead(B, &b);
450:   MatDenseGetArrayWrite(X, &x);
451:   KSPGetPC(ksp, &pc);
452:   /* in HPDDM, if PCASM is used, a call to PCASMGetSubKSP() is made to know if there is a single subsolver */
453:   PetscObjectTypeCompare((PetscObject)pc, PCASM, &same_local_solves);
454:   if (same_local_solves) {
455:     osm = (PC_ASM*)pc->data;
456:     same_local_solves = osm->same_local_solves;
457:   }
458:   MatGetSize(X, NULL, &n);
459:   KSPSolve_HPDDM_Private(ksp, b, x, n);
460:   /* if the PetscBool same_local_solves is not reset after the solve, KSPView() is way too verbose */
461:   if (same_local_solves) osm->same_local_solves = PETSC_TRUE;
462:   MatDenseRestoreArrayWrite(X, &x);
463:   MatDenseRestoreArrayRead(B, &b);
464:   return(0);
465: }

467: static PetscErrorCode KSPSetMatSolveBlockSize_HPDDM(KSP ksp, PetscInt bs)
468: {
469:   KSP_HPDDM *data = (KSP_HPDDM*)ksp->data;

472:   if (bs > std::numeric_limits<int>::max() || bs < std::numeric_limits<int>::min()) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_OUTOFRANGE, "KSPMatSolve() block size %D not representable by an integer", bs);
473:   else data->icntl[1] = static_cast<int>(bs);
474:   return(0);
475: }

477: static PetscErrorCode KSPGetMatSolveBlockSize_HPDDM(KSP ksp, PetscInt *bs)
478: {
479:   KSP_HPDDM *data = (KSP_HPDDM*)ksp->data;

482:   *bs = static_cast<PetscInt>(data->icntl[1]);
483:   return(0);
484: }

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

489:    Input Parameters:
490: +     ksp - iterative context
491: -     type - any of gmres, bgmres, cg, bcg, gcrodr, bgcrodr, bfbcg, or preonly

493:    Level: intermediate

495:    Notes:
496:      Unlike KSPReset(), this function does not destroy any deflation space attached to the KSP.
497:      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().

499: .seealso:  KSPCreate(), KSPType (for list of available types), KSPHPDDMGetType()
500: @*/
501: PetscErrorCode KSPHPDDMSetType(KSP ksp, KSPHPDDMType type)
502: {

507:   PetscUseMethod(ksp, "KSPHPDDMSetType_C", (KSP, KSPHPDDMType), (ksp, type));
508:   return(0);
509: }

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

514:    Input Parameter:
515: .     ksp - iterative context

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

520:    Level: intermediate

522: .seealso:  KSPCreate(), KSPType (for list of available types), KSPHPDDMSetType()
523: @*/
524: PetscErrorCode KSPHPDDMGetType(KSP ksp, KSPHPDDMType *type)
525: {

531:   PetscUseMethod(ksp, "KSPHPDDMGetType_C", (KSP, KSPHPDDMType*), (ksp, type));
532:   return(0);
533: }

535: static PetscErrorCode KSPHPDDMSetType_HPDDM(KSP ksp, KSPHPDDMType type)
536: {
537:   KSP_HPDDM      *data = (KSP_HPDDM*)ksp->data;
538:   PetscInt       i;
539:   PetscBool      flg = PETSC_FALSE;

543:   for (i = 0; i < static_cast<PetscInt>(ALEN(KSPHPDDMTypes)); ++i) {
544:     PetscStrcmp(KSPHPDDMTypes[type], KSPHPDDMTypes[i], &flg);
545:     if (flg) break;
546:   }
547:   if (i == ALEN(KSPHPDDMTypes)) SETERRQ1(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown KSPHPDDMType %s", type);
548:   if (data->cntl[0] != static_cast<char>(PETSC_DECIDE) && data->cntl[0] != i) {
549:     KSPHPDDMReset_Private(ksp);
550:   }
551:   data->cntl[0] = i;
552:   return(0);
553: }

555: static PetscErrorCode KSPHPDDMGetType_HPDDM(KSP ksp, KSPHPDDMType *type)
556: {
557:   KSP_HPDDM *data = (KSP_HPDDM*)ksp->data;

560:   if (data->cntl[0] == static_cast<char>(PETSC_DECIDE)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ORDER, "KSPHPDDMType not set yet");
561:   *type = static_cast<KSPHPDDMType>(data->cntl[0]);
562:   return(0);
563: }

565: /*MC
566:      KSPHPDDM - Interface with the HPDDM library.

568:    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.

570:    Options Database Keys:
571: +   -ksp_gmres_restart <restart, default=30> - see KSPGMRES
572: .   -ksp_hpddm_type <type, default=gmres> - any of gmres, bgmres, cg, bcg, gcrodr, bgcrodr, bfbcg, or preonly
573: .   -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)
574: .   -ksp_hpddm_enlarge_krylov_subspace <p, default=1> - split the initial right-hand side into multiple vectors (only relevant with nonblock methods)
575: .   -ksp_hpddm_orthogonalization <type, default=cgs> - any of cgs or mgs, see KSPGMRES
576: .   -ksp_hpddm_qr <type, default=cholqr> - distributed QR factorizations with any of cholqr, cgs, or mgs (only relevant with block methods)
577: .   -ksp_hpddm_variant <type, default=left> - any of left, right, or flexible (this option is superseded by KSPSetPCSide())
578: .   -ksp_hpddm_recycle <n, default=0> - number of harmonic Ritz vectors to compute (only relevant with GCRODR or BGCRODR)
579: .   -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_
580: -   -ksp_hpddm_recycle_strategy <type, default=A> - generalized eigenvalue problem A or B to solve for recycling (only relevant with flexible GCRODR or BGCRODR)

582:    References:
583: +   1980 - The Block Conjugate Gradient Algorithm and Related Methods. O'Leary. Linear Algebra and its Applications.
584: .   2006 - Recycling Krylov Subspaces for Sequences of Linear Systems. Parks, de Sturler, Mackey, Johnson, and Maiti. SIAM Journal on Scientific Computing
585: .   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.
586: .   2016 - Block Iterative Methods and Recycling for Improved Scalability of Linear Solvers. Jolivet and Tournier. SC16.
587: -   2017 - A breakdown-free block conjugate gradient method. Ji and Li. BIT Numerical Mathematics.

589:    Level: intermediate

591: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPGMRES, KSPCG, KSPLGMRES, KSPDGMRES
592: M*/
593: PETSC_EXTERN PetscErrorCode KSPCreate_HPDDM(KSP ksp)
594: {
595:   KSP_HPDDM      *data;
596:   PetscInt       i;
597:   const char     *common[] = { KSPGMRES, KSPCG, KSPPREONLY };
598:   PetscBool      flg = PETSC_FALSE;

602:   PetscNewLog(ksp, &data);
603:   ksp->data = (void*)data;
604:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2);
605:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 1);
606:   ksp->ops->solve          = KSPSolve_HPDDM;
607:   ksp->ops->matsolve       = KSPMatSolve_HPDDM;
608:   ksp->ops->setup          = KSPSetUp_HPDDM;
609:   ksp->ops->setfromoptions = KSPSetFromOptions_HPDDM;
610:   ksp->ops->destroy        = KSPDestroy_HPDDM;
611:   ksp->ops->view           = KSPView_HPDDM;
612:   ksp->ops->reset          = KSPReset_HPDDM;
613:   KSPHPDDMReset_Private(ksp);
614:   for (i = 0; i < static_cast<PetscInt>(ALEN(common)); ++i) {
615:     PetscStrcmp(((PetscObject)ksp)->type_name, common[i], &flg);
616:     if (flg) break;
617:   }
618:   if (!i) data->cntl[0] = HPDDM_KRYLOV_METHOD_GMRES;
619:   else if (i == 1) data->cntl[0] = HPDDM_KRYLOV_METHOD_CG;
620:   else if (i == 2) data->cntl[0] = 7; /* cannot use HPDDM_KRYLOV_METHOD_NONE because HPDDM_KRYLOV_METHOD_RICHARDSON is not registered in PETSc */
621:   if (data->cntl[0] != static_cast<char>(PETSC_DECIDE)) {
622:     PetscInfo1(ksp, "Using the previously set KSPType %s\n", common[i]);
623:   }
624:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetDeflationSpace_C", KSPHPDDMSetDeflationSpace_HPDDM);
625:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetDeflationSpace_C", KSPHPDDMGetDeflationSpace_HPDDM);
626:   PetscObjectComposeFunction((PetscObject)ksp, "KSPSetMatSolveBlockSize_C", KSPSetMatSolveBlockSize_HPDDM);
627:   PetscObjectComposeFunction((PetscObject)ksp, "KSPGetMatSolveBlockSize_C", KSPGetMatSolveBlockSize_HPDDM);
628:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetType_C", KSPHPDDMSetType_HPDDM);
629:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetType_C", KSPHPDDMGetType_HPDDM);
630: #if defined(PETSC_HAVE_SLEPC) && defined(PETSC_USE_SHARED_LIBRARIES)
631:   if (!loadedDL) {
632:     HPDDMLoadDL_Private(&loadedDL);
633:   }
634: #endif
635:   return(0);
636: }