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: PetscInt i;
17: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
18: PetscInt nnewvecs, nvecsprev;
20: PetscFunctionBegin;
21: /* Allocate enough new vectors to add chunksize new vectors, reach nvecsneedtotal, or to reach mmax+1, whichever is smallest */
22: if (fcg->nvecs < PetscMin(fcg->mmax + 1, nvecsneeded)) {
23: nvecsprev = fcg->nvecs;
24: nnewvecs = PetscMin(PetscMax(nvecsneeded - fcg->nvecs, chunksize), fcg->mmax + 1 - fcg->nvecs);
25: PetscCall(KSPCreateVecs(ksp, nnewvecs, &fcg->pCvecs[fcg->nchunks], 0, NULL));
26: PetscCall(KSPCreateVecs(ksp, nnewvecs, &fcg->pPvecs[fcg->nchunks], 0, NULL));
27: fcg->nvecs += nnewvecs;
28: for (i = 0; i < nnewvecs; ++i) {
29: fcg->Cvecs[nvecsprev + i] = fcg->pCvecs[fcg->nchunks][i];
30: fcg->Pvecs[nvecsprev + i] = fcg->pPvecs[fcg->nchunks][i];
31: }
32: fcg->chunksizes[fcg->nchunks] = nnewvecs;
33: ++fcg->nchunks;
34: }
35: PetscFunctionReturn(PETSC_SUCCESS);
36: }
38: static PetscErrorCode KSPSetUp_FCG(KSP ksp)
39: {
40: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
41: PetscInt maxit = ksp->max_it;
42: const PetscInt nworkstd = 2;
44: PetscFunctionBegin;
46: /* Allocate "standard" work vectors (not including the basis and transformed basis vectors) */
47: PetscCall(KSPSetWorkVecs(ksp, nworkstd));
49: /* Allocated space for pointers to additional work vectors
50: note that mmax is the number of previous directions, so we add 1 for the current direction,
51: and an extra 1 for the prealloc (which might be empty) */
52: 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));
54: /* If the requested number of preallocated vectors is greater than mmax reduce nprealloc */
55: 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));
57: /* Preallocate additional work vectors */
58: PetscCall(KSPAllocateVectors_FCG(ksp, fcg->nprealloc, fcg->nprealloc));
59: /*
60: If user requested computations of eigenvalues then allocate work
61: work space needed
62: */
63: if (ksp->calc_sings) {
64: /* get space to store tridiagonal matrix for Lanczos */
65: PetscCall(PetscMalloc4(maxit, &fcg->e, maxit, &fcg->d, maxit, &fcg->ee, maxit, &fcg->dd));
67: ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_CG;
68: ksp->ops->computeeigenvalues = KSPComputeEigenvalues_CG;
69: }
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: static PetscErrorCode KSPSolve_FCG(KSP ksp)
74: {
75: PetscInt i, k, idx, mi;
76: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
77: PetscScalar alpha = 0.0, beta = 0.0, dpi, s;
78: PetscReal dp = 0.0;
79: Vec B, R, Z, X, Pcurr, Ccurr;
80: Mat Amat, Pmat;
81: PetscInt eigs = ksp->calc_sings; /* Variables for eigen estimation - START*/
82: PetscInt stored_max_it = ksp->max_it;
83: PetscScalar alphaold = 0, betaold = 1.0, *e = NULL, *d = NULL; /* Variables for eigen estimation - FINISH */
85: PetscFunctionBegin;
87: #define VecXDot(x, y, a) (((fcg->type) == (KSP_CG_HERMITIAN)) ? VecDot(x, y, a) : VecTDot(x, y, a))
88: #define VecXMDot(a, b, c, d) (((fcg->type) == (KSP_CG_HERMITIAN)) ? VecMDot(a, b, c, d) : VecMTDot(a, b, c, d))
90: X = ksp->vec_sol;
91: B = ksp->vec_rhs;
92: R = ksp->work[0];
93: Z = ksp->work[1];
95: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
96: if (eigs) {
97: e = fcg->e;
98: d = fcg->d;
99: e[0] = 0.0;
100: }
101: /* Compute initial residual needed for convergence check*/
102: ksp->its = 0;
103: if (!ksp->guess_zero) {
104: PetscCall(KSP_MatMult(ksp, Amat, X, R));
105: PetscCall(VecAYPX(R, -1.0, B)); /* r <- b - Ax */
106: } else {
107: PetscCall(VecCopy(B, R)); /* r <- b (x is 0) */
108: }
109: switch (ksp->normtype) {
110: case KSP_NORM_PRECONDITIONED:
111: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
112: PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- dqrt(z'*z) = sqrt(e'*A'*B'*B*A*e) */
113: KSPCheckNorm(ksp, dp);
114: break;
115: case KSP_NORM_UNPRECONDITIONED:
116: PetscCall(VecNorm(R, NORM_2, &dp)); /* dp <- sqrt(r'*r) = sqrt(e'*A'*A*e) */
117: KSPCheckNorm(ksp, dp);
118: break;
119: case KSP_NORM_NATURAL:
120: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
121: PetscCall(VecXDot(R, Z, &s));
122: KSPCheckDot(ksp, s);
123: dp = PetscSqrtReal(PetscAbsScalar(s)); /* dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e) */
124: break;
125: case KSP_NORM_NONE:
126: dp = 0.0;
127: break;
128: default:
129: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
130: }
132: /* Initial Convergence Check */
133: PetscCall(KSPLogResidualHistory(ksp, dp));
134: PetscCall(KSPMonitor(ksp, 0, dp));
135: ksp->rnorm = dp;
136: if (ksp->normtype == KSP_NORM_NONE) {
137: PetscCall(KSPConvergedSkip(ksp, 0, dp, &ksp->reason, ksp->cnvP));
138: } else {
139: PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP));
140: }
141: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
143: /* Apply PC if not already done for convergence check */
144: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->normtype == KSP_NORM_NONE) { PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */ }
146: i = 0;
147: do {
148: ksp->its = i + 1;
150: /* If needbe, allocate a new chunk of vectors in P and C */
151: PetscCall(KSPAllocateVectors_FCG(ksp, i + 1, fcg->vecb));
153: /* Note that we wrap around and start clobbering old vectors */
154: idx = i % (fcg->mmax + 1);
155: Pcurr = fcg->Pvecs[idx];
156: Ccurr = fcg->Cvecs[idx];
158: /* number of old directions to orthogonalize against */
159: switch (fcg->truncstrat) {
160: case KSP_FCD_TRUNC_TYPE_STANDARD:
161: mi = fcg->mmax;
162: break;
163: case KSP_FCD_TRUNC_TYPE_NOTAY:
164: mi = ((i - 1) % fcg->mmax) + 1;
165: break;
166: default:
167: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized Truncation Strategy");
168: }
170: /* Compute a new column of P (Currently does not support modified G-S or iterative refinement)*/
171: PetscCall(VecCopy(Z, Pcurr));
173: {
174: PetscInt l, ndots;
176: l = PetscMax(0, i - mi);
177: ndots = i - l;
178: if (ndots) {
179: PetscInt j;
180: Vec *Pold, *Cold;
181: PetscScalar *dots;
183: PetscCall(PetscMalloc3(ndots, &dots, ndots, &Cold, ndots, &Pold));
184: for (k = l, j = 0; j < ndots; ++k, ++j) {
185: idx = k % (fcg->mmax + 1);
186: Cold[j] = fcg->Cvecs[idx];
187: Pold[j] = fcg->Pvecs[idx];
188: }
189: PetscCall(VecXMDot(Z, ndots, Cold, dots));
190: for (k = 0; k < ndots; ++k) dots[k] = -dots[k];
191: PetscCall(VecMAXPY(Pcurr, ndots, dots, Pold));
192: PetscCall(PetscFree3(dots, Cold, Pold));
193: }
194: }
196: /* Update X and R */
197: betaold = beta;
198: PetscCall(VecXDot(Pcurr, R, &beta)); /* beta <- pi'*r */
199: KSPCheckDot(ksp, beta);
200: PetscCall(KSP_MatMult(ksp, Amat, Pcurr, Ccurr)); /* w <- A*pi (stored in ci) */
201: PetscCall(VecXDot(Pcurr, Ccurr, &dpi)); /* dpi <- pi'*w */
202: alphaold = alpha;
203: alpha = beta / dpi; /* alpha <- beta/dpi */
204: PetscCall(VecAXPY(X, alpha, Pcurr)); /* x <- x + alpha * pi */
205: PetscCall(VecAXPY(R, -alpha, Ccurr)); /* r <- r - alpha * wi */
207: /* Compute norm for convergence check */
208: switch (ksp->normtype) {
209: case KSP_NORM_PRECONDITIONED:
210: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
211: PetscCall(VecNorm(Z, NORM_2, &dp)); /* dp <- sqrt(z'*z) = sqrt(e'*A'*B'*B*A*e) */
212: KSPCheckNorm(ksp, dp);
213: break;
214: case KSP_NORM_UNPRECONDITIONED:
215: PetscCall(VecNorm(R, NORM_2, &dp)); /* dp <- sqrt(r'*r) = sqrt(e'*A'*A*e) */
216: KSPCheckNorm(ksp, dp);
217: break;
218: case KSP_NORM_NATURAL:
219: PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */
220: PetscCall(VecXDot(R, Z, &s));
221: KSPCheckDot(ksp, s);
222: dp = PetscSqrtReal(PetscAbsScalar(s)); /* dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e) */
223: break;
224: case KSP_NORM_NONE:
225: dp = 0.0;
226: break;
227: default:
228: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
229: }
231: /* Check for convergence */
232: ksp->rnorm = dp;
233: PetscCall(KSPLogResidualHistory(ksp, dp));
234: PetscCall(KSPMonitor(ksp, i + 1, dp));
235: PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
236: if (ksp->reason) break;
238: /* Apply PC if not already done for convergence check */
239: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->normtype == KSP_NORM_NONE) { PetscCall(KSP_PCApply(ksp, R, Z)); /* z <- Br */ }
241: /* Compute current C (which is W/dpi) */
242: PetscCall(VecScale(Ccurr, 1.0 / dpi)); /* w <- ci/dpi */
244: if (eigs) {
245: if (i > 0) {
246: PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Can not change maxit AND calculate eigenvalues");
247: e[i] = PetscSqrtReal(PetscAbsScalar(beta / betaold)) / alphaold;
248: d[i] = PetscSqrtReal(PetscAbsScalar(beta / betaold)) * e[i] + 1.0 / alpha;
249: } else {
250: d[i] = PetscSqrtReal(PetscAbsScalar(beta)) * e[i] + 1.0 / alpha;
251: }
252: fcg->ned = ksp->its - 1;
253: }
254: ++i;
255: } while (i < ksp->max_it);
256: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: static PetscErrorCode KSPDestroy_FCG(KSP ksp)
261: {
262: PetscInt i;
263: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
265: PetscFunctionBegin;
267: /* Destroy "standard" work vecs */
268: PetscCall(VecDestroyVecs(ksp->nwork, &ksp->work));
270: /* Destroy P and C vectors and the arrays that manage pointers to them */
271: if (fcg->nvecs) {
272: for (i = 0; i < fcg->nchunks; ++i) {
273: PetscCall(VecDestroyVecs(fcg->chunksizes[i], &fcg->pPvecs[i]));
274: PetscCall(VecDestroyVecs(fcg->chunksizes[i], &fcg->pCvecs[i]));
275: }
276: }
277: PetscCall(PetscFree5(fcg->Pvecs, fcg->Cvecs, fcg->pPvecs, fcg->pCvecs, fcg->chunksizes));
278: /* free space used for singular value calculations */
279: if (ksp->calc_sings) PetscCall(PetscFree4(fcg->e, fcg->d, fcg->ee, fcg->dd));
280: PetscCall(KSPDestroyDefault(ksp));
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: static PetscErrorCode KSPView_FCG(KSP ksp, PetscViewer viewer)
285: {
286: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
287: PetscBool iascii, isstring;
288: const char *truncstr;
290: PetscFunctionBegin;
291: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
292: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
294: if (fcg->truncstrat == KSP_FCD_TRUNC_TYPE_STANDARD) truncstr = "Using standard truncation strategy";
295: else if (fcg->truncstrat == KSP_FCD_TRUNC_TYPE_NOTAY) truncstr = "Using Notay's truncation strategy";
296: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Undefined FCG truncation strategy");
298: if (iascii) {
299: PetscCall(PetscViewerASCIIPrintf(viewer, " m_max=%" PetscInt_FMT "\n", fcg->mmax));
300: PetscCall(PetscViewerASCIIPrintf(viewer, " preallocated %" PetscInt_FMT " directions\n", PetscMin(fcg->nprealloc, fcg->mmax + 1)));
301: PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", truncstr));
302: } else if (isstring) {
303: PetscCall(PetscViewerStringSPrintf(viewer, "m_max %" PetscInt_FMT " nprealloc %" PetscInt_FMT " %s", fcg->mmax, fcg->nprealloc, truncstr));
304: }
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }
308: /*@
309: KSPFCGSetMmax - set the maximum number of previous directions `KSPFCG` will store for orthogonalization
311: Logically Collective
313: Input Parameters:
314: + ksp - the Krylov space context
315: - mmax - the maximum number of previous directions to orthogonalize against
317: Options Database Key:
318: . -ksp_fcg_mmax <N> - maximum number of search directions
320: Level: intermediate
322: Note:
323: mmax + 1 directions are stored (mmax previous ones along with a current one)
324: and whether all are used in each iteration also depends on the truncation strategy
325: (see KSPFCGSetTruncationType())
327: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGGetNprealloc()`, `KSPFCGetMmax()`
328: @*/
329: PetscErrorCode KSPFCGSetMmax(KSP ksp, PetscInt mmax)
330: {
331: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
333: PetscFunctionBegin;
336: fcg->mmax = mmax;
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: /*@
341: KSPFCGGetMmax - get the maximum number of previous directions `KSPFCG` will store
343: Not Collective
345: Input Parameter:
346: . ksp - the Krylov space context
348: Output Parameter:
349: . mmax - the maximum number of previous directions allowed for orthogonalization
351: Level: intermediate
353: Note:
354: FCG stores mmax+1 directions at most (mmax previous ones, and one current one)
356: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGGetNprealloc()`, `KSPFCGSetMmax()`
357: @*/
358: PetscErrorCode KSPFCGGetMmax(KSP ksp, PetscInt *mmax)
359: {
360: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
362: PetscFunctionBegin;
364: *mmax = fcg->mmax;
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: /*@
369: KSPFCGSetNprealloc - set the number of directions to preallocate with `KSPFCG`
371: Logically Collective
373: Input Parameters:
374: + ksp - the Krylov space context
375: - nprealloc - the number of vectors to preallocate
377: Options Database Key:
378: . -ksp_fcg_nprealloc <N> - number of directions to preallocate
380: Level: advanced
382: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGGetNprealloc()`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`
383: @*/
384: PetscErrorCode KSPFCGSetNprealloc(KSP ksp, PetscInt nprealloc)
385: {
386: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
388: PetscFunctionBegin;
391: PetscCheck(nprealloc <= fcg->mmax + 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot preallocate more than m_max+1 vectors");
392: fcg->nprealloc = nprealloc;
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: /*@
397: KSPFCGGetNprealloc - get the number of directions preallocate by `KSPFCG`
399: Not Collective
401: Input Parameter:
402: . ksp - the Krylov space context
404: Output Parameter:
405: . nprealloc - the number of directions preallocated
407: Level: advanced
409: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGSetNprealloc()`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`
410: @*/
411: PetscErrorCode KSPFCGGetNprealloc(KSP ksp, PetscInt *nprealloc)
412: {
413: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
415: PetscFunctionBegin;
417: *nprealloc = fcg->nprealloc;
418: PetscFunctionReturn(PETSC_SUCCESS);
419: }
421: /*@
422: KSPFCGSetTruncationType - specify how many of its stored previous directions `KSPFCG` uses during orthoganalization
424: Logically Collective
426: Input Parameters:
427: + ksp - the Krylov space context
428: - truncstrat - the choice of strategy
429: .vb
430: KSP_FCD_TRUNC_TYPE_STANDARD uses all (up to mmax) stored directions
431: KSP_FCD_TRUNC_TYPE_NOTAY uses the last max(1,mod(i,mmax)) stored directions at iteration i=0,1,..
432: .ve
434: Options Database Key:
435: . -ksp_fcg_truncation_type <standard, notay> - specify how many of its stored previous directions `KSPFCG` uses during orthoganalization
437: Level: intermediate
439: .seealso: [](ch_ksp), `KSPFCDTruncationType`, `KSPFCGGetTruncationType`, `KSPFCGSetNprealloc()`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`
440: @*/
441: PetscErrorCode KSPFCGSetTruncationType(KSP ksp, KSPFCDTruncationType truncstrat)
442: {
443: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
445: PetscFunctionBegin;
448: fcg->truncstrat = truncstrat;
449: PetscFunctionReturn(PETSC_SUCCESS);
450: }
452: /*@
453: KSPFCGGetTruncationType - get the truncation strategy employed by `KSPFCG`
455: Not Collective
457: Input Parameter:
458: . ksp - the Krylov space context
460: Output Parameter:
461: . truncstrat - the strategy type
463: Level: intermediate
465: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGSetTruncationType`, `KSPFCDTruncationType`, `KSPFCGSetTruncationType()`
466: @*/
467: PetscErrorCode KSPFCGGetTruncationType(KSP ksp, KSPFCDTruncationType *truncstrat)
468: {
469: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
471: PetscFunctionBegin;
473: *truncstrat = fcg->truncstrat;
474: PetscFunctionReturn(PETSC_SUCCESS);
475: }
477: static PetscErrorCode KSPSetFromOptions_FCG(KSP ksp, PetscOptionItems *PetscOptionsObject)
478: {
479: KSP_FCG *fcg = (KSP_FCG *)ksp->data;
480: PetscInt mmax, nprealloc;
481: PetscBool flg;
483: PetscFunctionBegin;
484: PetscOptionsHeadBegin(PetscOptionsObject, "KSP FCG Options");
485: PetscCall(PetscOptionsInt("-ksp_fcg_mmax", "Maximum number of search directions to store", "KSPFCGSetMmax", fcg->mmax, &mmax, &flg));
486: if (flg) PetscCall(KSPFCGSetMmax(ksp, mmax));
487: PetscCall(PetscOptionsInt("-ksp_fcg_nprealloc", "Number of directions to preallocate", "KSPFCGSetNprealloc", fcg->nprealloc, &nprealloc, &flg));
488: if (flg) PetscCall(KSPFCGSetNprealloc(ksp, nprealloc));
489: PetscCall(PetscOptionsEnum("-ksp_fcg_truncation_type", "Truncation approach for directions", "KSPFCGSetTruncationType", KSPFCDTruncationTypes, (PetscEnum)fcg->truncstrat, (PetscEnum *)&fcg->truncstrat, NULL));
490: PetscOptionsHeadEnd();
491: PetscFunctionReturn(PETSC_SUCCESS);
492: }
494: /*MC
495: KSPFCG - Implements the Flexible Conjugate Gradient method (FCG). Unlike most `KSP` methods this allows the preconditioner to be nonlinear. [](sec_flexibleksp)
497: Options Database Keys:
498: + -ksp_fcg_mmax <N> - maximum number of search directions
499: . -ksp_fcg_nprealloc <N> - number of directions to preallocate
500: - -ksp_fcg_truncation_type <standard,notay> - truncation approach for directions
502: Level: beginner
504: Note:
505: Supports left preconditioning only.
507: Contributed by:
508: Patrick Sanan
510: References:
511: + * - Notay, Y."Flexible Conjugate Gradients", SIAM J. Sci. Comput. 22:4, 2000
512: - * - Axelsson, O. and Vassilevski, P. S. "A Black Box Generalized Conjugate Gradient Solver with Inner Iterations and Variable step Preconditioning",
513: SIAM J. Matrix Anal. Appl. 12:4, 1991
515: .seealso: [](ch_ksp), [](sec_flexibleksp), `KSPGCR`, `KSPFGMRES`, `KSPCG`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`, `KSPFCGSetNprealloc()`, `KSPFCGGetNprealloc()`, `KSPFCGSetTruncationType()`, `KSPFCGGetTruncationType()`,
516: `KSPFCGGetTruncationType`
517: M*/
518: PETSC_EXTERN PetscErrorCode KSPCreate_FCG(KSP ksp)
519: {
520: KSP_FCG *fcg;
522: PetscFunctionBegin;
523: PetscCall(PetscNew(&fcg));
524: #if !defined(PETSC_USE_COMPLEX)
525: fcg->type = KSP_CG_SYMMETRIC;
526: #else
527: fcg->type = KSP_CG_HERMITIAN;
528: #endif
529: fcg->mmax = KSPFCG_DEFAULT_MMAX;
530: fcg->nprealloc = KSPFCG_DEFAULT_NPREALLOC;
531: fcg->nvecs = 0;
532: fcg->vecb = KSPFCG_DEFAULT_VECB;
533: fcg->nchunks = 0;
534: fcg->truncstrat = KSPFCG_DEFAULT_TRUNCSTRAT;
536: ksp->data = (void *)fcg;
538: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
539: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 1));
540: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 1));
541: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
543: ksp->ops->setup = KSPSetUp_FCG;
544: ksp->ops->solve = KSPSolve_FCG;
545: ksp->ops->destroy = KSPDestroy_FCG;
546: ksp->ops->view = KSPView_FCG;
547: ksp->ops->setfromoptions = KSPSetFromOptions_FCG;
548: ksp->ops->buildsolution = KSPBuildSolutionDefault;
549: ksp->ops->buildresidual = KSPBuildResidualDefault;
550: PetscFunctionReturn(PETSC_SUCCESS);
551: }