Actual source code: cheby.c


  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;

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

 14: /*
 15:  * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
 16:  */
 17: static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest,PetscReal *emin,PetscReal *emax)
 18: {
 19:   PetscInt       n,neig;
 20:   PetscReal      *re,*im,min,max;

 22:   KSPGetIterationNumber(kspest,&n);
 23:   PetscMalloc2(n,&re,n,&im);
 24:   KSPComputeEigenvalues(kspest,n,re,im,&neig);
 25:   min  = PETSC_MAX_REAL;
 26:   max  = PETSC_MIN_REAL;
 27:   for (n=0; n<neig; n++) {
 28:     min = PetscMin(min,re[n]);
 29:     max = PetscMax(max,re[n]);
 30:   }
 31:   PetscFree2(re,im);
 32:   *emax = max;
 33:   *emin = min;
 34:   return 0;
 35: }

 37: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
 38: {
 39:   KSP_Chebyshev    *cheb = (KSP_Chebyshev*)ksp->data;
 40:   PetscBool        flg;
 41:   Mat              Pmat,Amat;
 42:   PetscObjectId    amatid,    pmatid;
 43:   PetscObjectState amatstate, pmatstate;
 44:   KSPSetWorkVecs(ksp,3);
 45:   if (cheb->emin == 0. || cheb->emax == 0.) { // User did not specify eigenvalues
 46:     PC pc;
 47:     KSPGetPC(ksp,&pc);
 48:     PetscObjectTypeCompare((PetscObject)pc,PCJACOBI,&flg);
 49:     if (!flg) { // Provided estimates are only relevant for Jacobi
 50:       cheb->emax_provided = 0;
 51:       cheb->emin_provided = 0;
 52:     }
 53:     if (!cheb->kspest) { /* We need to estimate eigenvalues */
 54:       KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
 55:     }
 56:   }
 57:   if (cheb->kspest) {
 58:     KSPGetOperators(ksp,&Amat,&Pmat);
 59:     MatGetOption(Pmat, MAT_SPD, &flg);
 60:     if (flg) {
 61:       const char *prefix;
 62:       KSPGetOptionsPrefix(cheb->kspest,&prefix);
 63:       PetscOptionsHasName(NULL,prefix,"-ksp_type",&flg);
 64:       if (!flg) {
 65:         KSPSetType(cheb->kspest, KSPCG);
 66:       }
 67:     }
 68:     PetscObjectGetId((PetscObject)Amat,&amatid);
 69:     PetscObjectGetId((PetscObject)Pmat,&pmatid);
 70:     PetscObjectStateGet((PetscObject)Amat,&amatstate);
 71:     PetscObjectStateGet((PetscObject)Pmat,&pmatstate);
 72:     if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
 73:       PetscReal          max=0.0,min=0.0;
 74:       Vec                B;
 75:       KSPConvergedReason reason;
 76:       KSPSetPC(cheb->kspest,ksp->pc);
 77:       if (cheb->usenoisy) {
 78:         B = ksp->work[1];
 79:         KSPSetNoisy_Private(B);
 80:       } else {
 81:         PetscBool change;

 84:         PCPreSolveChangeRHS(ksp->pc,&change);
 85:         if (change) {
 86:           B = ksp->work[1];
 87:           VecCopy(ksp->vec_rhs,B);
 88:         } else B = ksp->vec_rhs;
 89:       }
 90:       KSPSolve(cheb->kspest,B,ksp->work[0]);
 91:       KSPGetConvergedReason(cheb->kspest,&reason);
 92:       if (reason == KSP_DIVERGED_ITS) {
 93:         PetscInfo(ksp,"Eigen estimator ran for prescribed number of iterations\n");
 94:       } else if (reason == KSP_DIVERGED_PC_FAILED) {
 95:         PetscInt       its;
 96:         PCFailedReason pcreason;

 98:         KSPGetIterationNumber(cheb->kspest,&its);
 99:         if (ksp->normtype == KSP_NORM_NONE) {
100:           PetscInt  sendbuf,recvbuf;
101:           PCGetFailedReasonRank(ksp->pc,&pcreason);
102:           sendbuf = (PetscInt)pcreason;
103:           MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));
104:           PCSetFailedReason(ksp->pc,(PCFailedReason) recvbuf);
105:         }
106:         PCGetFailedReason(ksp->pc,&pcreason);
107:         ksp->reason = KSP_DIVERGED_PC_FAILED;
108:         PetscInfo(ksp,"Eigen estimator failed: %s %s at iteration %D",KSPConvergedReasons[reason],PCFailedReasons[pcreason],its);
109:         return 0;
110:       } else if (reason == KSP_CONVERGED_RTOL || reason == KSP_CONVERGED_ATOL) {
111:         PetscInfo(ksp,"Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n");
112:       } else if (reason < 0) {
113:         PetscInfo(ksp,"Eigen estimator failed %s, using estimates anyway\n",KSPConvergedReasons[reason]);
114:       }

116:       KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest,&min,&max);
117:       KSPSetPC(cheb->kspest,NULL);

119:       cheb->emin_computed = min;
120:       cheb->emax_computed = max;

122:       cheb->amatid    = amatid;
123:       cheb->pmatid    = pmatid;
124:       cheb->amatstate = amatstate;
125:       cheb->pmatstate = pmatstate;
126:     }
127:   }
128:   return 0;
129: }

131: static PetscErrorCode KSPChebyshevGetEigenvalues_Chebyshev(KSP ksp,PetscReal *emax,PetscReal *emin)
132: {
133:   KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;

135:   *emax = 0;
136:   *emin = 0;
137:   if (cheb->emax != 0.) {
138:     *emax = cheb->emax;
139:   } else if (cheb->emax_computed != 0.) {
140:     *emax = cheb->tform[2] * cheb->emin_computed + cheb->tform[3] * cheb->emax_computed;
141:   } else if (cheb->emax_provided != 0.) {
142:     *emax = cheb->tform[2] * cheb->emin_provided + cheb->tform[3] * cheb->emax_provided;
143:   }
144:   if (cheb->emin != 0.) {
145:     *emin = cheb->emin;
146:   } else if (cheb->emin_computed != 0.) {
147:     *emin = cheb->tform[0] * cheb->emin_computed + cheb->tform[1] * cheb->emax_computed;
148:   } else if (cheb->emin_provided != 0.) {
149:     *emin = cheb->tform[0] * cheb->emin_provided + cheb->tform[1] * cheb->emax_provided;
150:   }
151:   return 0;
152: }

154: static PetscErrorCode KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp,PetscReal emax,PetscReal emin)
155: {
156:   KSP_Chebyshev  *chebyshevP = (KSP_Chebyshev*)ksp->data;

160:   chebyshevP->emax = emax;
161:   chebyshevP->emin = emin;

163:   KSPChebyshevEstEigSet(ksp,0.,0.,0.,0.); /* Destroy any estimation setup */
164:   return 0;
165: }

167: static PetscErrorCode KSPChebyshevEstEigSet_Chebyshev(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
168: {
169:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

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

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

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

200: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp,PetscBool use)
201: {
202:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

204:   cheb->usenoisy = use;
205:   return 0;
206: }

208: /*@
209:    KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues
210:    of the preconditioned problem.

212:    Logically Collective on ksp

214:    Input Parameters:
215: +  ksp - the Krylov space context
216: -  emax, emin - the eigenvalue estimates

218:   Options Database:
219: .  -ksp_chebyshev_eigenvalues emin,emax - extreme eigenvalues

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

225:    When KSPCHEBYSHEV is used as a smoother, one often wants to target a portion of the spectrum rather than the entire
226:    spectrum. This function takes the range of target eigenvalues for Chebyshev, which will often slightly over-estimate
227:    the largest eigenvalue of the actual operator (for safety) and greatly overestimate the smallest eigenvalue to
228:    improve the smoothing properties of Chebyshev iteration on the higher frequencies in the spectrum.

230:    Level: intermediate

232: .seealso: KSPChebyshevEstEigSet()
233: @*/
234: PetscErrorCode  KSPChebyshevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
235: {
239:   PetscTryMethod(ksp,"KSPChebyshevSetEigenvalues_C",(KSP,PetscReal,PetscReal),(ksp,emax,emin));
240:   return 0;
241: }

243: /*@
244:    KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev

246:    Logically Collective on ksp

248:    Input Parameters:
249: +  ksp - the Krylov space context
250: .  a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
251: .  b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
252: .  c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
253: -  d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)

255:   Options Database:
256: .  -ksp_chebyshev_esteig a,b,c,d - estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds

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

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

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

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

273:    Level: intermediate

275: @*/
276: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
277: {
283:   PetscTryMethod(ksp,"KSPChebyshevEstEigSet_C",(KSP,PetscReal,PetscReal,PetscReal,PetscReal),(ksp,a,b,c,d));
284:   return 0;
285: }

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

290:    Logically Collective

292:    Input Parameters:
293: +  ksp - linear solver context
294: -  use - PETSC_TRUE to use noisy

296:    Options Database:
297: .  -ksp_chebyshev_esteig_noisy <true,false> - Use noisy right hand side for estimate

299:   Notes:
300:     This alledgely works better for multigrid smoothers

302:   Level: intermediate

304: .seealso: KSPChebyshevEstEigSet()
305: @*/
306: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp,PetscBool use)
307: {
309:   PetscTryMethod(ksp,"KSPChebyshevEstEigSetUseNoisy_C",(KSP,PetscBool),(ksp,use));
310:   return 0;
311: }

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

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

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

324:   Level: intermediate

326: .seealso: KSPChebyshevEstEigSet()
327: @*/
328: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
329: {
332:   *kspest = NULL;
333:   PetscTryMethod(ksp,"KSPChebyshevEstEigGetKSP_C",(KSP,KSP*),(ksp,kspest));
334:   return 0;
335: }

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

341:   *kspest = cheb->kspest;
342:   return 0;
343: }

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

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

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

381:   if (cheb->kspest) {
382:     PetscOptionsBool("-ksp_chebyshev_esteig_noisy","Use noisy right hand side for estimate","KSPChebyshevEstEigSetUseNoisy",cheb->usenoisy,&cheb->usenoisy,NULL);
383:     KSPSetFromOptions(cheb->kspest);
384:   }
385:   PetscOptionsTail();
386:   return 0;
387: }

389: static PetscErrorCode KSPSolve_Chebyshev(KSP ksp)
390: {
391:   PetscInt       k,kp1,km1,ktmp,i;
392:   PetscScalar    alpha,omegaprod,mu,omega,Gamma,c[3],scale;
393:   PetscReal      rnorm = 0.0,emax,emin;
394:   Vec            sol_orig,b,p[3],r;
395:   Mat            Amat,Pmat;
396:   PetscBool      diagonalscale;

398:   PCGetDiagonalScale(ksp->pc,&diagonalscale);

401:   PCGetOperators(ksp->pc,&Amat,&Pmat);
402:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
403:   ksp->its = 0;
404:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
405:   /* These three point to the three active solutions, we
406:      rotate these three at each solution update */
407:   km1      = 0; k = 1; kp1 = 2;
408:   sol_orig = ksp->vec_sol; /* ksp->vec_sol will be assigned to rotating vector p[k], thus save its address */
409:   b        = ksp->vec_rhs;
410:   p[km1]   = sol_orig;
411:   p[k]     = ksp->work[0];
412:   p[kp1]   = ksp->work[1];
413:   r        = ksp->work[2];

415:   KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin);
416:   /* use scale*B as our preconditioner */
417:   scale = 2.0/(emax + emin);

419:   /*   -alpha <=  scale*lambda(B^{-1}A) <= alpha   */
420:   alpha     = 1.0 - scale * emin;
421:   Gamma     = 1.0;
422:   mu        = 1.0/alpha;
423:   omegaprod = 2.0/alpha;

425:   c[km1] = 1.0;
426:   c[k]   = mu;

428:   if (!ksp->guess_zero) {
429:     KSP_MatMult(ksp,Amat,sol_orig,r);     /*  r = b - A*p[km1] */
430:     VecAYPX(r,-1.0,b);
431:   } else {
432:     VecCopy(b,r);
433:   }

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

468:   for (i=1; i<ksp->max_it; i++) {
469:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
470:     ksp->its++;
471:     PetscObjectSAWsGrantAccess((PetscObject)ksp);

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

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

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

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

550:   /* make sure solution is in vector x */
551:   ksp->vec_sol = sol_orig;
552:   if (k) {
553:     VecCopy(p[k],sol_orig);
554:   }
555:   if (ksp->reason == KSP_CONVERGED_ITS) {
556:     KSPLogErrorHistory(ksp);
557:   }
558:   return 0;
559: }

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

566:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
567:   if (iascii) {
568:     PetscReal emax,emin;
569:     KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin);
570:     PetscViewerASCIIPrintf(viewer,"  eigenvalue targets used: min %g, max %g\n",(double)emin,(double)emax);
571:     if (cheb->kspest) {
572:       PetscViewerASCIIPrintf(viewer,"  eigenvalues estimated via %s: min %g, max %g\n",((PetscObject)(cheb->kspest))->type_name,(double)cheb->emin_computed,(double)cheb->emax_computed);
573:       PetscViewerASCIIPrintf(viewer,"  eigenvalues estimated using %s with transform: [%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]);
574:       PetscViewerASCIIPushTab(viewer);
575:       KSPView(cheb->kspest,viewer);
576:       PetscViewerASCIIPopTab(viewer);
577:       if (cheb->usenoisy) {
578:         PetscViewerASCIIPrintf(viewer,"  estimating eigenvalues using noisy right hand side\n");
579:       }
580:     } else if (cheb->emax_provided != 0.) {
581:       PetscViewerASCIIPrintf(viewer, "  eigenvalues provided (min %g, max %g) with transform: [%g %g; %g %g]\n", (double)cheb->emin_provided, (double)cheb->emax_provided, (double)cheb->tform[0],(double)cheb->tform[1],(double)cheb->tform[2],(double)cheb->tform[3]);
582:     }
583:   }
584:   return 0;
585: }

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

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

600: /*MC
601:      KSPCHEBYSHEV - The preconditioned Chebyshev iterative method

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

611:    Level: beginner

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

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

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

625: M*/

627: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
628: {
629:   KSP_Chebyshev  *chebyshevP;

631:   PetscNewLog(ksp,&chebyshevP);

633:   ksp->data = (void*)chebyshevP;
634:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
635:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
636:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
637:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);

639:   chebyshevP->emin = 0.;
640:   chebyshevP->emax = 0.;

642:   chebyshevP->tform[0] = 0.0;
643:   chebyshevP->tform[1] = 0.1;
644:   chebyshevP->tform[2] = 0;
645:   chebyshevP->tform[3] = 1.1;
646:   chebyshevP->eststeps = 10;
647:   chebyshevP->usenoisy = PETSC_TRUE;
648:   ksp->setupnewmatrix = PETSC_TRUE;

650:   ksp->ops->setup          = KSPSetUp_Chebyshev;
651:   ksp->ops->solve          = KSPSolve_Chebyshev;
652:   ksp->ops->destroy        = KSPDestroy_Chebyshev;
653:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
654:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
655:   ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
656:   ksp->ops->view           = KSPView_Chebyshev;
657:   ksp->ops->reset          = KSPReset_Chebyshev;

659:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",KSPChebyshevSetEigenvalues_Chebyshev);
660:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",KSPChebyshevEstEigSet_Chebyshev);
661:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",KSPChebyshevEstEigSetUseNoisy_Chebyshev);
662:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",KSPChebyshevEstEigGetKSP_Chebyshev);
663:   return 0;
664: }