Actual source code: cheby.c

petsc-master 2020-09-18
Report Typos and Errors

  2: #include <../src/ksp/ksp/impls/cheby/chebyshevimpl.h>

  4: static PetscErrorCode KSPReset_Chebyshev(KSP ksp)
  5: {
  6:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

 10:   if (cheb->kspest) {
 11:     KSPReset(cheb->kspest);
 12:   }
 13:   return(0);
 14: }

 16: PETSC_STATIC_INLINE PetscScalar chebyhash(PetscInt xx)
 17: {
 18:   unsigned int x = xx;
 19:   x = ((x >> 16) ^ x) * 0x45d9f3b;
 20:   x = ((x >> 16) ^ x) * 0x45d9f3b;
 21:   x = ((x >> 16) ^ x);
 22:   return (PetscScalar)((PetscInt64)x-2147483648)*5.e-10; /* center around zero, scaled about -1. to 1.*/
 23: }

 25: /*
 26:  * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
 27:  */
 28: static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest,PetscReal *emin,PetscReal *emax)
 29: {
 31:   PetscInt       n,neig;
 32:   PetscReal      *re,*im,min,max;

 35:   KSPGetIterationNumber(kspest,&n);
 36:   PetscMalloc2(n,&re,n,&im);
 37:   KSPComputeEigenvalues(kspest,n,re,im,&neig);
 38:   min  = PETSC_MAX_REAL;
 39:   max  = PETSC_MIN_REAL;
 40:   for (n=0; n<neig; n++) {
 41:     min = PetscMin(min,re[n]);
 42:     max = PetscMax(max,re[n]);
 43:   }
 44:   PetscFree2(re,im);
 45:   *emax = max;
 46:   *emin = min;
 47:   return(0);
 48: }

 50: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
 51: {
 52:   KSP_Chebyshev    *cheb = (KSP_Chebyshev*)ksp->data;
 53:   PetscErrorCode   ierr;
 54:   PetscBool        flg;
 55:   Mat              Pmat,Amat;
 56:   PetscObjectId    amatid,    pmatid;
 57:   PetscObjectState amatstate, pmatstate;
 59:   KSPSetWorkVecs(ksp,3);
 60:   if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) { /* We need to estimate eigenvalues */
 61:     KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
 62:   }
 63:   if (cheb->kspest) {
 64:     KSPGetOperators(ksp,&Amat,&Pmat);
 65:     MatGetOption(Pmat, MAT_SPD, &flg);
 66:     if (flg) {
 67:       const char *prefix;
 68:       KSPGetOptionsPrefix(cheb->kspest,&prefix);
 69:       PetscOptionsHasName(NULL,prefix,"-ksp_type",&flg);
 70:       if (!flg) {
 71:         KSPSetType(cheb->kspest, KSPCG);
 72:       }
 73:     }
 74:     PetscObjectGetId((PetscObject)Amat,&amatid);
 75:     PetscObjectGetId((PetscObject)Pmat,&pmatid);
 76:     PetscObjectStateGet((PetscObject)Amat,&amatstate);
 77:     PetscObjectStateGet((PetscObject)Pmat,&pmatstate);
 78:     if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
 79:       PetscReal          max=0.0,min=0.0;
 80:       Vec                B;
 81:       KSPConvergedReason reason;
 82:       KSPSetPC(cheb->kspest,ksp->pc);
 83:       if (cheb->usenoisy) {
 84:         PetscInt       n,i,istart;
 85:         PetscScalar    *xx;

 87:         B    = ksp->work[1];
 88:         VecGetOwnershipRange(B,&istart,NULL);
 89:         VecGetLocalSize(B,&n);
 90:         VecGetArrayWrite(B,&xx);
 91:         for (i=0; i<n; i++) xx[i] = chebyhash(i+istart);
 92:         VecRestoreArrayWrite(B,&xx);
 93:       } else {
 94:         PetscBool change;

 96:         PCPreSolveChangeRHS(ksp->pc,&change);
 97:         if (change) {
 98:           B = ksp->work[1];
 99:           VecCopy(ksp->vec_rhs,B);
100:         } else B = ksp->vec_rhs;
101:       }
102:       KSPSolve(cheb->kspest,B,ksp->work[0]);
103:       KSPGetConvergedReason(cheb->kspest,&reason);
104:       if (reason == KSP_DIVERGED_ITS) {
105:         PetscInfo(ksp,"Eigen estimator ran for prescribed number of iterations\n");
106:       } else if (reason == KSP_DIVERGED_PC_FAILED) {
107:         PetscInt       its;
108:         PCFailedReason pcreason;

110:         KSPGetIterationNumber(cheb->kspest,&its);
111:         if (ksp->normtype == KSP_NORM_NONE) {
112:           PetscInt  sendbuf,recvbuf;
113:           PCGetFailedReasonRank(ksp->pc,&pcreason);
114:           sendbuf = (PetscInt)pcreason;
115:           MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));
116:           PCSetFailedReason(ksp->pc,(PCFailedReason) recvbuf);
117:         }
118:         PCGetFailedReason(ksp->pc,&pcreason);
119:         ksp->reason = KSP_DIVERGED_PC_FAILED;
120:         VecSetInf(ksp->vec_sol);
121:         PetscInfo3(ksp,"Eigen estimator failed: %s %s at iteration %D",KSPConvergedReasons[reason],PCFailedReasons[pcreason],its);
122:         return(0);
123:       } else if (reason == KSP_CONVERGED_RTOL || reason == KSP_CONVERGED_ATOL) {
124:         PetscInfo(ksp,"Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n");
125:       } else if (reason < 0) {
126:         PetscInfo1(ksp,"Eigen estimator failed %s, using estimates anyway\n",KSPConvergedReasons[reason]);
127:       }

129:       KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest,&min,&max);
130:       KSPSetPC(cheb->kspest,NULL);

132:       cheb->emin_computed = min;
133:       cheb->emax_computed = max;
134:       cheb->emin = cheb->tform[0]*min + cheb->tform[1]*max;
135:       cheb->emax = cheb->tform[2]*min + cheb->tform[3]*max;

137:       cheb->amatid    = amatid;
138:       cheb->pmatid    = pmatid;
139:       cheb->amatstate = amatstate;
140:       cheb->pmatstate = pmatstate;
141:     }
142:   }
143:   return(0);
144: }

146: static PetscErrorCode KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp,PetscReal emax,PetscReal emin)
147: {
148:   KSP_Chebyshev  *chebyshevP = (KSP_Chebyshev*)ksp->data;

152:   if (emax <= emin) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Maximum eigenvalue must be larger than minimum: max %g min %g",(double)emax,(double)emin);
153:   if (emax*emin <= 0.0) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Both eigenvalues must be of the same sign: max %g min %g",(double)emax,(double)emin);
154:   chebyshevP->emax = emax;
155:   chebyshevP->emin = emin;

157:   KSPChebyshevEstEigSet(ksp,0.,0.,0.,0.); /* Destroy any estimation setup */
158:   return(0);
159: }

161: static PetscErrorCode KSPChebyshevEstEigSet_Chebyshev(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
162: {
163:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

167:   if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
168:     if (!cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
169:       KSPCreate(PetscObjectComm((PetscObject)ksp),&cheb->kspest);
170:       KSPSetErrorIfNotConverged(cheb->kspest,ksp->errorifnotconverged);
171:       PetscObjectIncrementTabLevel((PetscObject)cheb->kspest,(PetscObject)ksp,1);
172:       /* use PetscObjectSet/AppendOptionsPrefix() instead of KSPSet/AppendOptionsPrefix() so that the PC prefix is not changed */
173:       PetscObjectSetOptionsPrefix((PetscObject)cheb->kspest,((PetscObject)ksp)->prefix);
174:       PetscObjectAppendOptionsPrefix((PetscObject)cheb->kspest,"esteig_");
175:       KSPSetSkipPCSetFromOptions(cheb->kspest,PETSC_TRUE);

177:       KSPSetComputeEigenvalues(cheb->kspest,PETSC_TRUE);

179:       /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
180:       KSPSetTolerances(cheb->kspest,1.e-12,PETSC_DEFAULT,PETSC_DEFAULT,cheb->eststeps);
181:     }
182:     if (a >= 0) cheb->tform[0] = a;
183:     if (b >= 0) cheb->tform[1] = b;
184:     if (c >= 0) cheb->tform[2] = c;
185:     if (d >= 0) cheb->tform[3] = d;
186:     cheb->amatid    = 0;
187:     cheb->pmatid    = 0;
188:     cheb->amatstate = -1;
189:     cheb->pmatstate = -1;
190:   } else {
191:     KSPDestroy(&cheb->kspest);
192:   }
193:   return(0);
194: }

196: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp,PetscBool use)
197: {
198:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

201:   cheb->usenoisy = use;
202:   return(0);
203: }

205: /*@
206:    KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues
207:    of the preconditioned problem.

209:    Logically Collective on ksp

211:    Input Parameters:
212: +  ksp - the Krylov space context
213: -  emax, emin - the eigenvalue estimates

215:   Options Database:
216: .  -ksp_chebyshev_eigenvalues emin,emax

218:    Note: Call KSPChebyshevEstEigSet() or use the option -ksp_chebyshev_esteig a,b,c,d to have the KSP
219:          estimate the eigenvalues and use these estimated values automatically

221:    Level: intermediate

223: @*/
224: PetscErrorCode  KSPChebyshevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
225: {

232:   PetscTryMethod(ksp,"KSPChebyshevSetEigenvalues_C",(KSP,PetscReal,PetscReal),(ksp,emax,emin));
233:   return(0);
234: }

236: /*@
237:    KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev

239:    Logically Collective on ksp

241:    Input Parameters:
242: +  ksp - the Krylov space context
243: .  a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
244: .  b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
245: .  c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
246: -  d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)

248:   Options Database:
249: .  -ksp_chebyshev_esteig a,b,c,d

251:    Notes:
252:    The Chebyshev bounds are set using
253: .vb
254:    minbound = a*minest + b*maxest
255:    maxbound = c*minest + d*maxest
256: .ve
257:    The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
258:    The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.

260:    If 0.0 is passed for all transform arguments (a,b,c,d), eigenvalue estimation is disabled.

262:    The default transform is (0,0.1; 0,1.1) which targets the "upper" part of the spectrum, as desirable for use with multigrid.

264:    The eigenvalues are estimated using the Lanczo (KSPCG) or Arnoldi (KSPGMRES) process using a noisy right hand side vector.

266:    Level: intermediate

268: @*/
269: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
270: {

279:   PetscTryMethod(ksp,"KSPChebyshevEstEigSet_C",(KSP,PetscReal,PetscReal,PetscReal,PetscReal),(ksp,a,b,c,d));
280:   return(0);
281: }

283: /*@
284:    KSPChebyshevEstEigSetUseNoisy - use a noisy right hand side in order to do the estimate instead of the given right hand side

286:    Logically Collective

288:    Input Arguments:
289: +  ksp - linear solver context
290: -  use - PETSC_TRUE to use noisy

292:    Options Database:
293: .  -ksp_chebyshev_esteig_noisy <true,false>

295:   Notes:
296:     This alledgely works better for multigrid smoothers

298:   Level: intermediate

300: .seealso: KSPChebyshevEstEigSet()
301: @*/
302: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp,PetscBool use)
303: {

308:   PetscTryMethod(ksp,"KSPChebyshevEstEigSetUseNoisy_C",(KSP,PetscBool),(ksp,use));
309:   return(0);
310: }

312: /*@
313:   KSPChebyshevEstEigGetKSP - Get the Krylov method context used to estimate eigenvalues for the Chebyshev method.  If
314:   a Krylov method is not being used for this purpose, NULL is returned.  The reference count of the returned KSP is
315:   not incremented: it should not be destroyed by the user.

317:   Input Parameters:
318: . ksp - the Krylov space context

320:   Output Parameters:
321: . kspest the eigenvalue estimation Krylov space context

323:   Level: intermediate

325: .seealso: KSPChebyshevEstEigSet()
326: @*/
327: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
328: {

334:   *kspest = NULL;
335:   PetscTryMethod(ksp,"KSPChebyshevEstEigGetKSP_C",(KSP,KSP*),(ksp,kspest));
336:   return(0);
337: }

339: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
340: {
341:   KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;

344:   *kspest = cheb->kspest;
345:   return(0);
346: }

348: static PetscErrorCode KSPSetFromOptions_Chebyshev(PetscOptionItems *PetscOptionsObject,KSP ksp)
349: {
350:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;
352:   PetscInt       neigarg = 2, nestarg = 4;
353:   PetscReal      eminmax[2] = {0., 0.};
354:   PetscReal      tform[4] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
355:   PetscBool      flgeig, flgest;

358:   PetscOptionsHead(PetscOptionsObject,"KSP Chebyshev Options");
359:   PetscOptionsInt("-ksp_chebyshev_esteig_steps","Number of est steps in Chebyshev","",cheb->eststeps,&cheb->eststeps,NULL);
360:   PetscOptionsRealArray("-ksp_chebyshev_eigenvalues","extreme eigenvalues","KSPChebyshevSetEigenvalues",eminmax,&neigarg,&flgeig);
361:   if (flgeig) {
362:     if (neigarg != 2) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
363:     KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]);
364:   }
365:   PetscOptionsRealArray("-ksp_chebyshev_esteig","estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds","KSPChebyshevEstEigSet",tform,&nestarg,&flgest);
366:   if (flgest) {
367:     switch (nestarg) {
368:     case 0:
369:       KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
370:       break;
371:     case 2:                     /* Base everything on the max eigenvalues */
372:       KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,tform[0],PETSC_DECIDE,tform[1]);
373:       break;
374:     case 4:                     /* Use the full 2x2 linear transformation */
375:       KSPChebyshevEstEigSet(ksp,tform[0],tform[1],tform[2],tform[3]);
376:       break;
377:     default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
378:     }
379:   }

381:   /* We need to estimate eigenvalues; need to set this here so that KSPSetFromOptions() is called on the estimator */
382:   if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) {
383:    KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
384:   }

386:   if (cheb->kspest) {
387:     PetscOptionsBool("-ksp_chebyshev_esteig_noisy","Use noisy right hand side for estimate","KSPChebyshevEstEigSetUseNoisy",cheb->usenoisy,&cheb->usenoisy,NULL);
388:     KSPSetFromOptions(cheb->kspest);
389:   }
390:   PetscOptionsTail();
391:   return(0);
392: }

394: static PetscErrorCode KSPSolve_Chebyshev(KSP ksp)
395: {
396:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;
398:   PetscInt       k,kp1,km1,ktmp,i;
399:   PetscScalar    alpha,omegaprod,mu,omega,Gamma,c[3],scale;
400:   PetscReal      rnorm = 0.0;
401:   Vec            sol_orig,b,p[3],r;
402:   Mat            Amat,Pmat;
403:   PetscBool      diagonalscale;

406:   PCGetDiagonalScale(ksp->pc,&diagonalscale);
407:   if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);

409:   PCGetOperators(ksp->pc,&Amat,&Pmat);
410:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
411:   ksp->its = 0;
412:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
413:   /* These three point to the three active solutions, we
414:      rotate these three at each solution update */
415:   km1      = 0; k = 1; kp1 = 2;
416:   sol_orig = ksp->vec_sol; /* ksp->vec_sol will be asigned to rotating vector p[k], thus save its address */
417:   b        = ksp->vec_rhs;
418:   p[km1]   = sol_orig;
419:   p[k]     = ksp->work[0];
420:   p[kp1]   = ksp->work[1];
421:   r        = ksp->work[2];

423:   /* use scale*B as our preconditioner */
424:   scale = 2.0/(cheb->emax + cheb->emin);

426:   /*   -alpha <=  scale*lambda(B^{-1}A) <= alpha   */
427:   alpha     = 1.0 - scale*(cheb->emin);
428:   Gamma     = 1.0;
429:   mu        = 1.0/alpha;
430:   omegaprod = 2.0/alpha;

432:   c[km1] = 1.0;
433:   c[k]   = mu;

435:   if (!ksp->guess_zero) {
436:     KSP_MatMult(ksp,Amat,sol_orig,r);     /*  r = b - A*p[km1] */
437:     VecAYPX(r,-1.0,b);
438:   } else {
439:     VecCopy(b,r);
440:   }

442:   /* calculate residual norm if requested, we have done one iteration */
443:   if (ksp->normtype) {
444:     switch (ksp->normtype) {
445:     case KSP_NORM_PRECONDITIONED:
446:       KSP_PCApply(ksp,r,p[k]);  /* p[k] = B^{-1}r */
447:       VecNorm(p[k],NORM_2,&rnorm);
448:       break;
449:     case KSP_NORM_UNPRECONDITIONED:
450:     case KSP_NORM_NATURAL:
451:       VecNorm(r,NORM_2,&rnorm);
452:       break;
453:     default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
454:     }
455:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
456:     ksp->rnorm   = rnorm;
457:     PetscObjectSAWsGrantAccess((PetscObject)ksp);
458:     KSPLogResidualHistory(ksp,rnorm);
459:     KSPMonitor(ksp,0,rnorm);
460:     (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
461:   } else ksp->reason = KSP_CONVERGED_ITERATING;
462:   if (ksp->reason || ksp->max_it==0) {
463:     if (ksp->max_it==0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
464:     return(0);
465:   }
466:   if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
467:     KSP_PCApply(ksp,r,p[k]);  /* p[k] = B^{-1}r */
468:   }
469:   VecAYPX(p[k],scale,p[km1]);  /* p[k] = scale B^{-1}r + p[km1] */
470:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
471:   ksp->its = 1;
472:   PetscObjectSAWsGrantAccess((PetscObject)ksp);

474:   for (i=1; i<ksp->max_it; i++) {
475:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
476:     ksp->its++;
477:     PetscObjectSAWsGrantAccess((PetscObject)ksp);

479:     KSP_MatMult(ksp,Amat,p[k],r);          /*  r = b - Ap[k]    */
480:     VecAYPX(r,-1.0,b);
481:     /* calculate residual norm if requested */
482:     if (ksp->normtype) {
483:       switch (ksp->normtype) {
484:       case KSP_NORM_PRECONDITIONED:
485:         KSP_PCApply(ksp,r,p[kp1]);             /*  p[kp1] = B^{-1}r  */
486:         VecNorm(p[kp1],NORM_2,&rnorm);
487:         break;
488:       case KSP_NORM_UNPRECONDITIONED:
489:       case KSP_NORM_NATURAL:
490:         VecNorm(r,NORM_2,&rnorm);
491:         break;
492:       default:
493:         rnorm = 0.0;
494:         break;
495:       }
496:       KSPCheckNorm(ksp,rnorm);
497:       PetscObjectSAWsTakeAccess((PetscObject)ksp);
498:       ksp->rnorm   = rnorm;
499:       PetscObjectSAWsGrantAccess((PetscObject)ksp);
500:       KSPLogResidualHistory(ksp,rnorm);
501:       KSPMonitor(ksp,i,rnorm);
502:       (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
503:       if (ksp->reason) break;
504:       if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
505:         KSP_PCApply(ksp,r,p[kp1]);             /*  p[kp1] = B^{-1}r  */
506:       }
507:     } else {
508:       KSP_PCApply(ksp,r,p[kp1]);             /*  p[kp1] = B^{-1}r  */
509:     }
510:     ksp->vec_sol = p[k];

512:     c[kp1] = 2.0*mu*c[k] - c[km1];
513:     omega  = omegaprod*c[k]/c[kp1];

515:     /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
516:     VecAXPBYPCZ(p[kp1],1.0-omega,omega,omega*Gamma*scale,p[km1],p[k]);

518:     ktmp = km1;
519:     km1  = k;
520:     k    = kp1;
521:     kp1  = ktmp;
522:   }
523:   if (!ksp->reason) {
524:     if (ksp->normtype) {
525:       KSP_MatMult(ksp,Amat,p[k],r);       /*  r = b - Ap[k]    */
526:       VecAYPX(r,-1.0,b);
527:       switch (ksp->normtype) {
528:       case KSP_NORM_PRECONDITIONED:
529:         KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
530:         VecNorm(p[kp1],NORM_2,&rnorm);
531:         break;
532:       case KSP_NORM_UNPRECONDITIONED:
533:       case KSP_NORM_NATURAL:
534:         VecNorm(r,NORM_2,&rnorm);
535:         break;
536:       default:
537:         rnorm = 0.0;
538:         break;
539:       }
540:       KSPCheckNorm(ksp,rnorm);
541:       PetscObjectSAWsTakeAccess((PetscObject)ksp);
542:       ksp->rnorm   = rnorm;
543:       PetscObjectSAWsGrantAccess((PetscObject)ksp);
544:       KSPLogResidualHistory(ksp,rnorm);
545:       KSPMonitor(ksp,i,rnorm);
546:     }
547:     if (ksp->its >= ksp->max_it) {
548:       if (ksp->normtype != KSP_NORM_NONE) {
549:         (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
550:         if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
551:       } else ksp->reason = KSP_CONVERGED_ITS;
552:     }
553:   }

555:   /* make sure solution is in vector x */
556:   ksp->vec_sol = sol_orig;
557:   if (k) {
558:     VecCopy(p[k],sol_orig);
559:   }
560:   return(0);
561: }

563: static  PetscErrorCode KSPView_Chebyshev(KSP ksp,PetscViewer viewer)
564: {
565:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;
567:   PetscBool      iascii;

570:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
571:   if (iascii) {
572:     PetscViewerASCIIPrintf(viewer,"  eigenvalue estimates used:  min = %g, max = %g\n",(double)cheb->emin,(double)cheb->emax);
573:     if (cheb->kspest) {
574:       PetscViewerASCIIPrintf(viewer,"  eigenvalues estimate via %s min %g, max %g\n",((PetscObject)(cheb->kspest))->type_name,(double)cheb->emin_computed,(double)cheb->emax_computed);
575:       PetscViewerASCIIPrintf(viewer,"  eigenvalues estimated using %s with translations  [%g %g; %g %g]\n",((PetscObject) cheb->kspest)->type_name,(double)cheb->tform[0],(double)cheb->tform[1],(double)cheb->tform[2],(double)cheb->tform[3]);
576:       PetscViewerASCIIPushTab(viewer);
577:       KSPView(cheb->kspest,viewer);
578:       PetscViewerASCIIPopTab(viewer);
579:       if (cheb->usenoisy) {
580:         PetscViewerASCIIPrintf(viewer,"  estimating eigenvalues using noisy right hand side\n");
581:       }
582:     }
583:   }
584:   return(0);
585: }

587: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
588: {
589:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

593:   KSPDestroy(&cheb->kspest);
594:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",NULL);
595:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",NULL);
596:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",NULL);
597:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",NULL);
598:   KSPDestroyDefault(ksp);
599:   return(0);
600: }

602: /*MC
603:      KSPCHEBYSHEV - The preconditioned Chebyshev iterative method

605:    Options Database Keys:
606: +   -ksp_chebyshev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
607:                   of the preconditioned operator. If these are accurate you will get much faster convergence.
608: .   -ksp_chebyshev_esteig <a,b,c,d> - estimate eigenvalues using a Krylov method, then use this
609:                          transform for Chebyshev eigenvalue bounds (KSPChebyshevEstEigSet())
610: .   -ksp_chebyshev_esteig_steps - number of estimation steps
611: -   -ksp_chebyshev_esteig_noisy - use noisy number generator to create right hand side for eigenvalue estimator

613:    Level: beginner

615:    Notes:
616:     The Chebyshev method requires both the matrix and preconditioner to
617:           be symmetric positive (semi) definite.
618:           Only support for left preconditioning.

620:           Chebyshev is configured as a smoother by default, targetting the "upper" part of the spectrum.
621:           The user should call KSPChebyshevSetEigenvalues() if they have eigenvalue estimates.

623: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
624:            KSPChebyshevSetEigenvalues(), KSPChebyshevEstEigSet(), KSPChebyshevEstEigSetUseNoisy()
625:            KSPRICHARDSON, KSPCG, PCMG

627: M*/

629: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
630: {
632:   KSP_Chebyshev  *chebyshevP;

635:   PetscNewLog(ksp,&chebyshevP);

637:   ksp->data = (void*)chebyshevP;
638:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
639:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
640:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
641:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);

643:   chebyshevP->emin = 0.;
644:   chebyshevP->emax = 0.;

646:   chebyshevP->tform[0] = 0.0;
647:   chebyshevP->tform[1] = 0.1;
648:   chebyshevP->tform[2] = 0;
649:   chebyshevP->tform[3] = 1.1;
650:   chebyshevP->eststeps = 10;
651:   chebyshevP->usenoisy = PETSC_TRUE;
652:   ksp->setupnewmatrix = PETSC_TRUE;

654:   ksp->ops->setup          = KSPSetUp_Chebyshev;
655:   ksp->ops->solve          = KSPSolve_Chebyshev;
656:   ksp->ops->destroy        = KSPDestroy_Chebyshev;
657:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
658:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
659:   ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
660:   ksp->ops->view           = KSPView_Chebyshev;
661:   ksp->ops->reset          = KSPReset_Chebyshev;

663:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",KSPChebyshevSetEigenvalues_Chebyshev);
664:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",KSPChebyshevEstEigSet_Chebyshev);
665:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",KSPChebyshevEstEigSetUseNoisy_Chebyshev);
666:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",KSPChebyshevEstEigGetKSP_Chebyshev);
667:   return(0);
668: }