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;
     45:   /* Allocate "standard" work vectors (not including the basis and transformed basis vectors) */
     46:   PetscCall(KSPSetWorkVecs(ksp, nworkstd));

     48:   /* Allocated space for pointers to additional work vectors
     49:    note that mmax is the number of previous directions, so we add 1 for the current direction,
     50:    and an extra 1 for the prealloc (which might be empty) */
     51:   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));

     53:   /* If the requested number of preallocated vectors is greater than mmax reduce nprealloc */
     54:   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));

     56:   /* Preallocate additional work vectors */
     57:   PetscCall(KSPAllocateVectors_FCG(ksp, fcg->nprealloc, fcg->nprealloc));
     58:   /*
     59:   If user requested computations of eigenvalues then allocate work
     60:   work space needed
     61:   */
     62:   if (ksp->calc_sings) {
     63:     /* get space to store tridiagonal matrix for Lanczos */
     64:     PetscCall(PetscMalloc4(maxit, &fcg->e, maxit, &fcg->d, maxit, &fcg->ee, maxit, &fcg->dd));

     66:     ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_CG;
     67:     ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_CG;
     68:   }
     69:   PetscFunctionReturn(PETSC_SUCCESS);
     70: }

     72: static PetscErrorCode KSPSolve_FCG(KSP ksp)
     73: {
     74:   PetscInt    i, k, idx, mi;
     75:   KSP_FCG    *fcg   = (KSP_FCG *)ksp->data;
     76:   PetscScalar alpha = 0.0, beta = 0.0, dpi, dpiold = 0.0, s;
     77:   PetscReal   dp = 0.0;
     78:   Vec         B, R, Z, X, Pcurr, Ccurr;
     79:   Mat         Amat, Pmat;
     80:   PetscInt    eigs          = ksp->calc_sings; /* Variables for eigen estimation - START*/
     81:   PetscInt    stored_max_it = ksp->max_it;
     82:   PetscScalar alphaold = 0, betaold = 1.0, *e = NULL, *d = NULL; /* Variables for eigen estimation  - FINISH */

     84:   PetscFunctionBegin;
     85: #define VecXDot(x, y, a)     (fcg->type == KSP_CG_HERMITIAN ? VecDot(x, y, a) : VecTDot(x, y, a))
     86: #define VecXMDot(a, b, c, d) (fcg->type == KSP_CG_HERMITIAN ? VecMDot(a, b, c, d) : VecMTDot(a, b, c, d))

     88:   X = ksp->vec_sol;
     89:   B = ksp->vec_rhs;
     90:   R = ksp->work[0];
     91:   Z = ksp->work[1];

     93:   PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
     94:   if (eigs) {
     95:     e    = fcg->e;
     96:     d    = fcg->d;
     97:     e[0] = 0.0;
     98:   }
     99:   /* Compute initial residual needed for convergence check*/
    100:   ksp->its = 0;
    101:   if (!ksp->guess_zero) {
    102:     PetscCall(KSP_MatMult(ksp, Amat, X, R));
    103:     PetscCall(VecAYPX(R, -1.0, B)); /*   r <- b - Ax     */
    104:   } else {
    105:     PetscCall(VecCopy(B, R)); /*   r <- b (x is 0) */
    106:   }
    107:   switch (ksp->normtype) {
    108:   case KSP_NORM_PRECONDITIONED:
    109:     PetscCall(KSP_PCApply(ksp, R, Z));  /*   z <- Br         */
    110:     PetscCall(VecNorm(Z, NORM_2, &dp)); /*   dp <- dqrt(z'*z) = sqrt(e'*A'*B'*B*A*e)     */
    111:     KSPCheckNorm(ksp, dp);
    112:     break;
    113:   case KSP_NORM_UNPRECONDITIONED:
    114:     PetscCall(VecNorm(R, NORM_2, &dp)); /*   dp <- sqrt(r'*r) = sqrt(e'*A'*A*e)     */
    115:     KSPCheckNorm(ksp, dp);
    116:     break;
    117:   case KSP_NORM_NATURAL:
    118:     PetscCall(KSP_PCApply(ksp, R, Z)); /*   z <- Br         */
    119:     PetscCall(VecXDot(R, Z, &s));
    120:     KSPCheckDot(ksp, s);
    121:     dp = PetscSqrtReal(PetscAbsScalar(s)); /*   dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e)  */
    122:     break;
    123:   case KSP_NORM_NONE:
    124:     dp = 0.0;
    125:     break;
    126:   default:
    127:     SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
    128:   }

    130:   /* Initial Convergence Check */
    131:   PetscCall(KSPLogResidualHistory(ksp, dp));
    132:   PetscCall(KSPMonitor(ksp, 0, dp));
    133:   ksp->rnorm = dp;
    134:   if (ksp->normtype == KSP_NORM_NONE) {
    135:     PetscCall(KSPConvergedSkip(ksp, 0, dp, &ksp->reason, ksp->cnvP));
    136:   } else {
    137:     PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP));
    138:   }
    139:   if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);

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

    144:   i = 0;
    145:   do {
    146:     ksp->its = i + 1;

    148:     /*  If needbe, allocate a new chunk of vectors in P and C */
    149:     PetscCall(KSPAllocateVectors_FCG(ksp, i + 1, fcg->vecb));

    151:     /* Note that we wrap around and start clobbering old vectors */
    152:     idx   = i % (fcg->mmax + 1);
    153:     Pcurr = fcg->Pvecs[idx];
    154:     Ccurr = fcg->Cvecs[idx];

    156:     /* number of old directions to orthogonalize against */
    157:     switch (fcg->truncstrat) {
    158:     case KSP_FCD_TRUNC_TYPE_STANDARD:
    159:       mi = fcg->mmax;
    160:       break;
    161:     case KSP_FCD_TRUNC_TYPE_NOTAY:
    162:       mi = ((i - 1) % fcg->mmax) + 1;
    163:       break;
    164:     default:
    165:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unrecognized Truncation Strategy");
    166:     }

    168:     /* Compute a new column of P (Currently does not support modified G-S or iterative refinement)*/
    169:     PetscCall(VecCopy(Z, Pcurr));

    171:     {
    172:       PetscInt l, ndots;

    174:       l     = PetscMax(0, i - mi);
    175:       ndots = i - l;
    176:       if (ndots) {
    177:         PetscInt     j;
    178:         Vec         *Pold, *Cold;
    179:         PetscScalar *dots;

    181:         PetscCall(PetscMalloc3(ndots, &dots, ndots, &Cold, ndots, &Pold));
    182:         for (k = l, j = 0; j < ndots; ++k, ++j) {
    183:           idx     = k % (fcg->mmax + 1);
    184:           Cold[j] = fcg->Cvecs[idx];
    185:           Pold[j] = fcg->Pvecs[idx];
    186:         }
    187:         PetscCall(VecXMDot(Z, ndots, Cold, dots));
    188:         for (k = 0; k < ndots; ++k) dots[k] = -dots[k];
    189:         PetscCall(VecMAXPY(Pcurr, ndots, dots, Pold));
    190:         PetscCall(PetscFree3(dots, Cold, Pold));
    191:       }
    192:     }

    194:     /* Update X and R */
    195:     betaold = beta;
    196:     PetscCall(VecXDot(Pcurr, R, &beta)); /*  beta <- pi'*r       */
    197:     KSPCheckDot(ksp, beta);
    198:     if ((i > 0) && (PetscAbsScalar(beta * betaold) < 0.0)) {
    199:       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));
    200:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
    201:       PetscCall(PetscInfo(ksp, "diverging due to indefinite preconditioner\n"));
    202:       break;
    203:     }
    204:     PetscCall(KSP_MatMult(ksp, Amat, Pcurr, Ccurr)); /*  w <- A*pi (stored in ci)   */
    205:     dpiold = dpi;
    206:     PetscCall(VecXDot(Pcurr, Ccurr, &dpi)); /*  dpi <- pi'*w        */
    207:     if ((dpi == 0.0) || ((i > 0) && ((PetscSign(PetscRealPart(dpi)) * PetscSign(PetscRealPart(dpiold))) < 0.0))) {
    208:       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));
    209:       ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
    210:       PetscCall(PetscInfo(ksp, "diverging due to indefinite matrix\n"));
    211:       break;
    212:     }
    213:     alphaold = alpha;
    214:     alpha    = beta / dpi;                /*  alpha <- beta/dpi    */
    215:     PetscCall(VecAXPY(X, alpha, Pcurr));  /*  x <- x + alpha * pi  */
    216:     PetscCall(VecAXPY(R, -alpha, Ccurr)); /*  r <- r - alpha * wi  */

    218:     /* Compute norm for convergence check */
    219:     switch (ksp->normtype) {
    220:     case KSP_NORM_PRECONDITIONED:
    221:       PetscCall(KSP_PCApply(ksp, R, Z));  /*   z <- Br             */
    222:       PetscCall(VecNorm(Z, NORM_2, &dp)); /*   dp <- sqrt(z'*z) = sqrt(e'*A'*B'*B*A*e)  */
    223:       KSPCheckNorm(ksp, dp);
    224:       break;
    225:     case KSP_NORM_UNPRECONDITIONED:
    226:       PetscCall(VecNorm(R, NORM_2, &dp)); /*   dp <- sqrt(r'*r) = sqrt(e'*A'*A*e)   */
    227:       KSPCheckNorm(ksp, dp);
    228:       break;
    229:     case KSP_NORM_NATURAL:
    230:       PetscCall(KSP_PCApply(ksp, R, Z)); /*   z <- Br             */
    231:       PetscCall(VecXDot(R, Z, &s));
    232:       KSPCheckDot(ksp, s);
    233:       dp = PetscSqrtReal(PetscAbsScalar(s)); /*   dp <- sqrt(r'*z) = sqrt(e'*A'*B*A*e)  */
    234:       break;
    235:     case KSP_NORM_NONE:
    236:       dp = 0.0;
    237:       break;
    238:     default:
    239:       SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
    240:     }

    242:     if (eigs) {
    243:       if (i > 0) {
    244:         PetscCheck(ksp->max_it == stored_max_it, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Can not change maxit AND calculate eigenvalues");
    245:         e[i] = PetscSqrtReal(PetscAbsScalar(beta / betaold)) / alphaold;
    246:         d[i] = PetscSqrtReal(PetscAbsScalar(beta / betaold)) * e[i] + 1.0 / alpha;
    247:       } else {
    248:         d[i] = PetscSqrtReal(PetscAbsScalar(beta)) * e[i] + 1.0 / alpha;
    249:       }
    250:     }

    252:     /* Check for convergence */
    253:     ksp->rnorm = dp;
    254:     PetscCall(KSPLogResidualHistory(ksp, dp));
    255:     PetscCall(KSPMonitor(ksp, i + 1, dp));
    256:     PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
    257:     if (ksp->reason) break;

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

    262:     /* Compute current C (which is W/dpi) */
    263:     PetscCall(VecScale(Ccurr, 1.0 / dpi)); /*   w <- ci/dpi   */
    264:     ++i;
    265:   } while (i < ksp->max_it);
    266:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
    267:   PetscFunctionReturn(PETSC_SUCCESS);
    268: }

    270: static PetscErrorCode KSPDestroy_FCG(KSP ksp)
    271: {
    272:   PetscInt i;
    273:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

    275:   PetscFunctionBegin;
    276:   /* Destroy "standard" work vecs */
    277:   PetscCall(VecDestroyVecs(ksp->nwork, &ksp->work));

    279:   /* Destroy P and C vectors and the arrays that manage pointers to them */
    280:   if (fcg->nvecs) {
    281:     for (i = 0; i < fcg->nchunks; ++i) {
    282:       PetscCall(VecDestroyVecs(fcg->chunksizes[i], &fcg->pPvecs[i]));
    283:       PetscCall(VecDestroyVecs(fcg->chunksizes[i], &fcg->pCvecs[i]));
    284:     }
    285:   }
    286:   PetscCall(PetscFree5(fcg->Pvecs, fcg->Cvecs, fcg->pPvecs, fcg->pCvecs, fcg->chunksizes));
    287:   /* free space used for singular value calculations */
    288:   if (ksp->calc_sings) PetscCall(PetscFree4(fcg->e, fcg->d, fcg->ee, fcg->dd));
    289:   PetscCall(KSPDestroyDefault(ksp));
    290:   PetscFunctionReturn(PETSC_SUCCESS);
    291: }

    293: static PetscErrorCode KSPView_FCG(KSP ksp, PetscViewer viewer)
    294: {
    295:   KSP_FCG    *fcg = (KSP_FCG *)ksp->data;
    296:   PetscBool   iascii, isstring;
    297:   const char *truncstr;

    299:   PetscFunctionBegin;
    300:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    301:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));

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

    307:   if (iascii) {
    308:     PetscCall(PetscViewerASCIIPrintf(viewer, "  m_max=%" PetscInt_FMT "\n", fcg->mmax));
    309:     PetscCall(PetscViewerASCIIPrintf(viewer, "  preallocated %" PetscInt_FMT " directions\n", PetscMin(fcg->nprealloc, fcg->mmax + 1)));
    310:     PetscCall(PetscViewerASCIIPrintf(viewer, "  %s\n", truncstr));
    311:   } else if (isstring) {
    312:     PetscCall(PetscViewerStringSPrintf(viewer, "m_max %" PetscInt_FMT " nprealloc %" PetscInt_FMT " %s", fcg->mmax, fcg->nprealloc, truncstr));
    313:   }
    314:   PetscFunctionReturn(PETSC_SUCCESS);
    315: }

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

    320:   Logically Collective

    322:   Input Parameters:
    323: + ksp  - the Krylov space context
    324: - mmax - the maximum number of previous directions to orthogonalize against

    326:   Options Database Key:
    327: . -ksp_fcg_mmax <N> - maximum number of search directions

    329:   Level: intermediate

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

    335: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGGetNprealloc()`, `KSPFCGetMmax()`
    336: @*/
    337: PetscErrorCode KSPFCGSetMmax(KSP ksp, PetscInt mmax)
    338: {
    339:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

    341:   PetscFunctionBegin;
    344:   fcg->mmax = mmax;
    345:   PetscFunctionReturn(PETSC_SUCCESS);
    346: }

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

    351:   Not Collective

    353:   Input Parameter:
    354: . ksp - the Krylov space context

    356:   Output Parameter:
    357: . mmax - the maximum number of previous directions allowed for orthogonalization

    359:   Level: intermediate

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

    364: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGGetNprealloc()`, `KSPFCGSetMmax()`
    365: @*/
    366: PetscErrorCode KSPFCGGetMmax(KSP ksp, PetscInt *mmax)
    367: {
    368:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

    370:   PetscFunctionBegin;
    372:   *mmax = fcg->mmax;
    373:   PetscFunctionReturn(PETSC_SUCCESS);
    374: }

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

    379:   Logically Collective

    381:   Input Parameters:
    382: + ksp       - the Krylov space context
    383: - nprealloc - the number of vectors to preallocate

    385:   Options Database Key:
    386: . -ksp_fcg_nprealloc <N> - number of directions to preallocate

    388:   Level: advanced

    390: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGGetNprealloc()`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`
    391: @*/
    392: PetscErrorCode KSPFCGSetNprealloc(KSP ksp, PetscInt nprealloc)
    393: {
    394:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

    396:   PetscFunctionBegin;
    399:   PetscCheck(nprealloc <= fcg->mmax + 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot preallocate more than m_max+1 vectors");
    400:   fcg->nprealloc = nprealloc;
    401:   PetscFunctionReturn(PETSC_SUCCESS);
    402: }

    404: /*@
    405:   KSPFCGGetNprealloc - get the number of directions preallocate by `KSPFCG`

    407:   Not Collective

    409:   Input Parameter:
    410: . ksp - the Krylov space context

    412:   Output Parameter:
    413: . nprealloc - the number of directions preallocated

    415:   Level: advanced

    417: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGGetTruncationType()`, `KSPFCGSetNprealloc()`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`
    418: @*/
    419: PetscErrorCode KSPFCGGetNprealloc(KSP ksp, PetscInt *nprealloc)
    420: {
    421:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

    423:   PetscFunctionBegin;
    425:   *nprealloc = fcg->nprealloc;
    426:   PetscFunctionReturn(PETSC_SUCCESS);
    427: }

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

    432:   Logically Collective

    434:   Input Parameters:
    435: + ksp        - the Krylov space context
    436: - truncstrat - the choice of strategy
    437: .vb
    438:   KSP_FCD_TRUNC_TYPE_STANDARD uses all (up to mmax) stored directions
    439:   KSP_FCD_TRUNC_TYPE_NOTAY uses the last max(1,mod(i,mmax)) stored directions at iteration i=0,1,..
    440: .ve

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

    445:   Level: intermediate

    447: .seealso: [](ch_ksp), `KSPFCDTruncationType`, `KSPFCGGetTruncationType()`, `KSPFCGSetNprealloc()`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`,
    448:           `KSP_FCD_TRUNC_TYPE_STANDARD`, `KSP_FCD_TRUNC_TYPE_NOTAY`
    449: @*/
    450: PetscErrorCode KSPFCGSetTruncationType(KSP ksp, KSPFCDTruncationType truncstrat)
    451: {
    452:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

    454:   PetscFunctionBegin;
    457:   fcg->truncstrat = truncstrat;
    458:   PetscFunctionReturn(PETSC_SUCCESS);
    459: }

    461: /*@
    462:   KSPFCGGetTruncationType - get the truncation strategy employed by `KSPFCG`

    464:   Not Collective

    466:   Input Parameter:
    467: . ksp - the Krylov space context

    469:   Output Parameter:
    470: . truncstrat - the strategy type

    472:   Level: intermediate

    474: .seealso: [](ch_ksp), `KSPFCG`, `KSPFCGSetTruncationType()`, `KSPFCDTruncationType`, `KSP_FCD_TRUNC_TYPE_STANDARD`, `KSP_FCD_TRUNC_TYPE_NOTAY`
    475: @*/
    476: PetscErrorCode KSPFCGGetTruncationType(KSP ksp, KSPFCDTruncationType *truncstrat)
    477: {
    478:   KSP_FCG *fcg = (KSP_FCG *)ksp->data;

    480:   PetscFunctionBegin;
    482:   *truncstrat = fcg->truncstrat;
    483:   PetscFunctionReturn(PETSC_SUCCESS);
    484: }

    486: static PetscErrorCode KSPSetFromOptions_FCG(KSP ksp, PetscOptionItems *PetscOptionsObject)
    487: {
    488:   KSP_FCG  *fcg = (KSP_FCG *)ksp->data;
    489:   PetscInt  mmax, nprealloc;
    490:   PetscBool flg;

    492:   PetscFunctionBegin;
    493:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP FCG Options");
    494:   PetscCall(PetscOptionsInt("-ksp_fcg_mmax", "Maximum number of search directions to store", "KSPFCGSetMmax", fcg->mmax, &mmax, &flg));
    495:   if (flg) PetscCall(KSPFCGSetMmax(ksp, mmax));
    496:   PetscCall(PetscOptionsInt("-ksp_fcg_nprealloc", "Number of directions to preallocate", "KSPFCGSetNprealloc", fcg->nprealloc, &nprealloc, &flg));
    497:   if (flg) PetscCall(KSPFCGSetNprealloc(ksp, nprealloc));
    498:   PetscCall(PetscOptionsEnum("-ksp_fcg_truncation_type", "Truncation approach for directions", "KSPFCGSetTruncationType", KSPFCDTruncationTypes, (PetscEnum)fcg->truncstrat, (PetscEnum *)&fcg->truncstrat, NULL));
    499:   PetscOptionsHeadEnd();
    500:   PetscFunctionReturn(PETSC_SUCCESS);
    501: }

    503: /*MC
    504:   KSPFCG - Implements the Flexible Conjugate Gradient method (FCG) {cite}`flexiblecg`, {cite}`generalizedcg`.
    505:   Unlike most `KSP` methods this allows the preconditioner to be nonlinear. [](sec_flexibleksp)

    507:   Options Database Keys:
    508: + -ksp_fcg_mmax <N>                         - maximum number of search directions, similar to the restart in `KSPGMRES` and `KSPFGMRES`
    509: . -ksp_fcg_nprealloc <N>                    - number of directions to preallocate
    510: - -ksp_fcg_truncation_type <standard,notay> - truncation approach for directions

    512:   Level: beginner

    514:   Notes:
    515:   `KSPFGMRES` provides a flexible GMRES which can be used for matrices that are not symmetric positive-definite (SPD).

    517:   Supports left preconditioning only.

    519:   Contributed by:
    520:   Patrick Sanan

    522: .seealso: [](ch_ksp), [](sec_flexibleksp), `KSPGCR`, `KSPFGMRES`, `KSPCG`, `KSPFCGSetMmax()`, `KSPFCGGetMmax()`, `KSPFCGSetNprealloc()`, `KSPFCGGetNprealloc()`, `KSPFCGSetTruncationType()`, `KSPFCGGetTruncationType()`,
    523:            `KSPFCGGetTruncationType`
    524: M*/
    525: PETSC_EXTERN PetscErrorCode KSPCreate_FCG(KSP ksp)
    526: {
    527:   KSP_FCG *fcg;

    529:   PetscFunctionBegin;
    530:   PetscCall(PetscNew(&fcg));
    531: #if !defined(PETSC_USE_COMPLEX)
    532:   fcg->type = KSP_CG_SYMMETRIC;
    533: #else
    534:   fcg->type = KSP_CG_HERMITIAN;
    535: #endif
    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->destroy        = KSPDestroy_FCG;
    553:   ksp->ops->view           = KSPView_FCG;
    554:   ksp->ops->setfromoptions = KSPSetFromOptions_FCG;
    555:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
    556:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
    557:   PetscFunctionReturn(PETSC_SUCCESS);
    558: }