Actual source code: fcg.c
1: /*
2: This file implements the FCG (Flexible Conjugate Gradient) method
3: */
5: #include <../src/ksp/ksp/impls/fcg/fcgimpl.h>
6: extern PetscErrorCode KSPComputeExtremeSingularValues_CG(KSP, PetscReal *, PetscReal *);
7: extern PetscErrorCode KSPComputeEigenvalues_CG(KSP, PetscInt, PetscReal *, PetscReal *, PetscInt *);
9: #define KSPFCG_DEFAULT_MMAX 30 /* maximum number of search directions to keep */
10: #define KSPFCG_DEFAULT_NPREALLOC 10 /* number of search directions to preallocate */
11: #define KSPFCG_DEFAULT_VECB 5 /* number of search directions to allocate each time new direction vectors are needed */
12: #define KSPFCG_DEFAULT_TRUNCSTRAT KSP_FCD_TRUNC_TYPE_NOTAY
14: static PetscErrorCode KSPAllocateVectors_FCG(KSP ksp, PetscInt nvecsneeded, PetscInt chunksize)
15: {
16: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
17: PetscInt nnewvecs, nvecsprev;
19: PetscFunctionBegin;
20: /* Allocate enough new vectors to add chunksize new vectors, reach nvecsneedtotal, or to reach mmax+1, whichever is smallest */
21: if (fcg->nvecs < PetscMin(fcg->mmax + 1, nvecsneeded)) {
22: nvecsprev = fcg->nvecs;
23: nnewvecs = PetscMin(PetscMax(nvecsneeded - fcg->nvecs, chunksize), fcg->mmax + 1 - fcg->nvecs);
24: PetscCall(KSPCreateVecs(ksp, nnewvecs, &fcg->pCvecs[fcg->nchunks], 0, NULL));
25: PetscCall(KSPCreateVecs(ksp, nnewvecs, &fcg->pPvecs[fcg->nchunks], 0, NULL));
26: fcg->nvecs += nnewvecs;
27: for (PetscInt i = 0; i < nnewvecs; ++i) {
28: fcg->Cvecs[nvecsprev + i] = fcg->pCvecs[fcg->nchunks][i];
29: fcg->Pvecs[nvecsprev + i] = fcg->pPvecs[fcg->nchunks][i];
30: }
31: fcg->chunksizes[fcg->nchunks] = nnewvecs;
32: ++fcg->nchunks;
33: }
34: PetscFunctionReturn(PETSC_SUCCESS);
35: }
37: static PetscErrorCode KSPSetUp_FCG(KSP ksp)
38: {
39: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
40: PetscInt maxit = ksp->max_it;
41: const PetscInt nworkstd = 2;
43: PetscFunctionBegin;
44: /* Allocate "standard" work vectors (not including the basis and transformed basis vectors) */
45: PetscCall(KSPSetWorkVecs(ksp, nworkstd));
47: /* Allocated space for pointers to additional work vectors
48: note that mmax is the number of previous directions, so we add 1 for the current direction,
49: and an extra 1 for the prealloc (which might be empty) */
50: PetscCall(PetscMalloc5(fcg->mmax + 1, &fcg->Pvecs, fcg->mmax + 1, &fcg->Cvecs, fcg->mmax + 1, &fcg->pPvecs, fcg->mmax + 1, &fcg->pCvecs, fcg->mmax + 2, &fcg->chunksizes));
52: /* If the requested number of preallocated vectors is greater than mmax reduce nprealloc */
53: if (fcg->nprealloc > fcg->mmax + 1) PetscCall(PetscInfo(NULL, "Requested nprealloc=%" PetscInt_FMT " is greater than m_max+1=%" PetscInt_FMT ". Resetting nprealloc = m_max+1.\n", fcg->nprealloc, fcg->mmax + 1));
55: /* Preallocate additional work vectors */
56: PetscCall(KSPAllocateVectors_FCG(ksp, fcg->nprealloc, fcg->nprealloc));
57: /*
58: If user requested computations of eigenvalues then allocate work
59: work space needed
60: */
61: if (ksp->calc_sings) {
62: /* get space to store tridiagonal matrix for Lanczos */
63: PetscCall(PetscMalloc4(maxit, &fcg->e, maxit, &fcg->d, maxit, &fcg->ee, maxit, &fcg->dd));
65: ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_CG;
66: ksp->ops->computeeigenvalues = KSPComputeEigenvalues_CG;
67: }
68: PetscFunctionReturn(PETSC_SUCCESS);
69: }
71: static PetscErrorCode KSPSolve_FCG(KSP ksp)
72: {
73: PetscInt i, k, idx, mi;
74: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
75: PetscScalar alpha = 0.0, beta = 0.0, dpi = 0.0, dpiold, s;
76: PetscReal dp = 0.0;
77: Vec B, R, Z, X, Pcurr, Ccurr;
78: Mat Amat, Pmat;
79: PetscInt eigs = ksp->calc_sings; /* Variables for eigen estimation - START*/
80: PetscInt stored_max_it = ksp->max_it;
81: PetscScalar alphaold = 0, betaold = 1.0, *e = NULL, *d = NULL; /* Variables for eigen estimation - FINISH */
83: PetscFunctionBegin;
84: #define VecXDot(x, y, a) (fcg->type == KSP_CG_HERMITIAN ? VecDot(x, y, a) : VecTDot(x, y, a))
85: #define VecXMDot(a, b, c, d) (fcg->type == KSP_CG_HERMITIAN ? VecMDot(a, b, c, d) : VecMTDot(a, b, c, d))
87: X = ksp->vec_sol;
88: B = ksp->vec_rhs;
89: R = ksp->work[0];
90: Z = ksp->work[1];
92: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
93: if (eigs) {
94: e = fcg->e;
95: d = fcg->d;
96: e[0] = 0.0;
97: }
98: /* Compute initial residual needed for convergence check*/
99: ksp->its = 0;
100: if (!ksp->guess_zero) {
101: PetscCall(KSP_MatMult(ksp, Amat, X, R));
102: PetscCall(VecAYPX(R, -1.0, B)); /* r <- b - Ax */
103: } else {
104: PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */
105: }
106: switch (ksp->normtype) {
107: case KSP_NORM_PRECONDITIONED:
108: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
109: PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- dqrt(z'*z) = sqrt(e'*A'*B'*B*A*e) */
110: KSPCheckNorm(ksp, dp);
111: break;
112: case KSP_NORM_UNPRECONDITIONED:
113: PetscCall(VecNorm(R, NORM_2, &dp)); /* dp <- sqrt(r'*r) = sqrt(e'*A'*A*e) */
114: KSPCheckNorm(ksp, dp);
115: break;
116: case KSP_NORM_NATURAL:
117: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
118: PetscCall(VecXDot(R, Z, &s));
119: KSPCheckDot(ksp, s);
120: dp = PetscSqrtReal(PetscAbsScalar(s)); /* dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e) */
121: break;
122: case KSP_NORM_NONE:
123: dp = 0.0;
124: break;
125: default:
126: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
127: }
129: /* Initial Convergence Check */
130: PetscCall(KSPLogResidualHistory(ksp, dp));
131: PetscCall(KSPMonitor(ksp, 0, dp));
132: ksp->rnorm = dp;
133: if (ksp->normtype == KSP_NORM_NONE) {
134: PetscCall(KSPConvergedSkip(ksp, 0, dp, &ksp->reason, ksp->cnvP));
135: } else {
136: PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP));
137: }
138: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
140: /* Apply PC if not already done for convergence check */
141: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->normtype == KSP_NORM_NONE) PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
143: i = 0;
144: do {
145: ksp->its = i + 1;
147: /* If needbe, allocate a new chunk of vectors in P and C */
148: PetscCall(KSPAllocateVectors_FCG(ksp, i + 1, fcg->vecb));
150: /* Note that we wrap around and start clobbering old vectors */
151: idx = i % (fcg->mmax + 1);
152: Pcurr = fcg->Pvecs[idx];
153: Ccurr = fcg->Cvecs[idx];
155: /* number of old directions to orthogonalize against */
156: switch (fcg->truncstrat) {
157: case KSP_FCD_TRUNC_TYPE_STANDARD:
158: mi = fcg->mmax;
159: break;
160: case KSP_FCD_TRUNC_TYPE_NOTAY:
161: mi = ((i - 1) % fcg->mmax) + 1;
162: break;
163: default:
164: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized Truncation Strategy");
165: }
167: /* Compute a new column of P (Currently does not support modified G-S or iterative refinement)*/
168: PetscCall(VecCopy(Z, Pcurr));
170: {
171: PetscInt l, ndots;
173: l = PetscMax(0, i - mi);
174: ndots = i - l;
175: if (ndots) {
176: PetscInt j;
177: Vec *Pold, *Cold;
178: PetscScalar *dots;
180: PetscCall(PetscMalloc3(ndots, &dots, ndots, &Cold, ndots, &Pold));
181: for (k = l, j = 0; j < ndots; ++k, ++j) {
182: idx = k % (fcg->mmax + 1);
183: Cold[j] = fcg->Cvecs[idx];
184: Pold[j] = fcg->Pvecs[idx];
185: }
186: PetscCall(VecXMDot(Z, ndots, Cold, dots));
187: for (k = 0; k < ndots; ++k) dots[k] = -dots[k];
188: PetscCall(VecMAXPY(Pcurr, ndots, dots, Pold));
189: PetscCall(PetscFree3(dots, Cold, Pold));
190: }
191: }
193: /* Update X and R */
194: betaold = beta;
195: PetscCall(VecXDot(Pcurr, R, &beta)); /* beta <- pi'*r */
196: KSPCheckDot(ksp, beta);
197: if ((i > 0) && (PetscAbsScalar(beta * betaold) < 0.0)) {
198: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Diverged due to indefinite preconditioner, beta %g, betaold %g", (double)PetscRealPart(beta), (double)PetscRealPart(betaold));
199: ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
200: PetscCall(PetscInfo(ksp, "diverging due to indefinite preconditioner\n"));
201: break;
202: }
203: PetscCall(KSP_MatMult(ksp, Amat, Pcurr, Ccurr)); /* w <- A*pi (stored in ci) */
204: dpiold = dpi;
205: PetscCall(VecXDot(Pcurr, Ccurr, &dpi)); /* dpi <- pi'*w */
206: if ((dpi == 0.0) || ((i > 0) && ((PetscSign(PetscRealPart(dpi)) * PetscSign(PetscRealPart(dpiold))) < 0.0))) {
207: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "Diverged due to indefinite matrix, dpi %g, dpiold %g", (double)PetscRealPart(dpi), (double)PetscRealPart(dpiold));
208: ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
209: PetscCall(PetscInfo(ksp, "diverging due to indefinite matrix\n"));
210: break;
211: }
212: alphaold = alpha;
213: alpha = beta / dpi; /* alpha <- beta/dpi */
214: PetscCall(VecAXPY(X, alpha, Pcurr)); /* x <- x + alpha * pi */
215: PetscCall(VecAXPY(R, -alpha, Ccurr)); /* r <- r - alpha * wi */
217: /* Compute norm for convergence check */
218: switch (ksp->normtype) {
219: case KSP_NORM_PRECONDITIONED:
220: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
221: PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- sqrt(z'*z) = sqrt(e'*A'*B'*B*A*e) */
222: KSPCheckNorm(ksp, dp);
223: break;
224: case KSP_NORM_UNPRECONDITIONED:
225: PetscCall(VecNorm(R, NORM_2, &dp)); /* dp <- sqrt(r'*r) = sqrt(e'*A'*A*e) */
226: KSPCheckNorm(ksp, dp);
227: break;
228: case KSP_NORM_NATURAL:
229: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
230: PetscCall(VecXDot(R, Z, &s));
231: KSPCheckDot(ksp, s);
232: dp = PetscSqrtReal(PetscAbsScalar(s)); /* dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e) */
233: break;
234: case KSP_NORM_NONE:
235: dp = 0.0;
236: break;
237: default:
238: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
239: }
241: if (eigs) {
242: if (i > 0) {
243: PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Can not change maxit AND calculate eigenvalues");
244: e[i] = PetscSqrtReal(PetscAbsScalar(beta / betaold)) / alphaold;
245: d[i] = PetscSqrtReal(PetscAbsScalar(beta / betaold)) * e[i] + 1.0 / alpha;
246: } else {
247: d[i] = PetscSqrtReal(PetscAbsScalar(beta)) * e[i] + 1.0 / alpha;
248: }
249: }
251: /* Check for convergence */
252: ksp->rnorm = dp;
253: PetscCall(KSPLogResidualHistory(ksp, dp));
254: PetscCall(KSPMonitor(ksp, i + 1, dp));
255: PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
256: if (ksp->reason) break;
258: /* Apply PC if not already done for convergence check */
259: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->normtype == KSP_NORM_NONE) PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
261: /* Compute current C (which is W/dpi) */
262: PetscCall(VecScale(Ccurr, 1.0 / dpi)); /* w <- ci/dpi */
263: ++i;
264: } while (i < ksp->max_it);
265: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: static PetscErrorCode KSPReset_FCG(KSP ksp)
270: {
271: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
273: PetscFunctionBegin;
274: /* Destroy P and C vectors and the arrays that manage pointers to them */
275: if (fcg->nvecs) {
276: for (PetscInt i = 0; i < fcg->nchunks; ++i) {
277: PetscCall(VecDestroyVecs(fcg->chunksizes[i], &fcg->pPvecs[i]));
278: PetscCall(VecDestroyVecs(fcg->chunksizes[i], &fcg->pCvecs[i]));
279: }
280: fcg->nchunks = fcg->nvecs = 0;
281: }
282: PetscCall(PetscFree5(fcg->Pvecs, fcg->Cvecs, fcg->pPvecs, fcg->pCvecs, fcg->chunksizes));
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: static PetscErrorCode KSPDestroy_FCG(KSP ksp)
287: {
288: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
290: PetscFunctionBegin;
291: /* free space used for singular value calculations */
292: if (ksp->calc_sings) PetscCall(PetscFree4(fcg->e, fcg->d, fcg->ee, fcg->dd));
293: PetscCall(KSPDestroyDefault(ksp));
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: static PetscErrorCode KSPView_FCG(KSP ksp, PetscViewer viewer)
298: {
299: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
300: PetscBool isascii, isstring;
301: const char *truncstr;
303: PetscFunctionBegin;
304: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
305: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
307: if (fcg->truncstrat == KSP_FCD_TRUNC_TYPE_STANDARD) truncstr = "Using standard truncation strategy";
308: else if (fcg->truncstrat == KSP_FCD_TRUNC_TYPE_NOTAY) truncstr = "Using Notay's truncation strategy";
309: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Undefined FCG truncation strategy");
311: if (isascii) {
312: PetscCall(PetscViewerASCIIPrintf(viewer, " m_max=%" PetscInt_FMT "\n", fcg->mmax));
313: PetscCall(PetscViewerASCIIPrintf(viewer, " preallocated %" PetscInt_FMT " directions\n", PetscMin(fcg->nprealloc, fcg->mmax + 1)));
314: PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", truncstr));
315: } else if (isstring) {
316: PetscCall(PetscViewerStringSPrintf(viewer, "m_max %" PetscInt_FMT " nprealloc %" PetscInt_FMT " %s", fcg->mmax, fcg->nprealloc, truncstr));
317: }
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: /*@
322: KSPFCGSetMmax - set the maximum number of previous directions `KSPFCG` will store for orthogonalization
324: Logically Collective
326: Input Parameters:
327: + ksp - the Krylov space context
328: - mmax - the maximum number of previous directions to orthogonalize against
330: Options Database Key:
331: . -ksp_fcg_mmax N - maximum number of search directions
333: Level: intermediate
335: Note:
336: `mmax` + 1 directions are stored (`mmax` previous ones along with a current one)
337: and whether all are used in each iteration also depends on the truncation strategy, see `KSPFCGSetTruncationType()`
339: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCDTruncationType`, `KSPFCGSetTruncationType()`, `KSPFCGGetNprealloc()`, `KSPFCGGetMmax()`
340: @*/
341: PetscErrorCode KSPFCGSetMmax(KSP ksp, PetscInt mmax)
342: {
343: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
345: PetscFunctionBegin;
348: fcg->mmax = mmax;
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: /*@
353: KSPFCGGetMmax - get the maximum number of previous directions `KSPFCG` will store
355: Not Collective
357: Input Parameter:
358: . ksp - the Krylov space context
360: Output Parameter:
361: . mmax - the maximum number of previous directions allowed for orthogonalization
363: Level: intermediate
365: Note:
366: `KSPFCG` stores `mmax`+1 directions at most (`mmax` previous ones, and one current one)
368: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGGetNprealloc()`, `KSPFCGSetMmax()`
369: @*/
370: PetscErrorCode KSPFCGGetMmax(KSP ksp, PetscInt *mmax)
371: {
372: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
374: PetscFunctionBegin;
376: *mmax = fcg->mmax;
377: PetscFunctionReturn(PETSC_SUCCESS);
378: }
380: /*@
381: KSPFCGSetNprealloc - set the number of directions to preallocate with `KSPFCG`
383: Logically Collective
385: Input Parameters:
386: + ksp - the Krylov space context
387: - nprealloc - the number of vectors to preallocate
389: Options Database Key:
390: . -ksp_fcg_nprealloc N - number of directions to preallocate
392: Level: advanced
394: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGGetNprealloc()`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`
395: @*/
396: PetscErrorCode KSPFCGSetNprealloc(KSP ksp, PetscInt nprealloc)
397: {
398: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
400: PetscFunctionBegin;
403: PetscCheck(nprealloc <= fcg->mmax + 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot preallocate more than m_max+1 vectors");
404: fcg->nprealloc = nprealloc;
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: /*@
409: KSPFCGGetNprealloc - get the number of directions preallocate by `KSPFCG`
411: Not Collective
413: Input Parameter:
414: . ksp - the Krylov space context
416: Output Parameter:
417: . nprealloc - the number of directions preallocated
419: Level: advanced
421: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGSetNprealloc()`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`
422: @*/
423: PetscErrorCode KSPFCGGetNprealloc(KSP ksp, PetscInt *nprealloc)
424: {
425: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
427: PetscFunctionBegin;
429: *nprealloc = fcg->nprealloc;
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }
433: /*@
434: KSPFCGSetTruncationType - specify how many of its stored previous directions `KSPFCG` uses during orthogonalization
436: Logically Collective
438: Input Parameters:
439: + ksp - the Krylov space context
440: - truncstrat - the choice of strategy
441: .vb
442: KSP_FCD_TRUNC_TYPE_STANDARD uses all (up to `mmax`) stored directions
443: KSP_FCD_TRUNC_TYPE_NOTAY uses the last `max(1,mod(i,mmax))` stored directions at iteration i = 0, 1, ...
444: .ve
446: Options Database Key:
447: . -ksp_fcg_truncation_type (standard|notay) - specify how many of its stored previous directions `KSPFCG` uses during orthogonalization
449: Level: intermediate
451: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCDTruncationType`, `KSPFCGGetTruncationType()`, `KSPFCGSetNprealloc()`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`,
452: `KSP_FCD_TRUNC_TYPE_STANDARD`, `KSP_FCD_TRUNC_TYPE_NOTAY`
453: @*/
454: PetscErrorCode KSPFCGSetTruncationType(KSP ksp, KSPFCDTruncationType truncstrat)
455: {
456: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
458: PetscFunctionBegin;
461: fcg->truncstrat = truncstrat;
462: PetscFunctionReturn(PETSC_SUCCESS);
463: }
465: /*@
466: KSPFCGGetTruncationType - get the truncation strategy employed by `KSPFCG`
468: Not Collective
470: Input Parameter:
471: . ksp - the Krylov space context
473: Output Parameter:
474: . truncstrat - the strategy type
476: Level: intermediate
478: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGSetTruncationType()`, `KSPFCDTruncationType`, `KSP_FCD_TRUNC_TYPE_STANDARD`, `KSP_FCD_TRUNC_TYPE_NOTAY`
479: @*/
480: PetscErrorCode KSPFCGGetTruncationType(KSP ksp, KSPFCDTruncationType *truncstrat)
481: {
482: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
484: PetscFunctionBegin;
486: *truncstrat = fcg->truncstrat;
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: static PetscErrorCode KSPSetFromOptions_FCG(KSP ksp, PetscOptionItems PetscOptionsObject)
491: {
492: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
493: PetscInt mmax, nprealloc;
494: PetscBool flg;
496: PetscFunctionBegin;
497: PetscOptionsHeadBegin(PetscOptionsObject, "KSP FCG Options");
498: PetscCall(PetscOptionsInt("-ksp_fcg_mmax", "Maximum number of search directions to store", "KSPFCGSetMmax", fcg->mmax, &mmax, &flg));
499: if (flg) PetscCall(KSPFCGSetMmax(ksp, mmax));
500: PetscCall(PetscOptionsInt("-ksp_fcg_nprealloc", "Number of directions to preallocate", "KSPFCGSetNprealloc", fcg->nprealloc, &nprealloc, &flg));
501: if (flg) PetscCall(KSPFCGSetNprealloc(ksp, nprealloc));
502: PetscCall(PetscOptionsEnum("-ksp_fcg_truncation_type", "Truncation approach for directions", "KSPFCGSetTruncationType", KSPFCDTruncationTypes, (PetscEnum)fcg->truncstrat, (PetscEnum *)&fcg->truncstrat, NULL));
503: PetscOptionsHeadEnd();
504: PetscFunctionReturn(PETSC_SUCCESS);
505: }
507: /*MC
508: KSPFCG - Implements the Flexible Conjugate Gradient method (FCG) {cite}`flexiblecg`, {cite}`generalizedcg`.
509: Unlike most `KSP` methods this allows the preconditioner to be nonlinear. [](sec_flexibleksp)
511: Options Database Keys:
512: + -ksp_fcg_mmax N - maximum number of search directions, similar to the restart in `KSPGMRES` and `KSPFGMRES`
513: . -ksp_fcg_nprealloc N - number of directions to preallocate
514: - -ksp_fcg_truncation_type (standard|notay) - truncation approach for directions
516: Level: beginner
518: Notes:
519: `KSPFCG` requires the matrix to be symmetric positive-definite (SPD); for non-SPD problems use `KSPFGMRES` or `KSPGCR`.
521: Supports left preconditioning only.
523: Contributed by:
524: Patrick Sanan
526: .seealso: [](ch_ksp), [](sec_flexibleksp), `KSPGCR`, `KSPPIPEGCR`, `KSPPIPEFCG`, `KSPFGMRES`, `KSPCG`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`, `KSPFCGSetNprealloc()`, `KSPFCGGetNprealloc()`,
527: `KSPFCGSetTruncationType()`, `KSPFCGGetTruncationType()`, `KSPFCDTruncationType`
528: M*/
529: PETSC_EXTERN PetscErrorCode KSPCreate_FCG(KSP ksp)
530: {
531: KSP_FCG *fcg;
533: PetscFunctionBegin;
534: PetscCall(PetscNew(&fcg));
535: fcg->type = !PetscDefined(USE_COMPLEX) ? KSP_CG_SYMMETRIC : KSP_CG_HERMITIAN;
536: fcg->mmax = KSPFCG_DEFAULT_MMAX;
537: fcg->nprealloc = KSPFCG_DEFAULT_NPREALLOC;
538: fcg->nvecs = 0;
539: fcg->vecb = KSPFCG_DEFAULT_VECB;
540: fcg->nchunks = 0;
541: fcg->truncstrat = KSPFCG_DEFAULT_TRUNCSTRAT;
543: ksp->data = (void *)fcg;
545: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
546: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 1));
547: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 1));
548: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
550: ksp->ops->setup = KSPSetUp_FCG;
551: ksp->ops->solve = KSPSolve_FCG;
552: ksp->ops->reset = KSPReset_FCG;
553: ksp->ops->destroy = KSPDestroy_FCG;
554: ksp->ops->view = KSPView_FCG;
555: ksp->ops->setfromoptions = KSPSetFromOptions_FCG;
556: ksp->ops->buildsolution = KSPBuildSolutionDefault;
557: ksp->ops->buildresidual = KSPBuildResidualDefault;
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }