Actual source code: cheby.c
petsc-3.11.1 2019-04-17
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: KSPSetErrorIfNotConverged(cheb->kspest,ksp->errorifnotconverged);
52: PetscObjectIncrementTabLevel((PetscObject)cheb->kspest,(PetscObject)ksp,1);
53: KSPSetPC(cheb->kspest,ksp->pc);
54: /* use PetscObjectSet/AppendOptionsPrefix() instead of KSPSet/AppendOptionsPrefix() so that the PC prefix is not changed */
55: PetscObjectSetOptionsPrefix((PetscObject)cheb->kspest,((PetscObject)ksp)->prefix);
56: PetscObjectAppendOptionsPrefix((PetscObject)cheb->kspest,"esteig_");
57: KSPSetSkipPCSetFromOptions(cheb->kspest,PETSC_TRUE);
59: KSPSetComputeEigenvalues(cheb->kspest,PETSC_TRUE);
61: /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
62: KSPSetTolerances(cheb->kspest,1.e-12,PETSC_DEFAULT,PETSC_DEFAULT,cheb->eststeps);
63: }
64: if (a >= 0) cheb->tform[0] = a;
65: if (b >= 0) cheb->tform[1] = b;
66: if (c >= 0) cheb->tform[2] = c;
67: if (d >= 0) cheb->tform[3] = d;
68: cheb->amatid = 0;
69: cheb->pmatid = 0;
70: cheb->amatstate = -1;
71: cheb->pmatstate = -1;
72: } else {
73: KSPDestroy(&cheb->kspest);
74: }
75: return(0);
76: }
78: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp,PetscBool use)
79: {
80: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
83: cheb->usenoisy = use;
84: return(0);
85: }
87: /*@
88: KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues
89: of the preconditioned problem.
91: Logically Collective on KSP
93: Input Parameters:
94: + ksp - the Krylov space context
95: - emax, emin - the eigenvalue estimates
97: Options Database:
98: . -ksp_chebyshev_eigenvalues emin,emax
100: Note: Call KSPChebyshevEstEigSet() or use the option -ksp_chebyshev_esteig a,b,c,d to have the KSP
101: estimate the eigenvalues and use these estimated values automatically
103: Level: intermediate
105: .keywords: KSP, Chebyshev, set, eigenvalues
106: @*/
107: PetscErrorCode KSPChebyshevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
108: {
115: PetscTryMethod(ksp,"KSPChebyshevSetEigenvalues_C",(KSP,PetscReal,PetscReal),(ksp,emax,emin));
116: return(0);
117: }
119: /*@
120: KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev
122: Logically Collective on KSP
124: Input Parameters:
125: + ksp - the Krylov space context
126: . a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
127: . b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
128: . c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
129: - d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
131: Options Database:
132: . -ksp_chebyshev_esteig a,b,c,d
134: Notes:
135: The Chebyshev bounds are set using
136: .vb
137: minbound = a*minest + b*maxest
138: maxbound = c*minest + d*maxest
139: .ve
140: The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
141: The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.
143: If 0.0 is passed for all transform arguments (a,b,c,d), eigenvalue estimation is disabled.
145: The default transform is (0,0.1; 0,1.1) which targets the "upper" part of the spectrum, as desirable for use with multigrid.
147: The eigenvalues are estimated using the Lanczo (KSPCG) or Arnoldi (KSPGMRES) process using a noisy right hand side vector.
149: Level: intermediate
151: .keywords: KSP, Chebyshev, set, eigenvalues, PCMG
152: @*/
153: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
154: {
163: PetscTryMethod(ksp,"KSPChebyshevEstEigSet_C",(KSP,PetscReal,PetscReal,PetscReal,PetscReal),(ksp,a,b,c,d));
164: return(0);
165: }
167: /*@
168: KSPChebyshevEstEigSetUseNoisy - use a noisy right hand side in order to do the estimate instead of the given right hand side
170: Logically Collective
172: Input Arguments:
173: + ksp - linear solver context
174: - use - PETSC_TRUE to use noisy
176: Options Database:
177: + -ksp_chebyshev_esteig_noisy <true,false>
179: Notes:
180: This alledgely works better for multigrid smoothers
182: Level: intermediate
184: .seealso: KSPChebyshevEstEigSet()
185: @*/
186: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp,PetscBool use)
187: {
192: PetscTryMethod(ksp,"KSPChebyshevEstEigSetUseNoisy_C",(KSP,PetscBool),(ksp,use));
193: return(0);
194: }
196: /*@
197: KSPChebyshevEstEigGetKSP - Get the Krylov method context used to estimate eigenvalues for the Chebyshev method. If
198: a Krylov method is not being used for this purpose, NULL is returned. The reference count of the returned KSP is
199: not incremented: it should not be destroyed by the user.
201: Input Parameters:
202: . ksp - the Krylov space context
204: Output Parameters:
205: . kspest the eigenvalue estimation Krylov space context
207: Level: intermediate
209: .seealso: KSPChebyshevEstEigSet()
210: @*/
211: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
212: {
217: *kspest = NULL;
218: PetscTryMethod(ksp,"KSPChebyshevEstEigGetKSP_C",(KSP,KSP*),(ksp,kspest));
219: return(0);
220: }
222: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
223: {
224: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
227: *kspest = cheb->kspest;
228: return(0);
229: }
231: static PetscErrorCode KSPSetFromOptions_Chebyshev(PetscOptionItems *PetscOptionsObject,KSP ksp)
232: {
233: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
235: PetscInt neigarg = 2, nestarg = 4;
236: PetscReal eminmax[2] = {0., 0.};
237: PetscReal tform[4] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
238: PetscBool flgeig, flgest;
241: PetscOptionsHead(PetscOptionsObject,"KSP Chebyshev Options");
242: PetscOptionsInt("-ksp_chebyshev_esteig_steps","Number of est steps in Chebyshev","",cheb->eststeps,&cheb->eststeps,NULL);
243: PetscOptionsRealArray("-ksp_chebyshev_eigenvalues","extreme eigenvalues","KSPChebyshevSetEigenvalues",eminmax,&neigarg,&flgeig);
244: if (flgeig) {
245: if (neigarg != 2) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
246: KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]);
247: }
248: PetscOptionsRealArray("-ksp_chebyshev_esteig","estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds","KSPChebyshevEstEigSet",tform,&nestarg,&flgest);
249: if (flgest) {
250: switch (nestarg) {
251: case 0:
252: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
253: break;
254: case 2: /* Base everything on the max eigenvalues */
255: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,tform[0],PETSC_DECIDE,tform[1]);
256: break;
257: case 4: /* Use the full 2x2 linear transformation */
258: KSPChebyshevEstEigSet(ksp,tform[0],tform[1],tform[2],tform[3]);
259: break;
260: default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
261: }
262: }
264: /* We need to estimate eigenvalues; need to set this here so that KSPSetFromOptions() is called on the estimator */
265: if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) {
266: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
267: }
269: if (cheb->kspest) {
270: PetscOptionsBool("-ksp_chebyshev_esteig_noisy","Use noisy right hand side for estimate","KSPChebyshevEstEigSetUseNoisy",cheb->usenoisy,&cheb->usenoisy,NULL);
271: KSPSetFromOptions(cheb->kspest);
272: }
273: PetscOptionsTail();
274: return(0);
275: }
277: /*
278: * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
279: */
280: static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest,PetscReal *emin,PetscReal *emax)
281: {
283: PetscInt n,neig;
284: PetscReal *re,*im,min,max;
287: KSPGetIterationNumber(kspest,&n);
288: PetscMalloc2(n,&re,n,&im);
289: KSPComputeEigenvalues(kspest,n,re,im,&neig);
290: min = PETSC_MAX_REAL;
291: max = PETSC_MIN_REAL;
292: for (n=0; n<neig; n++) {
293: min = PetscMin(min,re[n]);
294: max = PetscMax(max,re[n]);
295: }
296: PetscFree2(re,im);
297: *emax = max;
298: *emin = min;
299: return(0);
300: }
302: PETSC_STATIC_INLINE PetscScalar chebyhash(PetscInt xx)
303: {
304: unsigned int x = xx;
305: x = ((x >> 16) ^ x) * 0x45d9f3b;
306: x = ((x >> 16) ^ x) * 0x45d9f3b;
307: x = ((x >> 16) ^ x);
308: return (PetscScalar)((PetscInt64)x-2147483648)*5.e-10; /* center around zero, scaled about -1. to 1.*/
309: }
311: static PetscErrorCode KSPSolve_Chebyshev(KSP ksp)
312: {
313: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
315: PetscInt k,kp1,km1,ktmp,i;
316: PetscScalar alpha,omegaprod,mu,omega,Gamma,c[3],scale;
317: PetscReal rnorm = 0.0;
318: Vec sol_orig,b,p[3],r;
319: Mat Amat,Pmat;
320: PetscBool diagonalscale;
323: PCGetDiagonalScale(ksp->pc,&diagonalscale);
324: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
326: PCGetOperators(ksp->pc,&Amat,&Pmat);
327: if (cheb->kspest) {
328: PetscObjectId amatid, pmatid;
329: PetscObjectState amatstate, pmatstate;
330: PCFailedReason pcreason;
331: PC pc;
333: PetscObjectGetId((PetscObject)Amat,&amatid);
334: PetscObjectGetId((PetscObject)Pmat,&pmatid);
335: PetscObjectStateGet((PetscObject)Amat,&amatstate);
336: PetscObjectStateGet((PetscObject)Pmat,&pmatstate);
337: if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
338: PetscReal max=0.0,min=0.0;
339: Vec B;
340: KSPConvergedReason reason;
342: if (cheb->usenoisy) {
343: B = ksp->work[1];
344: {
346: PetscInt n,i,istart;
347: PetscScalar *xx;
348: VecGetOwnershipRange(B,&istart,NULL);
349: VecGetLocalSize(B,&n);
350: VecGetArray(B,&xx);
351: for (i=0; i<n; i++) {
352: PetscScalar v = chebyhash(i+istart);
353: xx[i] = v;
354: }
355: VecRestoreArray(B,&xx);
356: }
357: } else {
358: PC pc;
359: PetscBool change;
361: KSPGetPC(cheb->kspest,&pc);
362: PCPreSolveChangeRHS(pc,&change);
363: if (change) {
364: B = ksp->work[1];
365: VecCopy(ksp->vec_rhs,B);
366: } else {
367: B = ksp->vec_rhs;
368: }
369: }
370: KSPSolve(cheb->kspest,B,ksp->work[0]);
371: KSPGetConvergedReason(cheb->kspest,&reason);
372: KSPGetPC(cheb->kspest,&pc);
373: PCGetFailedReason(pc,&pcreason);
374: if (reason < 0 || pcreason) {
375: if (reason == KSP_DIVERGED_ITS) {
376: PetscInfo(ksp,"Eigen estimator ran for prescribed number of iterations\n");
377: } else {
378: PetscInt its;
379: KSPGetIterationNumber(cheb->kspest,&its);
380: ksp->reason = KSP_DIVERGED_PC_FAILED;
381: VecSetInf(ksp->vec_sol);
382: PetscInfo3(ksp,"Eigen estimator failed: %s %s at iteration %D",KSPConvergedReasons[reason],PCFailedReasons[pcreason],its);
383: return(0);
384: }
385: } else if (reason==KSP_CONVERGED_RTOL ||reason==KSP_CONVERGED_ATOL) {
386: PetscInfo(ksp,"Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n");
387: } else {
388: PetscInfo1(ksp,"Eigen estimator did not converge by iteration: %s\n",KSPConvergedReasons[reason]);
389: }
391: KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest,&min,&max);
393: cheb->emin_computed = min;
394: cheb->emax_computed = max;
395: cheb->emin = cheb->tform[0]*min + cheb->tform[1]*max;
396: cheb->emax = cheb->tform[2]*min + cheb->tform[3]*max;
398: cheb->amatid = amatid;
399: cheb->pmatid = pmatid;
400: cheb->amatstate = amatstate;
401: cheb->pmatstate = pmatstate;
402: }
403: }
404: PetscObjectSAWsTakeAccess((PetscObject)ksp);
405: ksp->its = 0;
406: PetscObjectSAWsGrantAccess((PetscObject)ksp);
407: /* These three point to the three active solutions, we
408: rotate these three at each solution update */
409: km1 = 0; k = 1; kp1 = 2;
410: sol_orig = ksp->vec_sol; /* ksp->vec_sol will be asigned to rotating vector p[k], thus save its address */
411: b = ksp->vec_rhs;
412: p[km1] = sol_orig;
413: p[k] = ksp->work[0];
414: p[kp1] = ksp->work[1];
415: r = ksp->work[2];
417: /* use scale*B as our preconditioner */
418: scale = 2.0/(cheb->emax + cheb->emin);
420: /* -alpha <= scale*lambda(B^{-1}A) <= alpha */
421: alpha = 1.0 - scale*(cheb->emin);
422: Gamma = 1.0;
423: mu = 1.0/alpha;
424: omegaprod = 2.0/alpha;
426: c[km1] = 1.0;
427: c[k] = mu;
429: if (!ksp->guess_zero) {
430: KSP_MatMult(ksp,Amat,sol_orig,r); /* r = b - A*p[km1] */
431: VecAYPX(r,-1.0,b);
432: } else {
433: VecCopy(b,r);
434: }
436: /* calculate residual norm if requested, we have done one iteration */
437: if (ksp->normtype) {
438: switch (ksp->normtype) {
439: case KSP_NORM_PRECONDITIONED:
440: KSP_PCApply(ksp,r,p[k]); /* p[k] = B^{-1}r */
441: VecNorm(p[k],NORM_2,&rnorm);
442: break;
443: case KSP_NORM_UNPRECONDITIONED:
444: case KSP_NORM_NATURAL:
445: VecNorm(r,NORM_2,&rnorm);
446: break;
447: case KSP_NORM_NONE:
448: rnorm = 0.0;
449: break;
450: default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
451: }
452: PetscObjectSAWsTakeAccess((PetscObject)ksp);
453: ksp->rnorm = rnorm;
454: PetscObjectSAWsGrantAccess((PetscObject)ksp);
455: KSPLogResidualHistory(ksp,rnorm);
456: KSPMonitor(ksp,0,rnorm);
457: (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
458: }
459: else ksp->reason = KSP_CONVERGED_ITERATING;
460: if (ksp->reason || ksp->max_it==0) {
461: if (ksp->max_it==0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
462: return(0);
463: }
464: if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
465: KSP_PCApply(ksp,r,p[k]); /* p[k] = B^{-1}r */
466: }
467: VecAYPX(p[k],scale,p[km1]); /* p[k] = scale B^{-1}r + p[km1] */
468: PetscObjectSAWsTakeAccess((PetscObject)ksp);
469: ksp->its = 1;
470: PetscObjectSAWsGrantAccess((PetscObject)ksp);
472: for (i=1; i<ksp->max_it; i++) {
473: PetscObjectSAWsTakeAccess((PetscObject)ksp);
474: ksp->its++;
475: PetscObjectSAWsGrantAccess((PetscObject)ksp);
477: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
478: VecAYPX(r,-1.0,b);
479: /* calculate residual norm if requested */
480: if (ksp->normtype) {
481: switch (ksp->normtype) {
482: case KSP_NORM_PRECONDITIONED:
483: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
484: VecNorm(p[kp1],NORM_2,&rnorm);
485: break;
486: case KSP_NORM_UNPRECONDITIONED:
487: case KSP_NORM_NATURAL:
488: VecNorm(r,NORM_2,&rnorm);
489: break;
490: default:
491: rnorm = 0.0;
492: break;
493: }
494: KSPCheckNorm(ksp,rnorm);
495: PetscObjectSAWsTakeAccess((PetscObject)ksp);
496: ksp->rnorm = rnorm;
497: PetscObjectSAWsGrantAccess((PetscObject)ksp);
498: KSPLogResidualHistory(ksp,rnorm);
499: KSPMonitor(ksp,i,rnorm);
500: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
501: if (ksp->reason) break;
502: if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
503: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
504: }
505: } else {
506: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
507: }
508: ksp->vec_sol = p[k];
510: c[kp1] = 2.0*mu*c[k] - c[km1];
511: omega = omegaprod*c[k]/c[kp1];
513: /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
514: VecAXPBYPCZ(p[kp1],1.0-omega,omega,omega*Gamma*scale,p[km1],p[k]);
516: ktmp = km1;
517: km1 = k;
518: k = kp1;
519: kp1 = ktmp;
520: }
521: if (!ksp->reason) {
522: if (ksp->normtype) {
523: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
524: VecAYPX(r,-1.0,b);
525: switch (ksp->normtype) {
526: case KSP_NORM_PRECONDITIONED:
527: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
528: VecNorm(p[kp1],NORM_2,&rnorm);
529: break;
530: case KSP_NORM_UNPRECONDITIONED:
531: case KSP_NORM_NATURAL:
532: VecNorm(r,NORM_2,&rnorm);
533: break;
534: default:
535: rnorm = 0.0;
536: break;
537: }
538: KSPCheckNorm(ksp,rnorm);
539: PetscObjectSAWsTakeAccess((PetscObject)ksp);
540: ksp->rnorm = rnorm;
541: PetscObjectSAWsGrantAccess((PetscObject)ksp);
542: KSPLogResidualHistory(ksp,rnorm);
543: KSPMonitor(ksp,i,rnorm);
544: }
545: if (ksp->its >= ksp->max_it) {
546: if (ksp->normtype != KSP_NORM_NONE) {
547: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
548: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
549: } else ksp->reason = KSP_CONVERGED_ITS;
550: }
551: }
553: /* make sure solution is in vector x */
554: ksp->vec_sol = sol_orig;
555: if (k) {
556: VecCopy(p[k],sol_orig);
557: }
558: return(0);
559: }
561: static PetscErrorCode KSPView_Chebyshev(KSP ksp,PetscViewer viewer)
562: {
563: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
565: PetscBool iascii;
568: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
569: if (iascii) {
570: PetscViewerASCIIPrintf(viewer," eigenvalue estimates used: min = %g, max = %g\n",(double)cheb->emin,(double)cheb->emax);
571: if (cheb->kspest) {
572: PetscViewerASCIIPrintf(viewer," eigenvalues estimate via %s min %g, max %g\n",((PetscObject)(cheb->kspest))->type_name,(double)cheb->emin_computed,(double)cheb->emax_computed);
573: 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]);
574: PetscViewerASCIIPushTab(viewer);
575: KSPView(cheb->kspest,viewer);
576: PetscViewerASCIIPopTab(viewer);
577: if (cheb->usenoisy) {
578: PetscViewerASCIIPrintf(viewer," estimating eigenvalues using noisy right hand side\n");
579: }
580: }
581: }
582: return(0);
583: }
585: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
586: {
587: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
591: KSPDestroy(&cheb->kspest);
592: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",NULL);
593: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",NULL);
594: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",NULL);
595: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",NULL);
596: KSPDestroyDefault(ksp);
597: return(0);
598: }
600: /*MC
601: KSPCHEBYSHEV - The preconditioned Chebyshev iterative method
603: Options Database Keys:
604: + -ksp_chebyshev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
605: of the preconditioned operator. If these are accurate you will get much faster convergence.
606: . -ksp_chebyshev_esteig <a,b,c,d> - estimate eigenvalues using a Krylov method, then use this
607: transform for Chebyshev eigenvalue bounds (KSPChebyshevEstEigSet())
608: . -ksp_chebyshev_esteig_steps - number of estimation steps
609: - -ksp_chebyshev_esteig_noisy - use noisy number generator to create right hand side for eigenvalue estimator
611: Level: beginner
613: Notes:
614: The Chebyshev method requires both the matrix and preconditioner to
615: be symmetric positive (semi) definite.
616: Only support for left preconditioning.
618: Chebyshev is configured as a smoother by default, targetting the "upper" part of the spectrum.
619: The user should call KSPChebyshevSetEigenvalues() if they have eigenvalue estimates.
621: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
622: KSPChebyshevSetEigenvalues(), KSPChebyshevEstEigSet(), KSPChebyshevEstEigSetUseNoisy()
623: KSPRICHARDSON, KSPCG, PCMG
625: M*/
627: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
628: {
630: KSP_Chebyshev *chebyshevP;
633: PetscNewLog(ksp,&chebyshevP);
635: ksp->data = (void*)chebyshevP;
636: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
637: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
638: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
639: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
641: chebyshevP->emin = 0.;
642: chebyshevP->emax = 0.;
644: chebyshevP->tform[0] = 0.0;
645: chebyshevP->tform[1] = 0.1;
646: chebyshevP->tform[2] = 0;
647: chebyshevP->tform[3] = 1.1;
648: chebyshevP->eststeps = 10;
649: chebyshevP->usenoisy = PETSC_TRUE;
651: ksp->ops->setup = KSPSetUp_Chebyshev;
652: ksp->ops->solve = KSPSolve_Chebyshev;
653: ksp->ops->destroy = KSPDestroy_Chebyshev;
654: ksp->ops->buildsolution = KSPBuildSolutionDefault;
655: ksp->ops->buildresidual = KSPBuildResidualDefault;
656: ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
657: ksp->ops->view = KSPView_Chebyshev;
658: ksp->ops->reset = KSPReset_Chebyshev;
660: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",KSPChebyshevSetEigenvalues_Chebyshev);
661: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",KSPChebyshevEstEigSet_Chebyshev);
662: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",KSPChebyshevEstEigSetUseNoisy_Chebyshev);
663: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",KSPChebyshevEstEigGetKSP_Chebyshev);
664: return(0);
665: }