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:     if (eigs) {
232:       if (i > 0) {
233:         PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Can not change maxit AND calculate eigenvalues");
234:         e[i] = PetscSqrtReal(PetscAbsScalar(beta / betaold)) / alphaold;
235:         d[i] = PetscSqrtReal(PetscAbsScalar(beta / betaold)) * e[i] + 1.0 / alpha;
236:       } else {
237:         d[i] = PetscSqrtReal(PetscAbsScalar(beta)) * e[i] + 1.0 / alpha;
238:       }
239:     }

241:     /* Check for convergence */
242:     ksp->rnorm = dp;
243:     PetscCall(KSPLogResidualHistory(ksp, dp));
244:     PetscCall(KSPMonitor(ksp, i + 1, dp));
245:     PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
246:     if (ksp->reason) break;

248:     /* Apply PC if not already done for convergence check */
249:     if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->normtype == KSP_NORM_NONE) { PetscCall(KSP_PCApply(ksp, R, Z)); /*   z <- Br         */ }

251:     /* Compute current C (which is W/dpi) */
252:     PetscCall(VecScale(Ccurr, 1.0 / dpi)); /*   w <- ci/dpi   */
253:     ++i;
254:   } while (i < ksp->max_it);
255:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
256:   PetscFunctionReturn(PETSC_SUCCESS);
257: }

259: static PetscErrorCode KSPDestroy_FCG(KSP ksp)
260: {
261:   PetscInt i;
262:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

264:   PetscFunctionBegin;

266:   /* Destroy "standard" work vecs */
267:   PetscCall(VecDestroyVecs(ksp->nwork, &ksp->work));

269:   /* Destroy P and C vectors and the arrays that manage pointers to them */
270:   if (fcg->nvecs) {
271:     for (i = 0; i < fcg->nchunks; ++i) {
272:       PetscCall(VecDestroyVecs(fcg->chunksizes[i], &fcg->pPvecs[i]));
273:       PetscCall(VecDestroyVecs(fcg->chunksizes[i], &fcg->pCvecs[i]));
274:     }
275:   }
276:   PetscCall(PetscFree5(fcg->Pvecs, fcg->Cvecs, fcg->pPvecs, fcg->pCvecs, fcg->chunksizes));
277:   /* free space used for singular value calculations */
278:   if (ksp->calc_sings) PetscCall(PetscFree4(fcg->e, fcg->d, fcg->ee, fcg->dd));
279:   PetscCall(KSPDestroyDefault(ksp));
280:   PetscFunctionReturn(PETSC_SUCCESS);
281: }

283: static PetscErrorCode KSPView_FCG(KSP ksp, PetscViewer viewer)
284: {
285:   KSP_FCG    *fcg = (KSP_FCG *)ksp->data;
286:   PetscBool   iascii, isstring;
287:   const char *truncstr;

289:   PetscFunctionBegin;
290:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
291:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));

293:   if (fcg->truncstrat == KSP_FCD_TRUNC_TYPE_STANDARD) truncstr = "Using standard truncation strategy";
294:   else if (fcg->truncstrat == KSP_FCD_TRUNC_TYPE_NOTAY) truncstr = "Using Notay's truncation strategy";
295:   else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Undefined FCG truncation strategy");

297:   if (iascii) {
298:     PetscCall(PetscViewerASCIIPrintf(viewer, "  m_max=%" PetscInt_FMT "\n", fcg->mmax));
299:     PetscCall(PetscViewerASCIIPrintf(viewer, "  preallocated %" PetscInt_FMT " directions\n", PetscMin(fcg->nprealloc, fcg->mmax + 1)));
300:     PetscCall(PetscViewerASCIIPrintf(viewer, "  %s\n", truncstr));
301:   } else if (isstring) {
302:     PetscCall(PetscViewerStringSPrintf(viewer, "m_max %" PetscInt_FMT " nprealloc %" PetscInt_FMT " %s", fcg->mmax, fcg->nprealloc, truncstr));
303:   }
304:   PetscFunctionReturn(PETSC_SUCCESS);
305: }

307: /*@
308:   KSPFCGSetMmax - set the maximum number of previous directions `KSPFCG` will store for orthogonalization

310:   Logically Collective

312:   Input Parameters:
313: + ksp  - the Krylov space context
314: - mmax - the maximum number of previous directions to orthogonalize against

316:   Options Database Key:
317: . -ksp_fcg_mmax <N> - maximum number of search directions

319:   Level: intermediate

321:   Note:
322:   `mmax` + 1 directions are stored (`mmax` previous ones along with a current one)
323:   and whether all are used in each iteration also depends on the truncation strategy, see `KSPFCGSetTruncationType()`

325: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGGetNprealloc()`, `KSPFCGetMmax()`
326: @*/
327: PetscErrorCode KSPFCGSetMmax(KSP ksp, PetscInt mmax)
328: {
329:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

331:   PetscFunctionBegin;
334:   fcg->mmax = mmax;
335:   PetscFunctionReturn(PETSC_SUCCESS);
336: }

338: /*@
339:   KSPFCGGetMmax - get the maximum number of previous directions `KSPFCG` will store

341:   Not Collective

343:   Input Parameter:
344: . ksp - the Krylov space context

346:   Output Parameter:
347: . mmax - the maximum number of previous directions allowed for orthogonalization

349:   Level: intermediate

351:   Note:
352:   `KSPFCG` stores `mmax`+1 directions at most (`mmax` previous ones, and one current one)

354: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGGetNprealloc()`, `KSPFCGSetMmax()`
355: @*/
356: PetscErrorCode KSPFCGGetMmax(KSP ksp, PetscInt *mmax)
357: {
358:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

360:   PetscFunctionBegin;
362:   *mmax = fcg->mmax;
363:   PetscFunctionReturn(PETSC_SUCCESS);
364: }

366: /*@
367:   KSPFCGSetNprealloc - set the number of directions to preallocate with `KSPFCG`

369:   Logically Collective

371:   Input Parameters:
372: + ksp       - the Krylov space context
373: - nprealloc - the number of vectors to preallocate

375:   Options Database Key:
376: . -ksp_fcg_nprealloc <N> - number of directions to preallocate

378:   Level: advanced

380: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGGetNprealloc()`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`
381: @*/
382: PetscErrorCode KSPFCGSetNprealloc(KSP ksp, PetscInt nprealloc)
383: {
384:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

386:   PetscFunctionBegin;
389:   PetscCheck(nprealloc <= fcg->mmax + 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot preallocate more than m_max+1 vectors");
390:   fcg->nprealloc = nprealloc;
391:   PetscFunctionReturn(PETSC_SUCCESS);
392: }

394: /*@
395:   KSPFCGGetNprealloc - get the number of directions preallocate by `KSPFCG`

397:   Not Collective

399:   Input Parameter:
400: . ksp - the Krylov space context

402:   Output Parameter:
403: . nprealloc - the number of directions preallocated

405:   Level: advanced

407: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGSetNprealloc()`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`
408: @*/
409: PetscErrorCode KSPFCGGetNprealloc(KSP ksp, PetscInt *nprealloc)
410: {
411:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

413:   PetscFunctionBegin;
415:   *nprealloc = fcg->nprealloc;
416:   PetscFunctionReturn(PETSC_SUCCESS);
417: }

419: /*@
420:   KSPFCGSetTruncationType - specify how many of its stored previous directions `KSPFCG` uses during orthoganalization

422:   Logically Collective

424:   Input Parameters:
425: + ksp        - the Krylov space context
426: - truncstrat - the choice of strategy
427: .vb
428:   KSP_FCD_TRUNC_TYPE_STANDARD uses all (up to mmax) stored directions
429:   KSP_FCD_TRUNC_TYPE_NOTAY uses the last max(1,mod(i,mmax)) stored directions at iteration i=0,1,..
430: .ve

432:   Options Database Key:
433: . -ksp_fcg_truncation_type <standard, notay> - specify how many of its stored previous directions `KSPFCG` uses during orthoganalization

435:   Level: intermediate

437: .seealso: [](ch_ksp), `KSPFCDTruncationType`, `KSPFCGGetTruncationType()`, `KSPFCGSetNprealloc()`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`,
438:           `KSP_FCD_TRUNC_TYPE_STANDARD`, `KSP_FCD_TRUNC_TYPE_NOTAY`
439: @*/
440: PetscErrorCode KSPFCGSetTruncationType(KSP ksp, KSPFCDTruncationType truncstrat)
441: {
442:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

444:   PetscFunctionBegin;
447:   fcg->truncstrat = truncstrat;
448:   PetscFunctionReturn(PETSC_SUCCESS);
449: }

451: /*@
452:   KSPFCGGetTruncationType - get the truncation strategy employed by `KSPFCG`

454:   Not Collective

456:   Input Parameter:
457: . ksp - the Krylov space context

459:   Output Parameter:
460: . truncstrat - the strategy type

462:   Level: intermediate

464: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGSetTruncationType()`, `KSPFCDTruncationType`, `KSP_FCD_TRUNC_TYPE_STANDARD`, `KSP_FCD_TRUNC_TYPE_NOTAY`
465: @*/
466: PetscErrorCode KSPFCGGetTruncationType(KSP ksp, KSPFCDTruncationType *truncstrat)
467: {
468:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

470:   PetscFunctionBegin;
472:   *truncstrat = fcg->truncstrat;
473:   PetscFunctionReturn(PETSC_SUCCESS);
474: }

476: static PetscErrorCode KSPSetFromOptions_FCG(KSP ksp, PetscOptionItems *PetscOptionsObject)
477: {
478:   KSP_FCG  *fcg = (KSP_FCG *)ksp->data;
479:   PetscInt  mmax, nprealloc;
480:   PetscBool flg;

482:   PetscFunctionBegin;
483:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP FCG Options");
484:   PetscCall(PetscOptionsInt("-ksp_fcg_mmax", "Maximum number of search directions to store", "KSPFCGSetMmax", fcg->mmax, &mmax, &flg));
485:   if (flg) PetscCall(KSPFCGSetMmax(ksp, mmax));
486:   PetscCall(PetscOptionsInt("-ksp_fcg_nprealloc", "Number of directions to preallocate", "KSPFCGSetNprealloc", fcg->nprealloc, &nprealloc, &flg));
487:   if (flg) PetscCall(KSPFCGSetNprealloc(ksp, nprealloc));
488:   PetscCall(PetscOptionsEnum("-ksp_fcg_truncation_type", "Truncation approach for directions", "KSPFCGSetTruncationType", KSPFCDTruncationTypes, (PetscEnum)fcg->truncstrat, (PetscEnum *)&fcg->truncstrat, NULL));
489:   PetscOptionsHeadEnd();
490:   PetscFunctionReturn(PETSC_SUCCESS);
491: }

493: /*MC
494:   KSPFCG - Implements the Flexible Conjugate Gradient method (FCG) {cite}`flexiblecg`, {cite}`generalizedcg`.
495:   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:   Notes:
505:   Compare to `KSPFCG`

507:   Supports left preconditioning only.

509:   Contributed by:
510:   Patrick Sanan

512: .seealso: [](ch_ksp), [](sec_flexibleksp), `KSPGCR`, `KSPFGMRES`, `KSPCG`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`, `KSPFCGSetNprealloc()`, `KSPFCGGetNprealloc()`, `KSPFCGSetTruncationType()`, `KSPFCGGetTruncationType()`,
513:            `KSPFCGGetTruncationType`
514: M*/
515: PETSC_EXTERN PetscErrorCode KSPCreate_FCG(KSP ksp)
516: {
517:   KSP_FCG *fcg;

519:   PetscFunctionBegin;
520:   PetscCall(PetscNew(&fcg));
521: #if !defined(PETSC_USE_COMPLEX)
522:   fcg->type = KSP_CG_SYMMETRIC;
523: #else
524:   fcg->type = KSP_CG_HERMITIAN;
525: #endif
526:   fcg->mmax       = KSPFCG_DEFAULT_MMAX;
527:   fcg->nprealloc  = KSPFCG_DEFAULT_NPREALLOC;
528:   fcg->nvecs      = 0;
529:   fcg->vecb       = KSPFCG_DEFAULT_VECB;
530:   fcg->nchunks    = 0;
531:   fcg->truncstrat = KSPFCG_DEFAULT_TRUNCSTRAT;

533:   ksp->data = (void *)fcg;

535:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2));
536:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 1));
537:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 1));
538:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));

540:   ksp->ops->setup          = KSPSetUp_FCG;
541:   ksp->ops->solve          = KSPSolve_FCG;
542:   ksp->ops->destroy        = KSPDestroy_FCG;
543:   ksp->ops->view           = KSPView_FCG;
544:   ksp->ops->setfromoptions = KSPSetFromOptions_FCG;
545:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
546:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
547:   PetscFunctionReturn(PETSC_SUCCESS);
548: }