Actual source code: hpddm.cxx

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

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

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

 16: static PetscBool citeKSP = PETSC_FALSE;
 17: 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";

 19: static PetscErrorCode KSPSetFromOptions_HPDDM(PetscOptionItems *PetscOptionsObject, KSP ksp)
 20: {
 21:   KSP_HPDDM      *data = (KSP_HPDDM*)ksp->data;
 22:   PetscInt       i, j;

 26:   data->scntl[0] = ksp->max_it;
 27:   data->rcntl[0] = ksp->rtol;
 28:   PetscOptionsHead(PetscOptionsObject, "KSPHPDDM options, cf. https://github.com/hpddm/hpddm");
 29:   i = HPDDM_KRYLOV_METHOD_GMRES;
 30:   PetscOptionsEList("-ksp_hpddm_type", "Type of Krylov method", "KSPHPDDM", HPDDMType, ALEN(HPDDMType), HPDDMType[HPDDM_KRYLOV_METHOD_GMRES], &i, NULL);
 31:   data->cntl[5] = i;
 32:   if (data->cntl[5] == HPDDM_KRYLOV_METHOD_RICHARDSON) {
 33:     data->rcntl[0] = 1.0;
 34:     PetscOptionsReal("-ksp_richardson_scale", "Damping factor used in Richardson iterations", "KSPHPDDM", data->rcntl[0], data->rcntl, NULL);
 35:   } else {
 36:     i = HPDDM_VARIANT_LEFT;
 37:     if (ksp->pc_side_set == PC_SIDE_DEFAULT) {
 38:       PetscOptionsEList("-ksp_hpddm_variant", "Left, right, or variable preconditioning", "KSPHPDDM", HPDDMVariant, ALEN(HPDDMVariant), HPDDMVariant[HPDDM_VARIANT_LEFT], &i, NULL);
 39:     } else if (ksp->pc_side_set == PC_RIGHT) i = HPDDM_VARIANT_RIGHT;
 40:     else if (ksp->pc_side_set == PC_SYMMETRIC) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Symmetric preconditioning not implemented");
 41:     if (i != HPDDM_VARIANT_LEFT && (data->cntl[5] == HPDDM_KRYLOV_METHOD_BCG || data->cntl[5] == HPDDM_KRYLOV_METHOD_BFBCG)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Right and flexible preconditioned (BF)BCG not implemented");
 42:     data->cntl[1] = i;
 43:     if (i > 0) {
 44:       KSPSetPCSide(ksp, PC_RIGHT);
 45:     }
 46:     if (data->cntl[5] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[5] == HPDDM_KRYLOV_METHOD_BGCRODR || data->cntl[5] == HPDDM_KRYLOV_METHOD_BFBCG) {
 47:       data->rcntl[1] = -1.0;
 48:       PetscOptionsReal("-ksp_hpddm_deflation_tol", "Tolerance when deflating right-hand sides inside block methods", "KSPHPDDM", data->rcntl[1], data->rcntl + 1, NULL);
 49:       i = 1;
 50:       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);
 51:       data->scntl[1 + (data->cntl[5] != HPDDM_KRYLOV_METHOD_BFBCG)] = i;
 52:     } else data->scntl[2] = 0;
 53:     if (data->cntl[5] == HPDDM_KRYLOV_METHOD_GMRES || data->cntl[5] == HPDDM_KRYLOV_METHOD_BGMRES || data->cntl[5] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[5] == HPDDM_KRYLOV_METHOD_BGCRODR) {
 54:       i = HPDDM_ORTHOGONALIZATION_CGS;
 55:       PetscOptionsEList("-ksp_hpddm_orthogonalization", "Classical (faster) or Modified (more robust) Gram--Schmidt process", "KSPHPDDM", HPDDMOrthogonalization, ALEN(HPDDMOrthogonalization), HPDDMOrthogonalization[HPDDM_ORTHOGONALIZATION_CGS], &i, NULL);
 56:       j = HPDDM_QR_CHOLQR;
 57:       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);
 58:       data->cntl[2] = static_cast<char>(i) + (static_cast<char>(j) << 2);
 59:       i = PetscMin(30, ksp->max_it - 1);
 60:       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));
 61:       data->scntl[1] = i;
 62:     }
 63:     if (data->cntl[5] == HPDDM_KRYLOV_METHOD_BCG || data->cntl[5] == HPDDM_KRYLOV_METHOD_BFBCG) {
 64:       j = HPDDM_QR_CHOLQR;
 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[1] = j;
 67:     }
 68:     if (data->cntl[5] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[5] == HPDDM_KRYLOV_METHOD_BGCRODR) {
 69:       i = PetscMin(20, data->scntl[1] - 1);
 70:       PetscOptionsRangeInt("-ksp_hpddm_recycle", "Number of harmonic Ritz vectors to compute", "KSPHPDDM", i, &i, NULL, 1, data->scntl[1] - 1);
 71:       data->icntl[0] = i;
 72:       i = HPDDM_RECYCLE_TARGET_SM;
 73:       PetscOptionsEList("-ksp_hpddm_recycle_target", "Criterion to select harmonic Ritz vectors", "KSPHPDDM", HPDDMRecycleTarget, ALEN(HPDDMRecycleTarget), HPDDMRecycleTarget[HPDDM_RECYCLE_TARGET_SM], &i, NULL);
 74:       data->cntl[3] = i;
 75:       i = HPDDM_RECYCLE_STRATEGY_A;
 76:       PetscOptionsEList("-ksp_hpddm_recycle_strategy", "Generalized eigenvalue problem to solve for recycling", "KSPHPDDM", HPDDMRecycleStrategy, ALEN(HPDDMRecycleStrategy), HPDDMRecycleStrategy[HPDDM_RECYCLE_STRATEGY_A], &i, NULL);
 77:       data->cntl[4] = i;
 78:     }
 79:   }
 80:   PetscOptionsTail();
 81:   return(0);
 82: }

 84: static PetscErrorCode KSPView_HPDDM(KSP ksp, PetscViewer viewer)
 85: {
 86:   KSP_HPDDM            *data = (KSP_HPDDM*)ksp->data;
 87:   HPDDM::PETScOperator *op = data->op;
 88:   const PetscScalar    *array = op ? op->storage() : NULL;
 89:   PetscBool            ascii;
 90:   PetscErrorCode       ierr;

 93:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &ascii);
 94:   if (op && ascii) {
 95:     PetscViewerASCIIPrintf(viewer, "HPDDM type: %s\n", HPDDMType[static_cast<PetscInt>(data->cntl[5])]);
 96:     if (data->cntl[5] == HPDDM_KRYLOV_METHOD_GCRODR || data->cntl[5] == HPDDM_KRYLOV_METHOD_BGCRODR) {
 97:       PetscViewerASCIIPrintf(viewer, "deflation subspace attached? %s\n", PetscBools[array ? PETSC_TRUE : PETSC_FALSE]);
 98:       PetscViewerASCIIPrintf(viewer, "deflation target: %s\n", HPDDMRecycleTarget[static_cast<PetscInt>(data->cntl[3])]);
 99:     }
100:   }
101:   return(0);
102: }

104: static PetscErrorCode KSPSetUp_HPDDM(KSP ksp)
105: {
106:   KSP_HPDDM      *data = (KSP_HPDDM*)ksp->data;
107:   Mat            A;
108:   PetscInt       n, bs;
109:   PetscBool      match;

113:   KSPGetOperators(ksp, &A, NULL);
114:   MatGetLocalSize(A, &n, NULL);
115:   MatGetBlockSize(A, &bs);
116:   PetscObjectTypeCompareAny((PetscObject)A, &match, MATSEQBAIJ, MATMPIBAIJ, MATSEQSBAIJ, MATMPISBAIJ, "");
117:   /* for block formats, the actual size of the underlying arrays are needed */
118:   if (match) n *= bs;
119:   PetscObjectTypeCompareAny((PetscObject)A, &match, MATSEQKAIJ, MATMPIKAIJ, "");
120:   if (match) n /= bs;
121: #if defined(PETSC_PKG_HPDDM_VERSION_MAJOR)
122: #if PETSC_PKG_HPDDM_VERSION_LT(2, 0, 4)
123:   data->op = new HPDDM::PETScOperator(ksp, n, 1);
124: #else
125:   data->op = new HPDDM::PETScOperator(ksp, n);
126: #endif
127: #else
128:   data->op = new HPDDM::PETScOperator(ksp, n, 1);
129: #endif
130:   return(0);
131: }

133: static PetscErrorCode KSPReset_HPDDM(KSP ksp)
134: {
135:   KSP_HPDDM *data = (KSP_HPDDM*)ksp->data;
137:   if (data->op) {
138:     delete data->op;
139:     data->op = NULL;
140:   }
141:   return(0);
142: }

144: static PetscErrorCode KSPDestroy_HPDDM(KSP ksp)
145: {

149:   KSPReset_HPDDM(ksp);
150:   KSPDestroyDefault(ksp);
151:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetDeflationSpace_C", NULL);
152:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetDeflationSpace_C", NULL);
153:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMMatSolve_C", NULL);
154:   return(0);
155: }

157: PETSC_STATIC_INLINE PetscErrorCode KSPSolve_HPDDM_Private(KSP ksp, const PetscScalar *b, PetscScalar *x, PetscInt n)
158: {
159:   KSP_HPDDM      *data = (KSP_HPDDM*)ksp->data;

163:   static_cast<PetscInt>(HPDDM::IterativeMethod::solve(*data->op, b, x, n, PetscObjectComm((PetscObject)ksp)));
164:   /* big assumption from HPDDM: all PetscErrorCode are positive                                            */
165:   /* if a PETSc call fails inside HPDDM, -ierr is returned (always negative given the previous assumption) */
166:   /* if a KSPSolve succeeds, the number of iterations is returned instead (always positive or null)        */
167:   ksp->its = 0;
168:   if (ierr >= 0) ksp->its = ierr;
169:   else           return PetscError(PETSC_COMM_SELF, __LINE__, "KSPSolve_HPDDM_Private", __FILE__, -ierr, PETSC_ERROR_INITIAL, "PETSc error detected in HPDDM");
170:   return(0);
171: }

173: static PetscErrorCode KSPSolve_HPDDM(KSP ksp)
174: {
175:   KSP_HPDDM         *data = (KSP_HPDDM*)ksp->data;
176:   Mat               A, B;
177:   PetscScalar       *x, *bt = NULL, **ptr;
178:   const PetscScalar *b;
179:   PetscInt          i, j, n;
180:   PetscBool         flg;
181:   PetscErrorCode    ierr;

184:   PetscCitationsRegister(hpddmCitationKSP, &citeKSP);
185:   KSPGetOperators(ksp, &A, NULL);
186:   PetscObjectTypeCompareAny((PetscObject)A, &flg, MATSEQKAIJ, MATMPIKAIJ, "");
187:   VecGetArray(ksp->vec_sol, &x);
188:   VecGetArrayRead(ksp->vec_rhs, &b);
189:   if (!flg) {
190:     KSPSolve_HPDDM_Private(ksp, b, x, 1);
191:   } else {
192:       MatKAIJGetScaledIdentity(A, &flg);
193:       MatKAIJGetAIJ(A, &B);
194:       MatGetBlockSize(A, &n);
195:       MatGetLocalSize(B, &i, NULL);
196:       j = data->op->getDof();
197:       if (!flg) i *= n; /* S and T are not scaled identities, cannot use block methods */
198:       if (i != j) { /* switching between block and standard methods */
199:         delete data->op;
200: #if defined(PETSC_PKG_HPDDM_VERSION_MAJOR)
201: #if PETSC_PKG_HPDDM_VERSION_LT(2, 0, 4)
202:         data->op = new HPDDM::PETScOperator(ksp, i, 1);
203: #else
204:         data->op = new HPDDM::PETScOperator(ksp, i);
205: #endif
206: #else
207:         data->op = new HPDDM::PETScOperator(ksp, i, 1);
208: #endif
209:       }
210:       if (flg && n > 1) {
211:         PetscMalloc1(i * n, &bt);
212:         /* from row- to column-major to be consistent with HPDDM */
213:         HPDDM::Wrapper<PetscScalar>::omatcopy<'T'>(i, n, b, n, bt, i);
214:         ptr = const_cast<PetscScalar**>(&b);
215:         std::swap(*ptr, bt);
216:         HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(i, n, x, n, i);
217:       }
218:       KSPSolve_HPDDM_Private(ksp, b, x, flg ? n : 1);
219:       if (flg && n > 1) {
220:         std::swap(*ptr, bt);
221:         PetscFree(bt);
222:         /* from column- to row-major to be consistent with MatKAIJ format */
223:         HPDDM::Wrapper<PetscScalar>::imatcopy<'T'>(n, i, x, i, n);
224:       }
225:   }
226:   VecRestoreArrayRead(ksp->vec_rhs, &b);
227:   VecRestoreArray(ksp->vec_sol, &x);
228:   if (ksp->its < ksp->max_it) ksp->reason = KSP_CONVERGED_RTOL;
229:   else ksp->reason = KSP_DIVERGED_ITS;
230:   return(0);
231: }

233: /*@
234:      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).

236:    Input Parameters:
237: +     ksp - iterative context
238: -     U - deflation space to be used during KSPSolve()

240:    Level: intermediate

242: .seealso:  KSPCreate(), KSPType (for list of available types), KSPHPDDMGetDeflationSpace()
243: @*/
244: PetscErrorCode KSPHPDDMSetDeflationSpace(KSP ksp, Mat U)
245: {

252:   PetscUseMethod(ksp, "KSPHPDDMSetDeflationSpace_C", (KSP, Mat), (ksp, U));
253:   return(0);
254: }

256: /*@
257:      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.

259:    Input Parameter:
260: .     ksp - iterative context

262:    Output Parameter:
263: .     U - deflation space generated during KSPSolve()

265:    Level: intermediate

267: .seealso:  KSPCreate(), KSPType (for list of available types), KSPHPDDMSetDeflationSpace()
268: @*/
269: PetscErrorCode KSPHPDDMGetDeflationSpace(KSP ksp, Mat *U)
270: {

275:   PetscUseMethod(ksp, "KSPHPDDMGetDeflationSpace_C", (KSP, Mat*), (ksp, U));
276:   return(0);
277: }

279: static PetscErrorCode KSPHPDDMSetDeflationSpace_HPDDM(KSP ksp, Mat U)
280: {
281:   KSP_HPDDM            *data = (KSP_HPDDM*)ksp->data;
282:   HPDDM::PETScOperator *op = data->op;
283:   Mat                  A;
284:   const PetscScalar    *array;
285:   PetscScalar          *copy;
286:   PetscInt             m1, M1, m2, M2, n2, N2, ldu;
287:   PetscBool            match;
288:   PetscErrorCode       ierr;

291:   if (!op) {
292:     KSPSetUp(ksp);
293:     op = data->op;
294:   }
295:   KSPGetOperators(ksp, &A, NULL);
296:   MatGetLocalSize(A, &m1, NULL);
297:   MatGetLocalSize(U, &m2, &n2);
298:   MatGetSize(A, &M1, NULL);
299:   MatGetSize(U, &M2, &N2);
300:   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);
301:   PetscObjectTypeCompareAny((PetscObject)U, &match, MATSEQDENSE, MATMPIDENSE, "");
302:   if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided deflation space not stored in a dense Mat");
303:   MatDenseGetArrayRead(U, &array);
304:   copy = op->allocate(m2, 1, N2);
305:   if (!copy) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_POINTER, "Memory allocation error");
306:   MatDenseGetLDA(U, &ldu);
307:   HPDDM::Wrapper<PetscScalar>::omatcopy<'N'>(N2, m2, array, ldu, copy, m2);
308:   MatDenseRestoreArrayRead(U, &array);
309:   return(0);
310: }

312: static PetscErrorCode KSPHPDDMGetDeflationSpace_HPDDM(KSP ksp, Mat *U)
313: {
314:   KSP_HPDDM            *data = (KSP_HPDDM*)ksp->data;
315:   HPDDM::PETScOperator *op = data->op;
316:   Mat                  A;
317:   const PetscScalar    *array;
318:   PetscScalar          *copy;
319:   PetscInt             m1, M1, N2;
320:   PetscErrorCode       ierr;

323:   if (!op) {
324:     KSPSetUp(ksp);
325:     op = data->op;
326:   }
327:   array = op->storage();
328:   N2 = op->k();
329:   if (!array) *U = NULL;
330:   else {
331:     KSPGetOperators(ksp, &A, NULL);
332:     MatGetLocalSize(A, &m1, NULL);
333:     MatGetSize(A, &M1, NULL);
334:     MatCreateDense(PetscObjectComm((PetscObject)ksp), m1, PETSC_DECIDE, M1, N2, NULL, U);
335:     MatDenseGetArray(*U, &copy);
336:     PetscArraycpy(copy, array, m1 * N2);
337:     MatDenseRestoreArray(*U, &copy);
338:   }
339:   return(0);
340: }

342: /*@
343:      KSPHPDDMMatSolve - Solves a linear system with multiple right-hand sides stored as a MATDENSE. Unlike KSPSolve(), B and X must be different matrices.

345:    Input Parameters:
346: +     ksp - iterative context
347: -     B - block of right-hand sides

349:    Output Parameter:
350: .     X - block of solutions

352:    Level: intermediate

354: .seealso:  KSPSolve(), MatMatSolve(), MATDENSE, PCBJACOBI, PCASM
355: @*/
356: PetscErrorCode KSPHPDDMMatSolve(KSP ksp, Mat B, Mat X)
357: {

364:   PetscUseMethod(ksp, "KSPHPDDMMatSolve_C", (KSP, Mat, Mat), (ksp, B, X));
365:   return(0);
366: }

368: static PetscErrorCode KSPViewFinalMatResidual_Internal(KSP ksp, Mat B, Mat X, PetscViewer viewer, PetscViewerFormat format)
369: {
370:   Mat            A, R;
371:   PetscReal      *norms;
372:   PetscInt       i, N;
373:   PetscBool      flg;

377:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &flg);
378:   if (flg) {
379:     PCGetOperators(ksp->pc, &A, NULL);
380:     MatAssembled(X, &flg);
381:     if (!flg) {
382:         MatSetOption(X, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
383:         MatAssemblyBegin(X, MAT_FINAL_ASSEMBLY);
384:         MatAssemblyEnd(X, MAT_FINAL_ASSEMBLY);
385:     }
386:     /* A and X must be assembled */
387:     MatMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &R);
388:     MatAYPX(R, -1.0, B, SAME_NONZERO_PATTERN);
389:     MatGetSize(R, NULL, &N);
390:     PetscMalloc1(N, &norms);
391:     MatGetColumnNorms(R, NORM_2, norms);
392:     MatDestroy(&R);
393:     for (i = 0; i < N; ++i) {
394:       PetscViewerASCIIPrintf(viewer, "%s #%D %g\n", i == 0 ? "KSP final norm of residual" : "                          ", i, (double)norms[i]);
395:     }
396:     PetscFree(norms);
397:   }
398:   return(0);
399: }

401: static PetscErrorCode KSPHPDDMMatSolve_HPDDM(KSP ksp, Mat B, Mat X)
402: {
403:   KSP_HPDDM            *data = (KSP_HPDDM*)ksp->data;
404:   PC                   pc;
405:   PC_BJacobi           *bjacobi = NULL;
406:   PC_ASM               *osm = NULL;
407:   HPDDM::PETScOperator *op = data->op;
408:   Mat                  A;
409:   Vec                  cb, cx;
410:   const PetscScalar    *b;
411:   PetscScalar          *x;
412:   PetscInt             m1, M1, m2, M2, n1, N1, n2, N2, lda;
413:   PetscBool            match, same_local_solves = PETSC_FALSE;
414:   PetscErrorCode       ierr;

417:   if (!op) {
418:     KSPSetUp(ksp);
419:     op = data->op;
420:   }
421:   if (B == X) SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_IDN, "B and X must be different matrices");
422:   KSPGetOperators(ksp, &A, NULL);
423:   MatGetLocalSize(A, &m1, NULL);
424:   MatGetLocalSize(B, &m2, &n2);
425:   MatGetSize(A, &M1, NULL);
426:   MatGetSize(B, &M2, &N2);
427:   if (m1 != m2 || M1 != M2) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Cannot use a block of right-hand sides with (m2,M2) = (%D,%D) for a linear system with (m1,M1) = (%D,%D)", m2, M2, m1, M1);
428:   MatGetLocalSize(X, &m1, &n1);
429:   MatGetSize(X, &M1, &N1);
430:   if (m1 != m2 || M1 != M2 || n1 != n2 || N1 != N2) SETERRQ8(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Incompatible block of right-hand sides (m2,M2)x(n2,N2) = (%D,%D)x(%D,%D) and solutions (m1,M1)x(n1,N1) = (%D,%D)x(%D,%D)", m2, M2, n2, N2, m1, M1, n1, N1);
431:   PetscObjectTypeCompareAny((PetscObject)B, &match, MATSEQDENSE, MATMPIDENSE, "");
432:   if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of right-hand sides not stored in a dense Mat");
433:   PetscObjectTypeCompareAny((PetscObject)X, &match, MATSEQDENSE, MATMPIDENSE, "");
434:   if (!match) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Provided block of solutions not stored in a dense Mat");
435:   MatDenseGetLDA(B, &lda);
436:   if (m2 != lda) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Unhandled leading dimension lda = %D with m2 = %D", lda, m2);
437:   MatDenseGetLDA(X, &lda);
438:   if (m1 != lda) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Unhandled leading dimension lda = %D with m1 = %D", lda, m1);
439:   MatDenseGetArrayRead(B, &b);
440:   MatDenseGetArray(X, &x);
441:   if (N1 > 1) {
442:     KSPGetPC(ksp, &pc);
443:     /* in HPDDM, if BJacobi or ASM is used, a call to PC[BJacobi|ASM]GetSubKSP() is made   */
444:     /* to know if there is a single subsolver and if it has a MatMatSolve() implementation */
445:     PetscObjectTypeCompare((PetscObject)pc, PCBJACOBI, &same_local_solves);
446:     if (same_local_solves) {
447:       bjacobi = (PC_BJacobi*)pc->data;
448:       same_local_solves = bjacobi->same_local_solves;
449:     }
450:     if (!bjacobi) {
451:       PetscObjectTypeCompare((PetscObject)pc, PCASM, &same_local_solves);
452:       if (same_local_solves) {
453:         osm = (PC_ASM*)pc->data;
454:         same_local_solves = osm->same_local_solves;
455:       }
456:     }
457:     PetscLogEventBegin(KSP_Solve, ksp, 0, 0, 0);
458:     KSPSolve_HPDDM_Private(ksp, b, x, N1);
459:     PetscLogEventEnd(KSP_Solve, ksp, 0, 0, 0);
460:     /* if the PetscBool same_local_solves is not reset after the solve, KSPView() is way too verbose */
461:     if (same_local_solves) {
462:       if (bjacobi) bjacobi->same_local_solves = PETSC_TRUE;
463:       if (osm) osm->same_local_solves = PETSC_TRUE;
464:     }
465:     if (ksp->its < ksp->max_it) ksp->reason = KSP_CONVERGED_RTOL;
466:     else ksp->reason = KSP_DIVERGED_ITS;
467:     /* stripped-down version of KSPSolve(), which only handles -ksp_view -ksp_converged_reason -ksp_view_final_residual */
468:     if (ksp->viewReason) {
469:       KSPReasonView(ksp, ksp->viewerReason);
470:     }
471:     if (ksp->viewFinalRes) {
472:       KSPViewFinalMatResidual_Internal(ksp, B, X, ksp->viewerFinalRes, ksp->formatFinalRes);
473:     }
474:     if (ksp->view) {
475:       KSPView(ksp, ksp->viewer);
476:     }
477:   } else {
478:     VecCreateMPIWithArray(PetscObjectComm((PetscObject)ksp), 1, m1, M1, b, &cb);
479:     VecCreateMPIWithArray(PetscObjectComm((PetscObject)ksp), 1, m1, M1, x, &cx);
480:     KSPSolve(ksp, cb, cx);
481:     VecDestroy(&cb);
482:     VecDestroy(&cx);
483:   }
484:   MatDenseRestoreArray(X, &x);
485:   MatDenseRestoreArrayRead(B, &b);
486:   return(0);
487: }

489: /*MC
490:      KSPHPDDM - Interface with the HPDDM library.

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

494:    Options Database Keys:
495: +   -ksp_richardson_scale <scale, default=1.0> - see KSPRICHARDSON
496: .   -ksp_gmres_restart <restart, default=40> - see KSPGMRES
497: .   -ksp_hpddm_type <type, default=gmres> - any of gmres, bgmres, cg, bcg, gcrodr, bgcrodr, or bfbcg
498: .   -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)
499: .   -ksp_hpddm_enlarge_krylov_subspace <p, default=1> - split the initial right-hand side into multiple vectors (only relevant with nonblock methods)
500: .   -ksp_hpddm_orthogonalization <type, default=cgs> - any of cgs or mgs, see KSPGMRES
501: .   -ksp_hpddm_qr <type, default=cholqr> - distributed QR factorizations with any of cholqr, cgs, or mgs (only relevant with block methods)
502: .   -ksp_hpddm_variant <type, default=left> - any of left, right, or flexible (this option is superseded by KSPSetPCSide())
503: .   -ksp_hpddm_recycle <n, default=0> - number of harmonic Ritz vectors to compute (only relevant with GCRODR or BGCRODR)
504: .   -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)
505: -   -ksp_hpddm_recycle_strategy <type, default=A> - generalized eigenvalue problem A or B to solve for recycling (only relevant with flexible GCRODR or BGCRODR)

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

514:    Level: intermediate

516: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPGMRES, KSPCG, KSPLGMRES, KSPDGMRES
517: M*/
518: PETSC_EXTERN PetscErrorCode KSPCreate_HPDDM(KSP ksp)
519: {
520:   KSP_HPDDM      *data;

524:   PetscNewLog(ksp, &data);
525:   ksp->data = (void*)data;
526:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2);
527:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 1);
528:   ksp->ops->setup          = KSPSetUp_HPDDM;
529:   ksp->ops->solve          = KSPSolve_HPDDM;
530:   ksp->ops->reset          = KSPReset_HPDDM;
531:   ksp->ops->destroy        = KSPDestroy_HPDDM;
532:   ksp->ops->setfromoptions = KSPSetFromOptions_HPDDM;
533:   ksp->ops->view           = KSPView_HPDDM;
534:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMSetDeflationSpace_C", KSPHPDDMSetDeflationSpace_HPDDM);
535:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMGetDeflationSpace_C", KSPHPDDMGetDeflationSpace_HPDDM);
536:   PetscObjectComposeFunction((PetscObject)ksp, "KSPHPDDMMatSolve_C", KSPHPDDMMatSolve_HPDDM);
537:   return(0);
538: }