Actual source code: cheby.c

petsc-3.11.2 2019-05-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:   KSPReset(cheb->kspest);
 11:   return(0);
 12: }

 14: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
 15: {
 16:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

 20:   KSPSetWorkVecs(ksp,3);
 21:   if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) { /* We need to estimate eigenvalues */
 22:     KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
 23:   }
 24:   return(0);
 25: }

 27: static PetscErrorCode KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp,PetscReal emax,PetscReal emin)
 28: {
 29:   KSP_Chebyshev  *chebyshevP = (KSP_Chebyshev*)ksp->data;

 33:   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);
 34:   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);
 35:   chebyshevP->emax = emax;
 36:   chebyshevP->emin = emin;

 38:   KSPChebyshevEstEigSet(ksp,0.,0.,0.,0.); /* Destroy any estimation setup */
 39:   return(0);
 40: }

 42: static PetscErrorCode KSPChebyshevEstEigSet_Chebyshev(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
 43: {
 44:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

 48:   if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
 49:     if (!cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
 50:       KSPCreate(PetscObjectComm((PetscObject)ksp),&cheb->kspest);
 51:       KSPSetErrorIfNotConverged(cheb->kspest,ksp->errorifnotconverged);
 52:       PetscObjectIncrementTabLevel((PetscObject)cheb->kspest,(PetscObject)ksp,1);
 53:       KSPSetPC(cheb->kspest,ksp->pc);
 54:       /* use PetscObjectSet/AppendOptionsPrefix() instead of KSPSet/AppendOptionsPrefix() so that the PC prefix is not changed */
 55:       PetscObjectSetOptionsPrefix((PetscObject)cheb->kspest,((PetscObject)ksp)->prefix);
 56:       PetscObjectAppendOptionsPrefix((PetscObject)cheb->kspest,"esteig_");
 57:       KSPSetSkipPCSetFromOptions(cheb->kspest,PETSC_TRUE);

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

 61:       /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
 62:       KSPSetTolerances(cheb->kspest,1.e-12,PETSC_DEFAULT,PETSC_DEFAULT,cheb->eststeps);
 63:     }
 64:     if (a >= 0) cheb->tform[0] = a;
 65:     if (b >= 0) cheb->tform[1] = b;
 66:     if (c >= 0) cheb->tform[2] = c;
 67:     if (d >= 0) cheb->tform[3] = d;
 68:     cheb->amatid    = 0;
 69:     cheb->pmatid    = 0;
 70:     cheb->amatstate = -1;
 71:     cheb->pmatstate = -1;
 72:   } else {
 73:     KSPDestroy(&cheb->kspest);
 74:   }
 75:   return(0);
 76: }

 78: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp,PetscBool use)
 79: {
 80:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

 83:   cheb->usenoisy = use;
 84:   return(0);
 85: }

 87: /*@
 88:    KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues
 89:    of the preconditioned problem.

 91:    Logically Collective on KSP

 93:    Input Parameters:
 94: +  ksp - the Krylov space context
 95: -  emax, emin - the eigenvalue estimates

 97:   Options Database:
 98: .  -ksp_chebyshev_eigenvalues emin,emax

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

103:    Level: intermediate

105: .keywords: KSP, Chebyshev, set, eigenvalues
106: @*/
107: PetscErrorCode  KSPChebyshevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
108: {

115:   PetscTryMethod(ksp,"KSPChebyshevSetEigenvalues_C",(KSP,PetscReal,PetscReal),(ksp,emax,emin));
116:   return(0);
117: }

119: /*@
120:    KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev

122:    Logically Collective on KSP

124:    Input Parameters:
125: +  ksp - the Krylov space context
126: .  a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
127: .  b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
128: .  c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
129: -  d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)

131:   Options Database:
132: .  -ksp_chebyshev_esteig a,b,c,d

134:    Notes:
135:    The Chebyshev bounds are set using
136: .vb
137:    minbound = a*minest + b*maxest
138:    maxbound = c*minest + d*maxest
139: .ve
140:    The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
141:    The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.

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

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

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

149:    Level: intermediate

151: .keywords: KSP, Chebyshev, set, eigenvalues, PCMG
152: @*/
153: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
154: {

163:   PetscTryMethod(ksp,"KSPChebyshevEstEigSet_C",(KSP,PetscReal,PetscReal,PetscReal,PetscReal),(ksp,a,b,c,d));
164:   return(0);
165: }

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

170:    Logically Collective

172:    Input Arguments:
173: +  ksp - linear solver context
174: -  use - PETSC_TRUE to use noisy

176:    Options Database:
177: +  -ksp_chebyshev_esteig_noisy <true,false>

179:   Notes:
180:     This alledgely works better for multigrid smoothers

182:   Level: intermediate

184: .seealso: KSPChebyshevEstEigSet()
185: @*/
186: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp,PetscBool use)
187: {

192:   PetscTryMethod(ksp,"KSPChebyshevEstEigSetUseNoisy_C",(KSP,PetscBool),(ksp,use));
193:   return(0);
194: }

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

201:   Input Parameters:
202: . ksp - the Krylov space context

204:   Output Parameters:
205: . kspest the eigenvalue estimation Krylov space context

207:   Level: intermediate

209: .seealso: KSPChebyshevEstEigSet()
210: @*/
211: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
212: {

217:   *kspest = NULL;
218:   PetscTryMethod(ksp,"KSPChebyshevEstEigGetKSP_C",(KSP,KSP*),(ksp,kspest));
219:   return(0);
220: }

222: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
223: {
224:   KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;

227:   *kspest = cheb->kspest;
228:   return(0);
229: }

231: static PetscErrorCode KSPSetFromOptions_Chebyshev(PetscOptionItems *PetscOptionsObject,KSP ksp)
232: {
233:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;
235:   PetscInt       neigarg = 2, nestarg = 4;
236:   PetscReal      eminmax[2] = {0., 0.};
237:   PetscReal      tform[4] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
238:   PetscBool      flgeig, flgest;

241:   PetscOptionsHead(PetscOptionsObject,"KSP Chebyshev Options");
242:   PetscOptionsInt("-ksp_chebyshev_esteig_steps","Number of est steps in Chebyshev","",cheb->eststeps,&cheb->eststeps,NULL);
243:   PetscOptionsRealArray("-ksp_chebyshev_eigenvalues","extreme eigenvalues","KSPChebyshevSetEigenvalues",eminmax,&neigarg,&flgeig);
244:   if (flgeig) {
245:     if (neigarg != 2) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
246:     KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]);
247:   }
248:   PetscOptionsRealArray("-ksp_chebyshev_esteig","estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds","KSPChebyshevEstEigSet",tform,&nestarg,&flgest);
249:   if (flgest) {
250:     switch (nestarg) {
251:     case 0:
252:       KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
253:       break;
254:     case 2:                     /* Base everything on the max eigenvalues */
255:       KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,tform[0],PETSC_DECIDE,tform[1]);
256:       break;
257:     case 4:                     /* Use the full 2x2 linear transformation */
258:       KSPChebyshevEstEigSet(ksp,tform[0],tform[1],tform[2],tform[3]);
259:       break;
260:     default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
261:     }
262:   }

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

269:   if (cheb->kspest) {
270:     PetscOptionsBool("-ksp_chebyshev_esteig_noisy","Use noisy right hand side for estimate","KSPChebyshevEstEigSetUseNoisy",cheb->usenoisy,&cheb->usenoisy,NULL);
271:     KSPSetFromOptions(cheb->kspest);
272:   }
273:   PetscOptionsTail();
274:   return(0);
275: }

277: /*
278:  * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
279:  */
280: static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest,PetscReal *emin,PetscReal *emax)
281: {
283:   PetscInt       n,neig;
284:   PetscReal      *re,*im,min,max;

287:   KSPGetIterationNumber(kspest,&n);
288:   PetscMalloc2(n,&re,n,&im);
289:   KSPComputeEigenvalues(kspest,n,re,im,&neig);
290:   min  = PETSC_MAX_REAL;
291:   max  = PETSC_MIN_REAL;
292:   for (n=0; n<neig; n++) {
293:     min = PetscMin(min,re[n]);
294:     max = PetscMax(max,re[n]);
295:   }
296:   PetscFree2(re,im);
297:   *emax = max;
298:   *emin = min;
299:   return(0);
300: }

302: PETSC_STATIC_INLINE PetscScalar chebyhash(PetscInt xx)
303: {
304:   unsigned int x = xx;
305:   x = ((x >> 16) ^ x) * 0x45d9f3b;
306:   x = ((x >> 16) ^ x) * 0x45d9f3b;
307:   x = ((x >> 16) ^ x);
308:   return (PetscScalar)((PetscInt64)x-2147483648)*5.e-10; /* center around zero, scaled about -1. to 1.*/
309: }

311: static PetscErrorCode KSPSolve_Chebyshev(KSP ksp)
312: {
313:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;
315:   PetscInt       k,kp1,km1,ktmp,i;
316:   PetscScalar    alpha,omegaprod,mu,omega,Gamma,c[3],scale;
317:   PetscReal      rnorm = 0.0;
318:   Vec            sol_orig,b,p[3],r;
319:   Mat            Amat,Pmat;
320:   PetscBool      diagonalscale;

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

326:   PCGetOperators(ksp->pc,&Amat,&Pmat);
327:   if (cheb->kspest) {
328:     PetscObjectId    amatid,    pmatid;
329:     PetscObjectState amatstate, pmatstate;
330:     PCFailedReason   pcreason;
331:     PC               pc;

333:     PetscObjectGetId((PetscObject)Amat,&amatid);
334:     PetscObjectGetId((PetscObject)Pmat,&pmatid);
335:     PetscObjectStateGet((PetscObject)Amat,&amatstate);
336:     PetscObjectStateGet((PetscObject)Pmat,&pmatstate);
337:     if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
338:       PetscReal          max=0.0,min=0.0;
339:       Vec                B;
340:       KSPConvergedReason reason;

342:       if (cheb->usenoisy) {
343:         B  = ksp->work[1];
344:         {
346:           PetscInt       n,i,istart;
347:           PetscScalar    *xx;
348:           VecGetOwnershipRange(B,&istart,NULL);
349:           VecGetLocalSize(B,&n);
350:           VecGetArray(B,&xx);
351:           for (i=0; i<n; i++) {
352:             PetscScalar v = chebyhash(i+istart);
353:             xx[i] = v;
354:           }
355:           VecRestoreArray(B,&xx);
356:         }
357:       } else {
358:         PC        pc;
359:         PetscBool change;

361:         KSPGetPC(cheb->kspest,&pc);
362:         PCPreSolveChangeRHS(pc,&change);
363:         if (change) {
364:           B = ksp->work[1];
365:           VecCopy(ksp->vec_rhs,B);
366:         } else {
367:           B = ksp->vec_rhs;
368:         }
369:       }
370:       KSPSolve(cheb->kspest,B,ksp->work[0]);
371:       KSPGetConvergedReason(cheb->kspest,&reason);
372:       KSPGetPC(cheb->kspest,&pc);
373:       PCGetFailedReason(pc,&pcreason);
374:       if (reason < 0 || pcreason) {
375:         if (reason == KSP_DIVERGED_ITS) {
376:           PetscInfo(ksp,"Eigen estimator ran for prescribed number of iterations\n");
377:         } else {
378:           PetscInt its;
379:           KSPGetIterationNumber(cheb->kspest,&its);
380:           ksp->reason = KSP_DIVERGED_PC_FAILED;
381:           VecSetInf(ksp->vec_sol);
382:           PetscInfo3(ksp,"Eigen estimator failed: %s %s at iteration %D",KSPConvergedReasons[reason],PCFailedReasons[pcreason],its);
383:           return(0);
384:         }
385:       } else if (reason==KSP_CONVERGED_RTOL ||reason==KSP_CONVERGED_ATOL) {
386:         PetscInfo(ksp,"Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n");
387:       } else {
388:         PetscInfo1(ksp,"Eigen estimator did not converge by iteration: %s\n",KSPConvergedReasons[reason]);
389:       }

391:       KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest,&min,&max);

393:       cheb->emin_computed = min;
394:       cheb->emax_computed = max;
395:       cheb->emin = cheb->tform[0]*min + cheb->tform[1]*max;
396:       cheb->emax = cheb->tform[2]*min + cheb->tform[3]*max;

398:       cheb->amatid    = amatid;
399:       cheb->pmatid    = pmatid;
400:       cheb->amatstate = amatstate;
401:       cheb->pmatstate = pmatstate;
402:     }
403:   }
404:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
405:   ksp->its = 0;
406:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
407:   /* These three point to the three active solutions, we
408:      rotate these three at each solution update */
409:   km1      = 0; k = 1; kp1 = 2;
410:   sol_orig = ksp->vec_sol; /* ksp->vec_sol will be asigned to rotating vector p[k], thus save its address */
411:   b        = ksp->vec_rhs;
412:   p[km1]   = sol_orig;
413:   p[k]     = ksp->work[0];
414:   p[kp1]   = ksp->work[1];
415:   r        = ksp->work[2];

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

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

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

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

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

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

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

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

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

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

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

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

568:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
569:   if (iascii) {
570:     PetscViewerASCIIPrintf(viewer,"  eigenvalue estimates used:  min = %g, max = %g\n",(double)cheb->emin,(double)cheb->emax);
571:     if (cheb->kspest) {
572:       PetscViewerASCIIPrintf(viewer,"  eigenvalues estimate 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 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]);
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:     }
581:   }
582:   return(0);
583: }

585: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
586: {
587:   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: {
630:   KSP_Chebyshev  *chebyshevP;

633:   PetscNewLog(ksp,&chebyshevP);

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

641:   chebyshevP->emin = 0.;
642:   chebyshevP->emax = 0.;

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

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

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