Actual source code: cheby.c
1: #include "chebyshevimpl.h"
2: #include <../src/ksp/ksp/impls/cheby/chebyshevimpl.h>
4: static const char *const KSPChebyshevKinds[] = {"FIRST", "FOURTH", "OPT_FOURTH", "KSPChebyshevKinds", "KSP_CHEBYSHEV_", NULL};
6: static PetscErrorCode KSPReset_Chebyshev(KSP ksp)
7: {
8: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
10: PetscFunctionBegin;
11: if (cheb->kspest) PetscCall(KSPReset(cheb->kspest));
12: PetscFunctionReturn(PETSC_SUCCESS);
13: }
15: /*
16: Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
17: */
18: static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest, PetscReal *emin, PetscReal *emax)
19: {
20: PetscInt n, neig;
21: PetscReal *re, *im, min, max;
23: PetscFunctionBegin;
24: PetscCall(KSPGetIterationNumber(kspest, &n));
25: PetscCall(PetscMalloc2(n, &re, n, &im));
26: PetscCall(KSPComputeEigenvalues(kspest, n, re, im, &neig));
27: min = PETSC_MAX_REAL;
28: max = PETSC_MIN_REAL;
29: for (n = 0; n < neig; n++) {
30: min = PetscMin(min, re[n]);
31: max = PetscMax(max, re[n]);
32: }
33: PetscCall(PetscFree2(re, im));
34: *emax = max;
35: *emin = min;
36: PetscFunctionReturn(PETSC_SUCCESS);
37: }
39: static PetscErrorCode KSPChebyshevGetEigenvalues_Chebyshev(KSP ksp, PetscReal *emax, PetscReal *emin)
40: {
41: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
43: PetscFunctionBegin;
44: *emax = 0;
45: *emin = 0;
46: if (cheb->emax != 0.) {
47: *emax = cheb->emax;
48: } else if (cheb->emax_computed != 0.) {
49: *emax = cheb->tform[2] * cheb->emin_computed + cheb->tform[3] * cheb->emax_computed;
50: } else if (cheb->emax_provided != 0.) {
51: *emax = cheb->tform[2] * cheb->emin_provided + cheb->tform[3] * cheb->emax_provided;
52: }
53: if (cheb->emin != 0.) {
54: *emin = cheb->emin;
55: } else if (cheb->emin_computed != 0.) {
56: *emin = cheb->tform[0] * cheb->emin_computed + cheb->tform[1] * cheb->emax_computed;
57: } else if (cheb->emin_provided != 0.) {
58: *emin = cheb->tform[0] * cheb->emin_provided + cheb->tform[1] * cheb->emax_provided;
59: }
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: static PetscErrorCode KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp, PetscReal emax, PetscReal emin)
64: {
65: KSP_Chebyshev *chebyshevP = (KSP_Chebyshev *)ksp->data;
67: PetscFunctionBegin;
68: PetscCheck(emax > emin || (emax == 0 && emin == 0) || (emax == -1 && emin == -1), PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "Maximum eigenvalue must be larger than minimum: max %g min %g", (double)emax, (double)emin);
69: PetscCheck(emax * emin > 0.0 || (emax == 0 && emin == 0), PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "Both eigenvalues must be of the same sign: max %g min %g", (double)emax, (double)emin);
70: chebyshevP->emax = emax;
71: chebyshevP->emin = emin;
73: PetscCall(KSPChebyshevEstEigSet(ksp, 0., 0., 0., 0.)); /* Destroy any estimation setup */
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: static PetscErrorCode KSPChebyshevEstEigSet_Chebyshev(KSP ksp, PetscReal a, PetscReal b, PetscReal c, PetscReal d)
78: {
79: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
81: PetscFunctionBegin;
82: if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
83: if ((cheb->emin_provided == 0. || cheb->emax_provided == 0.) && !cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
84: PetscCall(KSPCreate(PetscObjectComm((PetscObject)ksp), &cheb->kspest));
85: PetscCall(KSPSetErrorIfNotConverged(cheb->kspest, ksp->errorifnotconverged));
86: PetscCall(PetscObjectIncrementTabLevel((PetscObject)cheb->kspest, (PetscObject)ksp, 1));
87: /* use PetscObjectSet/AppendOptionsPrefix() instead of KSPSet/AppendOptionsPrefix() so that the PC prefix is not changed */
88: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)cheb->kspest, ((PetscObject)ksp)->prefix));
89: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)cheb->kspest, "esteig_"));
90: PetscCall(KSPSetSkipPCSetFromOptions(cheb->kspest, PETSC_TRUE));
91: PetscCall(KSPSetComputeEigenvalues(cheb->kspest, PETSC_TRUE));
93: /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
94: PetscCall(KSPSetTolerances(cheb->kspest, 1.e-12, PETSC_DEFAULT, PETSC_DEFAULT, cheb->eststeps));
95: PetscCall(PetscInfo(ksp, "Created eigen estimator KSP\n"));
96: }
97: if (a >= 0) cheb->tform[0] = a;
98: if (b >= 0) cheb->tform[1] = b;
99: if (c >= 0) cheb->tform[2] = c;
100: if (d >= 0) cheb->tform[3] = d;
101: cheb->amatid = 0;
102: cheb->pmatid = 0;
103: cheb->amatstate = -1;
104: cheb->pmatstate = -1;
105: } else {
106: PetscCall(KSPDestroy(&cheb->kspest));
107: }
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp, PetscBool use)
112: {
113: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
115: PetscFunctionBegin;
116: cheb->usenoisy = use;
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: static PetscErrorCode KSPChebyshevSetKind_Chebyshev(KSP ksp, KSPChebyshevKind kind)
121: {
122: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
124: PetscFunctionBegin;
125: cheb->chebykind = kind;
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: static PetscErrorCode KSPChebyshevGetKind_Chebyshev(KSP ksp, KSPChebyshevKind *kind)
130: {
131: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
133: PetscFunctionBegin;
134: *kind = cheb->chebykind;
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
137: /*@
138: KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues of the preconditioned problem.
140: Logically Collective
142: Input Parameters:
143: + ksp - the Krylov space context
144: . emax - the eigenvalue maximum estimate
145: - emin - the eigenvalue minimum estimate
147: Options Database Key:
148: . -ksp_chebyshev_eigenvalues emin,emax - extreme eigenvalues
150: Level: intermediate
152: Notes:
153: Call `KSPChebyshevEstEigSet()` or use the option `-ksp_chebyshev_esteig a,b,c,d` to have the `KSP`
154: estimate the eigenvalues and use these estimated values automatically.
156: When `KSPCHEBYSHEV` is used as a smoother, one often wants to target a portion of the spectrum rather than the entire
157: spectrum. This function takes the range of target eigenvalues for Chebyshev, which will often slightly over-estimate
158: the largest eigenvalue of the actual operator (for safety) and greatly overestimate the smallest eigenvalue to
159: improve the smoothing properties of Chebyshev iteration on the higher frequencies in the spectrum.
161: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`,
162: @*/
163: PetscErrorCode KSPChebyshevSetEigenvalues(KSP ksp, PetscReal emax, PetscReal emin)
164: {
165: PetscFunctionBegin;
169: PetscTryMethod(ksp, "KSPChebyshevSetEigenvalues_C", (KSP, PetscReal, PetscReal), (ksp, emax, emin));
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: /*@
174: KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev
176: Logically Collective
178: Input Parameters:
179: + ksp - the Krylov space context
180: . a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or `PETSC_DECIDE`)
181: . b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or `PETSC_DECIDE`)
182: . c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or `PETSC_DECIDE`)
183: - d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or `PETSC_DECIDE`)
185: Options Database Key:
186: . -ksp_chebyshev_esteig a,b,c,d - estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds
188: Notes:
189: The Chebyshev bounds are set using
190: .vb
191: minbound = a*minest + b*maxest
192: maxbound = c*minest + d*maxest
193: .ve
194: The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
195: The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.
197: If 0.0 is passed for all transform arguments (a,b,c,d), eigenvalue estimation is disabled.
199: The default transform is (0,0.1; 0,1.1) which targets the "upper" part of the spectrum, as desirable for use with multigrid.
201: The eigenvalues are estimated using the Lanczo (`KSPCG`) or Arnoldi (`KSPGMRES`) process
203: Level: intermediate
205: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSetUseNoisy()`, `KSPChebyshevEstEigGetKSP()`
206: @*/
207: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp, PetscReal a, PetscReal b, PetscReal c, PetscReal d)
208: {
209: PetscFunctionBegin;
215: PetscTryMethod(ksp, "KSPChebyshevEstEigSet_C", (KSP, PetscReal, PetscReal, PetscReal, PetscReal), (ksp, a, b, c, d));
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: /*@
220: KSPChebyshevEstEigSetUseNoisy - use a noisy right hand side in order to do the estimate instead of the given right hand side
222: Logically Collective
224: Input Parameters:
225: + ksp - linear solver context
226: - use - `PETSC_TRUE` to use noisy
228: Options Database Key:
229: . -ksp_chebyshev_esteig_noisy <true,false> - Use noisy right hand side for estimate
231: Note:
232: This allegedly works better for multigrid smoothers
234: Level: intermediate
236: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigGetKSP()`
237: @*/
238: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp, PetscBool use)
239: {
240: PetscFunctionBegin;
243: PetscTryMethod(ksp, "KSPChebyshevEstEigSetUseNoisy_C", (KSP, PetscBool), (ksp, use));
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: /*@
248: KSPChebyshevEstEigGetKSP - Get the Krylov method context used to estimate eigenvalues for the Chebyshev method.
250: Input Parameter:
251: . ksp - the Krylov space context
253: Output Parameter:
254: . kspest - the eigenvalue estimation Krylov space context
256: Level: advanced
258: Notes:
259: If a Krylov method is not being used for this purpose, `NULL` is returned. The reference count of the returned `KSP` is
260: not incremented: it should not be destroyed by the user.
262: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`
263: @*/
264: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
265: {
266: PetscFunctionBegin;
268: PetscAssertPointer(kspest, 2);
269: *kspest = NULL;
270: PetscTryMethod(ksp, "KSPChebyshevEstEigGetKSP_C", (KSP, KSP *), (ksp, kspest));
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: /*@
275: KSPChebyshevSetKind - set the kind of Chebyshev polynomial to use
277: Logically Collective
279: Input Parameters:
280: + ksp - Linear solver context
281: - kind - The kind of Chebyshev polynomial to use, see `KSPChebyshevKind`
283: Options Database Key:
284: . -ksp_chebyshev_kind <kind> - which kind of Chebyshev polynomial to use
286: Level: intermediate
288: Note:
289: When using multigrid methods for problems with a poor quality coarse space (e.g., due to anisotropy or aggressive
290: coarsening), it is necessary for the smoother to handle smaller eigenvalues. With first-kind Chebyshev smoothing, this
291: requires using higher degree Chebyhev polynomials and reducing the lower end of the target spectrum, at which point
292: the whole target spectrum experiences about the same damping. Fourth kind Chebyshev polynomials (and the "optimized"
293: fourth kind) avoid the ad-hoc choice of lower bound and extend smoothing to smaller eigenvalues while preferentially
294: smoothing higher modes faster as needed to minimize the energy norm of the error.
296: References:
297: + * - Malachi Phillips and Paul Fischer, Optimal Chebyshev Smoothers and One-sided V-cycles, https://arxiv.org/abs/2210.03179.
298: - * - James Lottes, Optimal Polynomial Smoothers for Multigrid V-cycles, https://arxiv.org/abs/2202.08830.
300: .seealso: [](ch_ksp), `KSPCHEBYSHEV` `KSPChebyshevKind`, `KSPChebyshevGetKind()`
301: @*/
302: PetscErrorCode KSPChebyshevSetKind(KSP ksp, KSPChebyshevKind kind)
303: {
304: PetscFunctionBegin;
307: PetscUseMethod(ksp, "KSPChebyshevSetKind_C", (KSP, KSPChebyshevKind), (ksp, kind));
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: /*@
312: KSPChebyshevGetKind - get the kind of Chebyshev polynomial to use
314: Logically Collective
316: Input Parameters:
317: + ksp - Linear solver context
318: - kind - The kind of Chebyshev polynomial used
320: Level: intermediate
322: .seealso: [](ch_ksp), `KSPCHEBYSHEV` `KSPChebyshevKind`, `KSPChebyshevSetKind()`
323: @*/
324: PetscErrorCode KSPChebyshevGetKind(KSP ksp, KSPChebyshevKind *kind)
325: {
326: PetscFunctionBegin;
328: PetscUseMethod(ksp, "KSPChebyshevGetKind_C", (KSP, KSPChebyshevKind *), (ksp, kind));
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
333: {
334: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
336: PetscFunctionBegin;
337: *kspest = cheb->kspest;
338: PetscFunctionReturn(PETSC_SUCCESS);
339: }
341: static PetscErrorCode KSPSetFromOptions_Chebyshev(KSP ksp, PetscOptionItems *PetscOptionsObject)
342: {
343: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
344: PetscInt neigarg = 2, nestarg = 4;
345: PetscReal eminmax[2] = {0., 0.};
346: PetscReal tform[4] = {PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE};
347: PetscBool flgeig, flgest;
349: PetscFunctionBegin;
350: PetscOptionsHeadBegin(PetscOptionsObject, "KSP Chebyshev Options");
351: PetscCall(PetscOptionsInt("-ksp_chebyshev_esteig_steps", "Number of est steps in Chebyshev", "", cheb->eststeps, &cheb->eststeps, NULL));
352: PetscCall(PetscOptionsRealArray("-ksp_chebyshev_eigenvalues", "extreme eigenvalues", "KSPChebyshevSetEigenvalues", eminmax, &neigarg, &flgeig));
353: if (flgeig) {
354: PetscCheck(neigarg == 2, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
355: PetscCall(KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]));
356: }
357: PetscCall(PetscOptionsRealArray("-ksp_chebyshev_esteig", "estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds", "KSPChebyshevEstEigSet", tform, &nestarg, &flgest));
358: if (flgest) {
359: switch (nestarg) {
360: case 0:
361: PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
362: break;
363: case 2: /* Base everything on the max eigenvalues */
364: PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, tform[0], PETSC_DECIDE, tform[1]));
365: break;
366: case 4: /* Use the full 2x2 linear transformation */
367: PetscCall(KSPChebyshevEstEigSet(ksp, tform[0], tform[1], tform[2], tform[3]));
368: break;
369: default:
370: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
371: }
372: }
374: cheb->chebykind = KSP_CHEBYSHEV_FIRST; /* Default to 1st-kind Chebyshev polynomial */
375: PetscCall(PetscOptionsEnum("-ksp_chebyshev_kind", "Type of Chebyshev polynomial", "KSPChebyshevKind", KSPChebyshevKinds, (PetscEnum)cheb->chebykind, (PetscEnum *)&cheb->chebykind, NULL));
377: /* We need to estimate eigenvalues; need to set this here so that KSPSetFromOptions() is called on the estimator */
378: if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
380: if (cheb->kspest) {
381: PetscCall(PetscOptionsBool("-ksp_chebyshev_esteig_noisy", "Use noisy right hand side for estimate", "KSPChebyshevEstEigSetUseNoisy", cheb->usenoisy, &cheb->usenoisy, NULL));
382: PetscCall(KSPSetFromOptions(cheb->kspest));
383: }
384: PetscOptionsHeadEnd();
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: static PetscErrorCode KSPSolve_Chebyshev_FirstKind(KSP ksp)
389: {
390: PetscInt k, kp1, km1, ktmp, i;
391: PetscScalar alpha, omegaprod, mu, omega, Gamma, c[3], scale;
392: PetscReal rnorm = 0.0, emax, emin;
393: Vec sol_orig, b, p[3], r;
394: Mat Amat, Pmat;
395: PetscBool diagonalscale;
397: PetscFunctionBegin;
398: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
399: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
401: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
402: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
403: ksp->its = 0;
404: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
405: /* These three point to the three active solutions, we
406: rotate these three at each solution update */
407: km1 = 0;
408: k = 1;
409: kp1 = 2;
410: sol_orig = ksp->vec_sol; /* ksp->vec_sol will be assigned 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: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
418: /* use scale*B as our preconditioner */
419: scale = 2.0 / (emax + emin);
421: /* -alpha <= scale*lambda(B^{-1}A) <= alpha */
422: alpha = 1.0 - scale * emin;
423: Gamma = 1.0;
424: mu = 1.0 / alpha;
425: omegaprod = 2.0 / alpha;
427: c[km1] = 1.0;
428: c[k] = mu;
430: if (!ksp->guess_zero) {
431: PetscCall(KSP_MatMult(ksp, Amat, sol_orig, r)); /* r = b - A*p[km1] */
432: PetscCall(VecAYPX(r, -1.0, b));
433: } else {
434: PetscCall(VecCopy(b, r));
435: }
437: /* calculate residual norm if requested, we have done one iteration */
438: if (ksp->normtype) {
439: switch (ksp->normtype) {
440: case KSP_NORM_PRECONDITIONED:
441: PetscCall(KSP_PCApply(ksp, r, p[k])); /* p[k] = B^{-1}r */
442: PetscCall(VecNorm(p[k], NORM_2, &rnorm));
443: break;
444: case KSP_NORM_UNPRECONDITIONED:
445: case KSP_NORM_NATURAL:
446: PetscCall(VecNorm(r, NORM_2, &rnorm));
447: break;
448: default:
449: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
450: }
451: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
452: ksp->rnorm = rnorm;
453: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
454: PetscCall(KSPLogResidualHistory(ksp, rnorm));
455: PetscCall(KSPLogErrorHistory(ksp));
456: PetscCall(KSPMonitor(ksp, 0, rnorm));
457: PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
458: } else ksp->reason = KSP_CONVERGED_ITERATING;
459: if (ksp->reason || ksp->max_it == 0) {
460: if (ksp->max_it == 0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
461: PetscFunctionReturn(PETSC_SUCCESS);
462: }
463: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, p[k])); /* p[k] = B^{-1}r */ }
464: PetscCall(VecAYPX(p[k], scale, p[km1])); /* p[k] = scale B^{-1}r + p[km1] */
465: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
466: ksp->its = 1;
467: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
469: for (i = 1; i < ksp->max_it; i++) {
470: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
471: ksp->its++;
472: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
474: PetscCall(KSP_MatMult(ksp, Amat, p[k], r)); /* r = b - Ap[k] */
475: PetscCall(VecAYPX(r, -1.0, b));
476: /* calculate residual norm if requested */
477: if (ksp->normtype) {
478: switch (ksp->normtype) {
479: case KSP_NORM_PRECONDITIONED:
480: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
481: PetscCall(VecNorm(p[kp1], NORM_2, &rnorm));
482: break;
483: case KSP_NORM_UNPRECONDITIONED:
484: case KSP_NORM_NATURAL:
485: PetscCall(VecNorm(r, NORM_2, &rnorm));
486: break;
487: default:
488: rnorm = 0.0;
489: break;
490: }
491: KSPCheckNorm(ksp, rnorm);
492: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
493: ksp->rnorm = rnorm;
494: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
495: PetscCall(KSPLogResidualHistory(ksp, rnorm));
496: PetscCall(KSPMonitor(ksp, i, rnorm));
497: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
498: if (ksp->reason) break;
499: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */ }
500: } else {
501: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
502: }
503: ksp->vec_sol = p[k];
504: PetscCall(KSPLogErrorHistory(ksp));
506: c[kp1] = 2.0 * mu * c[k] - c[km1];
507: omega = omegaprod * c[k] / c[kp1];
509: /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
510: PetscCall(VecAXPBYPCZ(p[kp1], 1.0 - omega, omega, omega * Gamma * scale, p[km1], p[k]));
512: ktmp = km1;
513: km1 = k;
514: k = kp1;
515: kp1 = ktmp;
516: }
517: if (!ksp->reason) {
518: if (ksp->normtype) {
519: PetscCall(KSP_MatMult(ksp, Amat, p[k], r)); /* r = b - Ap[k] */
520: PetscCall(VecAYPX(r, -1.0, b));
521: switch (ksp->normtype) {
522: case KSP_NORM_PRECONDITIONED:
523: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
524: PetscCall(VecNorm(p[kp1], NORM_2, &rnorm));
525: break;
526: case KSP_NORM_UNPRECONDITIONED:
527: case KSP_NORM_NATURAL:
528: PetscCall(VecNorm(r, NORM_2, &rnorm));
529: break;
530: default:
531: rnorm = 0.0;
532: break;
533: }
534: KSPCheckNorm(ksp, rnorm);
535: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
536: ksp->rnorm = rnorm;
537: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
538: PetscCall(KSPLogResidualHistory(ksp, rnorm));
539: PetscCall(KSPMonitor(ksp, i, rnorm));
540: }
541: if (ksp->its >= ksp->max_it) {
542: if (ksp->normtype != KSP_NORM_NONE) {
543: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
544: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
545: } else ksp->reason = KSP_CONVERGED_ITS;
546: }
547: }
549: /* make sure solution is in vector x */
550: ksp->vec_sol = sol_orig;
551: if (k) PetscCall(VecCopy(p[k], sol_orig));
552: if (ksp->reason == KSP_CONVERGED_ITS) PetscCall(KSPLogErrorHistory(ksp));
553: PetscFunctionReturn(PETSC_SUCCESS);
554: }
556: static PetscErrorCode KSPSolve_Chebyshev_FourthKind(KSP ksp)
557: {
558: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
559: PetscInt i;
560: PetscScalar scale, rScale, dScale;
561: PetscReal rnorm = 0.0, emax, emin;
562: Vec x, b, d, r, Br;
563: Mat Amat, Pmat;
564: PetscBool diagonalscale;
565: PetscReal *betas = cheb->betas;
567: PetscFunctionBegin;
568: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
569: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
571: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
572: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
573: ksp->its = 0;
574: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
576: x = ksp->vec_sol;
577: b = ksp->vec_rhs;
578: r = ksp->work[0];
579: d = ksp->work[1];
580: Br = ksp->work[2];
582: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
583: /* use scale*B as our preconditioner */
584: scale = 1.0 / emax;
586: if (!ksp->guess_zero) {
587: PetscCall(KSP_MatMult(ksp, Amat, x, r)); /* r = b - A*x */
588: PetscCall(VecAYPX(r, -1.0, b));
589: } else {
590: PetscCall(VecCopy(b, r));
591: }
593: /* calculate residual norm if requested, we have done one iteration */
594: if (ksp->normtype) {
595: switch (ksp->normtype) {
596: case KSP_NORM_PRECONDITIONED:
597: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
598: PetscCall(VecNorm(Br, NORM_2, &rnorm));
599: break;
600: case KSP_NORM_UNPRECONDITIONED:
601: case KSP_NORM_NATURAL:
602: PetscCall(VecNorm(r, NORM_2, &rnorm));
603: break;
604: default:
605: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
606: }
607: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
608: ksp->rnorm = rnorm;
609: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
610: PetscCall(KSPLogResidualHistory(ksp, rnorm));
611: PetscCall(KSPLogErrorHistory(ksp));
612: PetscCall(KSPMonitor(ksp, 0, rnorm));
613: PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
614: } else ksp->reason = KSP_CONVERGED_ITERATING;
615: if (ksp->reason || ksp->max_it == 0) {
616: if (ksp->max_it == 0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
617: PetscFunctionReturn(PETSC_SUCCESS);
618: }
619: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */ }
620: PetscCall(VecAXPBY(d, 4.0 / 3.0 * scale, 0.0, Br)); /* d = 4/3 * scale B^{-1}r */
621: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
622: ksp->its = 1;
623: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
625: for (i = 1; i < ksp->max_it; i++) {
626: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
627: ksp->its++;
628: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
630: PetscCall(VecAXPBY(x, betas[i - 1], 1.0, d)); /* x = x + \beta_k d */
632: PetscCall(KSP_MatMult(ksp, Amat, d, Br)); /* r = r - Ad */
633: PetscCall(VecAXPBY(r, -1.0, 1.0, Br));
635: /* calculate residual norm if requested */
636: if (ksp->normtype) {
637: switch (ksp->normtype) {
638: case KSP_NORM_PRECONDITIONED:
639: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
640: PetscCall(VecNorm(Br, NORM_2, &rnorm));
641: break;
642: case KSP_NORM_UNPRECONDITIONED:
643: case KSP_NORM_NATURAL:
644: PetscCall(VecNorm(r, NORM_2, &rnorm));
645: break;
646: default:
647: rnorm = 0.0;
648: break;
649: }
650: KSPCheckNorm(ksp, rnorm);
651: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
652: ksp->rnorm = rnorm;
653: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
654: PetscCall(KSPLogResidualHistory(ksp, rnorm));
655: PetscCall(KSPMonitor(ksp, i, rnorm));
656: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
657: if (ksp->reason) break;
658: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */ }
659: } else {
660: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
661: }
662: PetscCall(KSPLogErrorHistory(ksp));
664: rScale = scale * (8.0 * i + 4.0) / (2.0 * i + 3.0);
665: dScale = (2.0 * i - 1.0) / (2.0 * i + 3.0);
667: /* d_k+1 = \dfrac{2k-1}{2k+3} d_k + \dfrac{8k+4}{2k+3} \dfrac{1}{\rho(SA)} Br */
668: PetscCall(VecAXPBY(d, rScale, dScale, Br));
669: }
671: /* on last pass, update solution vector */
672: PetscCall(VecAXPBY(x, betas[ksp->max_it - 1], 1.0, d)); /* x = x + d */
674: if (!ksp->reason) {
675: if (ksp->normtype) {
676: PetscCall(KSP_MatMult(ksp, Amat, x, r)); /* r = b - Ax */
677: PetscCall(VecAYPX(r, -1.0, b));
678: switch (ksp->normtype) {
679: case KSP_NORM_PRECONDITIONED:
680: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br= B^{-1}r */
681: PetscCall(VecNorm(Br, NORM_2, &rnorm));
682: break;
683: case KSP_NORM_UNPRECONDITIONED:
684: case KSP_NORM_NATURAL:
685: PetscCall(VecNorm(r, NORM_2, &rnorm));
686: break;
687: default:
688: rnorm = 0.0;
689: break;
690: }
691: KSPCheckNorm(ksp, rnorm);
692: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
693: ksp->rnorm = rnorm;
694: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
695: PetscCall(KSPLogResidualHistory(ksp, rnorm));
696: PetscCall(KSPMonitor(ksp, i, rnorm));
697: }
698: if (ksp->its >= ksp->max_it) {
699: if (ksp->normtype != KSP_NORM_NONE) {
700: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
701: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
702: } else ksp->reason = KSP_CONVERGED_ITS;
703: }
704: }
706: if (ksp->reason == KSP_CONVERGED_ITS) PetscCall(KSPLogErrorHistory(ksp));
707: PetscFunctionReturn(PETSC_SUCCESS);
708: }
710: static PetscErrorCode KSPView_Chebyshev(KSP ksp, PetscViewer viewer)
711: {
712: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
713: PetscBool iascii;
715: PetscFunctionBegin;
716: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
717: if (iascii) {
718: switch (cheb->chebykind) {
719: case KSP_CHEBYSHEV_FIRST:
720: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of first kind\n"));
721: break;
722: case KSP_CHEBYSHEV_FOURTH:
723: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of fourth kind\n"));
724: break;
725: case KSP_CHEBYSHEV_OPT_FOURTH:
726: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of opt. fourth kind\n"));
727: break;
728: }
729: PetscReal emax, emin;
730: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
731: PetscCall(PetscViewerASCIIPrintf(viewer, " eigenvalue targets used: min %g, max %g\n", (double)emin, (double)emax));
732: if (cheb->kspest) {
733: PetscCall(PetscViewerASCIIPrintf(viewer, " eigenvalues estimated via %s: min %g, max %g\n", ((PetscObject)(cheb->kspest))->type_name, (double)cheb->emin_computed, (double)cheb->emax_computed));
734: PetscCall(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]));
735: PetscCall(PetscViewerASCIIPushTab(viewer));
736: PetscCall(KSPView(cheb->kspest, viewer));
737: PetscCall(PetscViewerASCIIPopTab(viewer));
738: if (cheb->usenoisy) PetscCall(PetscViewerASCIIPrintf(viewer, " estimating eigenvalues using noisy right hand side\n"));
739: } else if (cheb->emax_provided != 0.) {
740: 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],
741: (double)cheb->tform[3]));
742: }
743: }
744: PetscFunctionReturn(PETSC_SUCCESS);
745: }
747: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
748: {
749: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
750: PetscBool isset, flg;
751: Mat Pmat, Amat;
752: PetscObjectId amatid, pmatid;
753: PetscObjectState amatstate, pmatstate;
755: PetscFunctionBegin;
756: switch (cheb->chebykind) {
757: case KSP_CHEBYSHEV_FIRST:
758: ksp->ops->solve = KSPSolve_Chebyshev_FirstKind;
759: break;
760: case KSP_CHEBYSHEV_FOURTH:
761: case KSP_CHEBYSHEV_OPT_FOURTH:
762: ksp->ops->solve = KSPSolve_Chebyshev_FourthKind;
763: break;
764: }
766: if (ksp->max_it > cheb->num_betas_alloc) {
767: PetscCall(PetscFree(cheb->betas));
768: PetscCall(PetscMalloc1(ksp->max_it, &cheb->betas));
769: cheb->num_betas_alloc = ksp->max_it;
770: }
772: // coefficients for 4th-kind Chebyshev
773: for (PetscInt i = 0; i < ksp->max_it; i++) cheb->betas[i] = 1.0;
775: // coefficients for optimized 4th-kind Chebyshev
776: if (cheb->chebykind == KSP_CHEBYSHEV_OPT_FOURTH) PetscCall(KSPChebyshevGetBetas_Private(ksp));
778: PetscCall(KSPSetWorkVecs(ksp, 3));
779: if (cheb->emin == 0. || cheb->emax == 0.) { // User did not specify eigenvalues
780: PC pc;
782: PetscCall(KSPGetPC(ksp, &pc));
783: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &flg));
784: if (!flg) { // Provided estimates are only relevant for Jacobi
785: cheb->emax_provided = 0;
786: cheb->emin_provided = 0;
787: }
788: /* We need to estimate eigenvalues */
789: if (!cheb->kspest) PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
790: }
791: if (cheb->kspest) {
792: PetscCall(KSPGetOperators(ksp, &Amat, &Pmat));
793: PetscCall(MatIsSPDKnown(Pmat, &isset, &flg));
794: if (isset && flg) {
795: const char *prefix;
797: PetscCall(KSPGetOptionsPrefix(cheb->kspest, &prefix));
798: PetscCall(PetscOptionsHasName(NULL, prefix, "-ksp_type", &flg));
799: if (!flg) PetscCall(KSPSetType(cheb->kspest, KSPCG));
800: }
801: PetscCall(PetscObjectGetId((PetscObject)Amat, &amatid));
802: PetscCall(PetscObjectGetId((PetscObject)Pmat, &pmatid));
803: PetscCall(PetscObjectStateGet((PetscObject)Amat, &amatstate));
804: PetscCall(PetscObjectStateGet((PetscObject)Pmat, &pmatstate));
805: if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
806: PetscReal max = 0.0, min = 0.0;
807: Vec B;
808: KSPConvergedReason reason;
810: PetscCall(KSPSetPC(cheb->kspest, ksp->pc));
811: if (cheb->usenoisy) {
812: B = ksp->work[1];
813: PetscCall(KSPSetNoisy_Private(B));
814: } else {
815: PetscBool change;
817: PetscCheck(ksp->vec_rhs, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Chebyshev must use a noisy right hand side to estimate the eigenvalues when no right hand side is available");
818: PetscCall(PCPreSolveChangeRHS(ksp->pc, &change));
819: if (change) {
820: B = ksp->work[1];
821: PetscCall(VecCopy(ksp->vec_rhs, B));
822: } else B = ksp->vec_rhs;
823: }
824: if (ksp->setfromoptionscalled && !cheb->kspest->setfromoptionscalled) PetscCall(KSPSetFromOptions(cheb->kspest));
825: PetscCall(KSPSolve(cheb->kspest, B, ksp->work[0]));
826: PetscCall(KSPGetConvergedReason(cheb->kspest, &reason));
827: if (reason == KSP_DIVERGED_ITS) {
828: PetscCall(PetscInfo(ksp, "Eigen estimator ran for prescribed number of iterations\n"));
829: } else if (reason == KSP_DIVERGED_PC_FAILED) {
830: PetscInt its;
831: PCFailedReason pcreason;
833: PetscCall(KSPGetIterationNumber(cheb->kspest, &its));
834: if (ksp->normtype == KSP_NORM_NONE) PetscCall(PCReduceFailedReason(ksp->pc));
835: PetscCall(PCGetFailedReason(ksp->pc, &pcreason));
836: ksp->reason = KSP_DIVERGED_PC_FAILED;
837: PetscCall(PetscInfo(ksp, "Eigen estimator failed: %s %s at iteration %" PetscInt_FMT "\n", KSPConvergedReasons[reason], PCFailedReasons[pcreason], its));
838: PetscFunctionReturn(PETSC_SUCCESS);
839: } else if (reason == KSP_CONVERGED_RTOL || reason == KSP_CONVERGED_ATOL) {
840: PetscCall(PetscInfo(ksp, "Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n"));
841: } else if (reason < 0) {
842: PetscCall(PetscInfo(ksp, "Eigen estimator failed %s, using estimates anyway\n", KSPConvergedReasons[reason]));
843: }
845: PetscCall(KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest, &min, &max));
846: PetscCall(KSPSetPC(cheb->kspest, NULL));
848: cheb->emin_computed = min;
849: cheb->emax_computed = max;
851: cheb->amatid = amatid;
852: cheb->pmatid = pmatid;
853: cheb->amatstate = amatstate;
854: cheb->pmatstate = pmatstate;
855: }
856: }
857: PetscFunctionReturn(PETSC_SUCCESS);
858: }
860: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
861: {
862: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
864: PetscFunctionBegin;
865: PetscCall(PetscFree(cheb->betas));
866: PetscCall(KSPDestroy(&cheb->kspest));
867: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", NULL));
868: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", NULL));
869: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", NULL));
870: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetKind_C", NULL));
871: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevGetKind_C", NULL));
872: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", NULL));
873: PetscCall(KSPDestroyDefault(ksp));
874: PetscFunctionReturn(PETSC_SUCCESS);
875: }
877: /*MC
878: KSPCHEBYSHEV - The preconditioned Chebyshev iterative method
880: Options Database Keys:
881: + -ksp_chebyshev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
882: of the preconditioned operator. If these are accurate you will get much faster convergence.
883: . -ksp_chebyshev_esteig <a,b,c,d> - estimate eigenvalues using a Krylov method, then use this
884: transform for Chebyshev eigenvalue bounds (`KSPChebyshevEstEigSet()`)
885: . -ksp_chebyshev_esteig_steps - number of estimation steps
886: - -ksp_chebyshev_esteig_noisy - use noisy number generator to create right hand side for eigenvalue estimator
888: Level: beginner
890: Notes:
891: The Chebyshev method requires both the matrix and preconditioner to be symmetric positive (semi) definite, but it can work
892: as a smoother in other situations
894: Only support for left preconditioning.
896: Chebyshev is configured as a smoother by default, targeting the "upper" part of the spectrum.
898: The user should call `KSPChebyshevSetEigenvalues()` to get eigenvalue estimates.
900: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
901: `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigSetUseNoisy()`
902: `KSPRICHARDSON`, `KSPCG`, `PCMG`
903: M*/
905: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
906: {
907: KSP_Chebyshev *chebyshevP;
909: PetscFunctionBegin;
910: PetscCall(PetscNew(&chebyshevP));
912: ksp->data = (void *)chebyshevP;
913: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
914: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
915: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
916: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
918: chebyshevP->emin = 0.;
919: chebyshevP->emax = 0.;
921: chebyshevP->tform[0] = 0.0;
922: chebyshevP->tform[1] = 0.1;
923: chebyshevP->tform[2] = 0;
924: chebyshevP->tform[3] = 1.1;
925: chebyshevP->eststeps = 10;
926: chebyshevP->usenoisy = PETSC_TRUE;
927: ksp->setupnewmatrix = PETSC_TRUE;
929: ksp->ops->setup = KSPSetUp_Chebyshev;
930: ksp->ops->destroy = KSPDestroy_Chebyshev;
931: ksp->ops->buildsolution = KSPBuildSolutionDefault;
932: ksp->ops->buildresidual = KSPBuildResidualDefault;
933: ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
934: ksp->ops->view = KSPView_Chebyshev;
935: ksp->ops->reset = KSPReset_Chebyshev;
937: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", KSPChebyshevSetEigenvalues_Chebyshev));
938: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", KSPChebyshevEstEigSet_Chebyshev));
939: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", KSPChebyshevEstEigSetUseNoisy_Chebyshev));
940: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetKind_C", KSPChebyshevSetKind_Chebyshev));
941: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevGetKind_C", KSPChebyshevGetKind_Chebyshev));
942: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", KSPChebyshevEstEigGetKSP_Chebyshev));
943: PetscFunctionReturn(PETSC_SUCCESS);
944: }