Actual source code: cheby.c
petsc-3.7.5 2017-01-01
2: #include <petsc/private/kspimpl.h> /*I "petscksp.h" I*/
3: #include <../src/ksp/ksp/impls/cheby/chebyshevimpl.h>
7: static PetscErrorCode KSPReset_Chebyshev(KSP ksp)
8: {
9: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
13: KSPReset(cheb->kspest);
14: return(0);
15: }
19: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
20: {
21: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
25: KSPSetWorkVecs(ksp,3);
26: if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) { /* We need to estimate eigenvalues */
27: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
28: }
29: return(0);
30: }
34: static PetscErrorCode KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp,PetscReal emax,PetscReal emin)
35: {
36: KSP_Chebyshev *chebyshevP = (KSP_Chebyshev*)ksp->data;
40: 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);
41: 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);
42: chebyshevP->emax = emax;
43: chebyshevP->emin = emin;
45: KSPChebyshevEstEigSet(ksp,0.,0.,0.,0.); /* Destroy any estimation setup */
46: return(0);
47: }
51: static PetscErrorCode KSPChebyshevEstEigSet_Chebyshev(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
52: {
53: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
57: if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
58: if (!cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
59: KSPCreate(PetscObjectComm((PetscObject)ksp),&cheb->kspest);
60: PetscObjectIncrementTabLevel((PetscObject)cheb->kspest,(PetscObject)ksp,1);
61: KSPSetOptionsPrefix(cheb->kspest,((PetscObject)ksp)->prefix);
62: KSPAppendOptionsPrefix(cheb->kspest,"esteig_");
63: KSPSetSkipPCSetFromOptions(cheb->kspest,PETSC_TRUE);
65: KSPSetPC(cheb->kspest,ksp->pc);
67: KSPSetComputeEigenvalues(cheb->kspest,PETSC_TRUE);
69: /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
70: KSPSetTolerances(cheb->kspest,1.e-12,PETSC_DEFAULT,PETSC_DEFAULT,cheb->eststeps);
71: }
72: if (a >= 0) cheb->tform[0] = a;
73: if (b >= 0) cheb->tform[1] = b;
74: if (c >= 0) cheb->tform[2] = c;
75: if (d >= 0) cheb->tform[3] = d;
76: cheb->amatid = 0;
77: cheb->pmatid = 0;
78: cheb->amatstate = -1;
79: cheb->pmatstate = -1;
80: } else {
81: KSPDestroy(&cheb->kspest);
82: }
83: return(0);
84: }
88: static PetscErrorCode KSPChebyshevEstEigSetRandom_Chebyshev(KSP ksp,PetscRandom random)
89: {
90: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
94: if (random) {PetscObjectReference((PetscObject)random);}
95: PetscRandomDestroy(&cheb->random);
96: cheb->random = random;
97: return(0);
98: }
102: static PetscErrorCode KSPChebyshevEstEigSetUseRandom_Chebyshev(KSP ksp,PetscBool use)
103: {
104: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
107: cheb->userandom = use;
108: return(0);
109: }
113: /*@
114: KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues
115: of the preconditioned problem.
117: Logically Collective on KSP
119: Input Parameters:
120: + ksp - the Krylov space context
121: - emax, emin - the eigenvalue estimates
123: Options Database:
124: . -ksp_chebyshev_eigenvalues emin,emax
126: Note: Call KSPChebyshevEstEigSet() or use the option -ksp_chebyshev_esteig a,b,c,d to have the KSP
127: estimate the eigenvalues and use these estimated values automatically
129: Level: intermediate
131: .keywords: KSP, Chebyshev, set, eigenvalues
132: @*/
133: PetscErrorCode KSPChebyshevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
134: {
141: PetscTryMethod(ksp,"KSPChebyshevSetEigenvalues_C",(KSP,PetscReal,PetscReal),(ksp,emax,emin));
142: return(0);
143: }
147: /*@
148: KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev
150: Logically Collective on KSP
152: Input Parameters:
153: + ksp - the Krylov space context
154: . a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
155: . b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
156: . c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
157: - d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
159: Options Database:
160: . -ksp_chebyshev_esteig a,b,c,d
162: Notes:
163: The Chebyshev bounds are set using
164: .vb
165: minbound = a*minest + b*maxest
166: maxbound = c*minest + d*maxest
167: .ve
168: The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
169: The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.
171: If 0.0 is passed for all transform arguments (a,b,c,d), eigenvalue estimation is disabled.
173: The default transform is (0,0.1; 0,1.1) which targets the "upper" part of the spectrum, as desirable for use with multigrid.
175: The eigenvalues are estimated using the Lanczo (KSPCG) or Arnoldi (KSPGMRES) process using a random right hand side vector.
177: Level: intermediate
179: .keywords: KSP, Chebyshev, set, eigenvalues, PCMG
180: @*/
181: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
182: {
191: PetscTryMethod(ksp,"KSPChebyshevEstEigSet_C",(KSP,PetscReal,PetscReal,PetscReal,PetscReal),(ksp,a,b,c,d));
192: return(0);
193: }
197: /*@
198: KSPChebyshevEstEigSetUseRandom - use a random right hand side in order to do the estimate instead of the given right hand side
200: Logically Collective
202: Input Arguments:
203: + ksp - linear solver context
204: - use - PETSC_TRUE to use random
206: Options Database:
207: + -ksp_chebyshev_esteig_random <true,false>
209: Notes: This alledgely works better for multigrid smoothers
211: Use KSPChebyshevEstEigSetRandom() to provide the random number generator to be used. Otherwise it creates a default
213: Level: intermediate
215: .seealso: KSPChebyshevEstEigSet(), PetscRandomCreate(), KSPChebyshevEstEigSetRandom()
216: @*/
217: PetscErrorCode KSPChebyshevEstEigSetUseRandom(KSP ksp,PetscBool use)
218: {
223: PetscTryMethod(ksp,"KSPChebyshevEstEigSetUseRandom_C",(KSP,PetscBool),(ksp,use));
224: return(0);
225: }
229: /*@
230: KSPChebyshevEstEigSetRandom - set random context for estimating eigenvalues
232: Logically Collective
234: Input Arguments:
235: + ksp - linear solver context
236: - random - random number context or NULL to use default
238: Level: intermediate
240: .seealso: KSPChebyshevEstEigSet(), PetscRandomCreate()
241: @*/
242: PetscErrorCode KSPChebyshevEstEigSetRandom(KSP ksp,PetscRandom random)
243: {
249: PetscTryMethod(ksp,"KSPChebyshevEstEigSetRandom_C",(KSP,PetscRandom),(ksp,random));
250: return(0);
251: }
255: /*@
256: KSPChebyshevEstEigGetKSP - Get the Krylov method context used to estimate eigenvalues for the Chebyshev method. If
257: a Krylov method is not being used for this purpose, NULL is returned. The reference count of the returned KSP is
258: not incremented: it should not be destroyed by the user.
260: Input Parameters:
261: . ksp - the Krylov space context
263: Output Parameters:
264: . kspest the eigenvalue estimation Krylov space context
266: Level: intermediate
268: .seealso: KSPChebyshevEstEigSet()
269: @*/
270: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
271: {
276: *kspest = NULL;
277: PetscTryMethod(ksp,"KSPChebyshevEstEigGetKSP_C",(KSP,KSP*),(ksp,kspest));
278: return(0);
279: }
283: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
284: {
285: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
288: *kspest = cheb->kspest;
289: return(0);
290: }
294: static PetscErrorCode KSPSetFromOptions_Chebyshev(PetscOptionItems *PetscOptionsObject,KSP ksp)
295: {
296: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
298: PetscInt neigarg = 2, nestarg = 4;
299: PetscReal eminmax[2] = {0., 0.};
300: PetscReal tform[4] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
301: PetscBool flgeig, flgest;
304: PetscOptionsHead(PetscOptionsObject,"KSP Chebyshev Options");
305: PetscOptionsInt("-ksp_chebyshev_esteig_steps","Number of est steps in Chebyshev","",cheb->eststeps,&cheb->eststeps,NULL);
306: PetscOptionsRealArray("-ksp_chebyshev_eigenvalues","extreme eigenvalues","KSPChebyshevSetEigenvalues",eminmax,&neigarg,&flgeig);
307: if (flgeig) {
308: if (neigarg != 2) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
309: KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]);
310: }
311: PetscOptionsRealArray("-ksp_chebyshev_esteig","estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds","KSPChebyshevEstEigSet",tform,&nestarg,&flgest);
312: if (flgest) {
313: switch (nestarg) {
314: case 0:
315: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
316: break;
317: case 2: /* Base everything on the max eigenvalues */
318: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,tform[0],PETSC_DECIDE,tform[1]);
319: break;
320: case 4: /* Use the full 2x2 linear transformation */
321: KSPChebyshevEstEigSet(ksp,tform[0],tform[1],tform[2],tform[3]);
322: break;
323: default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
324: }
325: }
327: /* We need to estimate eigenvalues; need to set this here so that KSPSetFromOptions() is called on the estimator */
328: if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) {
329: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
330: }
332: if (cheb->kspest) {
333: PetscOptionsBool("-ksp_chebyshev_esteig_random","Use random right hand side for estimate","KSPChebyshevEstEigSetUseRandom",cheb->userandom,&cheb->userandom,NULL);
334: if (cheb->userandom) {
335: const char *ksppre;
336: if (!cheb->random) {
337: PetscRandomCreate(PetscObjectComm((PetscObject)ksp),&cheb->random);
338: }
339: KSPGetOptionsPrefix(cheb->kspest, &ksppre);
340: PetscObjectSetOptionsPrefix((PetscObject)cheb->random,ksppre);
341: PetscRandomSetFromOptions(cheb->random);
342: }
343: KSPSetFromOptions(cheb->kspest);
344: }
345: PetscOptionsTail();
346: return(0);
347: }
351: /*
352: * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
353: */
354: static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest,PetscReal *emin,PetscReal *emax)
355: {
357: PetscInt n,neig;
358: PetscReal *re,*im,min,max;
361: KSPGetIterationNumber(kspest,&n);
362: PetscMalloc2(n,&re,n,&im);
363: KSPComputeEigenvalues(kspest,n,re,im,&neig);
364: min = PETSC_MAX_REAL;
365: max = PETSC_MIN_REAL;
366: for (n=0; n<neig; n++) {
367: min = PetscMin(min,re[n]);
368: max = PetscMax(max,re[n]);
369: }
370: PetscFree2(re,im);
371: *emax = max;
372: *emin = min;
373: return(0);
374: }
378: static PetscErrorCode KSPSolve_Chebyshev(KSP ksp)
379: {
380: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
382: PetscInt k,kp1,km1,maxit,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: if (cheb->kspest) {
395: PetscObjectId amatid, pmatid;
396: PetscObjectState amatstate, pmatstate;
398: PetscObjectGetId((PetscObject)Amat,&amatid);
399: PetscObjectGetId((PetscObject)Pmat,&pmatid);
400: PetscObjectStateGet((PetscObject)Amat,&amatstate);
401: PetscObjectStateGet((PetscObject)Pmat,&pmatstate);
402: if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
403: PetscReal max=0.0,min=0.0;
404: Vec B;
405: KSPConvergedReason reason;
407: if (cheb->userandom) {
408: B = ksp->work[1];
409: if (!cheb->random) {
410: PetscRandomCreate(PetscObjectComm((PetscObject)B),&cheb->random);
411: }
412: VecSetRandom(B,cheb->random);
413: } else {
414: PC pc;
415: PetscBool change;
417: KSPGetPC(cheb->kspest,&pc);
418: PCPreSolveChangeRHS(pc,&change);
419: if (change) {
420: B = ksp->work[1];
421: VecCopy(ksp->vec_rhs,B);
422: } else {
423: B = ksp->vec_rhs;
424: }
425: }
426: KSPSolve(cheb->kspest,B,ksp->work[0]);
428: KSPGetConvergedReason(cheb->kspest,&reason);
429: if (reason < 0) {
430: if (reason == KSP_DIVERGED_ITS) {
431: PetscInfo(ksp,"Eigen estimator ran for prescribed number of iterations\n");
432: } else {
433: PetscInt its;
434: KSPGetIterationNumber(cheb->kspest,&its);
435: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
436: PetscInfo(ksp,"Eigen estimator KSP_DIVERGED_PCSETUP_FAILED\n");
437: ksp->reason = KSP_DIVERGED_PCSETUP_FAILED;
438: VecSetInf(ksp->vec_sol);
439: } else {
440: SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Eigen estimator failed: %s at iteration %D",KSPConvergedReasons[reason],its);
441: }
442: }
443: } else if (reason==KSP_CONVERGED_RTOL ||reason==KSP_CONVERGED_ATOL) {
444: PetscInfo(ksp,"Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n");
445: } else {
446: PetscInfo1(ksp,"Eigen estimator did not converge by iteration: %s\n",KSPConvergedReasons[reason]);
447: }
449: KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest,&min,&max);
451: cheb->emin = cheb->tform[0]*min + cheb->tform[1]*max;
452: cheb->emax = cheb->tform[2]*min + cheb->tform[3]*max;
454: cheb->amatid = amatid;
455: cheb->pmatid = pmatid;
456: cheb->amatstate = amatstate;
457: cheb->pmatstate = pmatstate;
458: }
459: }
461: ksp->its = 0;
462: maxit = ksp->max_it;
464: /* These three point to the three active solutions, we
465: rotate these three at each solution update */
466: km1 = 0; k = 1; kp1 = 2;
467: sol_orig = ksp->vec_sol; /* ksp->vec_sol will be asigned to rotating vector p[k], thus save its address */
468: b = ksp->vec_rhs;
469: p[km1] = sol_orig;
470: p[k] = ksp->work[0];
471: p[kp1] = ksp->work[1];
472: r = ksp->work[2];
474: /* use scale*B as our preconditioner */
475: scale = 2.0/(cheb->emax + cheb->emin);
477: /* -alpha <= scale*lambda(B^{-1}A) <= alpha */
478: alpha = 1.0 - scale*(cheb->emin);
479: Gamma = 1.0;
480: mu = 1.0/alpha;
481: omegaprod = 2.0/alpha;
483: c[km1] = 1.0;
484: c[k] = mu;
486: if (!ksp->guess_zero) {
487: KSP_MatMult(ksp,Amat,p[km1],r); /* r = b - A*p[km1] */
488: VecAYPX(r,-1.0,b);
489: } else {
490: VecCopy(b,r);
491: }
493: KSP_PCApply(ksp,r,p[k]); /* p[k] = scale B^{-1}r + p[km1] */
494: VecAYPX(p[k],scale,p[km1]);
496: for (i=0; i<maxit; i++) {
497: PetscObjectSAWsTakeAccess((PetscObject)ksp);
499: ksp->its++;
500: PetscObjectSAWsGrantAccess((PetscObject)ksp);
501: c[kp1] = 2.0*mu*c[k] - c[km1];
502: omega = omegaprod*c[k]/c[kp1];
504: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
505: VecAYPX(r,-1.0,b);
506: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
507: ksp->vec_sol = p[k];
509: /* calculate residual norm if requested */
510: if (ksp->normtype != KSP_NORM_NONE || ksp->numbermonitors) {
511: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
512: VecNorm(r,NORM_2,&rnorm);
513: } else {
514: VecNorm(p[kp1],NORM_2,&rnorm);
515: }
516: PetscObjectSAWsTakeAccess((PetscObject)ksp);
517: ksp->rnorm = rnorm;
518: PetscObjectSAWsGrantAccess((PetscObject)ksp);
519: KSPLogResidualHistory(ksp,rnorm);
520: KSPMonitor(ksp,i,rnorm);
521: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
522: if (ksp->reason) break;
523: }
525: /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
526: VecAXPBYPCZ(p[kp1],1.0-omega,omega,omega*Gamma*scale,p[km1],p[k]);
528: ktmp = km1;
529: km1 = k;
530: k = kp1;
531: kp1 = ktmp;
532: }
533: if (!ksp->reason) {
534: if (ksp->normtype != KSP_NORM_NONE) {
535: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
536: VecAYPX(r,-1.0,b);
537: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
538: VecNorm(r,NORM_2,&rnorm);
539: } else {
540: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
541: VecNorm(p[kp1],NORM_2,&rnorm);
542: }
543: PetscObjectSAWsTakeAccess((PetscObject)ksp);
544: ksp->rnorm = rnorm;
545: PetscObjectSAWsGrantAccess((PetscObject)ksp);
546: ksp->vec_sol = p[k];
547: KSPLogResidualHistory(ksp,rnorm);
548: KSPMonitor(ksp,i,rnorm);
549: }
550: if (ksp->its >= ksp->max_it) {
551: if (ksp->normtype != KSP_NORM_NONE) {
552: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
553: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
554: } else ksp->reason = KSP_CONVERGED_ITS;
555: }
556: }
558: /* make sure solution is in vector x */
559: ksp->vec_sol = sol_orig;
560: if (k) {
561: VecCopy(p[k],sol_orig);
562: }
563: return(0);
564: }
568: static PetscErrorCode KSPView_Chebyshev(KSP ksp,PetscViewer viewer)
569: {
570: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
572: PetscBool iascii;
575: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
576: if (iascii) {
577: PetscViewerASCIIPrintf(viewer," Chebyshev: eigenvalue estimates: min = %g, max = %g\n",(double)cheb->emin,(double)cheb->emax);
578: if (cheb->kspest) {
579: 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]);
580: PetscViewerASCIIPushTab(viewer);
581: KSPView(cheb->kspest,viewer);
582: PetscViewerASCIIPopTab(viewer);
583: if (cheb->userandom) {
584: PetscViewerASCIIPrintf(viewer," Chebyshev: estimating eigenvalues using random right hand side\n");
585: PetscViewerASCIIPushTab(viewer);
586: PetscRandomView(cheb->random,viewer);
587: PetscViewerASCIIPopTab(viewer);
588: }
589: }
590: }
591: return(0);
592: }
596: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
597: {
598: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
602: KSPDestroy(&cheb->kspest);
603: PetscRandomDestroy(&cheb->random);
604: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",NULL);
605: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",NULL);
606: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetRandom_C",NULL);
607: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",NULL);
608: KSPDestroyDefault(ksp);
609: return(0);
610: }
612: /*MC
613: KSPCHEBYSHEV - The preconditioned Chebyshev iterative method
615: Options Database Keys:
616: + -ksp_chebyshev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
617: of the preconditioned operator. If these are accurate you will get much faster convergence.
618: . -ksp_chebyshev_esteig <a,b,c,d> - estimate eigenvalues using a Krylov method, then use this
619: transform for Chebyshev eigenvalue bounds (KSPChebyshevEstEigSet())
620: . -ksp_chebyshev_esteig_steps - number of estimation steps
621: - -ksp_chebyshev_esteig_random - use random number generator to create right hand side for eigenvalue estimator
623: Level: beginner
625: Notes: The Chebyshev method requires both the matrix and preconditioner to
626: be symmetric positive (semi) definite.
627: Only support for left preconditioning.
629: Chebyshev is configured as a smoother by default, targetting the "upper" part of the spectrum.
630: The user should call KSPChebyshevSetEigenvalues() if they have eigenvalue estimates.
632: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
633: KSPChebyshevSetEigenvalues(), KSPChebyshevEstEigSet(), KSPChebyshevEstEigSetUseRandom(), KSPChebyshevEstEigSetRandom(),
634: KSPRICHARDSON, KSPCG, PCMG
636: M*/
640: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
641: {
643: KSP_Chebyshev *chebyshevP;
646: PetscNewLog(ksp,&chebyshevP);
648: ksp->data = (void*)chebyshevP;
649: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
650: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
652: chebyshevP->emin = 0.;
653: chebyshevP->emax = 0.;
655: chebyshevP->tform[0] = 0.0;
656: chebyshevP->tform[1] = 0.1;
657: chebyshevP->tform[2] = 0;
658: chebyshevP->tform[3] = 1.1;
659: chebyshevP->eststeps = 10;
660: chebyshevP->userandom = PETSC_FALSE;
662: ksp->ops->setup = KSPSetUp_Chebyshev;
663: ksp->ops->solve = KSPSolve_Chebyshev;
664: ksp->ops->destroy = KSPDestroy_Chebyshev;
665: ksp->ops->buildsolution = KSPBuildSolutionDefault;
666: ksp->ops->buildresidual = KSPBuildResidualDefault;
667: ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
668: ksp->ops->view = KSPView_Chebyshev;
669: ksp->ops->reset = KSPReset_Chebyshev;
671: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",KSPChebyshevSetEigenvalues_Chebyshev);
672: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",KSPChebyshevEstEigSet_Chebyshev);
673: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetRandom_C",KSPChebyshevEstEigSetRandom_Chebyshev);
674: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseRandom_C",KSPChebyshevEstEigSetUseRandom_Chebyshev);
675: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",KSPChebyshevEstEigGetKSP_Chebyshev);
676: return(0);
677: }