Actual source code: cheby.c
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;
8: if (cheb->kspest) {
9: KSPReset(cheb->kspest);
10: }
11: return 0;
12: }
14: /*
15: * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
16: */
17: static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest,PetscReal *emin,PetscReal *emax)
18: {
19: PetscInt n,neig;
20: PetscReal *re,*im,min,max;
22: KSPGetIterationNumber(kspest,&n);
23: PetscMalloc2(n,&re,n,&im);
24: KSPComputeEigenvalues(kspest,n,re,im,&neig);
25: min = PETSC_MAX_REAL;
26: max = PETSC_MIN_REAL;
27: for (n=0; n<neig; n++) {
28: min = PetscMin(min,re[n]);
29: max = PetscMax(max,re[n]);
30: }
31: PetscFree2(re,im);
32: *emax = max;
33: *emin = min;
34: return 0;
35: }
37: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
38: {
39: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
40: PetscBool flg;
41: Mat Pmat,Amat;
42: PetscObjectId amatid, pmatid;
43: PetscObjectState amatstate, pmatstate;
44: KSPSetWorkVecs(ksp,3);
45: if (cheb->emin == 0. || cheb->emax == 0.) { // User did not specify eigenvalues
46: PC pc;
47: KSPGetPC(ksp,&pc);
48: PetscObjectTypeCompare((PetscObject)pc,PCJACOBI,&flg);
49: if (!flg) { // Provided estimates are only relevant for Jacobi
50: cheb->emax_provided = 0;
51: cheb->emin_provided = 0;
52: }
53: if (!cheb->kspest) { /* We need to estimate eigenvalues */
54: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
55: }
56: }
57: if (cheb->kspest) {
58: KSPGetOperators(ksp,&Amat,&Pmat);
59: MatGetOption(Pmat, MAT_SPD, &flg);
60: if (flg) {
61: const char *prefix;
62: KSPGetOptionsPrefix(cheb->kspest,&prefix);
63: PetscOptionsHasName(NULL,prefix,"-ksp_type",&flg);
64: if (!flg) {
65: KSPSetType(cheb->kspest, KSPCG);
66: }
67: }
68: PetscObjectGetId((PetscObject)Amat,&amatid);
69: PetscObjectGetId((PetscObject)Pmat,&pmatid);
70: PetscObjectStateGet((PetscObject)Amat,&amatstate);
71: PetscObjectStateGet((PetscObject)Pmat,&pmatstate);
72: if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
73: PetscReal max=0.0,min=0.0;
74: Vec B;
75: KSPConvergedReason reason;
76: KSPSetPC(cheb->kspest,ksp->pc);
77: if (cheb->usenoisy) {
78: B = ksp->work[1];
79: KSPSetNoisy_Private(B);
80: } else {
81: PetscBool change;
84: PCPreSolveChangeRHS(ksp->pc,&change);
85: if (change) {
86: B = ksp->work[1];
87: VecCopy(ksp->vec_rhs,B);
88: } else B = ksp->vec_rhs;
89: }
90: KSPSolve(cheb->kspest,B,ksp->work[0]);
91: KSPGetConvergedReason(cheb->kspest,&reason);
92: if (reason == KSP_DIVERGED_ITS) {
93: PetscInfo(ksp,"Eigen estimator ran for prescribed number of iterations\n");
94: } else if (reason == KSP_DIVERGED_PC_FAILED) {
95: PetscInt its;
96: PCFailedReason pcreason;
98: KSPGetIterationNumber(cheb->kspest,&its);
99: if (ksp->normtype == KSP_NORM_NONE) {
100: PetscInt sendbuf,recvbuf;
101: PCGetFailedReasonRank(ksp->pc,&pcreason);
102: sendbuf = (PetscInt)pcreason;
103: MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));
104: PCSetFailedReason(ksp->pc,(PCFailedReason) recvbuf);
105: }
106: PCGetFailedReason(ksp->pc,&pcreason);
107: ksp->reason = KSP_DIVERGED_PC_FAILED;
108: PetscInfo(ksp,"Eigen estimator failed: %s %s at iteration %D",KSPConvergedReasons[reason],PCFailedReasons[pcreason],its);
109: return 0;
110: } else if (reason == KSP_CONVERGED_RTOL || reason == KSP_CONVERGED_ATOL) {
111: PetscInfo(ksp,"Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n");
112: } else if (reason < 0) {
113: PetscInfo(ksp,"Eigen estimator failed %s, using estimates anyway\n",KSPConvergedReasons[reason]);
114: }
116: KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest,&min,&max);
117: KSPSetPC(cheb->kspest,NULL);
119: cheb->emin_computed = min;
120: cheb->emax_computed = max;
122: cheb->amatid = amatid;
123: cheb->pmatid = pmatid;
124: cheb->amatstate = amatstate;
125: cheb->pmatstate = pmatstate;
126: }
127: }
128: return 0;
129: }
131: static PetscErrorCode KSPChebyshevGetEigenvalues_Chebyshev(KSP ksp,PetscReal *emax,PetscReal *emin)
132: {
133: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
135: *emax = 0;
136: *emin = 0;
137: if (cheb->emax != 0.) {
138: *emax = cheb->emax;
139: } else if (cheb->emax_computed != 0.) {
140: *emax = cheb->tform[2] * cheb->emin_computed + cheb->tform[3] * cheb->emax_computed;
141: } else if (cheb->emax_provided != 0.) {
142: *emax = cheb->tform[2] * cheb->emin_provided + cheb->tform[3] * cheb->emax_provided;
143: }
144: if (cheb->emin != 0.) {
145: *emin = cheb->emin;
146: } else if (cheb->emin_computed != 0.) {
147: *emin = cheb->tform[0] * cheb->emin_computed + cheb->tform[1] * cheb->emax_computed;
148: } else if (cheb->emin_provided != 0.) {
149: *emin = cheb->tform[0] * cheb->emin_provided + cheb->tform[1] * cheb->emax_provided;
150: }
151: return 0;
152: }
154: static PetscErrorCode KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp,PetscReal emax,PetscReal emin)
155: {
156: KSP_Chebyshev *chebyshevP = (KSP_Chebyshev*)ksp->data;
160: chebyshevP->emax = emax;
161: chebyshevP->emin = emin;
163: KSPChebyshevEstEigSet(ksp,0.,0.,0.,0.); /* Destroy any estimation setup */
164: return 0;
165: }
167: static PetscErrorCode KSPChebyshevEstEigSet_Chebyshev(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
168: {
169: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
171: if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
172: if ((cheb->emin_provided == 0. || cheb->emax_provided == 0.) && !cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
173: KSPCreate(PetscObjectComm((PetscObject)ksp),&cheb->kspest);
174: KSPSetErrorIfNotConverged(cheb->kspest,ksp->errorifnotconverged);
175: PetscObjectIncrementTabLevel((PetscObject)cheb->kspest,(PetscObject)ksp,1);
176: /* use PetscObjectSet/AppendOptionsPrefix() instead of KSPSet/AppendOptionsPrefix() so that the PC prefix is not changed */
177: PetscObjectSetOptionsPrefix((PetscObject)cheb->kspest,((PetscObject)ksp)->prefix);
178: PetscObjectAppendOptionsPrefix((PetscObject)cheb->kspest,"esteig_");
179: KSPSetSkipPCSetFromOptions(cheb->kspest,PETSC_TRUE);
181: KSPSetComputeEigenvalues(cheb->kspest,PETSC_TRUE);
183: /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
184: KSPSetTolerances(cheb->kspest,1.e-12,PETSC_DEFAULT,PETSC_DEFAULT,cheb->eststeps);
185: }
186: if (a >= 0) cheb->tform[0] = a;
187: if (b >= 0) cheb->tform[1] = b;
188: if (c >= 0) cheb->tform[2] = c;
189: if (d >= 0) cheb->tform[3] = d;
190: cheb->amatid = 0;
191: cheb->pmatid = 0;
192: cheb->amatstate = -1;
193: cheb->pmatstate = -1;
194: } else {
195: KSPDestroy(&cheb->kspest);
196: }
197: return 0;
198: }
200: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp,PetscBool use)
201: {
202: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
204: cheb->usenoisy = use;
205: return 0;
206: }
208: /*@
209: KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues
210: of the preconditioned problem.
212: Logically Collective on ksp
214: Input Parameters:
215: + ksp - the Krylov space context
216: - emax, emin - the eigenvalue estimates
218: Options Database:
219: . -ksp_chebyshev_eigenvalues emin,emax - extreme eigenvalues
221: Note:
222: Call KSPChebyshevEstEigSet() or use the option -ksp_chebyshev_esteig a,b,c,d to have the KSP
223: estimate the eigenvalues and use these estimated values automatically.
225: When KSPCHEBYSHEV is used as a smoother, one often wants to target a portion of the spectrum rather than the entire
226: spectrum. This function takes the range of target eigenvalues for Chebyshev, which will often slightly over-estimate
227: the largest eigenvalue of the actual operator (for safety) and greatly overestimate the smallest eigenvalue to
228: improve the smoothing properties of Chebyshev iteration on the higher frequencies in the spectrum.
230: Level: intermediate
232: .seealso: KSPChebyshevEstEigSet()
233: @*/
234: PetscErrorCode KSPChebyshevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
235: {
239: PetscTryMethod(ksp,"KSPChebyshevSetEigenvalues_C",(KSP,PetscReal,PetscReal),(ksp,emax,emin));
240: return 0;
241: }
243: /*@
244: KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev
246: Logically Collective on ksp
248: Input Parameters:
249: + ksp - the Krylov space context
250: . a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
251: . b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
252: . c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
253: - d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
255: Options Database:
256: . -ksp_chebyshev_esteig a,b,c,d - estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds
258: Notes:
259: The Chebyshev bounds are set using
260: .vb
261: minbound = a*minest + b*maxest
262: maxbound = c*minest + d*maxest
263: .ve
264: The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
265: The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.
267: If 0.0 is passed for all transform arguments (a,b,c,d), eigenvalue estimation is disabled.
269: The default transform is (0,0.1; 0,1.1) which targets the "upper" part of the spectrum, as desirable for use with multigrid.
271: The eigenvalues are estimated using the Lanczo (KSPCG) or Arnoldi (KSPGMRES) process using a noisy right hand side vector.
273: Level: intermediate
275: @*/
276: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
277: {
283: PetscTryMethod(ksp,"KSPChebyshevEstEigSet_C",(KSP,PetscReal,PetscReal,PetscReal,PetscReal),(ksp,a,b,c,d));
284: return 0;
285: }
287: /*@
288: KSPChebyshevEstEigSetUseNoisy - use a noisy right hand side in order to do the estimate instead of the given right hand side
290: Logically Collective
292: Input Parameters:
293: + ksp - linear solver context
294: - use - PETSC_TRUE to use noisy
296: Options Database:
297: . -ksp_chebyshev_esteig_noisy <true,false> - Use noisy right hand side for estimate
299: Notes:
300: This alledgely works better for multigrid smoothers
302: Level: intermediate
304: .seealso: KSPChebyshevEstEigSet()
305: @*/
306: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp,PetscBool use)
307: {
309: PetscTryMethod(ksp,"KSPChebyshevEstEigSetUseNoisy_C",(KSP,PetscBool),(ksp,use));
310: return 0;
311: }
313: /*@
314: KSPChebyshevEstEigGetKSP - Get the Krylov method context used to estimate eigenvalues for the Chebyshev method. If
315: a Krylov method is not being used for this purpose, NULL is returned. The reference count of the returned KSP is
316: not incremented: it should not be destroyed by the user.
318: Input Parameters:
319: . ksp - the Krylov space context
321: Output Parameters:
322: . kspest - the eigenvalue estimation Krylov space context
324: Level: intermediate
326: .seealso: KSPChebyshevEstEigSet()
327: @*/
328: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
329: {
332: *kspest = NULL;
333: PetscTryMethod(ksp,"KSPChebyshevEstEigGetKSP_C",(KSP,KSP*),(ksp,kspest));
334: return 0;
335: }
337: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
338: {
339: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
341: *kspest = cheb->kspest;
342: return 0;
343: }
345: static PetscErrorCode KSPSetFromOptions_Chebyshev(PetscOptionItems *PetscOptionsObject,KSP ksp)
346: {
347: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
348: PetscInt neigarg = 2, nestarg = 4;
349: PetscReal eminmax[2] = {0., 0.};
350: PetscReal tform[4] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
351: PetscBool flgeig, flgest;
353: PetscOptionsHead(PetscOptionsObject,"KSP Chebyshev Options");
354: PetscOptionsInt("-ksp_chebyshev_esteig_steps","Number of est steps in Chebyshev","",cheb->eststeps,&cheb->eststeps,NULL);
355: PetscOptionsRealArray("-ksp_chebyshev_eigenvalues","extreme eigenvalues","KSPChebyshevSetEigenvalues",eminmax,&neigarg,&flgeig);
356: if (flgeig) {
358: KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]);
359: }
360: PetscOptionsRealArray("-ksp_chebyshev_esteig","estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds","KSPChebyshevEstEigSet",tform,&nestarg,&flgest);
361: if (flgest) {
362: switch (nestarg) {
363: case 0:
364: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
365: break;
366: case 2: /* Base everything on the max eigenvalues */
367: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,tform[0],PETSC_DECIDE,tform[1]);
368: break;
369: case 4: /* Use the full 2x2 linear transformation */
370: KSPChebyshevEstEigSet(ksp,tform[0],tform[1],tform[2],tform[3]);
371: break;
372: default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
373: }
374: }
376: /* We need to estimate eigenvalues; need to set this here so that KSPSetFromOptions() is called on the estimator */
377: if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) {
378: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
379: }
381: if (cheb->kspest) {
382: PetscOptionsBool("-ksp_chebyshev_esteig_noisy","Use noisy right hand side for estimate","KSPChebyshevEstEigSetUseNoisy",cheb->usenoisy,&cheb->usenoisy,NULL);
383: KSPSetFromOptions(cheb->kspest);
384: }
385: PetscOptionsTail();
386: return 0;
387: }
389: static PetscErrorCode KSPSolve_Chebyshev(KSP ksp)
390: {
391: PetscInt k,kp1,km1,ktmp,i;
392: PetscScalar alpha,omegaprod,mu,omega,Gamma,c[3],scale;
393: PetscReal rnorm = 0.0,emax,emin;
394: Vec sol_orig,b,p[3],r;
395: Mat Amat,Pmat;
396: PetscBool diagonalscale;
398: PCGetDiagonalScale(ksp->pc,&diagonalscale);
401: PCGetOperators(ksp->pc,&Amat,&Pmat);
402: PetscObjectSAWsTakeAccess((PetscObject)ksp);
403: ksp->its = 0;
404: PetscObjectSAWsGrantAccess((PetscObject)ksp);
405: /* These three point to the three active solutions, we
406: rotate these three at each solution update */
407: km1 = 0; k = 1; kp1 = 2;
408: sol_orig = ksp->vec_sol; /* ksp->vec_sol will be assigned to rotating vector p[k], thus save its address */
409: b = ksp->vec_rhs;
410: p[km1] = sol_orig;
411: p[k] = ksp->work[0];
412: p[kp1] = ksp->work[1];
413: r = ksp->work[2];
415: KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin);
416: /* use scale*B as our preconditioner */
417: scale = 2.0/(emax + emin);
419: /* -alpha <= scale*lambda(B^{-1}A) <= alpha */
420: alpha = 1.0 - scale * 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,sol_orig,r); /* r = b - A*p[km1] */
430: VecAYPX(r,-1.0,b);
431: } else {
432: VecCopy(b,r);
433: }
435: /* calculate residual norm if requested, we have done one iteration */
436: if (ksp->normtype) {
437: switch (ksp->normtype) {
438: case KSP_NORM_PRECONDITIONED:
439: KSP_PCApply(ksp,r,p[k]); /* p[k] = B^{-1}r */
440: VecNorm(p[k],NORM_2,&rnorm);
441: break;
442: case KSP_NORM_UNPRECONDITIONED:
443: case KSP_NORM_NATURAL:
444: VecNorm(r,NORM_2,&rnorm);
445: break;
446: default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
447: }
448: PetscObjectSAWsTakeAccess((PetscObject)ksp);
449: ksp->rnorm = rnorm;
450: PetscObjectSAWsGrantAccess((PetscObject)ksp);
451: KSPLogResidualHistory(ksp,rnorm);
452: KSPLogErrorHistory(ksp);
453: KSPMonitor(ksp,0,rnorm);
454: (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
455: } else ksp->reason = KSP_CONVERGED_ITERATING;
456: if (ksp->reason || ksp->max_it==0) {
457: if (ksp->max_it==0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
458: return 0;
459: }
460: if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
461: KSP_PCApply(ksp,r,p[k]); /* p[k] = B^{-1}r */
462: }
463: VecAYPX(p[k],scale,p[km1]); /* p[k] = scale B^{-1}r + p[km1] */
464: PetscObjectSAWsTakeAccess((PetscObject)ksp);
465: ksp->its = 1;
466: PetscObjectSAWsGrantAccess((PetscObject)ksp);
468: for (i=1; i<ksp->max_it; i++) {
469: PetscObjectSAWsTakeAccess((PetscObject)ksp);
470: ksp->its++;
471: PetscObjectSAWsGrantAccess((PetscObject)ksp);
473: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
474: VecAYPX(r,-1.0,b);
475: /* calculate residual norm if requested */
476: if (ksp->normtype) {
477: switch (ksp->normtype) {
478: case KSP_NORM_PRECONDITIONED:
479: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
480: VecNorm(p[kp1],NORM_2,&rnorm);
481: break;
482: case KSP_NORM_UNPRECONDITIONED:
483: case KSP_NORM_NATURAL:
484: VecNorm(r,NORM_2,&rnorm);
485: break;
486: default:
487: rnorm = 0.0;
488: break;
489: }
490: KSPCheckNorm(ksp,rnorm);
491: PetscObjectSAWsTakeAccess((PetscObject)ksp);
492: ksp->rnorm = rnorm;
493: PetscObjectSAWsGrantAccess((PetscObject)ksp);
494: KSPLogResidualHistory(ksp,rnorm);
495: KSPMonitor(ksp,i,rnorm);
496: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
497: if (ksp->reason) break;
498: if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
499: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
500: }
501: } else {
502: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
503: }
504: ksp->vec_sol = p[k];
505: KSPLogErrorHistory(ksp);
507: c[kp1] = 2.0*mu*c[k] - c[km1];
508: omega = omegaprod*c[k]/c[kp1];
510: /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
511: VecAXPBYPCZ(p[kp1],1.0-omega,omega,omega*Gamma*scale,p[km1],p[k]);
513: ktmp = km1;
514: km1 = k;
515: k = kp1;
516: kp1 = ktmp;
517: }
518: if (!ksp->reason) {
519: if (ksp->normtype) {
520: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
521: VecAYPX(r,-1.0,b);
522: switch (ksp->normtype) {
523: case KSP_NORM_PRECONDITIONED:
524: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
525: VecNorm(p[kp1],NORM_2,&rnorm);
526: break;
527: case KSP_NORM_UNPRECONDITIONED:
528: case KSP_NORM_NATURAL:
529: VecNorm(r,NORM_2,&rnorm);
530: break;
531: default:
532: rnorm = 0.0;
533: break;
534: }
535: KSPCheckNorm(ksp,rnorm);
536: PetscObjectSAWsTakeAccess((PetscObject)ksp);
537: ksp->rnorm = rnorm;
538: PetscObjectSAWsGrantAccess((PetscObject)ksp);
539: KSPLogResidualHistory(ksp,rnorm);
540: KSPMonitor(ksp,i,rnorm);
541: }
542: if (ksp->its >= ksp->max_it) {
543: if (ksp->normtype != KSP_NORM_NONE) {
544: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
545: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
546: } else ksp->reason = KSP_CONVERGED_ITS;
547: }
548: }
550: /* make sure solution is in vector x */
551: ksp->vec_sol = sol_orig;
552: if (k) {
553: VecCopy(p[k],sol_orig);
554: }
555: if (ksp->reason == KSP_CONVERGED_ITS) {
556: KSPLogErrorHistory(ksp);
557: }
558: return 0;
559: }
561: static PetscErrorCode KSPView_Chebyshev(KSP ksp,PetscViewer viewer)
562: {
563: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
564: PetscBool iascii;
566: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
567: if (iascii) {
568: PetscReal emax,emin;
569: KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin);
570: PetscViewerASCIIPrintf(viewer," eigenvalue targets used: min %g, max %g\n",(double)emin,(double)emax);
571: if (cheb->kspest) {
572: PetscViewerASCIIPrintf(viewer," eigenvalues estimated 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 transform: [%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: } else if (cheb->emax_provided != 0.) {
581: PetscViewerASCIIPrintf(viewer, " eigenvalues provided (min %g, max %g) with transform: [%g %g; %g %g]\n", (double)cheb->emin_provided, (double)cheb->emax_provided, (double)cheb->tform[0],(double)cheb->tform[1],(double)cheb->tform[2],(double)cheb->tform[3]);
582: }
583: }
584: return 0;
585: }
587: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
588: {
589: 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: {
629: KSP_Chebyshev *chebyshevP;
631: PetscNewLog(ksp,&chebyshevP);
633: ksp->data = (void*)chebyshevP;
634: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
635: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
636: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
637: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
639: chebyshevP->emin = 0.;
640: chebyshevP->emax = 0.;
642: chebyshevP->tform[0] = 0.0;
643: chebyshevP->tform[1] = 0.1;
644: chebyshevP->tform[2] = 0;
645: chebyshevP->tform[3] = 1.1;
646: chebyshevP->eststeps = 10;
647: chebyshevP->usenoisy = PETSC_TRUE;
648: ksp->setupnewmatrix = PETSC_TRUE;
650: ksp->ops->setup = KSPSetUp_Chebyshev;
651: ksp->ops->solve = KSPSolve_Chebyshev;
652: ksp->ops->destroy = KSPDestroy_Chebyshev;
653: ksp->ops->buildsolution = KSPBuildSolutionDefault;
654: ksp->ops->buildresidual = KSPBuildResidualDefault;
655: ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
656: ksp->ops->view = KSPView_Chebyshev;
657: ksp->ops->reset = KSPReset_Chebyshev;
659: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",KSPChebyshevSetEigenvalues_Chebyshev);
660: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",KSPChebyshevEstEigSet_Chebyshev);
661: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",KSPChebyshevEstEigSetUseNoisy_Chebyshev);
662: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",KSPChebyshevEstEigGetKSP_Chebyshev);
663: return 0;
664: }