Actual source code: cheby.c

petsc-main 2021-04-20
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: /*
 17:  * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
 18:  */
 19: static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest,PetscReal *emin,PetscReal *emax)
 20: {
 22:   PetscInt       n,neig;
 23:   PetscReal      *re,*im,min,max;

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

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

 80:         if (!ksp->vec_rhs) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Chebyshev must use a noisy right hand side to estimate the eigenvalues when no right hand side is available");
 81:         PCPreSolveChangeRHS(ksp->pc,&change);
 82:         if (change) {
 83:           B = ksp->work[1];
 84:           VecCopy(ksp->vec_rhs,B);
 85:         } else B = ksp->vec_rhs;
 86:       }
 87:       KSPSolve(cheb->kspest,B,ksp->work[0]);
 88:       KSPGetConvergedReason(cheb->kspest,&reason);
 89:       if (reason == KSP_DIVERGED_ITS) {
 90:         PetscInfo(ksp,"Eigen estimator ran for prescribed number of iterations\n");
 91:       } else if (reason == KSP_DIVERGED_PC_FAILED) {
 92:         PetscInt       its;
 93:         PCFailedReason pcreason;

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

113:       KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest,&min,&max);
114:       KSPSetPC(cheb->kspest,NULL);

116:       cheb->emin_computed = min;
117:       cheb->emax_computed = max;
118:       cheb->emin = cheb->tform[0]*min + cheb->tform[1]*max;
119:       cheb->emax = cheb->tform[2]*min + cheb->tform[3]*max;

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

130: static PetscErrorCode KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp,PetscReal emax,PetscReal emin)
131: {
132:   KSP_Chebyshev  *chebyshevP = (KSP_Chebyshev*)ksp->data;

136:   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);
137:   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);
138:   chebyshevP->emax = emax;
139:   chebyshevP->emin = emin;

141:   KSPChebyshevEstEigSet(ksp,0.,0.,0.,0.); /* Destroy any estimation setup */
142:   return(0);
143: }

145: static PetscErrorCode KSPChebyshevEstEigSet_Chebyshev(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
146: {
147:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

151:   if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
152:     if (!cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
153:       KSPCreate(PetscObjectComm((PetscObject)ksp),&cheb->kspest);
154:       KSPSetErrorIfNotConverged(cheb->kspest,ksp->errorifnotconverged);
155:       PetscObjectIncrementTabLevel((PetscObject)cheb->kspest,(PetscObject)ksp,1);
156:       /* use PetscObjectSet/AppendOptionsPrefix() instead of KSPSet/AppendOptionsPrefix() so that the PC prefix is not changed */
157:       PetscObjectSetOptionsPrefix((PetscObject)cheb->kspest,((PetscObject)ksp)->prefix);
158:       PetscObjectAppendOptionsPrefix((PetscObject)cheb->kspest,"esteig_");
159:       KSPSetSkipPCSetFromOptions(cheb->kspest,PETSC_TRUE);

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

163:       /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
164:       KSPSetTolerances(cheb->kspest,1.e-12,PETSC_DEFAULT,PETSC_DEFAULT,cheb->eststeps);
165:     }
166:     if (a >= 0) cheb->tform[0] = a;
167:     if (b >= 0) cheb->tform[1] = b;
168:     if (c >= 0) cheb->tform[2] = c;
169:     if (d >= 0) cheb->tform[3] = d;
170:     cheb->amatid    = 0;
171:     cheb->pmatid    = 0;
172:     cheb->amatstate = -1;
173:     cheb->pmatstate = -1;
174:   } else {
175:     KSPDestroy(&cheb->kspest);
176:   }
177:   return(0);
178: }

180: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp,PetscBool use)
181: {
182:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

185:   cheb->usenoisy = use;
186:   return(0);
187: }

189: /*@
190:    KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues
191:    of the preconditioned problem.

193:    Logically Collective on ksp

195:    Input Parameters:
196: +  ksp - the Krylov space context
197: -  emax, emin - the eigenvalue estimates

199:   Options Database:
200: .  -ksp_chebyshev_eigenvalues emin,emax

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

205:    Level: intermediate

207: @*/
208: PetscErrorCode  KSPChebyshevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
209: {

216:   PetscTryMethod(ksp,"KSPChebyshevSetEigenvalues_C",(KSP,PetscReal,PetscReal),(ksp,emax,emin));
217:   return(0);
218: }

220: /*@
221:    KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev

223:    Logically Collective on ksp

225:    Input Parameters:
226: +  ksp - the Krylov space context
227: .  a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
228: .  b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
229: .  c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
230: -  d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)

232:   Options Database:
233: .  -ksp_chebyshev_esteig a,b,c,d

235:    Notes:
236:    The Chebyshev bounds are set using
237: .vb
238:    minbound = a*minest + b*maxest
239:    maxbound = c*minest + d*maxest
240: .ve
241:    The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
242:    The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.

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

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

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

250:    Level: intermediate

252: @*/
253: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
254: {

263:   PetscTryMethod(ksp,"KSPChebyshevEstEigSet_C",(KSP,PetscReal,PetscReal,PetscReal,PetscReal),(ksp,a,b,c,d));
264:   return(0);
265: }

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

270:    Logically Collective

272:    Input Arguments:
273: +  ksp - linear solver context
274: -  use - PETSC_TRUE to use noisy

276:    Options Database:
277: .  -ksp_chebyshev_esteig_noisy <true,false>

279:   Notes:
280:     This alledgely works better for multigrid smoothers

282:   Level: intermediate

284: .seealso: KSPChebyshevEstEigSet()
285: @*/
286: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp,PetscBool use)
287: {

292:   PetscTryMethod(ksp,"KSPChebyshevEstEigSetUseNoisy_C",(KSP,PetscBool),(ksp,use));
293:   return(0);
294: }

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

301:   Input Parameters:
302: . ksp - the Krylov space context

304:   Output Parameters:
305: . kspest the eigenvalue estimation Krylov space context

307:   Level: intermediate

309: .seealso: KSPChebyshevEstEigSet()
310: @*/
311: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
312: {

318:   *kspest = NULL;
319:   PetscTryMethod(ksp,"KSPChebyshevEstEigGetKSP_C",(KSP,KSP*),(ksp,kspest));
320:   return(0);
321: }

323: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
324: {
325:   KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;

328:   *kspest = cheb->kspest;
329:   return(0);
330: }

332: static PetscErrorCode KSPSetFromOptions_Chebyshev(PetscOptionItems *PetscOptionsObject,KSP ksp)
333: {
334:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;
336:   PetscInt       neigarg = 2, nestarg = 4;
337:   PetscReal      eminmax[2] = {0., 0.};
338:   PetscReal      tform[4] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
339:   PetscBool      flgeig, flgest;

342:   PetscOptionsHead(PetscOptionsObject,"KSP Chebyshev Options");
343:   PetscOptionsInt("-ksp_chebyshev_esteig_steps","Number of est steps in Chebyshev","",cheb->eststeps,&cheb->eststeps,NULL);
344:   PetscOptionsRealArray("-ksp_chebyshev_eigenvalues","extreme eigenvalues","KSPChebyshevSetEigenvalues",eminmax,&neigarg,&flgeig);
345:   if (flgeig) {
346:     if (neigarg != 2) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
347:     KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]);
348:   }
349:   PetscOptionsRealArray("-ksp_chebyshev_esteig","estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds","KSPChebyshevEstEigSet",tform,&nestarg,&flgest);
350:   if (flgest) {
351:     switch (nestarg) {
352:     case 0:
353:       KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
354:       break;
355:     case 2:                     /* Base everything on the max eigenvalues */
356:       KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,tform[0],PETSC_DECIDE,tform[1]);
357:       break;
358:     case 4:                     /* Use the full 2x2 linear transformation */
359:       KSPChebyshevEstEigSet(ksp,tform[0],tform[1],tform[2],tform[3]);
360:       break;
361:     default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
362:     }
363:   }

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

370:   if (cheb->kspest) {
371:     PetscOptionsBool("-ksp_chebyshev_esteig_noisy","Use noisy right hand side for estimate","KSPChebyshevEstEigSetUseNoisy",cheb->usenoisy,&cheb->usenoisy,NULL);
372:     KSPSetFromOptions(cheb->kspest);
373:   }
374:   PetscOptionsTail();
375:   return(0);
376: }

378: static PetscErrorCode KSPSolve_Chebyshev(KSP ksp)
379: {
380:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;
382:   PetscInt       k,kp1,km1,ktmp,i;
383:   PetscScalar    alpha,omegaprod,mu,omega,Gamma,c[3],scale;
384:   PetscReal      rnorm = 0.0;
385:   Vec            sol_orig,b,p[3],r;
386:   Mat            Amat,Pmat;
387:   PetscBool      diagonalscale;

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

393:   PCGetOperators(ksp->pc,&Amat,&Pmat);
394:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
395:   ksp->its = 0;
396:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
397:   /* These three point to the three active solutions, we
398:      rotate these three at each solution update */
399:   km1      = 0; k = 1; kp1 = 2;
400:   sol_orig = ksp->vec_sol; /* ksp->vec_sol will be assigned to rotating vector p[k], thus save its address */
401:   b        = ksp->vec_rhs;
402:   p[km1]   = sol_orig;
403:   p[k]     = ksp->work[0];
404:   p[kp1]   = ksp->work[1];
405:   r        = ksp->work[2];

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

410:   /*   -alpha <=  scale*lambda(B^{-1}A) <= alpha   */
411:   alpha     = 1.0 - scale*(cheb->emin);
412:   Gamma     = 1.0;
413:   mu        = 1.0/alpha;
414:   omegaprod = 2.0/alpha;

416:   c[km1] = 1.0;
417:   c[k]   = mu;

419:   if (!ksp->guess_zero) {
420:     KSP_MatMult(ksp,Amat,sol_orig,r);     /*  r = b - A*p[km1] */
421:     VecAYPX(r,-1.0,b);
422:   } else {
423:     VecCopy(b,r);
424:   }

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

459:   for (i=1; i<ksp->max_it; i++) {
460:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
461:     ksp->its++;
462:     PetscObjectSAWsGrantAccess((PetscObject)ksp);

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

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

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

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

541:   /* make sure solution is in vector x */
542:   ksp->vec_sol = sol_orig;
543:   if (k) {
544:     VecCopy(p[k],sol_orig);
545:   }
546:   if (ksp->reason == KSP_CONVERGED_ITS) {
547:     KSPLogErrorHistory(ksp);
548:   }
549:   return(0);
550: }

552: static  PetscErrorCode KSPView_Chebyshev(KSP ksp,PetscViewer viewer)
553: {
554:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;
556:   PetscBool      iascii;

559:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
560:   if (iascii) {
561:     PetscViewerASCIIPrintf(viewer,"  eigenvalue estimates used:  min = %g, max = %g\n",(double)cheb->emin,(double)cheb->emax);
562:     if (cheb->kspest) {
563:       PetscViewerASCIIPrintf(viewer,"  eigenvalues estimate via %s min %g, max %g\n",((PetscObject)(cheb->kspest))->type_name,(double)cheb->emin_computed,(double)cheb->emax_computed);
564:       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]);
565:       PetscViewerASCIIPushTab(viewer);
566:       KSPView(cheb->kspest,viewer);
567:       PetscViewerASCIIPopTab(viewer);
568:       if (cheb->usenoisy) {
569:         PetscViewerASCIIPrintf(viewer,"  estimating eigenvalues using noisy right hand side\n");
570:       }
571:     }
572:   }
573:   return(0);
574: }

576: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
577: {
578:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

582:   KSPDestroy(&cheb->kspest);
583:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",NULL);
584:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",NULL);
585:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",NULL);
586:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",NULL);
587:   KSPDestroyDefault(ksp);
588:   return(0);
589: }

591: /*MC
592:      KSPCHEBYSHEV - The preconditioned Chebyshev iterative method

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

602:    Level: beginner

604:    Notes:
605:     The Chebyshev method requires both the matrix and preconditioner to
606:           be symmetric positive (semi) definite.
607:           Only support for left preconditioning.

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

612: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
613:            KSPChebyshevSetEigenvalues(), KSPChebyshevEstEigSet(), KSPChebyshevEstEigSetUseNoisy()
614:            KSPRICHARDSON, KSPCG, PCMG

616: M*/

618: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
619: {
621:   KSP_Chebyshev  *chebyshevP;

624:   PetscNewLog(ksp,&chebyshevP);

626:   ksp->data = (void*)chebyshevP;
627:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
628:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
629:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
630:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);

632:   chebyshevP->emin = 0.;
633:   chebyshevP->emax = 0.;

635:   chebyshevP->tform[0] = 0.0;
636:   chebyshevP->tform[1] = 0.1;
637:   chebyshevP->tform[2] = 0;
638:   chebyshevP->tform[3] = 1.1;
639:   chebyshevP->eststeps = 10;
640:   chebyshevP->usenoisy = PETSC_TRUE;
641:   ksp->setupnewmatrix = PETSC_TRUE;

643:   ksp->ops->setup          = KSPSetUp_Chebyshev;
644:   ksp->ops->solve          = KSPSolve_Chebyshev;
645:   ksp->ops->destroy        = KSPDestroy_Chebyshev;
646:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
647:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
648:   ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
649:   ksp->ops->view           = KSPView_Chebyshev;
650:   ksp->ops->reset          = KSPReset_Chebyshev;

652:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",KSPChebyshevSetEigenvalues_Chebyshev);
653:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",KSPChebyshevEstEigSet_Chebyshev);
654:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",KSPChebyshevEstEigSetUseNoisy_Chebyshev);
655:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",KSPChebyshevEstEigGetKSP_Chebyshev);
656:   return(0);
657: }