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: }