Actual source code: fcg.c

      1: /*
      2:     This file implements the FCG (Flexible Conjugate Gradient) method

      4: */

      6: #include <../src/ksp/ksp/impls/fcg/fcgimpl.h>
      7: extern PetscErrorCode KSPComputeExtremeSingularValues_CG(KSP,PetscReal*,PetscReal*);
      8: extern PetscErrorCode KSPComputeEigenvalues_CG(KSP,PetscInt,PetscReal*,PetscReal*,PetscInt*);

     10: #define KSPFCG_DEFAULT_MMAX 30          /* maximum number of search directions to keep */
     11: #define KSPFCG_DEFAULT_NPREALLOC 10     /* number of search directions to preallocate */
     12: #define KSPFCG_DEFAULT_VECB 5           /* number of search directions to allocate each time new direction vectors are needed */
     13: #define KSPFCG_DEFAULT_TRUNCSTRAT KSP_FCD_TRUNC_TYPE_NOTAY

     15: static PetscErrorCode KSPAllocateVectors_FCG(KSP ksp, PetscInt nvecsneeded, PetscInt chunksize)
     16: {
     17:   PetscErrorCode  ierr;
     18:   PetscInt        i;
     19:   KSP_FCG         *fcg = (KSP_FCG*)ksp->data;
     20:   PetscInt        nnewvecs, nvecsprev;

     23:   /* Allocate enough new vectors to add chunksize new vectors, reach nvecsneedtotal, or to reach mmax+1, whichever is smallest */
     24:   if (fcg->nvecs < PetscMin(fcg->mmax+1,nvecsneeded)){
     25:     nvecsprev = fcg->nvecs;
     26:     nnewvecs = PetscMin(PetscMax(nvecsneeded-fcg->nvecs,chunksize),fcg->mmax+1-fcg->nvecs);
     27:     KSPCreateVecs(ksp,nnewvecs,&fcg->pCvecs[fcg->nchunks],0,NULL);
     28:     PetscLogObjectParents((PetscObject)ksp,nnewvecs,fcg->pCvecs[fcg->nchunks]);
     29:     KSPCreateVecs(ksp,nnewvecs,&fcg->pPvecs[fcg->nchunks],0,NULL);
     30:     PetscLogObjectParents((PetscObject)ksp,nnewvecs,fcg->pPvecs[fcg->nchunks]);
     31:     fcg->nvecs += nnewvecs;
     32:     for (i=0;i<nnewvecs;++i){
     33:       fcg->Cvecs[nvecsprev + i] = fcg->pCvecs[fcg->nchunks][i];
     34:       fcg->Pvecs[nvecsprev + i] = fcg->pPvecs[fcg->nchunks][i];
     35:     }
     36:     fcg->chunksizes[fcg->nchunks] = nnewvecs;
     37:     ++fcg->nchunks;
     38:   }
     39:   return(0);
     40: }

     42: static PetscErrorCode    KSPSetUp_FCG(KSP ksp)
     43: {
     45:   KSP_FCG        *fcg = (KSP_FCG*)ksp->data;
     46:   PetscInt       maxit = ksp->max_it;
     47:   const PetscInt nworkstd = 2;

     51:   /* Allocate "standard" work vectors (not including the basis and transformed basis vectors) */
     52:   KSPSetWorkVecs(ksp,nworkstd);

     54:   /* Allocated space for pointers to additional work vectors
     55:    note that mmax is the number of previous directions, so we add 1 for the current direction,
     56:    and an extra 1 for the prealloc (which might be empty) */
     57:   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);
     58:   PetscLogObjectMemory((PetscObject)ksp,2*(fcg->mmax+1)*sizeof(Vec*) + 2*(fcg->mmax + 1)*sizeof(Vec**) + (fcg->mmax + 2)*sizeof(PetscInt));

     60:   /* If the requested number of preallocated vectors is greater than mmax reduce nprealloc */
     61:   if (fcg->nprealloc > fcg->mmax+1){
     62:     PetscInfo2(NULL,"Requested nprealloc=%d is greater than m_max+1=%d. Resetting nprealloc = m_max+1.\n",fcg->nprealloc, fcg->mmax+1);
     63:   }

     65:   /* Preallocate additional work vectors */
     66:   KSPAllocateVectors_FCG(ksp,fcg->nprealloc,fcg->nprealloc);
     67:   /*
     68:   If user requested computations of eigenvalues then allocate work
     69:   work space needed
     70:   */
     71:   if (ksp->calc_sings) {
     72:     /* get space to store tridiagonal matrix for Lanczos */
     73:     PetscMalloc4(maxit,&fcg->e,maxit,&fcg->d,maxit,&fcg->ee,maxit,&fcg->dd);
     74:     PetscLogObjectMemory((PetscObject)ksp,2*(maxit+1)*(sizeof(PetscScalar)+sizeof(PetscReal)));

     76:     ksp->ops->computeextremesingularvalues = KSPComputeExtremeSingularValues_CG;
     77:     ksp->ops->computeeigenvalues           = KSPComputeEigenvalues_CG;
     78:   }
     79:   return(0);
     80: }

     82: static PetscErrorCode KSPSolve_FCG(KSP ksp)
     83: {
     85:   PetscInt       i,k,idx,mi;
     86:   KSP_FCG        *fcg = (KSP_FCG*)ksp->data;
     87:   PetscScalar    alpha=0.0,beta = 0.0,dpi,s;
     88:   PetscReal      dp=0.0;
     89:   Vec            B,R,Z,X,Pcurr,Ccurr;
     90:   Mat            Amat,Pmat;
     91:   PetscInt       eigs = ksp->calc_sings; /* Variables for eigen estimation - START*/
     92:   PetscInt       stored_max_it = ksp->max_it;
     93:   PetscScalar    alphaold = 0,betaold = 1.0,*e = NULL,*d = NULL;/* Variables for eigen estimation  - FINISH */

     97: #define VecXDot(x,y,a) (((fcg->type) == (KSP_CG_HERMITIAN)) ? VecDot(x,y,a) : VecTDot(x,y,a))
     98: #define VecXMDot(a,b,c,d) (((fcg->type) == (KSP_CG_HERMITIAN)) ? VecMDot(a,b,c,d) : VecMTDot(a,b,c,d))

    100:   X             = ksp->vec_sol;
    101:   B             = ksp->vec_rhs;
    102:   R             = ksp->work[0];
    103:   Z             = ksp->work[1];

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

    137:   /* Initial Convergence Check */
    138:   KSPLogResidualHistory(ksp,dp);
    139:   KSPMonitor(ksp,0,dp);
    140:   ksp->rnorm = dp;
    141:   if (ksp->normtype == KSP_NORM_NONE) {
    142:     KSPConvergedSkip(ksp,0,dp,&ksp->reason,ksp->cnvP);
    143:   } else {
    144:     (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);
    145:   }
    146:   if (ksp->reason) return(0);

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

    153:   i = 0;
    154:   do {
    155:     ksp->its = i+1;

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

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

    165:     /* number of old directions to orthogonalize against */
    166:     switch(fcg->truncstrat){
    167:       case KSP_FCD_TRUNC_TYPE_STANDARD:
    168:         mi = fcg->mmax;
    169:         break;
    170:       case KSP_FCD_TRUNC_TYPE_NOTAY:
    171:         mi = ((i-1) % fcg->mmax)+1;
    172:         break;
    173:       default:
    174:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unrecognized Truncation Strategy");
    175:     }

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

    180:     {
    181:       PetscInt l,ndots;

    183:       l = PetscMax(0,i-mi);
    184:       ndots = i-l;
    185:       if (ndots){
    186:         PetscInt    j;
    187:         Vec         *Pold,  *Cold;
    188:         PetscScalar *dots;

    190:         PetscMalloc3(ndots,&dots,ndots,&Cold,ndots,&Pold);
    191:         for (k=l,j=0;j<ndots;++k,++j){
    192:           idx = k % (fcg->mmax+1);
    193:           Cold[j] = fcg->Cvecs[idx];
    194:           Pold[j] = fcg->Pvecs[idx];
    195:         }
    196:         VecXMDot(Z,ndots,Cold,dots);
    197:         for (k=0;k<ndots;++k){
    198:           dots[k] = -dots[k];
    199:         }
    200:         VecMAXPY(Pcurr,ndots,dots,Pold);
    201:         PetscFree3(dots,Cold,Pold);
    202:       }
    203:     }

    205:     /* Update X and R */
    206:     betaold = beta;
    207:     VecXDot(Pcurr,R,&beta);                 /*  beta <- pi'*r       */
    208:     KSPCheckDot(ksp,beta);
    209:     KSP_MatMult(ksp,Amat,Pcurr,Ccurr);      /*  w <- A*pi (stored in ci)   */
    210:     VecXDot(Pcurr,Ccurr,&dpi);              /*  dpi <- pi'*w        */
    211:     alphaold = alpha;
    212:     alpha = beta / dpi;                                          /*  alpha <- beta/dpi    */
    213:     VecAXPY(X,alpha,Pcurr);                 /*  x <- x + alpha * pi  */
    214:     VecAXPY(R,-alpha,Ccurr);                /*  r <- r - alpha * wi  */

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

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

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

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

    254:     if (eigs) {
    255:       if (i > 0) {
    256:         if (ksp->max_it != stored_max_it) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Can not change maxit AND calculate eigenvalues");
    257:         e[i] = PetscSqrtReal(PetscAbsScalar(beta/betaold))/alphaold;
    258:         d[i] = PetscSqrtReal(PetscAbsScalar(beta/betaold))*e[i] + 1.0/alpha;
    259:       } else {
    260:         d[i] = PetscSqrtReal(PetscAbsScalar(beta))*e[i] + 1.0/alpha;
    261:       }
    262:       fcg->ned = ksp->its-1;
    263:     }
    264:     ++i;
    265:   } while (i<ksp->max_it);
    266:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
    267:   return(0);
    268: }

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

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

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

    297: static PetscErrorCode KSPView_FCG(KSP ksp,PetscViewer viewer)
    298: {
    299:   KSP_FCG        *fcg = (KSP_FCG*)ksp->data;
    301:   PetscBool      iascii,isstring;
    302:   const char     *truncstr;

    305:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
    306:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);

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

    312:   if (iascii) {
    313:     PetscViewerASCIIPrintf(viewer,"  m_max=%D\n",fcg->mmax);
    314:     PetscViewerASCIIPrintf(viewer,"  preallocated %D directions\n",PetscMin(fcg->nprealloc,fcg->mmax+1));
    315:     PetscViewerASCIIPrintf(viewer,"  %s\n",truncstr);
    316:   } else if (isstring) {
    317:     PetscViewerStringSPrintf(viewer,"m_max %D nprealloc %D %s",fcg->mmax,fcg->nprealloc,truncstr);
    318:   }
    319:   return(0);
    320: }

    322: /*@
    323:   KSPFCGSetMmax - set the maximum number of previous directions FCG will store for orthogonalization

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

    329:   Logically Collective on ksp

    331:   Input Parameters:
    332: +  ksp - the Krylov space context
    333: -  mmax - the maximum number of previous directions to orthogonalize againt

    335:   Level: intermediate

    337:   Options Database:
    338: . -ksp_fcg_mmax <N>

    340: .seealso: KSPFCG, KSPFCGGetTruncationType(), KSPFCGGetNprealloc()
    341: @*/
    342: PetscErrorCode KSPFCGSetMmax(KSP ksp,PetscInt mmax)
    343: {
    344:   KSP_FCG *fcg = (KSP_FCG*)ksp->data;

    349:   fcg->mmax = mmax;
    350:   return(0);
    351: }

    353: /*@
    354:   KSPFCGGetMmax - get the maximum number of previous directions FCG will store

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

    358:    Not Collective

    360:    Input Parameter:
    361: .  ksp - the Krylov space context

    363:    Output Parameter:
    364: .  mmax - the maximum number of previous directons allowed for orthogonalization

    366:   Options Database:
    367: . -ksp_fcg_mmax <N>

    369:    Level: intermediate

    371: .seealso: KSPFCG, KSPFCGGetTruncationType(), KSPFCGGetNprealloc(), KSPFCGSetMmax()
    372: @*/

    374: PetscErrorCode KSPFCGGetMmax(KSP ksp,PetscInt *mmax)
    375: {
    376:   KSP_FCG *fcg=(KSP_FCG*)ksp->data;

    380:   *mmax = fcg->mmax;
    381:   return(0);
    382: }

    384: /*@
    385:   KSPFCGSetNprealloc - set the number of directions to preallocate with FCG

    387:   Logically Collective on ksp

    389:   Input Parameters:
    390: +  ksp - the Krylov space context
    391: -  nprealloc - the number of vectors to preallocate

    393:   Level: advanced

    395:   Options Database:
    396: . -ksp_fcg_nprealloc <N> - number of directions to preallocate

    398: .seealso: KSPFCG, KSPFCGGetTruncationType(), KSPFCGGetNprealloc()
    399: @*/
    400: PetscErrorCode KSPFCGSetNprealloc(KSP ksp,PetscInt nprealloc)
    401: {
    402:   KSP_FCG *fcg=(KSP_FCG*)ksp->data;

    407:   if (nprealloc > fcg->mmax+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Cannot preallocate more than m_max+1 vectors");
    408:   fcg->nprealloc = nprealloc;
    409:   return(0);
    410: }

    412: /*@
    413:   KSPFCGGetNprealloc - get the number of directions preallocate by FCG

    415:    Not Collective

    417:    Input Parameter:
    418: .  ksp - the Krylov space context

    420:    Output Parameter:
    421: .  nprealloc - the number of directions preallocated

    423:    Level: advanced

    425: .seealso: KSPFCG, KSPFCGGetTruncationType(), KSPFCGSetNprealloc()
    426: @*/
    427: PetscErrorCode KSPFCGGetNprealloc(KSP ksp,PetscInt *nprealloc)
    428: {
    429:   KSP_FCG *fcg=(KSP_FCG*)ksp->data;

    433:   *nprealloc = fcg->nprealloc;
    434:   return(0);
    435: }

    437: /*@
    438:   KSPFCGSetTruncationType - specify how many of its stored previous directions FCG uses during orthoganalization

    440:   Logically Collective on ksp

    442:   KSP_FCD_TRUNC_TYPE_STANDARD uses all (up to mmax) stored directions
    443:   KSP_FCD_TRUNC_TYPE_NOTAY uses the last max(1,mod(i,mmax)) stored directions at iteration i=0,1,..

    445:   Input Parameters:
    446: +  ksp - the Krylov space context
    447: -  truncstrat - the choice of strategy

    449:   Level: intermediate

    451:   Options Database:
    452: . -ksp_fcg_truncation_type <standard, notay> - specify how many of its stored previous directions FCG uses during orthoganalization

    454:   .seealso: KSPFCDTruncationType, KSPFCGGetTruncationType
    455: @*/
    456: PetscErrorCode KSPFCGSetTruncationType(KSP ksp,KSPFCDTruncationType truncstrat)
    457: {
    458:   KSP_FCG *fcg=(KSP_FCG*)ksp->data;

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

    467: /*@
    468:   KSPFCGGetTruncationType - get the truncation strategy employed by FCG

    470:    Not Collective

    472:    Input Parameter:
    473: .  ksp - the Krylov space context

    475:    Output Parameter:
    476: .  truncstrat - the strategy type

    478:    Level: intermediate

    480: .seealso: KSPFCG, KSPFCGSetTruncationType, KSPFCDTruncationType
    481: @*/
    482: PetscErrorCode KSPFCGGetTruncationType(KSP ksp,KSPFCDTruncationType *truncstrat)
    483: {
    484:   KSP_FCG *fcg=(KSP_FCG*)ksp->data;

    488:   *truncstrat=fcg->truncstrat;
    489:   return(0);
    490: }

    492: static PetscErrorCode KSPSetFromOptions_FCG(PetscOptionItems *PetscOptionsObject,KSP ksp)
    493: {
    495:   KSP_FCG        *fcg=(KSP_FCG*)ksp->data;
    496:   PetscInt       mmax,nprealloc;
    497:   PetscBool      flg;

    500:   PetscOptionsHead(PetscOptionsObject,"KSP FCG Options");
    501:   PetscOptionsInt("-ksp_fcg_mmax","Maximum number of search directions to store","KSPFCGSetMmax",fcg->mmax,&mmax,&flg);
    502:   if (flg) {
    503:     KSPFCGSetMmax(ksp,mmax);
    504:   }
    505:   PetscOptionsInt("-ksp_fcg_nprealloc","Number of directions to preallocate","KSPFCGSetNprealloc",fcg->nprealloc,&nprealloc,&flg);
    506:   if (flg) {
    507:     KSPFCGSetNprealloc(ksp,nprealloc);
    508:   }
    509:   PetscOptionsEnum("-ksp_fcg_truncation_type","Truncation approach for directions","KSPFCGSetTruncationType",KSPFCDTruncationTypes,(PetscEnum)fcg->truncstrat,(PetscEnum*)&fcg->truncstrat,NULL);
    510:   PetscOptionsTail();
    511:   return(0);
    512: }

    514: /*MC
    515:       KSPFCG - Implements the Flexible Conjugate Gradient method (FCG)

    517:   Options Database Keys:
    518: +   -ksp_fcg_mmax <N>  - maximum number of search directions
    519: .   -ksp_fcg_nprealloc <N> - number of directions to preallocate
    520: -   -ksp_fcg_truncation_type <standard,notay> - truncation approach for directions

    522:     Contributed by Patrick Sanan

    524:    Notes:
    525:    Supports left preconditioning only.

    527:    Level: beginner

    529:   References:
    530: +    1. - Notay, Y."Flexible Conjugate Gradients", SIAM J. Sci. Comput. 22:4, 2000
    531: -    2. - Axelsson, O. and Vassilevski, P. S. "A Black Box Generalized Conjugate Gradient Solver with Inner Iterations and Variable step Preconditioning",
    532:     SIAM J. Matrix Anal. Appl. 12:4, 1991

    534:  .seealso : KSPGCR, KSPFGMRES, KSPCG, KSPFCGSetMmax(), KSPFCGGetMmax(), KSPFCGSetNprealloc(), KSPFCGGetNprealloc(), KSPFCGSetTruncationType(), KSPFCGGetTruncationType()

    536: M*/
    537: PETSC_EXTERN PetscErrorCode KSPCreate_FCG(KSP ksp)
    538: {
    540:   KSP_FCG        *fcg;

    543:   PetscNewLog(ksp,&fcg);
    544: #if !defined(PETSC_USE_COMPLEX)
    545:   fcg->type       = KSP_CG_SYMMETRIC;
    546: #else
    547:   fcg->type       = KSP_CG_HERMITIAN;
    548: #endif
    549:   fcg->mmax       = KSPFCG_DEFAULT_MMAX;
    550:   fcg->nprealloc  = KSPFCG_DEFAULT_NPREALLOC;
    551:   fcg->nvecs      = 0;
    552:   fcg->vecb       = KSPFCG_DEFAULT_VECB;
    553:   fcg->nchunks    = 0;
    554:   fcg->truncstrat = KSPFCG_DEFAULT_TRUNCSTRAT;

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

    558:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
    559:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,1);
    560:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);
    561:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);

    563:   ksp->ops->setup          = KSPSetUp_FCG;
    564:   ksp->ops->solve          = KSPSolve_FCG;
    565:   ksp->ops->destroy        = KSPDestroy_FCG;
    566:   ksp->ops->view           = KSPView_FCG;
    567:   ksp->ops->setfromoptions = KSPSetFromOptions_FCG;
    568:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
    569:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
    570:   return(0);
    571: }