Actual source code: cheby.c
petsc-3.4.5 2014-06-29
2: #include <petsc-private/kspimpl.h> /*I "petscksp.h" I*/
3: #include <../src/ksp/ksp/impls/cheby/chebyshevimpl.h>
7: PetscErrorCode KSPReset_Chebyshev(KSP ksp)
8: {
9: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
13: KSPReset(cheb->kspest);
14: return(0);
15: }
19: 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.) { /* We need to estimate eigenvalues */
27: KSPChebyshevSetEstimateEigenvalues(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
28: /* Enable runtime options for cheb->kspest:
29: KSPChebyshevSetEstimateEigenvalues() creates cheb->kspest, but does not call KSPSetFromOptions(cheb->kspest)! */
30: KSPSetFromOptions(cheb->kspest);
31: } else if (cheb->hybrid && !cheb->kspest) { /* We need to create cheb->kspest */
32: PetscBool nonzero;
33: KSPCreate(PetscObjectComm((PetscObject)ksp),&cheb->kspest);
34: PetscObjectIncrementTabLevel((PetscObject)cheb->kspest,(PetscObject)ksp,1);
35: KSPSetOptionsPrefix(cheb->kspest,((PetscObject)ksp)->prefix);
36: KSPAppendOptionsPrefix(cheb->kspest,"est_");
38: KSPGetPC(cheb->kspest,&cheb->pcnone);
39: PetscObjectReference((PetscObject)cheb->pcnone);
40: PCSetType(cheb->pcnone,PCNONE);
41: KSPSetPC(cheb->kspest,ksp->pc);
42:
43: KSPGetInitialGuessNonzero(ksp,&nonzero);
44: KSPSetInitialGuessNonzero(cheb->kspest,nonzero);
45: KSPSetComputeEigenvalues(cheb->kspest,PETSC_TRUE);
47: /* Estimate with a fixed number of iterations */
48: KSPSetConvergenceTest(cheb->kspest,KSPSkipConverged,0,0);
49: KSPSetNormType(cheb->kspest,KSP_NORM_NONE);
50: KSPSetTolerances(cheb->kspest,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,cheb->eststeps);
52: /* Enable runtime options for cheb->kspest */
53: KSPSetFromOptions(cheb->kspest);
54: }
55: return(0);
56: }
60: static PetscErrorCode KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp,PetscReal emax,PetscReal emin)
61: {
62: KSP_Chebyshev *chebyshevP = (KSP_Chebyshev*)ksp->data;
66: if (emax <= emin) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Maximum eigenvalue must be larger than minimum: max %g min %G",emax,emin);
67: 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",emax,emin);
68: chebyshevP->emax = emax;
69: chebyshevP->emin = emin;
71: KSPChebyshevSetEstimateEigenvalues(ksp,0.,0.,0.,0.); /* Destroy any estimation setup */
72: return(0);
73: }
77: static PetscErrorCode KSPChebyshevSetEstimateEigenvalues_Chebyshev(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
78: {
79: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
83: if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
84: if (!cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
85: PetscBool nonzero;
87: KSPCreate(PetscObjectComm((PetscObject)ksp),&cheb->kspest);
88: PetscObjectIncrementTabLevel((PetscObject)cheb->kspest,(PetscObject)ksp,1);
89: KSPSetOptionsPrefix(cheb->kspest,((PetscObject)ksp)->prefix);
90: KSPAppendOptionsPrefix(cheb->kspest,"est_");
92: KSPGetPC(cheb->kspest,&cheb->pcnone);
93: PetscObjectReference((PetscObject)cheb->pcnone);
94: PCSetType(cheb->pcnone,PCNONE);
95: KSPSetPC(cheb->kspest,ksp->pc);
96:
97: KSPGetInitialGuessNonzero(ksp,&nonzero);
98: KSPSetInitialGuessNonzero(cheb->kspest,nonzero);
99: KSPSetComputeEigenvalues(cheb->kspest,PETSC_TRUE);
101: /* Estimate with a fixed number of iterations */
102: KSPSetConvergenceTest(cheb->kspest,KSPSkipConverged,0,0);
103: KSPSetNormType(cheb->kspest,KSP_NORM_NONE);
104: KSPSetTolerances(cheb->kspest,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,cheb->eststeps);
105: }
106: if (a >= 0) cheb->tform[0] = a;
107: if (b >= 0) cheb->tform[1] = b;
108: if (c >= 0) cheb->tform[2] = c;
109: if (d >= 0) cheb->tform[3] = d;
110: cheb->estimate_current = PETSC_FALSE;
111: } else {
112: KSPDestroy(&cheb->kspest);
113: PCDestroy(&cheb->pcnone);
114: }
115: return(0);
116: }
120: static PetscErrorCode KSPChebyshevEstEigSetRandom_Chebyshev(KSP ksp,PetscRandom random)
121: {
122: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
126: if (random) {PetscObjectReference((PetscObject)random);}
127: PetscRandomDestroy(&cheb->random);
129: cheb->random = random;
130: return(0);
131: }
135: static PetscErrorCode KSPChebyshevSetNewMatrix_Chebyshev(KSP ksp)
136: {
137: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
140: cheb->estimate_current = PETSC_FALSE;
141: return(0);
142: }
146: /*@
147: KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues
148: of the preconditioned problem.
150: Logically Collective on KSP
152: Input Parameters:
153: + ksp - the Krylov space context
154: - emax, emin - the eigenvalue estimates
156: Options Database:
157: . -ksp_chebyshev_eigenvalues emin,emax
159: Note: If you run with the Krylov method of KSPCG with the option -ksp_monitor_singular_value it will
160: for that given matrix and preconditioner an estimate of the extreme eigenvalues.
162: Level: intermediate
164: .keywords: KSP, Chebyshev, set, eigenvalues
165: @*/
166: PetscErrorCode KSPChebyshevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
167: {
174: PetscTryMethod(ksp,"KSPChebyshevSetEigenvalues_C",(KSP,PetscReal,PetscReal),(ksp,emax,emin));
175: return(0);
176: }
180: /*@
181: KSPChebyshevSetEstimateEigenvalues - Automatically estimate the eigenvalues to use for Chebyshev
183: Logically Collective on KSP
185: Input Parameters:
186: + ksp - the Krylov space context
187: . a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
188: . b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
189: . c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
190: - d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
192: Options Database:
193: . -ksp_chebyshev_estimate_eigenvalues a,b,c,d
195: Notes:
196: The Chebyshev bounds are estimated using
197: .vb
198: minbound = a*minest + b*maxest
199: maxbound = c*minest + d*maxest
200: .ve
201: The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
202: The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.
204: If 0.0 is passed for all transform arguments (a,b,c,d), eigenvalue estimation is disabled.
206: The default transform is (0,0.1; 0,1.1) which targets the "upper" part of the spectrum, as desirable for use with multigrid.
208: Level: intermediate
210: .keywords: KSP, Chebyshev, set, eigenvalues, PCMG
211: @*/
212: PetscErrorCode KSPChebyshevSetEstimateEigenvalues(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
213: {
222: PetscTryMethod(ksp,"KSPChebyshevSetEstimateEigenvalues_C",(KSP,PetscReal,PetscReal,PetscReal,PetscReal),(ksp,a,b,c,d));
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 disable randomized RHS
237: Level: intermediate
239: .seealso: KSPChebyshevSetEstimateEigenvalues(), 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: KSPChebyshevSetNewMatrix - Indicates that the matrix has changed, causes eigenvalue estimates to be recomputed if appropriate.
257: Logically Collective on KSP
259: Input Parameter:
260: . ksp - the Krylov space context
262: Level: developer
264: .keywords: KSP, Chebyshev, set, eigenvalues
266: .seealso: KSPChebyshevSetEstimateEigenvalues()
267: @*/
268: PetscErrorCode KSPChebyshevSetNewMatrix(KSP ksp)
269: {
274: PetscTryMethod(ksp,"KSPChebyshevSetNewMatrix_C",(KSP),(ksp));
275: return(0);
276: }
280: PetscErrorCode KSPSetFromOptions_Chebyshev(KSP ksp)
281: {
282: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
284: PetscInt two = 2,four = 4;
285: PetscReal tform[4] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
286: PetscBool flg;
289: PetscOptionsHead("KSP Chebyshev Options");
291: /*
292: Use hybrid Chebyshev.
293: Ref: "A hybrid Chebyshev Krylov-subspace algorithm for solving nonsymmetric systems of linear equations",
294: Howard Elman and Y. Saad and P. E. Saylor, SIAM Journal on Scientific and Statistical Computing, 1986.
295: */
296: PetscOptionsBool("-ksp_chebyshev_hybrid","Use hybrid Chebyshev","",cheb->hybrid,&cheb->hybrid,NULL);
297: PetscOptionsInt("-ksp_chebyshev_hybrid_chebysteps","Number of Chebyshev steps in hybrid Chebyshev","",cheb->chebysteps,&cheb->chebysteps,NULL);
298: PetscOptionsInt("-ksp_chebyshev_hybrid_eststeps","Number of adaptive/est steps in hybrid Chebyshev","",cheb->eststeps,&cheb->eststeps,NULL);
299: PetscOptionsBool("-ksp_chebyshev_hybrid_purification","Use purification in hybrid Chebyshev","",cheb->purification,&cheb->purification,NULL);
301: PetscOptionsRealArray("-ksp_chebyshev_eigenvalues","extreme eigenvalues","KSPChebyshevSetEigenvalues",&cheb->emin,&two,0);
302: PetscOptionsRealArray("-ksp_chebyshev_estimate_eigenvalues","estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds","KSPChebyshevSetEstimateEigenvalues",tform,&four,&flg);
303: if (flg) {
304: switch (four) {
305: case 0:
306: KSPChebyshevSetEstimateEigenvalues(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
307: break;
308: case 2: /* Base everything on the max eigenvalues */
309: KSPChebyshevSetEstimateEigenvalues(ksp,PETSC_DECIDE,tform[0],PETSC_DECIDE,tform[1]);
310: break;
311: case 4: /* Use the full 2x2 linear transformation */
312: KSPChebyshevSetEstimateEigenvalues(ksp,tform[0],tform[1],tform[2],tform[3]);
313: break;
314: default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
315: }
316: }
318: if (cheb->kspest) {
319: PetscBool estrand = PETSC_FALSE;
320: PetscOptionsBool("-ksp_chebyshev_estimate_eigenvalues_random","Use Random right hand side for eigenvalue estimation","KSPChebyshevEstEigSetRandom",estrand,&estrand,NULL);
321: if (estrand) {
322: PetscRandom random;
323: PetscRandomCreate(PetscObjectComm((PetscObject)ksp),&random);
324: PetscObjectSetOptionsPrefix((PetscObject)random,((PetscObject)ksp)->prefix);
325: PetscObjectAppendOptionsPrefix((PetscObject)random,"ksp_chebyshev_estimate_eigenvalues_");
326: PetscRandomSetFromOptions(random);
327: KSPChebyshevEstEigSetRandom(ksp,random);
328: PetscRandomDestroy(&random);
329: }
330: }
332: if (cheb->kspest) {
333: /* Mask the PC so that PCSetFromOptions does not do anything */
334: KSPSetPC(cheb->kspest,cheb->pcnone);
335: KSPSetOptionsPrefix(cheb->kspest,((PetscObject)ksp)->prefix);
336: KSPAppendOptionsPrefix(cheb->kspest,"est_");
337: if (!((PetscObject)cheb->kspest)->type_name) {
338: KSPSetType(cheb->kspest,KSPGMRES);
339: }
340: KSPSetFromOptions(cheb->kspest);
341: KSPSetPC(cheb->kspest,ksp->pc);
342: }
343: PetscOptionsTail();
344: return(0);
345: }
349: /*
350: * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
351: */
352: static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest,PetscReal *emin,PetscReal *emax)
353: {
355: PetscInt n,neig;
356: PetscReal *re,*im,min,max;
359: KSPGetIterationNumber(kspest,&n);
360: PetscMalloc2(n,PetscReal,&re,n,PetscReal,&im);
361: KSPComputeEigenvalues(kspest,n,re,im,&neig);
362: min = PETSC_MAX_REAL;
363: max = PETSC_MIN_REAL;
364: for (n=0; n<neig; n++) {
365: min = PetscMin(min,re[n]);
366: max = PetscMax(max,re[n]);
367: }
368: PetscFree2(re,im);
369: *emax = max;
370: *emin = min;
371: return(0);
372: }
376: PetscErrorCode KSPSolve_Chebyshev(KSP ksp)
377: {
378: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
380: PetscInt k,kp1,km1,maxit,ktmp,i;
381: PetscScalar alpha,omegaprod,mu,omega,Gamma,c[3],scale;
382: PetscReal rnorm = 0.0;
383: Vec sol_orig,b,p[3],r;
384: Mat Amat,Pmat;
385: MatStructure pflag;
386: PetscBool diagonalscale,hybrid=cheb->hybrid;
387: PetscBool purification=cheb->purification;
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: if (cheb->kspest && !cheb->estimate_current) {
394: PetscReal max,min;
395: Vec X,B;
396:
397: if (hybrid && purification) {
398: X = ksp->vec_sol;
399: } else {
400: X = ksp->work[0];
401: }
403: if (cheb->random) {
404: B = ksp->work[1];
405: VecSetRandom(B,cheb->random);
406: } else {
407: B = ksp->vec_rhs;
408: }
409: KSPSolve(cheb->kspest,B,X);
410: if (hybrid) {
411: cheb->its = 0; /* initialize Chebyshev iteration associated to kspest */
412: KSPSetInitialGuessNonzero(cheb->kspest,PETSC_TRUE);
413: } else if (ksp->guess_zero) {
414: VecZeroEntries(X);
415: }
416: KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest,&min,&max);
418: cheb->emin = cheb->tform[0]*min + cheb->tform[1]*max;
419: cheb->emax = cheb->tform[2]*min + cheb->tform[3]*max;
421: cheb->estimate_current = PETSC_TRUE;
422: }
424: ksp->its = 0;
425: PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);
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 || (hybrid && cheb->its== 0)) {
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: PetscObjectAMSTakeAccess((PetscObject)ksp);
462: if (hybrid && cheb->its && (cheb->its%cheb->chebysteps==0)) {
463: /* Adaptive step: update eigenvalue estimate - does not seem to improve convergence */
464: PetscReal max,min;
465: Vec X;
467: if (purification) {
468: X = p[km1]; /* will be updated by adaptive steps */
469: } else {
470: X = p[kp1]; /* temp vector */
471: }
473: VecCopy(p[k],X); /* p[k] = previous p[kp1] */
474: KSPSolve(cheb->kspest,ksp->vec_rhs,X);
475: KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest,&min,&max);
477: cheb->emin = cheb->tform[0]*min + cheb->tform[1]*max;
478: cheb->emax = cheb->tform[2]*min + cheb->tform[3]*max;
479: cheb->estimate_current = PETSC_TRUE;
481: /* update parameters that are dependent on emax and emin */
482: scale = 2.0/(cheb->emax + cheb->emin);
483: alpha = 1.0 - scale*(cheb->emin);
484: mu = 1.0/alpha;
485: omegaprod = 2.0/alpha;
487: c[km1] = 1.0;
488: c[k] = mu;
489: if (purification) { /* update p[k] */
490: KSP_MatMult(ksp,Amat,X,r); /* r = b - A*X */
491: VecAYPX(r,-1.0,b);
493: KSP_PCApply(ksp,r,p[k]); /* p[k] = scale B^{-1}r + X */
494: VecAYPX(p[k],scale,X);
495: }
496: }
498: ksp->its++;
499: cheb->its++;
500: PetscObjectAMSGrantAccess((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: PetscObjectAMSTakeAccess((PetscObject)ksp);
517: ksp->rnorm = rnorm;
518: PetscObjectAMSGrantAccess((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: VecScale(p[kp1],omega*Gamma*scale);
527: VecAXPY(p[kp1],1.0-omega,p[km1]);
528: VecAXPY(p[kp1],omega,p[k]);
530: ktmp = km1;
531: km1 = k;
532: k = kp1;
533: kp1 = ktmp;
534: }
535: if (!ksp->reason) {
536: if (ksp->normtype != KSP_NORM_NONE) {
537: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
538: VecAYPX(r,-1.0,b);
539: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
540: VecNorm(r,NORM_2,&rnorm);
541: } else {
542: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
543: VecNorm(p[kp1],NORM_2,&rnorm);
544: }
545: PetscObjectAMSTakeAccess((PetscObject)ksp);
546: ksp->rnorm = rnorm;
547: PetscObjectAMSGrantAccess((PetscObject)ksp);
548: ksp->vec_sol = p[k];
549: KSPLogResidualHistory(ksp,rnorm);
550: KSPMonitor(ksp,i,rnorm);
551: }
552: if (ksp->its >= ksp->max_it) {
553: if (ksp->normtype != KSP_NORM_NONE) {
554: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
555: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
556: } else ksp->reason = KSP_CONVERGED_ITS;
557: }
558: }
560: /* make sure solution is in vector x */
561: ksp->vec_sol = sol_orig;
562: if (k) {
563: VecCopy(p[k],sol_orig);
564: }
565: return(0);
566: }
570: PetscErrorCode KSPView_Chebyshev(KSP ksp,PetscViewer viewer)
571: {
572: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
574: PetscBool iascii;
577: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
578: if (iascii) {
579: PetscViewerASCIIPrintf(viewer," Chebyshev: eigenvalue estimates: min = %G, max = %G\n",cheb->emin,cheb->emax);
580: if (cheb->kspest) {
581: PetscViewerASCIIPrintf(viewer," Chebyshev: estimated using: [%G %G; %G %G]\n",cheb->tform[0],cheb->tform[1],cheb->tform[2],cheb->tform[3]);
582: if (cheb->hybrid) { /* display info about hybrid options being used */
583: PetscViewerASCIIPrintf(viewer," Chebyshev: hybrid is used, eststeps %D, chebysteps %D, purification %D\n",cheb->eststeps,cheb->chebysteps,cheb->purification);
584: }
585: if (cheb->random) {
586: PetscViewerASCIIPrintf(viewer," Chebyshev: estimating eigenvalues using random right hand side\n");
587: PetscViewerASCIIPushTab(viewer);
588: PetscRandomView(cheb->random,viewer);
589: PetscViewerASCIIPopTab(viewer);
590: }
591: PetscViewerASCIIPushTab(viewer);
592: KSPView(cheb->kspest,viewer);
593: PetscViewerASCIIPopTab(viewer);
594: }
595: }
596: return(0);
597: }
601: PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
602: {
603: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
607: KSPDestroy(&cheb->kspest);
608: PCDestroy(&cheb->pcnone);
609: PetscRandomDestroy(&cheb->random);
610: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",NULL);
611: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEstimateEigenvalues_C",NULL);
612: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetRandom_C",NULL);
613: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetNewMatrix_C",NULL);
614: KSPDestroyDefault(ksp);
615: return(0);
616: }
618: /*MC
619: KSPCHEBYSHEV - The preconditioned Chebyshev iterative method
621: Options Database Keys:
622: + -ksp_chebyshev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
623: of the preconditioned operator. If these are accurate you will get much faster convergence.
624: - -ksp_chebyshev_estimate_eigenvalues <a,b,c,d> - estimate eigenvalues using a Krylov method, then use this
625: transform for Chebyshev eigenvalue bounds (KSPChebyshevSetEstimateEigenvalues)
628: Level: beginner
630: Notes: The Chebyshev method requires both the matrix and preconditioner to
631: be symmetric positive (semi) definite.
632: Only support for left preconditioning.
634: Chebyshev is configured as a smoother by default, targetting the "upper" part of the spectrum.
635: The user should call KSPChebyshevSetEigenvalues() if they have eigenvalue estimates.
637: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
638: KSPChebyshevSetEigenvalues(), KSPChebyshevSetEstimateEigenvalues(), KSPRICHARDSON, KSPCG, PCMG
640: M*/
644: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
645: {
647: KSP_Chebyshev *chebyshevP;
650: PetscNewLog(ksp,KSP_Chebyshev,&chebyshevP);
652: ksp->data = (void*)chebyshevP;
653: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
654: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,1);
656: chebyshevP->emin = 0.;
657: chebyshevP->emax = 0.;
659: chebyshevP->tform[0] = 0.0;
660: chebyshevP->tform[1] = 0.1;
661: chebyshevP->tform[2] = 0;
662: chebyshevP->tform[3] = 1.1;
664: chebyshevP->hybrid = PETSC_FALSE;
665: chebyshevP->chebysteps = 20000;
666: chebyshevP->eststeps = 10;
667: chebyshevP->its = 0;
668: chebyshevP->purification = PETSC_TRUE;
670: ksp->ops->setup = KSPSetUp_Chebyshev;
671: ksp->ops->solve = KSPSolve_Chebyshev;
672: ksp->ops->destroy = KSPDestroy_Chebyshev;
673: ksp->ops->buildsolution = KSPBuildSolutionDefault;
674: ksp->ops->buildresidual = KSPBuildResidualDefault;
675: ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
676: ksp->ops->view = KSPView_Chebyshev;
677: ksp->ops->reset = KSPReset_Chebyshev;
679: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",KSPChebyshevSetEigenvalues_Chebyshev);
680: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEstimateEigenvalues_C",KSPChebyshevSetEstimateEigenvalues_Chebyshev);
681: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetRandom_C",KSPChebyshevEstEigSetRandom_Chebyshev);
682: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetNewMatrix_C",KSPChebyshevSetNewMatrix_Chebyshev);
683: return(0);
684: }