Actual source code: fcg.c

petsc-3.6.4 2016-04-12
Report Typos and Errors
  1: /*
  2:     This file implements the FCG (Flexible Conjugate Gradient) method

  4: */

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

 10: const char *const KSPFCGTruncationTypes[]     = {"STANDARD","NOTAY","KSPFCGTrunctionTypes","KSP_FCG_TRUNC_TYPE_",0};

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

 19: static PetscErrorCode KSPAllocateVectors_FCG(KSP ksp, PetscInt nvecsneeded, PetscInt chunksize)
 20: {
 21:   PetscErrorCode  ierr;
 22:   PetscInt        i;
 23:   KSP_FCG         *fcg = (KSP_FCG*)ksp->data;
 24:   PetscInt        nnewvecs, nvecsprev;

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

 48: PetscErrorCode    KSPSetUp_FCG(KSP ksp)
 49: {
 51:   KSP_FCG        *fcg = (KSP_FCG*)ksp->data;
 52:   PetscInt       maxit = ksp->max_it;
 53:   const PetscInt nworkstd = 2;


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

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

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

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

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


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

103:   X             = ksp->vec_sol;
104:   B             = ksp->vec_rhs;
105:   R             = ksp->work[0];
106:   Z             = ksp->work[1];

108:   PCGetOperators(ksp->pc,&Amat,&Pmat);
109:   if (eigs) {e = fcg->e; d = fcg->d; e[0] = 0.0; }
110:   /* Compute initial residual needed for convergence check*/
111:   ksp->its = 0;
112:   if (!ksp->guess_zero) {
113:     KSP_MatMult(ksp,Amat,X,R);
114:     VecAYPX(R,-1.0,B);                    /*   r <- b - Ax     */
115:   } else {
116:     VecCopy(B,R);                         /*   r <- b (x is 0) */
117:   }
118:   switch (ksp->normtype) {
119:     case KSP_NORM_PRECONDITIONED:
120:       KSP_PCApply(ksp,R,Z);               /*   z <- Br         */
121:       VecNorm(Z,NORM_2,&dp);              /*   dp <- dqrt(z'*z) = sqrt(e'*A'*B'*B*A*e)     */
122:       break;
123:     case KSP_NORM_UNPRECONDITIONED:
124:       VecNorm(R,NORM_2,&dp);              /*   dp <- sqrt(r'*r) = sqrt(e'*A'*A*e)     */
125:       break;
126:     case KSP_NORM_NATURAL:
127:       KSP_PCApply(ksp,R,Z);               /*   z <- Br         */
128:       VecXDot(R,Z,&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:     /* Compute a new column of P (Currently does not support modified G-S or iterative refinement)*/
166:     switch(fcg->truncstrat){
167:       case KSP_FCG_TRUNC_TYPE_NOTAY :
168:         mi = PetscMax(1,i%(fcg->mmax+1));
169:         break;
170:       case KSP_FCG_TRUNC_TYPE_STANDARD :
171:         mi = fcg->mmax;
172:         break;
173:       default:
174:         SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unrecognized FCG Truncation Strategy");
175:     }
176:     VecCopy(Z,Pcurr);

178:     {
179:       PetscInt l,ndots;

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

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

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

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

233:     /* Check for convergence */
234:     ksp->rnorm = dp;
235:     KSPLogResidualHistory(ksp,dp);
236:     KSPMonitor(ksp,i+1,dp);
237:     (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);
238:     if (ksp->reason) break;

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

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

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

266: PetscErrorCode KSPDestroy_FCG(KSP ksp)
267: {
269:   PetscInt       i;
270:   KSP_FCG        *fcg = (KSP_FCG*)ksp->data;


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

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

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

303:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
304:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);

306:   if (fcg->truncstrat == KSP_FCG_TRUNC_TYPE_STANDARD) truncstr = "Using standard truncation strategy";
307:   else if (fcg->truncstrat == KSP_FCG_TRUNC_TYPE_NOTAY) truncstr = "Using Notay's truncation strategy";
308:   else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Undefined FCG truncation strategy");

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

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

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

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

360:    Not Collective

362:    Input Parameter:
363: .  ksp - the Krylov space context

365:    Output Parameter:
366: .  mmax - the maximum number of previous directons allowed for orthogonalization

368:   Options Database:
369: . -ksp_fcg_mmax <N>

371:    Level: intermediate

373: .keywords: KSP, FCG, truncation

375: .seealso: KSPFCG, KSPFCGGetTruncationType(), KSPFCGGetNprealloc(), KSPFCGSetMmax()
376: @*/

378: PetscErrorCode KSPFCGGetMmax(KSP ksp,PetscInt *mmax)
379: {
380:   KSP_FCG *fcg=(KSP_FCG*)ksp->data;

384:   *mmax = fcg->mmax;
385:   return(0);
386: }

390: /*@
391:   KSPFCGSetNprealloc - set the number of directions to preallocate with FCG

393:   Logically Collective on KSP

395:   Input Parameters:
396: +  ksp - the Krylov space context
397: -  nprealloc - the number of vectors to preallocate

399:   Level: advanced

401:   Options Database:
402: . -ksp_fcg_nprealloc <N>

404: .seealso: KSPFCG, KSPFCGGetTruncationType(), KSPFCGGetNprealloc()
405: @*/
406: PetscErrorCode KSPFCGSetNprealloc(KSP ksp,PetscInt nprealloc)
407: {
408:   KSP_FCG *fcg=(KSP_FCG*)ksp->data;

413:   if(nprealloc > fcg->mmax+1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Cannot preallocate more than m_max+1 vectors");
414:   fcg->nprealloc = nprealloc;
415:   return(0);
416: }

420: /*@
421:   KSPFCGGetNprealloc - get the number of directions preallocate by FCG

423:    Not Collective

425:    Input Parameter:
426: .  ksp - the Krylov space context

428:    Output Parameter:
429: .  nprealloc - the number of directions preallocated

431:   Options Database:
432: . -ksp_fcg_nprealloc <N>

434:    Level: advanced

436: .keywords: KSP, FCG, truncation

438: .seealso: KSPFCG, KSPFCGGetTruncationType(), KSPFCGSetNprealloc()
439: @*/
440: PetscErrorCode KSPFCGGetNprealloc(KSP ksp,PetscInt *nprealloc)
441: {
442:   KSP_FCG *fcg=(KSP_FCG*)ksp->data;

446:   *nprealloc = fcg->nprealloc;
447:   return(0);
448: }

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

455:   Logically Collective on KSP

457:   KSP_FCG_TRUNC_TYPE_STANDARD uses all (up to mmax) stored directions
458:   KSP_FCG_TRUNC_TYPE_NOTAY uses the last max(1,mod(i,mmax)) stored directions at iteration i=0,1..

460:   Input Parameters:
461: +  ksp - the Krylov space context
462: -  truncstrat - the choice of strategy

464:   Level: intermediate

466:   Options Database:
467: . -ksp_fcg_truncation_type <standard, notay>

469:   .seealso: KSPFCGTruncationType, KSPFCGGetTruncationType
470: @*/
471: PetscErrorCode KSPFCGSetTruncationType(KSP ksp,KSPFCGTruncationType truncstrat)
472: {
473:   KSP_FCG *fcg=(KSP_FCG*)ksp->data;

478:   fcg->truncstrat=truncstrat;
479:   return(0);
480: }

484: /*@
485:   KSPFCGGetTruncationType - get the truncation strategy employed by FCG

487:    Not Collective

489:    Input Parameter:
490: .  ksp - the Krylov space context

492:    Output Parameter:
493: .  truncstrat - the strategy type

495:   Options Database:
496: . -ksp_fcg_truncation_type <standard, notay>

498:    Level: intermediate

500: .keywords: KSP, FCG, truncation

502: .seealso: KSPFCG, KSPFCGSetTruncationType, KSPFCGTruncationType
503: @*/
504: PetscErrorCode KSPFCGGetTruncationType(KSP ksp,KSPFCGTruncationType *truncstrat)
505: {
506:   KSP_FCG *fcg=(KSP_FCG*)ksp->data;

510:   *truncstrat=fcg->truncstrat;
511:   return(0);
512: }

516: PetscErrorCode KSPSetFromOptions_FCG(PetscOptions *PetscOptionsObject,KSP ksp)
517: {
519:   KSP_FCG        *fcg=(KSP_FCG*)ksp->data;
520:   PetscInt       mmax,nprealloc;
521:   PetscBool      flg;

524:   PetscOptionsHead(PetscOptionsObject,"KSP FCG Options");
525:   PetscOptionsInt("-ksp_fcg_mmax","Maximum number of search directions to store","KSPFCGSetMmax",fcg->mmax,&mmax,&flg);
526:   if (flg) {
527:     KSPFCGSetMmax(ksp,mmax);
528:   }
529:   PetscOptionsInt("-ksp_fcg_nprealloc","Number of directions to preallocate","KSPFCGSetNprealloc",fcg->nprealloc,&nprealloc,&flg);
530:   if (flg) {
531:     KSPFCGSetNprealloc(ksp,nprealloc);
532:   }
533:   PetscOptionsEnum("-ksp_fcg_truncation_type","Truncation approach for directions","KSPFCGSetTruncationType",KSPFCGTruncationTypes,(PetscEnum)fcg->truncstrat,(PetscEnum*)&fcg->truncstrat,NULL);
534:   PetscOptionsTail();
535:   return(0);
536: }

538: /*MC
539:       KSPFCG - Implements the Flexible Conjugate Gradient method (FCG)

541:   Options Database Keys:
542: +   -ksp_fcg_mmax <N>
543: .   -ksp_fcg_nprealloc <N>
544: -   -ksp_fcg_truncation_type <standard,notay>

546:     Contributed by Patrick Sanan

548:    Notes:
549:    Supports left preconditioning only.

551:    Level: beginner

553:   References:
554:     1) Notay, Y."Flexible Conjugate Gradients", SIAM J. Sci. Comput. 22:4, pp 1444-1460, 2000

556:     2) Axelsson, O. and Vassilevski, P. S. "A Black Box Generalized Conjugate Gradient Solver with Inner Iterations and Variable-Step Preconditioning",
557:     SIAM J. Matrix Anal. Appl. 12:4, pp 625-44, 1991

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

561: M*/
564: PETSC_EXTERN PetscErrorCode KSPCreate_FCG(KSP ksp)
565: {
567:   KSP_FCG        *fcg;

570:   PetscNewLog(ksp,&fcg);
571: #if !defined(PETSC_USE_COMPLEX)
572:   fcg->type       = KSP_CG_SYMMETRIC;
573: #else
574:   fcg->type       = KSP_CG_HERMITIAN;
575: #endif
576:   fcg->mmax       = KSPFCG_DEFAULT_MMAX;
577:   fcg->nprealloc  = KSPFCG_DEFAULT_NPREALLOC;
578:   fcg->nvecs      = 0;
579:   fcg->vecb       = KSPFCG_DEFAULT_VECB;
580:   fcg->nchunks    = 0;
581:   fcg->truncstrat = KSPFCG_DEFAULT_TRUNCSTRAT;

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

585:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
586:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,1);
587:   KSPSetSupportedNorm(ksp,KSP_NORM_NATURAL,PC_LEFT,1);

589:   ksp->ops->setup          = KSPSetUp_FCG;
590:   ksp->ops->solve          = KSPSolve_FCG;
591:   ksp->ops->destroy        = KSPDestroy_FCG;
592:   ksp->ops->view           = KSPView_FCG;
593:   ksp->ops->setfromoptions = KSPSetFromOptions_FCG;
594:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
595:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
596:   return(0);
597: }