Actual source code: hpddm.cxx
petsc-3.13.6 2020-09-29
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, ©);
336: PetscArraycpy(copy, array, m1 * N2);
337: MatDenseRestoreArray(*U, ©);
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: }