Actual source code: cheby.c
petsc-3.8.4 2018-03-24
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: PetscObjectIncrementTabLevel((PetscObject)cheb->kspest,(PetscObject)ksp,1);
52: KSPSetPC(cheb->kspest,ksp->pc);
53: /* use PetscObjectSet/AppendOptionsPrefix() instead of KSPSet/AppendOptionsPrefix() so that the PC prefix is not changed */
54: PetscObjectSetOptionsPrefix((PetscObject)cheb->kspest,((PetscObject)ksp)->prefix);
55: PetscObjectAppendOptionsPrefix((PetscObject)cheb->kspest,"esteig_");
56: KSPSetSkipPCSetFromOptions(cheb->kspest,PETSC_TRUE);
58: KSPSetComputeEigenvalues(cheb->kspest,PETSC_TRUE);
60: /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
61: KSPSetTolerances(cheb->kspest,1.e-12,PETSC_DEFAULT,PETSC_DEFAULT,cheb->eststeps);
62: }
63: if (a >= 0) cheb->tform[0] = a;
64: if (b >= 0) cheb->tform[1] = b;
65: if (c >= 0) cheb->tform[2] = c;
66: if (d >= 0) cheb->tform[3] = d;
67: cheb->amatid = 0;
68: cheb->pmatid = 0;
69: cheb->amatstate = -1;
70: cheb->pmatstate = -1;
71: } else {
72: KSPDestroy(&cheb->kspest);
73: }
74: return(0);
75: }
77: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp,PetscBool use)
78: {
79: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
82: cheb->usenoisy = use;
83: return(0);
84: }
86: /*@
87: KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues
88: of the preconditioned problem.
90: Logically Collective on KSP
92: Input Parameters:
93: + ksp - the Krylov space context
94: - emax, emin - the eigenvalue estimates
96: Options Database:
97: . -ksp_chebyshev_eigenvalues emin,emax
99: Note: Call KSPChebyshevEstEigSet() or use the option -ksp_chebyshev_esteig a,b,c,d to have the KSP
100: estimate the eigenvalues and use these estimated values automatically
102: Level: intermediate
104: .keywords: KSP, Chebyshev, set, eigenvalues
105: @*/
106: PetscErrorCode KSPChebyshevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
107: {
114: PetscTryMethod(ksp,"KSPChebyshevSetEigenvalues_C",(KSP,PetscReal,PetscReal),(ksp,emax,emin));
115: return(0);
116: }
118: /*@
119: KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev
121: Logically Collective on KSP
123: Input Parameters:
124: + ksp - the Krylov space context
125: . a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
126: . b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
127: . c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
128: - d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
130: Options Database:
131: . -ksp_chebyshev_esteig a,b,c,d
133: Notes:
134: The Chebyshev bounds are set using
135: .vb
136: minbound = a*minest + b*maxest
137: maxbound = c*minest + d*maxest
138: .ve
139: The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
140: The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.
142: If 0.0 is passed for all transform arguments (a,b,c,d), eigenvalue estimation is disabled.
144: The default transform is (0,0.1; 0,1.1) which targets the "upper" part of the spectrum, as desirable for use with multigrid.
146: The eigenvalues are estimated using the Lanczo (KSPCG) or Arnoldi (KSPGMRES) process using a noisy right hand side vector.
148: Level: intermediate
150: .keywords: KSP, Chebyshev, set, eigenvalues, PCMG
151: @*/
152: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
153: {
162: PetscTryMethod(ksp,"KSPChebyshevEstEigSet_C",(KSP,PetscReal,PetscReal,PetscReal,PetscReal),(ksp,a,b,c,d));
163: return(0);
164: }
166: /*@
167: KSPChebyshevEstEigSetUseNoisy - use a noisy right hand side in order to do the estimate instead of the given right hand side
169: Logically Collective
171: Input Arguments:
172: + ksp - linear solver context
173: - use - PETSC_TRUE to use noisy
175: Options Database:
176: + -ksp_chebyshev_esteig_noisy <true,false>
178: Notes: This alledgely works better for multigrid smoothers
180: Level: intermediate
182: .seealso: KSPChebyshevEstEigSet()
183: @*/
184: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp,PetscBool use)
185: {
190: PetscTryMethod(ksp,"KSPChebyshevEstEigSetUseNoisy_C",(KSP,PetscBool),(ksp,use));
191: return(0);
192: }
194: /*@
195: KSPChebyshevEstEigGetKSP - Get the Krylov method context used to estimate eigenvalues for the Chebyshev method. If
196: a Krylov method is not being used for this purpose, NULL is returned. The reference count of the returned KSP is
197: not incremented: it should not be destroyed by the user.
199: Input Parameters:
200: . ksp - the Krylov space context
202: Output Parameters:
203: . kspest the eigenvalue estimation Krylov space context
205: Level: intermediate
207: .seealso: KSPChebyshevEstEigSet()
208: @*/
209: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
210: {
215: *kspest = NULL;
216: PetscTryMethod(ksp,"KSPChebyshevEstEigGetKSP_C",(KSP,KSP*),(ksp,kspest));
217: return(0);
218: }
220: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
221: {
222: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
225: *kspest = cheb->kspest;
226: return(0);
227: }
229: static PetscErrorCode KSPSetFromOptions_Chebyshev(PetscOptionItems *PetscOptionsObject,KSP ksp)
230: {
231: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
233: PetscInt neigarg = 2, nestarg = 4;
234: PetscReal eminmax[2] = {0., 0.};
235: PetscReal tform[4] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
236: PetscBool flgeig, flgest;
239: PetscOptionsHead(PetscOptionsObject,"KSP Chebyshev Options");
240: PetscOptionsInt("-ksp_chebyshev_esteig_steps","Number of est steps in Chebyshev","",cheb->eststeps,&cheb->eststeps,NULL);
241: PetscOptionsRealArray("-ksp_chebyshev_eigenvalues","extreme eigenvalues","KSPChebyshevSetEigenvalues",eminmax,&neigarg,&flgeig);
242: if (flgeig) {
243: if (neigarg != 2) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
244: KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]);
245: }
246: PetscOptionsRealArray("-ksp_chebyshev_esteig","estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds","KSPChebyshevEstEigSet",tform,&nestarg,&flgest);
247: if (flgest) {
248: switch (nestarg) {
249: case 0:
250: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
251: break;
252: case 2: /* Base everything on the max eigenvalues */
253: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,tform[0],PETSC_DECIDE,tform[1]);
254: break;
255: case 4: /* Use the full 2x2 linear transformation */
256: KSPChebyshevEstEigSet(ksp,tform[0],tform[1],tform[2],tform[3]);
257: break;
258: default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
259: }
260: }
262: /* We need to estimate eigenvalues; need to set this here so that KSPSetFromOptions() is called on the estimator */
263: if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) {
264: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
265: }
267: if (cheb->kspest) {
268: PetscOptionsBool("-ksp_chebyshev_esteig_noisy","Use noisy right hand side for estimate","KSPChebyshevEstEigSetUseNoisy",cheb->usenoisy,&cheb->usenoisy,NULL);
269: KSPSetFromOptions(cheb->kspest);
270: }
271: PetscOptionsTail();
272: return(0);
273: }
275: /*
276: * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
277: */
278: static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest,PetscReal *emin,PetscReal *emax)
279: {
281: PetscInt n,neig;
282: PetscReal *re,*im,min,max;
285: KSPGetIterationNumber(kspest,&n);
286: PetscMalloc2(n,&re,n,&im);
287: KSPComputeEigenvalues(kspest,n,re,im,&neig);
288: min = PETSC_MAX_REAL;
289: max = PETSC_MIN_REAL;
290: for (n=0; n<neig; n++) {
291: min = PetscMin(min,re[n]);
292: max = PetscMax(max,re[n]);
293: }
294: PetscFree2(re,im);
295: *emax = max;
296: *emin = min;
297: return(0);
298: }
300: PETSC_STATIC_INLINE PetscScalar chebyhash(PetscInt xx)
301: {
302: unsigned int x = xx;
303: x = ((x >> 16) ^ x) * 0x45d9f3b;
304: x = ((x >> 16) ^ x) * 0x45d9f3b;
305: x = ((x >> 16) ^ x);
306: return (PetscScalar)((PetscInt64)x-2147483648)*5.e-10; /* center around zero, scaled about -1. to 1.*/
307: }
309: static PetscErrorCode KSPSolve_Chebyshev(KSP ksp)
310: {
311: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
313: PetscInt k,kp1,km1,maxit,ktmp,i;
314: PetscScalar alpha,omegaprod,mu,omega,Gamma,c[3],scale;
315: PetscReal rnorm = 0.0;
316: Vec sol_orig,b,p[3],r;
317: Mat Amat,Pmat;
318: PetscBool diagonalscale;
321: PCGetDiagonalScale(ksp->pc,&diagonalscale);
322: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
324: PCGetOperators(ksp->pc,&Amat,&Pmat);
325: if (cheb->kspest) {
326: PetscObjectId amatid, pmatid;
327: PetscObjectState amatstate, pmatstate;
329: PetscObjectGetId((PetscObject)Amat,&amatid);
330: PetscObjectGetId((PetscObject)Pmat,&pmatid);
331: PetscObjectStateGet((PetscObject)Amat,&amatstate);
332: PetscObjectStateGet((PetscObject)Pmat,&pmatstate);
333: if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
334: PetscReal max=0.0,min=0.0;
335: Vec B;
336: KSPConvergedReason reason;
338: if (cheb->usenoisy) {
339: B = ksp->work[1];
340: {
342: PetscInt n,i,istart;
343: PetscScalar *xx;
344: VecGetOwnershipRange(B,&istart,NULL);
345: VecGetLocalSize(B,&n);
346: VecGetArray(B,&xx);
347: for (i=0; i<n; i++) {
348: PetscScalar v = chebyhash(i+istart);
349: xx[i] = v;
350: }
351: VecRestoreArray(B,&xx);
352: }
353: } else {
354: PC pc;
355: PetscBool change;
357: KSPGetPC(cheb->kspest,&pc);
358: PCPreSolveChangeRHS(pc,&change);
359: if (change) {
360: B = ksp->work[1];
361: VecCopy(ksp->vec_rhs,B);
362: } else {
363: B = ksp->vec_rhs;
364: }
365: }
366: KSPSolve(cheb->kspest,B,ksp->work[0]);
368: KSPGetConvergedReason(cheb->kspest,&reason);
369: if (reason < 0) {
370: if (reason == KSP_DIVERGED_ITS) {
371: PetscInfo(ksp,"Eigen estimator ran for prescribed number of iterations\n");
372: } else {
373: PetscInt its;
374: KSPGetIterationNumber(cheb->kspest,&its);
375: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
376: PetscInfo(ksp,"Eigen estimator KSP_DIVERGED_PCSETUP_FAILED\n");
377: ksp->reason = KSP_DIVERGED_PCSETUP_FAILED;
378: VecSetInf(ksp->vec_sol);
379: } else {
380: SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Eigen estimator failed: %s at iteration %D",KSPConvergedReasons[reason],its);
381: }
382: }
383: } else if (reason==KSP_CONVERGED_RTOL ||reason==KSP_CONVERGED_ATOL) {
384: PetscInfo(ksp,"Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n");
385: } else {
386: PetscInfo1(ksp,"Eigen estimator did not converge by iteration: %s\n",KSPConvergedReasons[reason]);
387: }
389: KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest,&min,&max);
391: cheb->emin_computed = min;
392: cheb->emax_computed = max;
393: cheb->emin = cheb->tform[0]*min + cheb->tform[1]*max;
394: cheb->emax = cheb->tform[2]*min + cheb->tform[3]*max;
396: cheb->amatid = amatid;
397: cheb->pmatid = pmatid;
398: cheb->amatstate = amatstate;
399: cheb->pmatstate = pmatstate;
400: }
401: }
403: ksp->its = 0;
404: maxit = ksp->max_it;
406: /* These three point to the three active solutions, we
407: rotate these three at each solution update */
408: km1 = 0; k = 1; kp1 = 2;
409: sol_orig = ksp->vec_sol; /* ksp->vec_sol will be asigned to rotating vector p[k], thus save its address */
410: b = ksp->vec_rhs;
411: p[km1] = sol_orig;
412: p[k] = ksp->work[0];
413: p[kp1] = ksp->work[1];
414: r = ksp->work[2];
416: /* use scale*B as our preconditioner */
417: scale = 2.0/(cheb->emax + cheb->emin);
419: /* -alpha <= scale*lambda(B^{-1}A) <= alpha */
420: alpha = 1.0 - scale*(cheb->emin);
421: Gamma = 1.0;
422: mu = 1.0/alpha;
423: omegaprod = 2.0/alpha;
425: c[km1] = 1.0;
426: c[k] = mu;
428: if (!ksp->guess_zero) {
429: KSP_MatMult(ksp,Amat,p[km1],r); /* r = b - A*p[km1] */
430: VecAYPX(r,-1.0,b);
431: } else {
432: VecCopy(b,r);
433: }
435: KSP_PCApply(ksp,r,p[k]); /* p[k] = scale B^{-1}r + p[km1] */
436: VecAYPX(p[k],scale,p[km1]);
438: for (i=0; i<maxit; i++) {
439: PetscObjectSAWsTakeAccess((PetscObject)ksp);
441: ksp->its++;
442: PetscObjectSAWsGrantAccess((PetscObject)ksp);
443: c[kp1] = 2.0*mu*c[k] - c[km1];
444: omega = omegaprod*c[k]/c[kp1];
446: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
447: VecAYPX(r,-1.0,b);
448: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
449: ksp->vec_sol = p[k];
451: /* calculate residual norm if requested */
452: if (ksp->normtype != KSP_NORM_NONE || ksp->numbermonitors) {
453: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
454: VecNorm(r,NORM_2,&rnorm);
455: } else {
456: VecNorm(p[kp1],NORM_2,&rnorm);
457: }
458: PetscObjectSAWsTakeAccess((PetscObject)ksp);
459: ksp->rnorm = rnorm;
460: PetscObjectSAWsGrantAccess((PetscObject)ksp);
461: KSPLogResidualHistory(ksp,rnorm);
462: KSPMonitor(ksp,i,rnorm);
463: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
464: if (ksp->reason) break;
465: }
467: /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
468: VecAXPBYPCZ(p[kp1],1.0-omega,omega,omega*Gamma*scale,p[km1],p[k]);
470: ktmp = km1;
471: km1 = k;
472: k = kp1;
473: kp1 = ktmp;
474: }
475: if (!ksp->reason) {
476: if (ksp->normtype != KSP_NORM_NONE) {
477: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
478: VecAYPX(r,-1.0,b);
479: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
480: VecNorm(r,NORM_2,&rnorm);
481: } else {
482: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
483: VecNorm(p[kp1],NORM_2,&rnorm);
484: }
485: PetscObjectSAWsTakeAccess((PetscObject)ksp);
486: ksp->rnorm = rnorm;
487: PetscObjectSAWsGrantAccess((PetscObject)ksp);
488: ksp->vec_sol = p[k];
489: KSPLogResidualHistory(ksp,rnorm);
490: KSPMonitor(ksp,i,rnorm);
491: }
492: if (ksp->its >= ksp->max_it) {
493: if (ksp->normtype != KSP_NORM_NONE) {
494: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
495: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
496: } else ksp->reason = KSP_CONVERGED_ITS;
497: }
498: }
500: /* make sure solution is in vector x */
501: ksp->vec_sol = sol_orig;
502: if (k) {
503: VecCopy(p[k],sol_orig);
504: }
505: return(0);
506: }
508: static PetscErrorCode KSPView_Chebyshev(KSP ksp,PetscViewer viewer)
509: {
510: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
512: PetscBool iascii;
515: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
516: if (iascii) {
517: PetscViewerASCIIPrintf(viewer," eigenvalue estimates used: min = %g, max = %g\n",(double)cheb->emin,(double)cheb->emax);
518: if (cheb->kspest) {
519: PetscViewerASCIIPrintf(viewer," eigenvalues estimate via %s min %g, max %g\n",((PetscObject)(cheb->kspest))->type_name,(double)cheb->emin_computed,(double)cheb->emax_computed);
520: 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]);
521: PetscViewerASCIIPushTab(viewer);
522: KSPView(cheb->kspest,viewer);
523: PetscViewerASCIIPopTab(viewer);
524: if (cheb->usenoisy) {
525: PetscViewerASCIIPrintf(viewer," estimating eigenvalues using noisy right hand side\n");
526: }
527: }
528: }
529: return(0);
530: }
532: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
533: {
534: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
538: KSPDestroy(&cheb->kspest);
539: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",NULL);
540: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",NULL);
541: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",NULL);
542: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",NULL);
543: KSPDestroyDefault(ksp);
544: return(0);
545: }
547: /*MC
548: KSPCHEBYSHEV - The preconditioned Chebyshev iterative method
550: Options Database Keys:
551: + -ksp_chebyshev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
552: of the preconditioned operator. If these are accurate you will get much faster convergence.
553: . -ksp_chebyshev_esteig <a,b,c,d> - estimate eigenvalues using a Krylov method, then use this
554: transform for Chebyshev eigenvalue bounds (KSPChebyshevEstEigSet())
555: . -ksp_chebyshev_esteig_steps - number of estimation steps
556: - -ksp_chebyshev_esteig_noisy - use noisy number generator to create right hand side for eigenvalue estimator
558: Level: beginner
560: Notes: The Chebyshev method requires both the matrix and preconditioner to
561: be symmetric positive (semi) definite.
562: Only support for left preconditioning.
564: Chebyshev is configured as a smoother by default, targetting the "upper" part of the spectrum.
565: The user should call KSPChebyshevSetEigenvalues() if they have eigenvalue estimates.
567: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
568: KSPChebyshevSetEigenvalues(), KSPChebyshevEstEigSet(), KSPChebyshevEstEigSetUseNoisy()
569: KSPRICHARDSON, KSPCG, PCMG
571: M*/
573: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
574: {
576: KSP_Chebyshev *chebyshevP;
579: PetscNewLog(ksp,&chebyshevP);
581: ksp->data = (void*)chebyshevP;
582: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
583: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
584: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
585: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
587: chebyshevP->emin = 0.;
588: chebyshevP->emax = 0.;
590: chebyshevP->tform[0] = 0.0;
591: chebyshevP->tform[1] = 0.1;
592: chebyshevP->tform[2] = 0;
593: chebyshevP->tform[3] = 1.1;
594: chebyshevP->eststeps = 10;
595: chebyshevP->usenoisy = PETSC_TRUE;
597: ksp->ops->setup = KSPSetUp_Chebyshev;
598: ksp->ops->solve = KSPSolve_Chebyshev;
599: ksp->ops->destroy = KSPDestroy_Chebyshev;
600: ksp->ops->buildsolution = KSPBuildSolutionDefault;
601: ksp->ops->buildresidual = KSPBuildResidualDefault;
602: ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
603: ksp->ops->view = KSPView_Chebyshev;
604: ksp->ops->reset = KSPReset_Chebyshev;
606: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",KSPChebyshevSetEigenvalues_Chebyshev);
607: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",KSPChebyshevEstEigSet_Chebyshev);
608: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",KSPChebyshevEstEigSetUseNoisy_Chebyshev);
609: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",KSPChebyshevEstEigGetKSP_Chebyshev);
610: return(0);
611: }