Actual source code: cheby.c

petsc-master 2016-06-29
Report Typos and Errors
  2: #include <../src/ksp/ksp/impls/cheby/chebyshevimpl.h>    /*I "petscksp.h" I*/

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

 12:   KSPReset(cheb->kspest);
 13:   return(0);
 14: }

 18: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
 19: {
 20:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

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

 33: static PetscErrorCode KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp,PetscReal emax,PetscReal emin)
 34: {
 35:   KSP_Chebyshev  *chebyshevP = (KSP_Chebyshev*)ksp->data;

 39:   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);
 40:   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);
 41:   chebyshevP->emax = emax;
 42:   chebyshevP->emin = emin;

 44:   KSPChebyshevEstEigSet(ksp,0.,0.,0.,0.); /* Destroy any estimation setup */
 45:   return(0);
 46: }

 50: static PetscErrorCode KSPChebyshevEstEigSet_Chebyshev(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
 51: {
 52:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

 56:   if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
 57:     if (!cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
 58:       KSPCreate(PetscObjectComm((PetscObject)ksp),&cheb->kspest);
 59:       PetscObjectIncrementTabLevel((PetscObject)cheb->kspest,(PetscObject)ksp,1);
 60:       KSPSetOptionsPrefix(cheb->kspest,((PetscObject)ksp)->prefix);
 61:       KSPAppendOptionsPrefix(cheb->kspest,"esteig_");
 62:       KSPSetSkipPCSetFromOptions(cheb->kspest,PETSC_TRUE);

 64:       KSPSetPC(cheb->kspest,ksp->pc);

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

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

 87: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp,PetscBool use)
 88: {
 89:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

 92:   cheb->usenoisy = use;
 93:   return(0);
 94: }

 98: /*@
 99:    KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues
100:    of the preconditioned problem.

102:    Logically Collective on KSP

104:    Input Parameters:
105: +  ksp - the Krylov space context
106: -  emax, emin - the eigenvalue estimates

108:   Options Database:
109: .  -ksp_chebyshev_eigenvalues emin,emax

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

114:    Level: intermediate

116: .keywords: KSP, Chebyshev, set, eigenvalues
117: @*/
118: PetscErrorCode  KSPChebyshevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
119: {

126:   PetscTryMethod(ksp,"KSPChebyshevSetEigenvalues_C",(KSP,PetscReal,PetscReal),(ksp,emax,emin));
127:   return(0);
128: }

132: /*@
133:    KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev

135:    Logically Collective on KSP

137:    Input Parameters:
138: +  ksp - the Krylov space context
139: .  a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
140: .  b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
141: .  c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
142: -  d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)

144:   Options Database:
145: .  -ksp_chebyshev_esteig a,b,c,d

147:    Notes:
148:    The Chebyshev bounds are set using
149: .vb
150:    minbound = a*minest + b*maxest
151:    maxbound = c*minest + d*maxest
152: .ve
153:    The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
154:    The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.

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

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

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

162:    Level: intermediate

164: .keywords: KSP, Chebyshev, set, eigenvalues, PCMG
165: @*/
166: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
167: {

176:   PetscTryMethod(ksp,"KSPChebyshevEstEigSet_C",(KSP,PetscReal,PetscReal,PetscReal,PetscReal),(ksp,a,b,c,d));
177:   return(0);
178: }

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

185:    Logically Collective

187:    Input Arguments:
188: +  ksp - linear solver context
189: -  use - PETSC_TRUE to use noisy

191:    Options Database:
192: +  -ksp_chebyshev_esteig_noisy <true,false>

194:   Notes: This alledgely works better for multigrid smoothers

196:   Level: intermediate

198: .seealso: KSPChebyshevEstEigSet()
199: @*/
200: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp,PetscBool use)
201: {

206:   PetscTryMethod(ksp,"KSPChebyshevEstEigSetUseNoisy_C",(KSP,PetscBool),(ksp,use));
207:   return(0);
208: }

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

217:   Input Parameters:
218: . ksp - the Krylov space context

220:   Output Parameters:
221: . kspest the eigenvalue estimation Krylov space context

223:   Level: intermediate

225: .seealso: KSPChebyshevEstEigSet()
226: @*/
227: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
228: {

233:   *kspest = NULL;
234:   PetscTryMethod(ksp,"KSPChebyshevEstEigGetKSP_C",(KSP,KSP*),(ksp,kspest));
235:   return(0);
236: }

240: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
241: {
242:   KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;

245:   *kspest = cheb->kspest;
246:   return(0);
247: }

251: static PetscErrorCode KSPSetFromOptions_Chebyshev(PetscOptionItems *PetscOptionsObject,KSP ksp)
252: {
253:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;
255:   PetscInt       neigarg = 2, nestarg = 4;
256:   PetscReal      eminmax[2] = {0., 0.};
257:   PetscReal      tform[4] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
258:   PetscBool      flgeig, flgest;

261:   PetscOptionsHead(PetscOptionsObject,"KSP Chebyshev Options");
262:   PetscOptionsInt("-ksp_chebyshev_esteig_steps","Number of est steps in Chebyshev","",cheb->eststeps,&cheb->eststeps,NULL);
263:   PetscOptionsRealArray("-ksp_chebyshev_eigenvalues","extreme eigenvalues","KSPChebyshevSetEigenvalues",eminmax,&neigarg,&flgeig);
264:   if (flgeig) {
265:     if (neigarg != 2) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
266:     KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]);
267:   }
268:   PetscOptionsRealArray("-ksp_chebyshev_esteig","estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds","KSPChebyshevEstEigSet",tform,&nestarg,&flgest);
269:   if (flgest) {
270:     switch (nestarg) {
271:     case 0:
272:       KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
273:       break;
274:     case 2:                     /* Base everything on the max eigenvalues */
275:       KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,tform[0],PETSC_DECIDE,tform[1]);
276:       break;
277:     case 4:                     /* Use the full 2x2 linear transformation */
278:       KSPChebyshevEstEigSet(ksp,tform[0],tform[1],tform[2],tform[3]);
279:       break;
280:     default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
281:     }
282:   }

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

289:   if (cheb->kspest) {
290:     PetscOptionsBool("-ksp_chebyshev_esteig_noisy","Use noisy right hand side for estimate","KSPChebyshevEstEigSetUseNoisy",cheb->usenoisy,&cheb->usenoisy,NULL);
291:     KSPSetFromOptions(cheb->kspest);
292:   }
293:   PetscOptionsTail();
294:   return(0);
295: }

299: /*
300:  * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
301:  */
302: static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest,PetscReal *emin,PetscReal *emax)
303: {
305:   PetscInt       n,neig;
306:   PetscReal      *re,*im,min,max;

309:   KSPGetIterationNumber(kspest,&n);
310:   PetscMalloc2(n,&re,n,&im);
311:   KSPComputeEigenvalues(kspest,n,re,im,&neig);
312:   min  = PETSC_MAX_REAL;
313:   max  = PETSC_MIN_REAL;
314:   for (n=0; n<neig; n++) {
315:     min = PetscMin(min,re[n]);
316:     max = PetscMax(max,re[n]);
317:   }
318:   PetscFree2(re,im);
319:   *emax = max;
320:   *emin = min;
321:   return(0);
322: }

324: static PetscScalar chebyhash(PetscInt xx) {
325:   unsigned int x = xx;
326:   x = ((x >> 16) ^ x) * 0x45d9f3b;
327:   x = ((x >> 16) ^ x) * 0x45d9f3b;
328:   x = ((x >> 16) ^ x);
329:   return (PetscScalar)((PetscInt64)x-2147483648)*5.e-10; /* center around zero, scaled about -1. to 1.*/
330: }
333: static PetscErrorCode KSPSolve_Chebyshev(KSP ksp)
334: {
335:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;
337:   PetscInt       k,kp1,km1,maxit,ktmp,i;
338:   PetscScalar    alpha,omegaprod,mu,omega,Gamma,c[3],scale;
339:   PetscReal      rnorm = 0.0;
340:   Vec            sol_orig,b,p[3],r;
341:   Mat            Amat,Pmat;
342:   PetscBool      diagonalscale;

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

348:   PCGetOperators(ksp->pc,&Amat,&Pmat);
349:   if (cheb->kspest) {
350:     PetscObjectId    amatid,    pmatid;
351:     PetscObjectState amatstate, pmatstate;

353:     PetscObjectGetId((PetscObject)Amat,&amatid);
354:     PetscObjectGetId((PetscObject)Pmat,&pmatid);
355:     PetscObjectStateGet((PetscObject)Amat,&amatstate);
356:     PetscObjectStateGet((PetscObject)Pmat,&pmatstate);
357:     if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
358:       PetscReal          max=0.0,min=0.0;
359:       Vec                B;
360:       KSPConvergedReason reason;

362:       if (cheb->usenoisy) {
363:         B  = ksp->work[1];
364:         {
366:           PetscInt       n,i,istart;
367:           PetscScalar    *xx;
368:           VecGetOwnershipRange(B,&istart,NULL);
369:           VecGetLocalSize(B,&n);
370:           VecGetArray(B,&xx);
371:           for (i=0; i<n; i++) {
372:             PetscScalar v = chebyhash(i+istart);
373:             xx[i] = v;
374:           }
375:           VecRestoreArray(B,&xx);
376:         }
377:       } else {
378:         PC        pc;
379:         PetscBool change;

381:         KSPGetPC(cheb->kspest,&pc);
382:         PCPreSolveChangeRHS(pc,&change);
383:         if (change) {
384:           B = ksp->work[1];
385:           VecCopy(ksp->vec_rhs,B);
386:         } else {
387:           B = ksp->vec_rhs;
388:         }
389:       }
390:       KSPSolve(cheb->kspest,B,ksp->work[0]);

392:       KSPGetConvergedReason(cheb->kspest,&reason);
393:       if (reason < 0) {
394:         if (reason == KSP_DIVERGED_ITS) {
395:           PetscInfo(ksp,"Eigen estimator ran for prescribed number of iterations\n");
396:         } else {
397:           PetscInt its;
398:           KSPGetIterationNumber(cheb->kspest,&its);
399:           if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
400:             PetscInfo(ksp,"Eigen estimator KSP_DIVERGED_PCSETUP_FAILED\n");
401:             ksp->reason = KSP_DIVERGED_PCSETUP_FAILED;
402:             VecSetInf(ksp->vec_sol);
403:           } else {
404:             SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Eigen estimator failed: %s at iteration %D",KSPConvergedReasons[reason],its);
405:           }
406:         }
407:       } else if (reason==KSP_CONVERGED_RTOL ||reason==KSP_CONVERGED_ATOL) {
408:         PetscInfo(ksp,"Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n");
409:       } else {
410:         PetscInfo1(ksp,"Eigen estimator did not converge by iteration: %s\n",KSPConvergedReasons[reason]);
411:       }

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

415:       cheb->emin = cheb->tform[0]*min + cheb->tform[1]*max;
416:       cheb->emax = cheb->tform[2]*min + cheb->tform[3]*max;

418:       cheb->amatid    = amatid;
419:       cheb->pmatid    = pmatid;
420:       cheb->amatstate = amatstate;
421:       cheb->pmatstate = pmatstate;
422:     }
423:   }

425:   ksp->its = 0;
426:   maxit    = ksp->max_it;

428:   /* These three point to the three active solutions, we
429:      rotate these three at each solution update */
430:   km1      = 0; k = 1; kp1 = 2;
431:   sol_orig = ksp->vec_sol; /* ksp->vec_sol will be asigned to rotating vector p[k], thus save its address */
432:   b        = ksp->vec_rhs;
433:   p[km1]   = sol_orig;
434:   p[k]     = ksp->work[0];
435:   p[kp1]   = ksp->work[1];
436:   r        = ksp->work[2];

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

441:   /*   -alpha <=  scale*lambda(B^{-1}A) <= alpha   */
442:   alpha     = 1.0 - scale*(cheb->emin);
443:   Gamma     = 1.0;
444:   mu        = 1.0/alpha;
445:   omegaprod = 2.0/alpha;

447:   c[km1] = 1.0;
448:   c[k]   = mu;

450:   if (!ksp->guess_zero) {
451:     KSP_MatMult(ksp,Amat,p[km1],r);     /*  r = b - A*p[km1] */
452:     VecAYPX(r,-1.0,b);
453:   } else {
454:     VecCopy(b,r);
455:   }

457:   KSP_PCApply(ksp,r,p[k]);  /* p[k] = scale B^{-1}r + p[km1] */
458:   VecAYPX(p[k],scale,p[km1]);

460:   for (i=0; i<maxit; i++) {
461:     PetscObjectSAWsTakeAccess((PetscObject)ksp);

463:     ksp->its++;
464:     PetscObjectSAWsGrantAccess((PetscObject)ksp);
465:     c[kp1] = 2.0*mu*c[k] - c[km1];
466:     omega  = omegaprod*c[k]/c[kp1];

468:     KSP_MatMult(ksp,Amat,p[k],r);          /*  r = b - Ap[k]    */
469:     VecAYPX(r,-1.0,b);
470:     KSP_PCApply(ksp,r,p[kp1]);             /*  p[kp1] = B^{-1}r  */
471:     ksp->vec_sol = p[k];

473:     /* calculate residual norm if requested */
474:     if (ksp->normtype != KSP_NORM_NONE || ksp->numbermonitors) {
475:       if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
476:         VecNorm(r,NORM_2,&rnorm);
477:       } else {
478:         VecNorm(p[kp1],NORM_2,&rnorm);
479:       }
480:       PetscObjectSAWsTakeAccess((PetscObject)ksp);
481:       ksp->rnorm   = rnorm;
482:       PetscObjectSAWsGrantAccess((PetscObject)ksp);
483:       KSPLogResidualHistory(ksp,rnorm);
484:       KSPMonitor(ksp,i,rnorm);
485:       (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
486:       if (ksp->reason) break;
487:     }

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

492:     ktmp = km1;
493:     km1  = k;
494:     k    = kp1;
495:     kp1  = ktmp;
496:   }
497:   if (!ksp->reason) {
498:     if (ksp->normtype != KSP_NORM_NONE) {
499:       KSP_MatMult(ksp,Amat,p[k],r);       /*  r = b - Ap[k]    */
500:       VecAYPX(r,-1.0,b);
501:       if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
502:         VecNorm(r,NORM_2,&rnorm);
503:       } else {
504:         KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
505:         VecNorm(p[kp1],NORM_2,&rnorm);
506:       }
507:       PetscObjectSAWsTakeAccess((PetscObject)ksp);
508:       ksp->rnorm   = rnorm;
509:       PetscObjectSAWsGrantAccess((PetscObject)ksp);
510:       ksp->vec_sol = p[k];
511:       KSPLogResidualHistory(ksp,rnorm);
512:       KSPMonitor(ksp,i,rnorm);
513:     }
514:     if (ksp->its >= ksp->max_it) {
515:       if (ksp->normtype != KSP_NORM_NONE) {
516:         (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
517:         if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
518:       } else ksp->reason = KSP_CONVERGED_ITS;
519:     }
520:   }

522:   /* make sure solution is in vector x */
523:   ksp->vec_sol = sol_orig;
524:   if (k) {
525:     VecCopy(p[k],sol_orig);
526:   }
527:   return(0);
528: }

532: static  PetscErrorCode KSPView_Chebyshev(KSP ksp,PetscViewer viewer)
533: {
534:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;
536:   PetscBool      iascii;

539:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
540:   if (iascii) {
541:     PetscViewerASCIIPrintf(viewer,"  Chebyshev: eigenvalue estimates:  min = %g, max = %g\n",(double)cheb->emin,(double)cheb->emax);
542:     if (cheb->kspest) {
543:       PetscViewerASCIIPrintf(viewer,"  Chebyshev: 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]);
544:       PetscViewerASCIIPushTab(viewer);
545:       KSPView(cheb->kspest,viewer);
546:       PetscViewerASCIIPopTab(viewer);
547:       if (cheb->usenoisy) {
548:         PetscViewerASCIIPrintf(viewer,"  Chebyshev: estimating eigenvalues using noisy right hand side\n");
549:       }
550:     }
551:   }
552:   return(0);
553: }

557: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
558: {
559:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

563:   KSPDestroy(&cheb->kspest);
564:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",NULL);
565:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",NULL);
566:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",NULL);
567:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",NULL);
568:   KSPDestroyDefault(ksp);
569:   return(0);
570: }

572: /*MC
573:      KSPCHEBYSHEV - The preconditioned Chebyshev iterative method

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

583:    Level: beginner

585:    Notes: The Chebyshev method requires both the matrix and preconditioner to
586:           be symmetric positive (semi) definite.
587:           Only support for left preconditioning.

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

592: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
593:            KSPChebyshevSetEigenvalues(), KSPChebyshevEstEigSet(), KSPChebyshevEstEigSetUseNoisy()
594:            KSPRICHARDSON, KSPCG, PCMG

596: M*/

600: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
601: {
603:   KSP_Chebyshev  *chebyshevP;

606:   PetscNewLog(ksp,&chebyshevP);

608:   ksp->data = (void*)chebyshevP;
609:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
610:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);

612:   chebyshevP->emin = 0.;
613:   chebyshevP->emax = 0.;

615:   chebyshevP->tform[0] = 0.0;
616:   chebyshevP->tform[1] = 0.1;
617:   chebyshevP->tform[2] = 0;
618:   chebyshevP->tform[3] = 1.1;
619:   chebyshevP->eststeps = 10;
620:   chebyshevP->usenoisy = PETSC_TRUE;

622:   ksp->ops->setup          = KSPSetUp_Chebyshev;
623:   ksp->ops->solve          = KSPSolve_Chebyshev;
624:   ksp->ops->destroy        = KSPDestroy_Chebyshev;
625:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
626:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
627:   ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
628:   ksp->ops->view           = KSPView_Chebyshev;
629:   ksp->ops->reset          = KSPReset_Chebyshev;

631:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",KSPChebyshevSetEigenvalues_Chebyshev);
632:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",KSPChebyshevEstEigSet_Chebyshev);
633:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",KSPChebyshevEstEigSetUseNoisy_Chebyshev);
634:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",KSPChebyshevEstEigGetKSP_Chebyshev);
635:   return(0);
636: }