Actual source code: cheby.c

petsc-master 2016-05-28
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 KSPChebyshevEstEigSetRandom_Chebyshev(KSP ksp,PetscRandom random)
 88: {
 89:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

 93:   if (random) {PetscObjectReference((PetscObject)random);}
 94:   PetscRandomDestroy(&cheb->random);
 95:   cheb->random = random;
 96:   return(0);
 97: }

101: static PetscErrorCode KSPChebyshevEstEigSetUseRandom_Chebyshev(KSP ksp,PetscBool use)
102: {
103:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

106:   cheb->userandom = use;
107:   return(0);
108: }

112: /*@
113:    KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues
114:    of the preconditioned problem.

116:    Logically Collective on KSP

118:    Input Parameters:
119: +  ksp - the Krylov space context
120: -  emax, emin - the eigenvalue estimates

122:   Options Database:
123: .  -ksp_chebyshev_eigenvalues emin,emax

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

128:    Level: intermediate

130: .keywords: KSP, Chebyshev, set, eigenvalues
131: @*/
132: PetscErrorCode  KSPChebyshevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
133: {

140:   PetscTryMethod(ksp,"KSPChebyshevSetEigenvalues_C",(KSP,PetscReal,PetscReal),(ksp,emax,emin));
141:   return(0);
142: }

146: /*@
147:    KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev

149:    Logically Collective on KSP

151:    Input Parameters:
152: +  ksp - the Krylov space context
153: .  a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
154: .  b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
155: .  c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
156: -  d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)

158:   Options Database:
159: .  -ksp_chebyshev_esteig a,b,c,d

161:    Notes:
162:    The Chebyshev bounds are set using
163: .vb
164:    minbound = a*minest + b*maxest
165:    maxbound = c*minest + d*maxest
166: .ve
167:    The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
168:    The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.

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

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

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

176:    Level: intermediate

178: .keywords: KSP, Chebyshev, set, eigenvalues, PCMG
179: @*/
180: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
181: {

190:   PetscTryMethod(ksp,"KSPChebyshevEstEigSet_C",(KSP,PetscReal,PetscReal,PetscReal,PetscReal),(ksp,a,b,c,d));
191:   return(0);
192: }

196: /*@
197:    KSPChebyshevEstEigSetUseRandom - use a random right hand side in order to do the estimate instead of the given right hand side

199:    Logically Collective

201:    Input Arguments:
202: +  ksp - linear solver context
203: -  use - PETSC_TRUE to use random

205:    Options Database:
206: +  -ksp_chebyshev_esteig_random <true,false>

208:   Notes: This alledgely works better for multigrid smoothers

210:   Use KSPChebyshevEstEigSetRandom() to provide the random number generator to be used. Otherwise it creates a default

212:   Level: intermediate

214: .seealso: KSPChebyshevEstEigSet(), PetscRandomCreate(), KSPChebyshevEstEigSetRandom()
215: @*/
216: PetscErrorCode KSPChebyshevEstEigSetUseRandom(KSP ksp,PetscBool use)
217: {

222:   PetscTryMethod(ksp,"KSPChebyshevEstEigSetUseRandom_C",(KSP,PetscBool),(ksp,use));
223:   return(0);
224: }

228: /*@
229:    KSPChebyshevEstEigSetRandom - set random context for estimating eigenvalues

231:    Logically Collective

233:    Input Arguments:
234: +  ksp - linear solver context
235: -  random - random number context or NULL to use default

237:   Level: intermediate

239: .seealso: KSPChebyshevEstEigSet(), PetscRandomCreate()
240: @*/
241: PetscErrorCode KSPChebyshevEstEigSetRandom(KSP ksp,PetscRandom random)
242: {

248:   PetscTryMethod(ksp,"KSPChebyshevEstEigSetRandom_C",(KSP,PetscRandom),(ksp,random));
249:   return(0);
250: }

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

259:   Input Parameters:
260: . ksp - the Krylov space context

262:   Output Parameters:
263: . kspest the eigenvalue estimation Krylov space context

265:   Level: intermediate

267: .seealso: KSPChebyshevEstEigSet()
268: @*/
269: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
270: {

275:   *kspest = NULL;
276:   PetscTryMethod(ksp,"KSPChebyshevEstEigGetKSP_C",(KSP,KSP*),(ksp,kspest));
277:   return(0);
278: }

282: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
283: {
284:   KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;

287:   *kspest = cheb->kspest;
288:   return(0);
289: }

293: static PetscErrorCode KSPSetFromOptions_Chebyshev(PetscOptionItems *PetscOptionsObject,KSP ksp)
294: {
295:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;
297:   PetscInt       neigarg = 2, nestarg = 4;
298:   PetscReal      eminmax[2] = {0., 0.};
299:   PetscReal      tform[4] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
300:   PetscBool      flgeig, flgest;

303:   PetscOptionsHead(PetscOptionsObject,"KSP Chebyshev Options");
304:   PetscOptionsInt("-ksp_chebyshev_esteig_steps","Number of est steps in Chebyshev","",cheb->eststeps,&cheb->eststeps,NULL);
305:   PetscOptionsRealArray("-ksp_chebyshev_eigenvalues","extreme eigenvalues","KSPChebyshevSetEigenvalues",eminmax,&neigarg,&flgeig);
306:   if (flgeig) {
307:     if (neigarg != 2) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
308:     KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]);
309:   }
310:   PetscOptionsRealArray("-ksp_chebyshev_esteig","estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds","KSPChebyshevEstEigSet",tform,&nestarg,&flgest);
311:   if (flgest) {
312:     switch (nestarg) {
313:     case 0:
314:       KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
315:       break;
316:     case 2:                     /* Base everything on the max eigenvalues */
317:       KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,tform[0],PETSC_DECIDE,tform[1]);
318:       break;
319:     case 4:                     /* Use the full 2x2 linear transformation */
320:       KSPChebyshevEstEigSet(ksp,tform[0],tform[1],tform[2],tform[3]);
321:       break;
322:     default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
323:     }
324:   }

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

331:   if (cheb->kspest) {
332:     PetscOptionsBool("-ksp_chebyshev_esteig_random","Use random right hand side for estimate","KSPChebyshevEstEigSetUseRandom",cheb->userandom,&cheb->userandom,NULL);
333:     if (cheb->userandom) {
334:       const char  *ksppre;
335:       if (!cheb->random) {
336:         PetscRandomCreate(PetscObjectComm((PetscObject)ksp),&cheb->random);
337:       }
338:       KSPGetOptionsPrefix(cheb->kspest, &ksppre);
339:       PetscObjectSetOptionsPrefix((PetscObject)cheb->random,ksppre);
340:       PetscRandomSetFromOptions(cheb->random);
341:     }
342:     KSPSetFromOptions(cheb->kspest);
343:   }
344:   PetscOptionsTail();
345:   return(0);
346: }

350: /*
351:  * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
352:  */
353: static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest,PetscReal *emin,PetscReal *emax)
354: {
356:   PetscInt       n,neig;
357:   PetscReal      *re,*im,min,max;

360:   KSPGetIterationNumber(kspest,&n);
361:   PetscMalloc2(n,&re,n,&im);
362:   KSPComputeEigenvalues(kspest,n,re,im,&neig);
363:   min  = PETSC_MAX_REAL;
364:   max  = PETSC_MIN_REAL;
365:   for (n=0; n<neig; n++) {
366:     min = PetscMin(min,re[n]);
367:     max = PetscMax(max,re[n]);
368:   }
369:   PetscFree2(re,im);
370:   *emax = max;
371:   *emin = min;
372:   return(0);
373: }

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

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

392:   PCGetOperators(ksp->pc,&Amat,&Pmat);
393:   if (cheb->kspest) {
394:     PetscObjectId    amatid,    pmatid;
395:     PetscObjectState amatstate, pmatstate;

397:     PetscObjectGetId((PetscObject)Amat,&amatid);
398:     PetscObjectGetId((PetscObject)Pmat,&pmatid);
399:     PetscObjectStateGet((PetscObject)Amat,&amatstate);
400:     PetscObjectStateGet((PetscObject)Pmat,&pmatstate);
401:     if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
402:       PetscReal          max=0.0,min=0.0;
403:       Vec                B;
404:       KSPConvergedReason reason;

406:       if (cheb->userandom) {
407:         B  = ksp->work[1];
408:         if (!cheb->random) {
409:           PetscRandomCreate(PetscObjectComm((PetscObject)B),&cheb->random);
410:         }
411:         VecSetRandom(B,cheb->random);
412:       } else {
413:         PC        pc;
414:         PetscBool change;

416:         KSPGetPC(cheb->kspest,&pc);
417:         PCPreSolveChangeRHS(pc,&change);
418:         if (change) {
419:           B = ksp->work[1];
420:           VecCopy(ksp->vec_rhs,B);
421:         } else {
422:           B = ksp->vec_rhs;
423:         }
424:       }
425:       KSPSolve(cheb->kspest,B,ksp->work[0]);

427:       KSPGetConvergedReason(cheb->kspest,&reason);
428:       if (reason < 0) {
429:         if (reason == KSP_DIVERGED_ITS) {
430:           PetscInfo(ksp,"Eigen estimator ran for prescribed number of iterations\n");
431:         } else {
432:           PetscInt its;
433:           KSPGetIterationNumber(cheb->kspest,&its);
434:           if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
435:             PetscInfo(ksp,"Eigen estimator KSP_DIVERGED_PCSETUP_FAILED\n");
436:             ksp->reason = KSP_DIVERGED_PCSETUP_FAILED;
437:             VecSetInf(ksp->vec_sol);
438:           } else {
439:             SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Eigen estimator failed: %s at iteration %D",KSPConvergedReasons[reason],its);
440:           }
441:         }
442:       } else if (reason==KSP_CONVERGED_RTOL ||reason==KSP_CONVERGED_ATOL) {
443:         PetscInfo(ksp,"Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n");
444:       } else {
445:         PetscInfo1(ksp,"Eigen estimator did not converge by iteration: %s\n",KSPConvergedReasons[reason]);
446:       }

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

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

453:       cheb->amatid    = amatid;
454:       cheb->pmatid    = pmatid;
455:       cheb->amatstate = amatstate;
456:       cheb->pmatstate = pmatstate;
457:     }
458:   }

460:   ksp->its = 0;
461:   maxit    = ksp->max_it;

463:   /* These three point to the three active solutions, we
464:      rotate these three at each solution update */
465:   km1      = 0; k = 1; kp1 = 2;
466:   sol_orig = ksp->vec_sol; /* ksp->vec_sol will be asigned to rotating vector p[k], thus save its address */
467:   b        = ksp->vec_rhs;
468:   p[km1]   = sol_orig;
469:   p[k]     = ksp->work[0];
470:   p[kp1]   = ksp->work[1];
471:   r        = ksp->work[2];

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

476:   /*   -alpha <=  scale*lambda(B^{-1}A) <= alpha   */
477:   alpha     = 1.0 - scale*(cheb->emin);
478:   Gamma     = 1.0;
479:   mu        = 1.0/alpha;
480:   omegaprod = 2.0/alpha;

482:   c[km1] = 1.0;
483:   c[k]   = mu;

485:   if (!ksp->guess_zero) {
486:     KSP_MatMult(ksp,Amat,p[km1],r);     /*  r = b - A*p[km1] */
487:     VecAYPX(r,-1.0,b);
488:   } else {
489:     VecCopy(b,r);
490:   }

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

495:   for (i=0; i<maxit; i++) {
496:     PetscObjectSAWsTakeAccess((PetscObject)ksp);

498:     ksp->its++;
499:     PetscObjectSAWsGrantAccess((PetscObject)ksp);
500:     c[kp1] = 2.0*mu*c[k] - c[km1];
501:     omega  = omegaprod*c[k]/c[kp1];

503:     KSP_MatMult(ksp,Amat,p[k],r);          /*  r = b - Ap[k]    */
504:     VecAYPX(r,-1.0,b);
505:     KSP_PCApply(ksp,r,p[kp1]);             /*  p[kp1] = B^{-1}r  */
506:     ksp->vec_sol = p[k];

508:     /* calculate residual norm if requested */
509:     if (ksp->normtype != KSP_NORM_NONE || ksp->numbermonitors) {
510:       if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
511:         VecNorm(r,NORM_2,&rnorm);
512:       } else {
513:         VecNorm(p[kp1],NORM_2,&rnorm);
514:       }
515:       PetscObjectSAWsTakeAccess((PetscObject)ksp);
516:       ksp->rnorm   = rnorm;
517:       PetscObjectSAWsGrantAccess((PetscObject)ksp);
518:       KSPLogResidualHistory(ksp,rnorm);
519:       KSPMonitor(ksp,i,rnorm);
520:       (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
521:       if (ksp->reason) break;
522:     }

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

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

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

567: static  PetscErrorCode KSPView_Chebyshev(KSP ksp,PetscViewer viewer)
568: {
569:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;
571:   PetscBool      iascii;

574:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
575:   if (iascii) {
576:     PetscViewerASCIIPrintf(viewer,"  Chebyshev: eigenvalue estimates:  min = %g, max = %g\n",(double)cheb->emin,(double)cheb->emax);
577:     if (cheb->kspest) {
578:       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]);
579:       PetscViewerASCIIPushTab(viewer);
580:       KSPView(cheb->kspest,viewer);
581:       PetscViewerASCIIPopTab(viewer);
582:       if (cheb->userandom) {
583:         PetscViewerASCIIPrintf(viewer,"  Chebyshev: estimating eigenvalues using random right hand side\n");
584:         PetscViewerASCIIPushTab(viewer);
585:         PetscRandomView(cheb->random,viewer);
586:         PetscViewerASCIIPopTab(viewer);
587:       }
588:     }
589:   }
590:   return(0);
591: }

595: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
596: {
597:   KSP_Chebyshev  *cheb = (KSP_Chebyshev*)ksp->data;

601:   KSPDestroy(&cheb->kspest);
602:   PetscRandomDestroy(&cheb->random);
603:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",NULL);
604:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",NULL);
605:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetRandom_C",NULL);
606:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",NULL);
607:   KSPDestroyDefault(ksp);
608:   return(0);
609: }

611: /*MC
612:      KSPCHEBYSHEV - The preconditioned Chebyshev iterative method

614:    Options Database Keys:
615: +   -ksp_chebyshev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
616:                   of the preconditioned operator. If these are accurate you will get much faster convergence.
617: .   -ksp_chebyshev_esteig <a,b,c,d> - estimate eigenvalues using a Krylov method, then use this
618:                          transform for Chebyshev eigenvalue bounds (KSPChebyshevEstEigSet())
619: .   -ksp_chebyshev_esteig_steps - number of estimation steps
620: -   -ksp_chebyshev_esteig_random - use random number generator to create right hand side for eigenvalue estimator

622:    Level: beginner

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

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

631: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
632:            KSPChebyshevSetEigenvalues(), KSPChebyshevEstEigSet(), KSPChebyshevEstEigSetUseRandom(), KSPChebyshevEstEigSetRandom(),
633:            KSPRICHARDSON, KSPCG, PCMG

635: M*/

639: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
640: {
642:   KSP_Chebyshev  *chebyshevP;

645:   PetscNewLog(ksp,&chebyshevP);

647:   ksp->data = (void*)chebyshevP;
648:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
649:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);

651:   chebyshevP->emin = 0.;
652:   chebyshevP->emax = 0.;

654:   chebyshevP->tform[0] = 0.0;
655:   chebyshevP->tform[1] = 0.1;
656:   chebyshevP->tform[2] = 0;
657:   chebyshevP->tform[3] = 1.1;
658:   chebyshevP->eststeps = 10;
659:   chebyshevP->userandom = PETSC_FALSE;

661:   ksp->ops->setup          = KSPSetUp_Chebyshev;
662:   ksp->ops->solve          = KSPSolve_Chebyshev;
663:   ksp->ops->destroy        = KSPDestroy_Chebyshev;
664:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
665:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
666:   ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
667:   ksp->ops->view           = KSPView_Chebyshev;
668:   ksp->ops->reset          = KSPReset_Chebyshev;

670:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",KSPChebyshevSetEigenvalues_Chebyshev);
671:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",KSPChebyshevEstEigSet_Chebyshev);
672:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetRandom_C",KSPChebyshevEstEigSetRandom_Chebyshev);
673:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseRandom_C",KSPChebyshevEstEigSetUseRandom_Chebyshev);
674:   PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",KSPChebyshevEstEigGetKSP_Chebyshev);
675:   return(0);
676: }