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) KSPReset(cheb->kspest);
9: return 0;
10: }
12: /*
13: * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
14: */
15: static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest, PetscReal *emin, PetscReal *emax)
16: {
17: PetscInt n, neig;
18: PetscReal *re, *im, min, max;
20: KSPGetIterationNumber(kspest, &n);
21: PetscMalloc2(n, &re, n, &im);
22: KSPComputeEigenvalues(kspest, n, re, im, &neig);
23: min = PETSC_MAX_REAL;
24: max = PETSC_MIN_REAL;
25: for (n = 0; n < neig; n++) {
26: min = PetscMin(min, re[n]);
27: max = PetscMax(max, re[n]);
28: }
29: PetscFree2(re, im);
30: *emax = max;
31: *emin = min;
32: return 0;
33: }
35: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
36: {
37: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
38: PetscBool isset, flg;
39: Mat Pmat, Amat;
40: PetscObjectId amatid, pmatid;
41: PetscObjectState amatstate, pmatstate;
43: KSPSetWorkVecs(ksp, 3);
44: if (cheb->emin == 0. || cheb->emax == 0.) { // User did not specify eigenvalues
45: PC pc;
46: KSPGetPC(ksp, &pc);
47: PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &flg);
48: if (!flg) { // Provided estimates are only relevant for Jacobi
49: cheb->emax_provided = 0;
50: cheb->emin_provided = 0;
51: }
52: if (!cheb->kspest) { /* We need to estimate eigenvalues */
53: KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);
54: }
55: }
56: if (cheb->kspest) {
57: KSPGetOperators(ksp, &Amat, &Pmat);
58: MatIsSPDKnown(Pmat, &isset, &flg);
59: if (isset && flg) {
60: const char *prefix;
61: KSPGetOptionsPrefix(cheb->kspest, &prefix);
62: PetscOptionsHasName(NULL, prefix, "-ksp_type", &flg);
63: if (!flg) KSPSetType(cheb->kspest, KSPCG);
64: }
65: PetscObjectGetId((PetscObject)Amat, &amatid);
66: PetscObjectGetId((PetscObject)Pmat, &pmatid);
67: PetscObjectStateGet((PetscObject)Amat, &amatstate);
68: PetscObjectStateGet((PetscObject)Pmat, &pmatstate);
69: if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
70: PetscReal max = 0.0, min = 0.0;
71: Vec B;
72: KSPConvergedReason reason;
73: KSPSetPC(cheb->kspest, ksp->pc);
74: if (cheb->usenoisy) {
75: B = ksp->work[1];
76: KSPSetNoisy_Private(B);
77: } else {
78: PetscBool change;
81: PCPreSolveChangeRHS(ksp->pc, &change);
82: if (change) {
83: B = ksp->work[1];
84: VecCopy(ksp->vec_rhs, B);
85: } else B = ksp->vec_rhs;
86: }
87: KSPSolve(cheb->kspest, B, ksp->work[0]);
88: KSPGetConvergedReason(cheb->kspest, &reason);
89: if (reason == KSP_DIVERGED_ITS) {
90: PetscInfo(ksp, "Eigen estimator ran for prescribed number of iterations\n");
91: } else if (reason == KSP_DIVERGED_PC_FAILED) {
92: PetscInt its;
93: PCFailedReason pcreason;
95: KSPGetIterationNumber(cheb->kspest, &its);
96: if (ksp->normtype == KSP_NORM_NONE) {
97: PetscInt sendbuf, recvbuf;
98: PCGetFailedReasonRank(ksp->pc, &pcreason);
99: sendbuf = (PetscInt)pcreason;
100: MPI_Allreduce(&sendbuf, &recvbuf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)ksp));
101: PCSetFailedReason(ksp->pc, (PCFailedReason)recvbuf);
102: }
103: PCGetFailedReason(ksp->pc, &pcreason);
104: ksp->reason = KSP_DIVERGED_PC_FAILED;
105: PetscInfo(ksp, "Eigen estimator failed: %s %s at iteration %" PetscInt_FMT, KSPConvergedReasons[reason], PCFailedReasons[pcreason], its);
106: return 0;
107: } else if (reason == KSP_CONVERGED_RTOL || reason == KSP_CONVERGED_ATOL) {
108: PetscInfo(ksp, "Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n");
109: } else if (reason < 0) {
110: PetscInfo(ksp, "Eigen estimator failed %s, using estimates anyway\n", KSPConvergedReasons[reason]);
111: }
113: KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest, &min, &max);
114: KSPSetPC(cheb->kspest, NULL);
116: cheb->emin_computed = min;
117: cheb->emax_computed = max;
119: cheb->amatid = amatid;
120: cheb->pmatid = pmatid;
121: cheb->amatstate = amatstate;
122: cheb->pmatstate = pmatstate;
123: }
124: }
125: return 0;
126: }
128: static PetscErrorCode KSPChebyshevGetEigenvalues_Chebyshev(KSP ksp, PetscReal *emax, PetscReal *emin)
129: {
130: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
132: *emax = 0;
133: *emin = 0;
134: if (cheb->emax != 0.) {
135: *emax = cheb->emax;
136: } else if (cheb->emax_computed != 0.) {
137: *emax = cheb->tform[2] * cheb->emin_computed + cheb->tform[3] * cheb->emax_computed;
138: } else if (cheb->emax_provided != 0.) {
139: *emax = cheb->tform[2] * cheb->emin_provided + cheb->tform[3] * cheb->emax_provided;
140: }
141: if (cheb->emin != 0.) {
142: *emin = cheb->emin;
143: } else if (cheb->emin_computed != 0.) {
144: *emin = cheb->tform[0] * cheb->emin_computed + cheb->tform[1] * cheb->emax_computed;
145: } else if (cheb->emin_provided != 0.) {
146: *emin = cheb->tform[0] * cheb->emin_provided + cheb->tform[1] * cheb->emax_provided;
147: }
148: return 0;
149: }
151: static PetscErrorCode KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp, PetscReal emax, PetscReal emin)
152: {
153: KSP_Chebyshev *chebyshevP = (KSP_Chebyshev *)ksp->data;
157: chebyshevP->emax = emax;
158: chebyshevP->emin = emin;
160: KSPChebyshevEstEigSet(ksp, 0., 0., 0., 0.); /* Destroy any estimation setup */
161: return 0;
162: }
164: static PetscErrorCode KSPChebyshevEstEigSet_Chebyshev(KSP ksp, PetscReal a, PetscReal b, PetscReal c, PetscReal d)
165: {
166: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
168: if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
169: if ((cheb->emin_provided == 0. || cheb->emax_provided == 0.) && !cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
170: KSPCreate(PetscObjectComm((PetscObject)ksp), &cheb->kspest);
171: KSPSetErrorIfNotConverged(cheb->kspest, ksp->errorifnotconverged);
172: PetscObjectIncrementTabLevel((PetscObject)cheb->kspest, (PetscObject)ksp, 1);
173: /* use PetscObjectSet/AppendOptionsPrefix() instead of KSPSet/AppendOptionsPrefix() so that the PC prefix is not changed */
174: PetscObjectSetOptionsPrefix((PetscObject)cheb->kspest, ((PetscObject)ksp)->prefix);
175: PetscObjectAppendOptionsPrefix((PetscObject)cheb->kspest, "esteig_");
176: KSPSetSkipPCSetFromOptions(cheb->kspest, PETSC_TRUE);
178: KSPSetComputeEigenvalues(cheb->kspest, PETSC_TRUE);
180: /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
181: KSPSetTolerances(cheb->kspest, 1.e-12, PETSC_DEFAULT, PETSC_DEFAULT, cheb->eststeps);
182: }
183: if (a >= 0) cheb->tform[0] = a;
184: if (b >= 0) cheb->tform[1] = b;
185: if (c >= 0) cheb->tform[2] = c;
186: if (d >= 0) cheb->tform[3] = d;
187: cheb->amatid = 0;
188: cheb->pmatid = 0;
189: cheb->amatstate = -1;
190: cheb->pmatstate = -1;
191: } else {
192: KSPDestroy(&cheb->kspest);
193: }
194: return 0;
195: }
197: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp, PetscBool use)
198: {
199: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
201: cheb->usenoisy = use;
202: return 0;
203: }
205: /*@
206: KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues
207: of the preconditioned problem.
209: Logically Collective
211: Input Parameters:
212: + ksp - the Krylov space context
213: - emax, emin - the eigenvalue estimates
215: Options Database Key:
216: . -ksp_chebyshev_eigenvalues emin,emax - extreme eigenvalues
218: Notes:
219: Call `KSPChebyshevEstEigSet()` or use the option -ksp_chebyshev_esteig a,b,c,d to have the KSP
220: estimate the eigenvalues and use these estimated values automatically.
222: When `KSPCHEBYSHEV` is used as a smoother, one often wants to target a portion of the spectrum rather than the entire
223: spectrum. This function takes the range of target eigenvalues for Chebyshev, which will often slightly over-estimate
224: the largest eigenvalue of the actual operator (for safety) and greatly overestimate the smallest eigenvalue to
225: improve the smoothing properties of Chebyshev iteration on the higher frequencies in the spectrum.
227: Level: intermediate
229: .seealso: [](chapter_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`,
230: @*/
231: PetscErrorCode KSPChebyshevSetEigenvalues(KSP ksp, PetscReal emax, PetscReal emin)
232: {
236: PetscTryMethod(ksp, "KSPChebyshevSetEigenvalues_C", (KSP, PetscReal, PetscReal), (ksp, emax, emin));
237: return 0;
238: }
240: /*@
241: KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev
243: Logically Collective
245: Input Parameters:
246: + ksp - the Krylov space context
247: . a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
248: . b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
249: . c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
250: - d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
252: Options Database Key:
253: . -ksp_chebyshev_esteig a,b,c,d - estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds
255: Notes:
256: The Chebyshev bounds are set using
257: .vb
258: minbound = a*minest + b*maxest
259: maxbound = c*minest + d*maxest
260: .ve
261: The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
262: The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.
264: If 0.0 is passed for all transform arguments (a,b,c,d), eigenvalue estimation is disabled.
266: The default transform is (0,0.1; 0,1.1) which targets the "upper" part of the spectrum, as desirable for use with multigrid.
268: The eigenvalues are estimated using the Lanczo (`KSPCG`) or Arnoldi (`KSPGMRES`) process using a noisy right hand side vector.
270: Level: intermediate
272: .seealso: [](chapter_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigSetUseNoisy()`, `KSPChebyshevEstEigGetKSP()`
273: @*/
274: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp, PetscReal a, PetscReal b, PetscReal c, PetscReal d)
275: {
281: PetscTryMethod(ksp, "KSPChebyshevEstEigSet_C", (KSP, PetscReal, PetscReal, PetscReal, PetscReal), (ksp, a, b, c, d));
282: return 0;
283: }
285: /*@
286: KSPChebyshevEstEigSetUseNoisy - use a noisy right hand side in order to do the estimate instead of the given right hand side
288: Logically Collective
290: Input Parameters:
291: + ksp - linear solver context
292: - use - `PETSC_TRUE` to use noisy
294: Options Database Key:
295: . -ksp_chebyshev_esteig_noisy <true,false> - Use noisy right hand side for estimate
297: Note:
298: This allegedly works better for multigrid smoothers
300: Level: intermediate
302: .seealso: [](chapter_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigGetKSP()`
303: @*/
304: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp, PetscBool use)
305: {
307: PetscTryMethod(ksp, "KSPChebyshevEstEigSetUseNoisy_C", (KSP, PetscBool), (ksp, use));
308: return 0;
309: }
311: /*@
312: KSPChebyshevEstEigGetKSP - Get the Krylov method context used to estimate eigenvalues for the Chebyshev method. If
313: a Krylov method is not being used for this purpose, NULL is returned. The reference count of the returned `KSP` is
314: not incremented: it should not be destroyed by the user.
316: Input Parameters:
317: . ksp - the Krylov space context
319: Output Parameters:
320: . kspest - the eigenvalue estimation Krylov space context
322: Level: advanced
324: .seealso: [](chapter_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`
325: @*/
326: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
327: {
330: *kspest = NULL;
331: PetscTryMethod(ksp, "KSPChebyshevEstEigGetKSP_C", (KSP, KSP *), (ksp, kspest));
332: return 0;
333: }
335: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
336: {
337: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
339: *kspest = cheb->kspest;
340: return 0;
341: }
343: static PetscErrorCode KSPSetFromOptions_Chebyshev(KSP ksp, PetscOptionItems *PetscOptionsObject)
344: {
345: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
346: PetscInt neigarg = 2, nestarg = 4;
347: PetscReal eminmax[2] = {0., 0.};
348: PetscReal tform[4] = {PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE};
349: PetscBool flgeig, flgest;
351: PetscOptionsHeadBegin(PetscOptionsObject, "KSP Chebyshev Options");
352: PetscOptionsInt("-ksp_chebyshev_esteig_steps", "Number of est steps in Chebyshev", "", cheb->eststeps, &cheb->eststeps, NULL);
353: PetscOptionsRealArray("-ksp_chebyshev_eigenvalues", "extreme eigenvalues", "KSPChebyshevSetEigenvalues", eminmax, &neigarg, &flgeig);
354: if (flgeig) {
356: KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]);
357: }
358: PetscOptionsRealArray("-ksp_chebyshev_esteig", "estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds", "KSPChebyshevEstEigSet", tform, &nestarg, &flgest);
359: if (flgest) {
360: switch (nestarg) {
361: case 0:
362: KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);
363: break;
364: case 2: /* Base everything on the max eigenvalues */
365: KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, tform[0], PETSC_DECIDE, tform[1]);
366: break;
367: case 4: /* Use the full 2x2 linear transformation */
368: KSPChebyshevEstEigSet(ksp, tform[0], tform[1], tform[2], tform[3]);
369: break;
370: default:
371: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
372: }
373: }
375: /* We need to estimate eigenvalues; need to set this here so that KSPSetFromOptions() is called on the estimator */
376: if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);
378: if (cheb->kspest) {
379: PetscOptionsBool("-ksp_chebyshev_esteig_noisy", "Use noisy right hand side for estimate", "KSPChebyshevEstEigSetUseNoisy", cheb->usenoisy, &cheb->usenoisy, NULL);
380: KSPSetFromOptions(cheb->kspest);
381: }
382: PetscOptionsHeadEnd();
383: return 0;
384: }
386: static PetscErrorCode KSPSolve_Chebyshev(KSP ksp)
387: {
388: PetscInt k, kp1, km1, ktmp, i;
389: PetscScalar alpha, omegaprod, mu, omega, Gamma, c[3], scale;
390: PetscReal rnorm = 0.0, emax, emin;
391: Vec sol_orig, b, p[3], r;
392: Mat Amat, Pmat;
393: PetscBool diagonalscale;
395: PCGetDiagonalScale(ksp->pc, &diagonalscale);
398: PCGetOperators(ksp->pc, &Amat, &Pmat);
399: PetscObjectSAWsTakeAccess((PetscObject)ksp);
400: ksp->its = 0;
401: PetscObjectSAWsGrantAccess((PetscObject)ksp);
402: /* These three point to the three active solutions, we
403: rotate these three at each solution update */
404: km1 = 0;
405: k = 1;
406: kp1 = 2;
407: sol_orig = ksp->vec_sol; /* ksp->vec_sol will be assigned to rotating vector p[k], thus save its address */
408: b = ksp->vec_rhs;
409: p[km1] = sol_orig;
410: p[k] = ksp->work[0];
411: p[kp1] = ksp->work[1];
412: r = ksp->work[2];
414: KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin);
415: /* use scale*B as our preconditioner */
416: scale = 2.0 / (emax + emin);
418: /* -alpha <= scale*lambda(B^{-1}A) <= alpha */
419: alpha = 1.0 - scale * emin;
420: Gamma = 1.0;
421: mu = 1.0 / alpha;
422: omegaprod = 2.0 / alpha;
424: c[km1] = 1.0;
425: c[k] = mu;
427: if (!ksp->guess_zero) {
428: KSP_MatMult(ksp, Amat, sol_orig, r); /* r = b - A*p[km1] */
429: VecAYPX(r, -1.0, b);
430: } else {
431: VecCopy(b, r);
432: }
434: /* calculate residual norm if requested, we have done one iteration */
435: if (ksp->normtype) {
436: switch (ksp->normtype) {
437: case KSP_NORM_PRECONDITIONED:
438: KSP_PCApply(ksp, r, p[k]); /* p[k] = B^{-1}r */
439: VecNorm(p[k], NORM_2, &rnorm);
440: break;
441: case KSP_NORM_UNPRECONDITIONED:
442: case KSP_NORM_NATURAL:
443: VecNorm(r, NORM_2, &rnorm);
444: break;
445: default:
446: 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) { KSP_PCApply(ksp, r, p[k]); /* p[k] = B^{-1}r */ }
461: VecAYPX(p[k], scale, p[km1]); /* p[k] = scale B^{-1}r + p[km1] */
462: PetscObjectSAWsTakeAccess((PetscObject)ksp);
463: ksp->its = 1;
464: PetscObjectSAWsGrantAccess((PetscObject)ksp);
466: for (i = 1; i < ksp->max_it; i++) {
467: PetscObjectSAWsTakeAccess((PetscObject)ksp);
468: ksp->its++;
469: PetscObjectSAWsGrantAccess((PetscObject)ksp);
471: KSP_MatMult(ksp, Amat, p[k], r); /* r = b - Ap[k] */
472: VecAYPX(r, -1.0, b);
473: /* calculate residual norm if requested */
474: if (ksp->normtype) {
475: switch (ksp->normtype) {
476: case KSP_NORM_PRECONDITIONED:
477: KSP_PCApply(ksp, r, p[kp1]); /* p[kp1] = B^{-1}r */
478: VecNorm(p[kp1], NORM_2, &rnorm);
479: break;
480: case KSP_NORM_UNPRECONDITIONED:
481: case KSP_NORM_NATURAL:
482: VecNorm(r, NORM_2, &rnorm);
483: break;
484: default:
485: rnorm = 0.0;
486: break;
487: }
488: KSPCheckNorm(ksp, rnorm);
489: PetscObjectSAWsTakeAccess((PetscObject)ksp);
490: ksp->rnorm = rnorm;
491: PetscObjectSAWsGrantAccess((PetscObject)ksp);
492: KSPLogResidualHistory(ksp, rnorm);
493: KSPMonitor(ksp, i, rnorm);
494: (*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP);
495: if (ksp->reason) break;
496: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { KSP_PCApply(ksp, r, p[kp1]); /* p[kp1] = B^{-1}r */ }
497: } else {
498: KSP_PCApply(ksp, r, p[kp1]); /* p[kp1] = B^{-1}r */
499: }
500: ksp->vec_sol = p[k];
501: KSPLogErrorHistory(ksp);
503: c[kp1] = 2.0 * mu * c[k] - c[km1];
504: omega = omegaprod * c[k] / c[kp1];
506: /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
507: VecAXPBYPCZ(p[kp1], 1.0 - omega, omega, omega * Gamma * scale, p[km1], p[k]);
509: ktmp = km1;
510: km1 = k;
511: k = kp1;
512: kp1 = ktmp;
513: }
514: if (!ksp->reason) {
515: if (ksp->normtype) {
516: KSP_MatMult(ksp, Amat, p[k], r); /* r = b - Ap[k] */
517: VecAYPX(r, -1.0, b);
518: switch (ksp->normtype) {
519: case KSP_NORM_PRECONDITIONED:
520: KSP_PCApply(ksp, r, p[kp1]); /* p[kp1] = B^{-1}r */
521: VecNorm(p[kp1], NORM_2, &rnorm);
522: break;
523: case KSP_NORM_UNPRECONDITIONED:
524: case KSP_NORM_NATURAL:
525: VecNorm(r, NORM_2, &rnorm);
526: break;
527: default:
528: rnorm = 0.0;
529: break;
530: }
531: KSPCheckNorm(ksp, rnorm);
532: PetscObjectSAWsTakeAccess((PetscObject)ksp);
533: ksp->rnorm = rnorm;
534: PetscObjectSAWsGrantAccess((PetscObject)ksp);
535: KSPLogResidualHistory(ksp, rnorm);
536: KSPMonitor(ksp, i, rnorm);
537: }
538: if (ksp->its >= ksp->max_it) {
539: if (ksp->normtype != KSP_NORM_NONE) {
540: (*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP);
541: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
542: } else ksp->reason = KSP_CONVERGED_ITS;
543: }
544: }
546: /* make sure solution is in vector x */
547: ksp->vec_sol = sol_orig;
548: if (k) VecCopy(p[k], sol_orig);
549: if (ksp->reason == KSP_CONVERGED_ITS) KSPLogErrorHistory(ksp);
550: return 0;
551: }
553: static PetscErrorCode KSPView_Chebyshev(KSP ksp, PetscViewer viewer)
554: {
555: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
556: PetscBool iascii;
558: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
559: if (iascii) {
560: PetscReal emax, emin;
561: KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin);
562: PetscViewerASCIIPrintf(viewer, " eigenvalue targets used: min %g, max %g\n", (double)emin, (double)emax);
563: if (cheb->kspest) {
564: PetscViewerASCIIPrintf(viewer, " eigenvalues estimated via %s: min %g, max %g\n", ((PetscObject)(cheb->kspest))->type_name, (double)cheb->emin_computed, (double)cheb->emax_computed);
565: 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]);
566: PetscViewerASCIIPushTab(viewer);
567: KSPView(cheb->kspest, viewer);
568: PetscViewerASCIIPopTab(viewer);
569: if (cheb->usenoisy) PetscViewerASCIIPrintf(viewer, " estimating eigenvalues using noisy right hand side\n");
570: } else if (cheb->emax_provided != 0.) {
571: PetscCall(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],
572: (double)cheb->tform[3]));
573: }
574: }
575: return 0;
576: }
578: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
579: {
580: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
582: KSPDestroy(&cheb->kspest);
583: PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", NULL);
584: PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", NULL);
585: PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", NULL);
586: PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", NULL);
587: KSPDestroyDefault(ksp);
588: return 0;
589: }
591: /*MC
592: KSPCHEBYSHEV - The preconditioned Chebyshev iterative method
594: Options Database Keys:
595: + -ksp_chebyshev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
596: of the preconditioned operator. If these are accurate you will get much faster convergence.
597: . -ksp_chebyshev_esteig <a,b,c,d> - estimate eigenvalues using a Krylov method, then use this
598: transform for Chebyshev eigenvalue bounds (`KSPChebyshevEstEigSet()`)
599: . -ksp_chebyshev_esteig_steps - number of estimation steps
600: - -ksp_chebyshev_esteig_noisy - use noisy number generator to create right hand side for eigenvalue estimator
602: Level: beginner
604: Notes:
605: The Chebyshev method requires both the matrix and preconditioner to be symmetric positive (semi) definite, but it can work as a smoother in other situations
607: Only support for left preconditioning.
609: Chebyshev is configured as a smoother by default, targetting the "upper" part of the spectrum.
611: The user should call `KSPChebyshevSetEigenvalues()` to get eigenvalue estimates.
613: .seealso: [](chapter_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
614: `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigSetUseNoisy()`
615: `KSPRICHARDSON`, `KSPCG`, `PCMG`
616: M*/
618: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
619: {
620: KSP_Chebyshev *chebyshevP;
622: PetscNew(&chebyshevP);
624: ksp->data = (void *)chebyshevP;
625: KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3);
626: KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2);
627: KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);
628: KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1);
630: chebyshevP->emin = 0.;
631: chebyshevP->emax = 0.;
633: chebyshevP->tform[0] = 0.0;
634: chebyshevP->tform[1] = 0.1;
635: chebyshevP->tform[2] = 0;
636: chebyshevP->tform[3] = 1.1;
637: chebyshevP->eststeps = 10;
638: chebyshevP->usenoisy = PETSC_TRUE;
639: ksp->setupnewmatrix = PETSC_TRUE;
641: ksp->ops->setup = KSPSetUp_Chebyshev;
642: ksp->ops->solve = KSPSolve_Chebyshev;
643: ksp->ops->destroy = KSPDestroy_Chebyshev;
644: ksp->ops->buildsolution = KSPBuildSolutionDefault;
645: ksp->ops->buildresidual = KSPBuildResidualDefault;
646: ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
647: ksp->ops->view = KSPView_Chebyshev;
648: ksp->ops->reset = KSPReset_Chebyshev;
650: PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", KSPChebyshevSetEigenvalues_Chebyshev);
651: PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", KSPChebyshevEstEigSet_Chebyshev);
652: PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", KSPChebyshevEstEigSetUseNoisy_Chebyshev);
653: PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", KSPChebyshevEstEigGetKSP_Chebyshev);
654: return 0;
655: }