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:   /* 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:     KSPCreateVecs(ksp, nnewvecs, &fcg->pCvecs[fcg->nchunks], 0, NULL);
 25:     KSPCreateVecs(ksp, nnewvecs, &fcg->pPvecs[fcg->nchunks], 0, NULL);
 26:     fcg->nvecs += nnewvecs;
 27:     for (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:   return 0;
 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;


 44:   /* Allocate "standard" work vectors (not including the basis and transformed basis vectors) */
 45:   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:   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) 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:   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:     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:   return 0;
 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, 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 */


 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:   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:     KSP_MatMult(ksp, Amat, X, R);
102:     VecAYPX(R, -1.0, B); /*   r <- b - Ax     */
103:   } else {
104:     VecCopy(B, R); /*   r <- b (x is 0) */
105:   }
106:   switch (ksp->normtype) {
107:   case KSP_NORM_PRECONDITIONED:
108:     KSP_PCApply(ksp, R, Z);  /*   z <- Br         */
109:     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:     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:     KSP_PCApply(ksp, R, Z); /*   z <- Br         */
118:     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:   KSPLogResidualHistory(ksp, dp);
131:   KSPMonitor(ksp, 0, dp);
132:   ksp->rnorm = dp;
133:   if (ksp->normtype == KSP_NORM_NONE) {
134:     KSPConvergedSkip(ksp, 0, dp, &ksp->reason, ksp->cnvP);
135:   } else {
136:     (*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP);
137:   }
138:   if (ksp->reason) return 0;

140:   /* Apply PC if not already done for convergence check */
141:   if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->normtype == KSP_NORM_NONE) { 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:     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:     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:         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:         VecXMDot(Z, ndots, Cold, dots);
187:         for (k = 0; k < ndots; ++k) dots[k] = -dots[k];
188:         VecMAXPY(Pcurr, ndots, dots, Pold);
189:         PetscFree3(dots, Cold, Pold);
190:       }
191:     }

193:     /* Update X and R */
194:     betaold = beta;
195:     VecXDot(Pcurr, R, &beta); /*  beta <- pi'*r       */
196:     KSPCheckDot(ksp, beta);
197:     KSP_MatMult(ksp, Amat, Pcurr, Ccurr); /*  w <- A*pi (stored in ci)   */
198:     VecXDot(Pcurr, Ccurr, &dpi);          /*  dpi <- pi'*w        */
199:     alphaold = alpha;
200:     alpha    = beta / dpi;                /*  alpha <- beta/dpi    */
201:     VecAXPY(X, alpha, Pcurr);  /*  x <- x + alpha * pi  */
202:     VecAXPY(R, -alpha, Ccurr); /*  r <- r - alpha * wi  */

204:     /* Compute norm for convergence check */
205:     switch (ksp->normtype) {
206:     case KSP_NORM_PRECONDITIONED:
207:       KSP_PCApply(ksp, R, Z);  /*   z <- Br             */
208:       VecNorm(Z, NORM_2, &dp); /*   dp <- sqrt(z'*z) = sqrt(e'*A'*B'*B*A*e)  */
209:       KSPCheckNorm(ksp, dp);
210:       break;
211:     case KSP_NORM_UNPRECONDITIONED:
212:       VecNorm(R, NORM_2, &dp); /*   dp <- sqrt(r'*r) = sqrt(e'*A'*A*e)   */
213:       KSPCheckNorm(ksp, dp);
214:       break;
215:     case KSP_NORM_NATURAL:
216:       KSP_PCApply(ksp, R, Z); /*   z <- Br             */
217:       VecXDot(R, Z, &s);
218:       KSPCheckDot(ksp, s);
219:       dp = PetscSqrtReal(PetscAbsScalar(s)); /*   dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e)  */
220:       break;
221:     case KSP_NORM_NONE:
222:       dp = 0.0;
223:       break;
224:     default:
225:       SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
226:     }

228:     /* Check for convergence */
229:     ksp->rnorm = dp;
230:     KSPLogResidualHistory(ksp, dp);
231:     KSPMonitor(ksp, i + 1, dp);
232:     (*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP);
233:     if (ksp->reason) break;

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

238:     /* Compute current C (which is W/dpi) */
239:     VecScale(Ccurr, 1.0 / dpi); /*   w <- ci/dpi   */

241:     if (eigs) {
242:       if (i > 0) {
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:       fcg->ned = ksp->its - 1;
250:     }
251:     ++i;
252:   } while (i < ksp->max_it);
253:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
254:   return 0;
255: }

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


263:   /* Destroy "standard" work vecs */
264:   VecDestroyVecs(ksp->nwork, &ksp->work);

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

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

286:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
287:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring);

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

293:   if (iascii) {
294:     PetscViewerASCIIPrintf(viewer, "  m_max=%" PetscInt_FMT "\n", fcg->mmax);
295:     PetscViewerASCIIPrintf(viewer, "  preallocated %" PetscInt_FMT " directions\n", PetscMin(fcg->nprealloc, fcg->mmax + 1));
296:     PetscViewerASCIIPrintf(viewer, "  %s\n", truncstr);
297:   } else if (isstring) {
298:     PetscViewerStringSPrintf(viewer, "m_max %" PetscInt_FMT " nprealloc %" PetscInt_FMT " %s", fcg->mmax, fcg->nprealloc, truncstr);
299:   }
300:   return 0;
301: }

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

306:   Logically Collective

308:   Input Parameters:
309: +  ksp - the Krylov space context
310: -  mmax - the maximum number of previous directions to orthogonalize against

312:   Options Database Key:
313: .   -ksp_fcg_mmax <N>  - maximum number of search directions

315:   Level: intermediate

317:   Note:
318:    mmax + 1 directions are stored (mmax previous ones along with a current one)
319:   and whether all are used in each iteration also depends on the truncation strategy
320:   (see KSPFCGSetTruncationType())

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

330:   fcg->mmax = mmax;
331:   return 0;
332: }

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

337:    Not Collective

339:    Input Parameter:
340: .  ksp - the Krylov space context

342:    Output Parameter:
343: .  mmax - the maximum number of previous directions allowed for orthogonalization

345:    Level: intermediate

347:   Note:
348:   FCG stores mmax+1 directions at most (mmax previous ones, and one current one)

350: .seealso: [](chapter_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGGetNprealloc()`, `KSPFCGSetMmax()`
351: @*/

353: PetscErrorCode KSPFCGGetMmax(KSP ksp, PetscInt *mmax)
354: {
355:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

358:   *mmax = fcg->mmax;
359:   return 0;
360: }

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

365:   Logically Collective

367:   Input Parameters:
368: +  ksp - the Krylov space context
369: -  nprealloc - the number of vectors to preallocate

371:   Options Database Key:
372: . -ksp_fcg_nprealloc <N> - number of directions to preallocate

374:   Level: advanced

376: .seealso: [](chapter_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGGetNprealloc()`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`
377: @*/
378: PetscErrorCode KSPFCGSetNprealloc(KSP ksp, PetscInt nprealloc)
379: {
380:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

385:   fcg->nprealloc = nprealloc;
386:   return 0;
387: }

389: /*@
390:   KSPFCGGetNprealloc - get the number of directions preallocate by `KSPFCG`

392:    Not Collective

394:    Input Parameter:
395: .  ksp - the Krylov space context

397:    Output Parameter:
398: .  nprealloc - the number of directions preallocated

400:    Level: advanced

402: .seealso: [](chapter_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGSetNprealloc()`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`
403: @*/
404: PetscErrorCode KSPFCGGetNprealloc(KSP ksp, PetscInt *nprealloc)
405: {
406:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

409:   *nprealloc = fcg->nprealloc;
410:   return 0;
411: }

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

416:   Logically Collective

418:   Input Parameters:
419: +  ksp - the Krylov space context
420: -  truncstrat - the choice of strategy
421: .vb
422:   KSP_FCD_TRUNC_TYPE_STANDARD uses all (up to mmax) stored directions
423:   KSP_FCD_TRUNC_TYPE_NOTAY uses the last max(1,mod(i,mmax)) stored directions at iteration i=0,1,..
424: .ve

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

429:   Level: intermediate

431: .seealso: [](chapter_ksp), `KSPFCDTruncationType`, `KSPFCGGetTruncationType`, `KSPFCGSetNprealloc()`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`
432: @*/
433: PetscErrorCode KSPFCGSetTruncationType(KSP ksp, KSPFCDTruncationType truncstrat)
434: {
435:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

439:   fcg->truncstrat = truncstrat;
440:   return 0;
441: }

443: /*@
444:   KSPFCGGetTruncationType - get the truncation strategy employed by `KSPFCG`

446:    Not Collective

448:    Input Parameter:
449: .  ksp - the Krylov space context

451:    Output Parameter:
452: .  truncstrat - the strategy type

454:    Level: intermediate

456: .seealso: [](chapter_ksp), `KSPFCG`, `KSPFCGSetTruncationType`, `KSPFCDTruncationType`, `KSPFCGSetTruncationType()`
457: @*/
458: PetscErrorCode KSPFCGGetTruncationType(KSP ksp, KSPFCDTruncationType *truncstrat)
459: {
460:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

463:   *truncstrat = fcg->truncstrat;
464:   return 0;
465: }

467: static PetscErrorCode KSPSetFromOptions_FCG(KSP ksp, PetscOptionItems *PetscOptionsObject)
468: {
469:   KSP_FCG  *fcg = (KSP_FCG *)ksp->data;
470:   PetscInt  mmax, nprealloc;
471:   PetscBool flg;

473:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP FCG Options");
474:   PetscOptionsInt("-ksp_fcg_mmax", "Maximum number of search directions to store", "KSPFCGSetMmax", fcg->mmax, &mmax, &flg);
475:   if (flg) KSPFCGSetMmax(ksp, mmax);
476:   PetscOptionsInt("-ksp_fcg_nprealloc", "Number of directions to preallocate", "KSPFCGSetNprealloc", fcg->nprealloc, &nprealloc, &flg);
477:   if (flg) KSPFCGSetNprealloc(ksp, nprealloc);
478:   PetscOptionsEnum("-ksp_fcg_truncation_type", "Truncation approach for directions", "KSPFCGSetTruncationType", KSPFCDTruncationTypes, (PetscEnum)fcg->truncstrat, (PetscEnum *)&fcg->truncstrat, NULL);
479:   PetscOptionsHeadEnd();
480:   return 0;
481: }

483: /*MC
484:       KSPFCG - Implements the Flexible Conjugate Gradient method (FCG). Unlike most `KSP` methods this allows the preconditioner to be nonlinear. [](sec_flexibleksp)

486:   Options Database Keys:
487: +   -ksp_fcg_mmax <N>  - maximum number of search directions
488: .   -ksp_fcg_nprealloc <N> - number of directions to preallocate
489: -   -ksp_fcg_truncation_type <standard,notay> - truncation approach for directions

491:   Level: beginner

493:    Note:
494:    Supports left preconditioning only.

496:   Contributed by:
497:   Patrick Sanan

499:   References:
500: + * - Notay, Y."Flexible Conjugate Gradients", SIAM J. Sci. Comput. 22:4, 2000
501: - * - Axelsson, O. and Vassilevski, P. S. "A Black Box Generalized Conjugate Gradient Solver with Inner Iterations and Variable step Preconditioning",
502:     SIAM J. Matrix Anal. Appl. 12:4, 1991

504: .seealso: [](chapter_ksp), [](sec_flexibleksp), `KSPGCR`, `KSPFGMRES`, `KSPCG`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`, `KSPFCGSetNprealloc()`, `KSPFCGGetNprealloc()`, `KSPFCGSetTruncationType()`, `KSPFCGGetTruncationType()`,
505:            `KSPFCGGetTruncationType`
506: M*/
507: PETSC_EXTERN PetscErrorCode KSPCreate_FCG(KSP ksp)
508: {
509:   KSP_FCG *fcg;

511:   PetscNew(&fcg);
512: #if !defined(PETSC_USE_COMPLEX)
513:   fcg->type = KSP_CG_SYMMETRIC;
514: #else
515:   fcg->type = KSP_CG_HERMITIAN;
516: #endif
517:   fcg->mmax       = KSPFCG_DEFAULT_MMAX;
518:   fcg->nprealloc  = KSPFCG_DEFAULT_NPREALLOC;
519:   fcg->nvecs      = 0;
520:   fcg->vecb       = KSPFCG_DEFAULT_VECB;
521:   fcg->nchunks    = 0;
522:   fcg->truncstrat = KSPFCG_DEFAULT_TRUNCSTRAT;

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

526:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 2);
527:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 1);
528:   KSPSetSupportedNorm(ksp, KSP_NORM_NATURAL, PC_LEFT, 1);
529:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);

531:   ksp->ops->setup          = KSPSetUp_FCG;
532:   ksp->ops->solve          = KSPSolve_FCG;
533:   ksp->ops->destroy        = KSPDestroy_FCG;
534:   ksp->ops->view           = KSPView_FCG;
535:   ksp->ops->setfromoptions = KSPSetFromOptions_FCG;
536:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
537:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
538:   return 0;
539: }