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: PetscCall(PetscInfo(kspest, " eigen estimate min/max = %g %g\n", (double)min, (double)max));
37: PetscFunctionReturn(PETSC_SUCCESS);
38: }
40: static PetscErrorCode KSPChebyshevGetEigenvalues_Chebyshev(KSP ksp, PetscReal *emax, PetscReal *emin)
41: {
42: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
44: PetscFunctionBegin;
45: *emax = 0;
46: *emin = 0;
47: if (cheb->emax != 0.) {
48: *emax = cheb->emax;
49: } else if (cheb->emax_computed != 0.) {
50: *emax = cheb->tform[2] * cheb->emin_computed + cheb->tform[3] * cheb->emax_computed;
51: } else if (cheb->emax_provided != 0.) {
52: *emax = cheb->tform[2] * cheb->emin_provided + cheb->tform[3] * cheb->emax_provided;
53: }
54: if (cheb->emin != 0.) {
55: *emin = cheb->emin;
56: } else if (cheb->emin_computed != 0.) {
57: *emin = cheb->tform[0] * cheb->emin_computed + cheb->tform[1] * cheb->emax_computed;
58: } else if (cheb->emin_provided != 0.) {
59: *emin = cheb->tform[0] * cheb->emin_provided + cheb->tform[1] * cheb->emax_provided;
60: }
61: PetscFunctionReturn(PETSC_SUCCESS);
62: }
64: static PetscErrorCode KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp, PetscReal emax, PetscReal emin)
65: {
66: KSP_Chebyshev *chebyshevP = (KSP_Chebyshev *)ksp->data;
68: PetscFunctionBegin;
69: 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);
70: 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);
71: chebyshevP->emax = emax;
72: chebyshevP->emin = emin;
74: PetscCall(KSPChebyshevEstEigSet(ksp, 0., 0., 0., 0.)); /* Destroy any estimation setup */
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: static PetscErrorCode KSPChebyshevEstEigSet_Chebyshev(KSP ksp, PetscReal a, PetscReal b, PetscReal c, PetscReal d)
79: {
80: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
81: PetscInt nestlevel;
83: PetscFunctionBegin;
84: if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
85: if ((cheb->emin_provided == 0. || cheb->emax_provided == 0.) && !cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
86: PetscCall(KSPCreate(PetscObjectComm((PetscObject)ksp), &cheb->kspest));
87: PetscCall(KSPGetNestLevel(ksp, &nestlevel));
88: PetscCall(KSPSetNestLevel(cheb->kspest, nestlevel + 1));
89: PetscCall(KSPSetErrorIfNotConverged(cheb->kspest, ksp->errorifnotconverged));
90: PetscCall(PetscObjectIncrementTabLevel((PetscObject)cheb->kspest, (PetscObject)ksp, 1));
91: /* use PetscObjectSet/AppendOptionsPrefix() instead of KSPSet/AppendOptionsPrefix() so that the PC prefix is not changed */
92: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)cheb->kspest, ((PetscObject)ksp)->prefix));
93: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)cheb->kspest, "esteig_"));
94: PetscCall(KSPSetSkipPCSetFromOptions(cheb->kspest, PETSC_TRUE));
95: PetscCall(KSPSetComputeEigenvalues(cheb->kspest, PETSC_TRUE));
97: /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
98: PetscCall(KSPSetTolerances(cheb->kspest, 1.e-12, PETSC_CURRENT, PETSC_CURRENT, cheb->eststeps));
99: PetscCall(PetscInfo(ksp, "Created eigen estimator KSP\n"));
100: }
101: if (a >= 0) cheb->tform[0] = a;
102: if (b >= 0) cheb->tform[1] = b;
103: if (c >= 0) cheb->tform[2] = c;
104: if (d >= 0) cheb->tform[3] = d;
105: PetscCall(MatStateInvalidate(cheb->amatstate));
106: PetscCall(MatStateInvalidate(cheb->pmatstate));
107: } else {
108: PetscCall(KSPDestroy(&cheb->kspest));
109: }
110: PetscFunctionReturn(PETSC_SUCCESS);
111: }
113: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp, PetscBool use)
114: {
115: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
117: PetscFunctionBegin;
118: cheb->usenoisy = use;
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: static PetscErrorCode KSPChebyshevSetKind_Chebyshev(KSP ksp, KSPChebyshevKind kind)
123: {
124: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
126: PetscFunctionBegin;
127: cheb->chebykind = kind;
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: static PetscErrorCode KSPChebyshevGetKind_Chebyshev(KSP ksp, KSPChebyshevKind *kind)
132: {
133: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
135: PetscFunctionBegin;
136: *kind = cheb->chebykind;
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
139: /*@
140: KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues of the preconditioned problem.
142: Logically Collective
144: Input Parameters:
145: + ksp - the Krylov space context
146: . emax - the eigenvalue maximum estimate
147: - emin - the eigenvalue minimum estimate
149: Options Database Key:
150: . -ksp_chebyshev_eigenvalues emin,emax - extreme eigenvalues
152: Level: intermediate
154: Notes:
155: Call `KSPChebyshevEstEigSet()` or use the option `-ksp_chebyshev_esteig a,b,c,d` to have the `KSP`
156: estimate the eigenvalues and use these estimated values automatically.
158: When `KSPCHEBYSHEV` is used as a smoother, one often wants to target a portion of the spectrum rather than the entire
159: spectrum. This function takes the range of target eigenvalues for Chebyshev, which will often slightly over-estimate
160: the largest eigenvalue of the actual operator (for safety) and greatly overestimate the smallest eigenvalue to
161: improve the smoothing properties of Chebyshev iteration on the higher frequencies in the spectrum.
163: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`
164: @*/
165: PetscErrorCode KSPChebyshevSetEigenvalues(KSP ksp, PetscReal emax, PetscReal emin)
166: {
167: PetscFunctionBegin;
171: PetscTryMethod(ksp, "KSPChebyshevSetEigenvalues_C", (KSP, PetscReal, PetscReal), (ksp, emax, emin));
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: /*@
176: KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev
178: Logically Collective
180: Input Parameters:
181: + ksp - the Krylov space context
182: . a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or `PETSC_DECIDE`)
183: . b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or `PETSC_DECIDE`)
184: . c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or `PETSC_DECIDE`)
185: - d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or `PETSC_DECIDE`)
187: Options Database Key:
188: . -ksp_chebyshev_esteig a,b,c,d - estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds
190: Notes:
191: The Chebyshev bounds are set using
192: .vb
193: minbound = a*minest + b*maxest
194: maxbound = c*minest + d*maxest
195: .ve
196: The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
197: The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.
199: If 0.0 is passed for all transform arguments (a,b,c,d), eigenvalue estimation is disabled.
201: The default transform is (0,0.1; 0,1.1) which targets the "upper" part of the spectrum, as desirable for use with multigrid.
203: The eigenvalues are estimated using the Lanczos (`KSPCG`) or Arnoldi (`KSPGMRES`) process depending on if the operator is
204: symmetric definite or not.
206: Level: intermediate
208: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSetUseNoisy()`, `KSPChebyshevEstEigGetKSP()`
209: @*/
210: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp, PetscReal a, PetscReal b, PetscReal c, PetscReal d)
211: {
212: PetscFunctionBegin;
218: PetscTryMethod(ksp, "KSPChebyshevEstEigSet_C", (KSP, PetscReal, PetscReal, PetscReal, PetscReal), (ksp, a, b, c, d));
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: /*@
223: KSPChebyshevEstEigSetUseNoisy - use a noisy random number generated right-hand side to estimate the extreme eigenvalues instead of the given right-hand side
225: Logically Collective
227: Input Parameters:
228: + ksp - linear solver context
229: - use - `PETSC_TRUE` to use noisy
231: Options Database Key:
232: . -ksp_chebyshev_esteig_noisy (true|false) - Use noisy right-hand side for estimate
234: Level: intermediate
236: Note:
237: This allegedly works better for multigrid smoothers
239: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigGetKSP()`
240: @*/
241: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp, PetscBool use)
242: {
243: PetscFunctionBegin;
246: PetscTryMethod(ksp, "KSPChebyshevEstEigSetUseNoisy_C", (KSP, PetscBool), (ksp, use));
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: /*@
251: KSPChebyshevEstEigGetKSP - Get the Krylov method context used to estimate the eigenvalues for the Chebyshev method.
253: Input Parameter:
254: . ksp - the Krylov space context
256: Output Parameter:
257: . kspest - the eigenvalue estimation Krylov space context
259: Level: advanced
261: Notes:
262: If a Krylov method is not being used for this purpose, `NULL` is returned. The reference count of the returned `KSP` is
263: not incremented: it should not be destroyed by the user.
265: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`
266: @*/
267: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
268: {
269: PetscFunctionBegin;
271: PetscAssertPointer(kspest, 2);
272: *kspest = NULL;
273: PetscTryMethod(ksp, "KSPChebyshevEstEigGetKSP_C", (KSP, KSP *), (ksp, kspest));
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: /*@
278: KSPChebyshevSetKind - set the kind of Chebyshev polynomial to use
280: Logically Collective
282: Input Parameters:
283: + ksp - Linear solver context
284: - kind - The kind of Chebyshev polynomial to use, see `KSPChebyshevKind`, one of `KSP_CHEBYSHEV_FIRST`, `KSP_CHEBYSHEV_FOURTH`, or `KSP_CHEBYSHEV_OPT_FOURTH`
286: Options Database Key:
287: . -ksp_chebyshev_kind (first|fourth|opt_fourth) - which kind of Chebyshev polynomial to use
289: Level: intermediate
291: Note:
292: When using multigrid methods for problems with a poor quality coarse space (e.g., due to anisotropy or aggressive
293: coarsening), it is necessary for the smoother to handle smaller eigenvalues. With first-kind Chebyshev smoothing, this
294: requires using higher degree Chebyhev polynomials and reducing the lower end of the target spectrum, at which point
295: the whole target spectrum experiences about the same damping. Fourth kind Chebyshev polynomials (and the "optimized"
296: fourth kind) avoid the ad-hoc choice of lower bound and extend smoothing to smaller eigenvalues while preferentially
297: smoothing higher modes faster as needed to minimize the energy norm of the error. {cite}`phillips2022optimal`, {cite}`lottes2023optimal`
299: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevKind`, `KSPChebyshevGetKind()`, `KSP_CHEBYSHEV_FIRST`, `KSP_CHEBYSHEV_FOURTH`, `KSP_CHEBYSHEV_OPT_FOURTH`
300: @*/
301: PetscErrorCode KSPChebyshevSetKind(KSP ksp, KSPChebyshevKind kind)
302: {
303: PetscFunctionBegin;
306: PetscUseMethod(ksp, "KSPChebyshevSetKind_C", (KSP, KSPChebyshevKind), (ksp, kind));
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: /*@
311: KSPChebyshevGetKind - get the kind of Chebyshev polynomial to use
313: Logically Collective
315: Input Parameters:
316: + ksp - Linear solver context
317: - kind - The kind of Chebyshev polynomial used
319: Level: intermediate
321: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevKind`, `KSPChebyshevSetKind()`, `KSP_CHEBYSHEV_FIRST`, `KSP_CHEBYSHEV_FOURTH`, `KSP_CHEBYSHEV_OPT_FOURTH`
322: @*/
323: PetscErrorCode KSPChebyshevGetKind(KSP ksp, KSPChebyshevKind *kind)
324: {
325: PetscFunctionBegin;
327: PetscUseMethod(ksp, "KSPChebyshevGetKind_C", (KSP, KSPChebyshevKind *), (ksp, kind));
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
332: {
333: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
335: PetscFunctionBegin;
336: *kspest = cheb->kspest;
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: static PetscErrorCode KSPSetFromOptions_Chebyshev(KSP ksp, PetscOptionItems PetscOptionsObject)
341: {
342: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
343: PetscInt neigarg = 2, nestarg = 4;
344: PetscReal eminmax[2] = {0., 0.};
345: PetscReal tform[4] = {PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE};
346: PetscBool flgeig, flgest;
348: PetscFunctionBegin;
349: PetscOptionsHeadBegin(PetscOptionsObject, "KSP Chebyshev Options");
350: PetscCall(PetscOptionsInt("-ksp_chebyshev_esteig_steps", "Number of est steps in Chebyshev", "", cheb->eststeps, &cheb->eststeps, NULL));
351: PetscCall(PetscOptionsRealArray("-ksp_chebyshev_eigenvalues", "extreme eigenvalues", "KSPChebyshevSetEigenvalues", eminmax, &neigarg, &flgeig));
352: if (flgeig) {
353: PetscCheck(neigarg == 2, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
354: PetscCall(KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]));
355: }
356: PetscCall(PetscOptionsRealArray("-ksp_chebyshev_esteig", "estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds", "KSPChebyshevEstEigSet", tform, &nestarg, &flgest));
357: if (flgest) {
358: switch (nestarg) {
359: case 0:
360: PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
361: break;
362: case 2: /* Base everything on the max eigenvalues */
363: PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, tform[0], PETSC_DECIDE, tform[1]));
364: break;
365: case 4: /* Use the full 2x2 linear transformation */
366: PetscCall(KSPChebyshevEstEigSet(ksp, tform[0], tform[1], tform[2], tform[3]));
367: break;
368: default:
369: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
370: }
371: }
373: cheb->chebykind = KSP_CHEBYSHEV_FIRST; /* Default to 1st-kind Chebyshev polynomial */
374: PetscCall(PetscOptionsEnum("-ksp_chebyshev_kind", "Type of Chebyshev polynomial", "KSPChebyshevKind", KSPChebyshevKinds, (PetscEnum)cheb->chebykind, (PetscEnum *)&cheb->chebykind, NULL));
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) PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
379: if (cheb->kspest) {
380: PetscCall(PetscOptionsBool("-ksp_chebyshev_esteig_noisy", "Use noisy random number generated right-hand side for estimate", "KSPChebyshevEstEigSetUseNoisy", cheb->usenoisy, &cheb->usenoisy, NULL));
381: PetscCall(KSPSetFromOptions(cheb->kspest));
382: }
383: PetscOptionsHeadEnd();
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: static PetscErrorCode KSPSolve_Chebyshev_FirstKind(KSP ksp)
388: {
389: PetscInt k, kp1, km1, ktmp, i;
390: PetscScalar alpha, omegaprod, mu, omega, Gamma, c[3], scale;
391: PetscReal rnorm = 0.0, emax, emin;
392: Vec sol_orig, b, p[3], r;
393: Mat Amat, Pmat;
395: PetscFunctionBegin;
396: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
397: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
398: ksp->its = 0;
399: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
400: /* These three point to the three active solutions, we
401: rotate these three at each solution update */
402: km1 = 0;
403: k = 1;
404: kp1 = 2;
405: sol_orig = ksp->vec_sol; /* ksp->vec_sol will be assigned to rotating vector p[k], thus save its address */
406: b = ksp->vec_rhs;
407: p[km1] = sol_orig;
408: p[k] = ksp->work[0];
409: p[kp1] = ksp->work[1];
410: r = ksp->work[2];
412: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
413: /* use scale*B as our preconditioner */
414: scale = 2.0 / (emax + emin);
416: /* -alpha <= scale*lambda(B^{-1}A) <= alpha */
417: alpha = 1.0 - scale * emin;
418: Gamma = 1.0;
419: mu = 1.0 / alpha;
420: omegaprod = 2.0 / alpha;
422: c[km1] = 1.0;
423: c[k] = mu;
425: if (!ksp->guess_zero) {
426: PetscCall(KSP_MatMult(ksp, Amat, sol_orig, r)); /* r = b - A*p[km1] */
427: PetscCall(VecAYPX(r, -1.0, b));
428: } else {
429: PetscCall(VecCopy(b, r));
430: }
432: /* calculate residual norm if requested, we have done one iteration */
433: if (ksp->normtype) {
434: switch (ksp->normtype) {
435: case KSP_NORM_PRECONDITIONED:
436: PetscCall(KSP_PCApply(ksp, r, p[k])); /* p[k] = B^{-1}r */
437: PetscCall(VecNorm(p[k], NORM_2, &rnorm));
438: break;
439: case KSP_NORM_UNPRECONDITIONED:
440: case KSP_NORM_NATURAL:
441: PetscCall(VecNorm(r, NORM_2, &rnorm));
442: break;
443: default:
444: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
445: }
446: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
447: ksp->rnorm = rnorm;
448: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
449: PetscCall(KSPLogResidualHistory(ksp, rnorm));
450: PetscCall(KSPLogErrorHistory(ksp));
451: PetscCall(KSPMonitor(ksp, 0, rnorm));
452: PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
453: } else ksp->reason = KSP_CONVERGED_ITERATING;
454: if (ksp->reason || ksp->max_it == 0) {
455: if (ksp->max_it == 0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
456: PetscFunctionReturn(PETSC_SUCCESS);
457: }
458: if (ksp->normtype != KSP_NORM_PRECONDITIONED) PetscCall(KSP_PCApply(ksp, r, p[k])); /* p[k] = B^{-1}r */
459: PetscCall(VecAYPX(p[k], scale, p[km1])); /* p[k] = scale B^{-1}r + p[km1] */
460: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
461: ksp->its = 1;
462: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
464: for (i = 1; i < ksp->max_it; i++) {
465: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
466: ksp->its++;
467: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
469: PetscCall(KSP_MatMult(ksp, Amat, p[k], r)); /* r = b - Ap[k] */
470: PetscCall(VecAYPX(r, -1.0, b));
471: /* calculate residual norm if requested */
472: if (ksp->normtype) {
473: switch (ksp->normtype) {
474: case KSP_NORM_PRECONDITIONED:
475: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
476: PetscCall(VecNorm(p[kp1], NORM_2, &rnorm));
477: break;
478: case KSP_NORM_UNPRECONDITIONED:
479: case KSP_NORM_NATURAL:
480: PetscCall(VecNorm(r, NORM_2, &rnorm));
481: break;
482: default:
483: rnorm = 0.0;
484: break;
485: }
486: KSPCheckNorm(ksp, rnorm);
487: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
488: ksp->rnorm = rnorm;
489: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
490: PetscCall(KSPLogResidualHistory(ksp, rnorm));
491: PetscCall(KSPMonitor(ksp, i, rnorm));
492: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
493: if (ksp->reason) break;
494: if (ksp->normtype != KSP_NORM_PRECONDITIONED) PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
495: } else {
496: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
497: }
498: ksp->vec_sol = p[k];
499: PetscCall(KSPLogErrorHistory(ksp));
501: c[kp1] = 2.0 * mu * c[k] - c[km1];
502: omega = omegaprod * c[k] / c[kp1];
504: /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
505: PetscCall(VecAXPBYPCZ(p[kp1], 1.0 - omega, omega, omega * Gamma * scale, p[km1], p[k]));
507: ktmp = km1;
508: km1 = k;
509: k = kp1;
510: kp1 = ktmp;
511: }
512: if (!ksp->reason) {
513: if (ksp->normtype) {
514: PetscCall(KSP_MatMult(ksp, Amat, p[k], r)); /* r = b - Ap[k] */
515: PetscCall(VecAYPX(r, -1.0, b));
516: switch (ksp->normtype) {
517: case KSP_NORM_PRECONDITIONED:
518: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
519: PetscCall(VecNorm(p[kp1], NORM_2, &rnorm));
520: break;
521: case KSP_NORM_UNPRECONDITIONED:
522: case KSP_NORM_NATURAL:
523: PetscCall(VecNorm(r, NORM_2, &rnorm));
524: break;
525: default:
526: rnorm = 0.0;
527: break;
528: }
529: KSPCheckNorm(ksp, rnorm);
530: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
531: ksp->rnorm = rnorm;
532: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
533: PetscCall(KSPLogResidualHistory(ksp, rnorm));
534: PetscCall(KSPMonitor(ksp, i, rnorm));
535: }
536: if (ksp->its >= ksp->max_it) {
537: if (ksp->normtype != KSP_NORM_NONE) {
538: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
539: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
540: } else ksp->reason = KSP_CONVERGED_ITS;
541: }
542: }
544: /* make sure solution is in vector x */
545: ksp->vec_sol = sol_orig;
546: if (k) PetscCall(VecCopy(p[k], sol_orig));
547: if (ksp->reason == KSP_CONVERGED_ITS) PetscCall(KSPLogErrorHistory(ksp));
548: PetscFunctionReturn(PETSC_SUCCESS);
549: }
551: static PetscErrorCode KSPSolve_Chebyshev_FourthKind(KSP ksp)
552: {
553: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
554: PetscInt i;
555: PetscScalar scale, rScale, dScale;
556: PetscReal rnorm = 0.0, emax, emin;
557: Vec x, b, d, r, Br;
558: Mat Amat, Pmat;
559: PetscReal *betas = cheb->betas;
561: PetscFunctionBegin;
562: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
563: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
564: ksp->its = 0;
565: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
567: x = ksp->vec_sol;
568: b = ksp->vec_rhs;
569: r = ksp->work[0];
570: d = ksp->work[1];
571: Br = ksp->work[2];
573: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
574: /* use scale*B as our preconditioner */
575: scale = 1.0 / emax;
577: if (!ksp->guess_zero) {
578: PetscCall(KSP_MatMult(ksp, Amat, x, r)); /* r = b - A*x */
579: PetscCall(VecAYPX(r, -1.0, b));
580: } else {
581: PetscCall(VecCopy(b, r));
582: }
584: /* calculate residual norm if requested, we have done one iteration */
585: if (ksp->normtype) {
586: switch (ksp->normtype) {
587: case KSP_NORM_PRECONDITIONED:
588: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
589: PetscCall(VecNorm(Br, NORM_2, &rnorm));
590: break;
591: case KSP_NORM_UNPRECONDITIONED:
592: case KSP_NORM_NATURAL:
593: PetscCall(VecNorm(r, NORM_2, &rnorm));
594: break;
595: default:
596: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
597: }
598: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
599: ksp->rnorm = rnorm;
600: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
601: PetscCall(KSPLogResidualHistory(ksp, rnorm));
602: PetscCall(KSPLogErrorHistory(ksp));
603: PetscCall(KSPMonitor(ksp, 0, rnorm));
604: PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
605: } else ksp->reason = KSP_CONVERGED_ITERATING;
606: if (ksp->reason || ksp->max_it == 0) {
607: if (ksp->max_it == 0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
608: PetscFunctionReturn(PETSC_SUCCESS);
609: }
610: if (ksp->normtype != KSP_NORM_PRECONDITIONED) PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
611: PetscCall(VecAXPBY(d, 4.0 / 3.0 * scale, 0.0, Br)); /* d = 4/3 * scale B^{-1}r */
612: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
613: ksp->its = 1;
614: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
616: for (i = 1; i < ksp->max_it; i++) {
617: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
618: ksp->its++;
619: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
621: PetscCall(VecAXPBY(x, betas[i - 1], 1.0, d)); /* x = x + \beta_k d */
623: PetscCall(KSP_MatMult(ksp, Amat, d, Br)); /* r = r - Ad */
624: PetscCall(VecAXPBY(r, -1.0, 1.0, Br));
626: /* calculate residual norm if requested */
627: if (ksp->normtype) {
628: switch (ksp->normtype) {
629: case KSP_NORM_PRECONDITIONED:
630: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
631: PetscCall(VecNorm(Br, NORM_2, &rnorm));
632: break;
633: case KSP_NORM_UNPRECONDITIONED:
634: case KSP_NORM_NATURAL:
635: PetscCall(VecNorm(r, NORM_2, &rnorm));
636: break;
637: default:
638: rnorm = 0.0;
639: break;
640: }
641: KSPCheckNorm(ksp, rnorm);
642: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
643: ksp->rnorm = rnorm;
644: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
645: PetscCall(KSPLogResidualHistory(ksp, rnorm));
646: PetscCall(KSPMonitor(ksp, i, rnorm));
647: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
648: if (ksp->reason) break;
649: if (ksp->normtype != KSP_NORM_PRECONDITIONED) PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
650: } else {
651: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
652: }
653: PetscCall(KSPLogErrorHistory(ksp));
655: rScale = scale * (8.0 * i + 4.0) / (2.0 * i + 3.0);
656: dScale = (2.0 * i - 1.0) / (2.0 * i + 3.0);
658: /* d_k+1 = \dfrac{2k-1}{2k+3} d_k + \dfrac{8k+4}{2k+3} \dfrac{1}{\rho(SA)} Br */
659: PetscCall(VecAXPBY(d, rScale, dScale, Br));
660: }
662: /* on last pass, update solution vector */
663: PetscCall(VecAXPBY(x, betas[ksp->max_it - 1], 1.0, d)); /* x = x + d */
665: if (!ksp->reason) {
666: if (ksp->normtype) {
667: PetscCall(KSP_MatMult(ksp, Amat, x, r)); /* r = b - Ax */
668: PetscCall(VecAYPX(r, -1.0, b));
669: switch (ksp->normtype) {
670: case KSP_NORM_PRECONDITIONED:
671: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br= B^{-1}r */
672: PetscCall(VecNorm(Br, NORM_2, &rnorm));
673: break;
674: case KSP_NORM_UNPRECONDITIONED:
675: case KSP_NORM_NATURAL:
676: PetscCall(VecNorm(r, NORM_2, &rnorm));
677: break;
678: default:
679: rnorm = 0.0;
680: break;
681: }
682: KSPCheckNorm(ksp, rnorm);
683: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
684: ksp->rnorm = rnorm;
685: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
686: PetscCall(KSPLogResidualHistory(ksp, rnorm));
687: PetscCall(KSPMonitor(ksp, i, rnorm));
688: }
689: if (ksp->its >= ksp->max_it) {
690: if (ksp->normtype != KSP_NORM_NONE) {
691: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
692: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
693: } else ksp->reason = KSP_CONVERGED_ITS;
694: }
695: }
697: if (ksp->reason == KSP_CONVERGED_ITS) PetscCall(KSPLogErrorHistory(ksp));
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: static PetscErrorCode KSPView_Chebyshev(KSP ksp, PetscViewer viewer)
702: {
703: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
704: PetscBool isascii;
706: PetscFunctionBegin;
707: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
708: if (isascii) {
709: switch (cheb->chebykind) {
710: case KSP_CHEBYSHEV_FIRST:
711: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of first kind\n"));
712: break;
713: case KSP_CHEBYSHEV_FOURTH:
714: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of fourth kind\n"));
715: break;
716: case KSP_CHEBYSHEV_OPT_FOURTH:
717: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of opt. fourth kind\n"));
718: break;
719: }
720: PetscReal emax, emin;
721: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
722: PetscCall(PetscViewerASCIIPrintf(viewer, " eigenvalue targets used: min %g, max %g\n", (double)emin, (double)emax));
723: if (cheb->kspest) {
724: 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));
725: 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]));
726: PetscCall(PetscViewerASCIIPushTab(viewer));
727: PetscCall(KSPView(cheb->kspest, viewer));
728: PetscCall(PetscViewerASCIIPopTab(viewer));
729: if (cheb->usenoisy) PetscCall(PetscViewerASCIIPrintf(viewer, " estimating eigenvalues using a noisy random number generated right-hand side\n"));
730: } else if (cheb->emax_provided != 0.) {
731: 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],
732: (double)cheb->tform[3]));
733: }
734: }
735: PetscFunctionReturn(PETSC_SUCCESS);
736: }
738: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
739: {
740: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
741: PetscBool isset, flg;
743: PetscFunctionBegin;
744: switch (cheb->chebykind) {
745: case KSP_CHEBYSHEV_FIRST:
746: ksp->ops->solve = KSPSolve_Chebyshev_FirstKind;
747: break;
748: case KSP_CHEBYSHEV_FOURTH:
749: case KSP_CHEBYSHEV_OPT_FOURTH:
750: ksp->ops->solve = KSPSolve_Chebyshev_FourthKind;
751: break;
752: }
754: if (ksp->max_it > cheb->num_betas_alloc) {
755: PetscCall(PetscFree(cheb->betas));
756: PetscCall(PetscMalloc1(ksp->max_it, &cheb->betas));
757: cheb->num_betas_alloc = ksp->max_it;
758: }
760: // coefficients for 4th-kind Chebyshev
761: for (PetscInt i = 0; i < ksp->max_it; i++) cheb->betas[i] = 1.0;
763: // coefficients for optimized 4th-kind Chebyshev
764: if (cheb->chebykind == KSP_CHEBYSHEV_OPT_FOURTH) PetscCall(KSPChebyshevGetBetas_Private(ksp));
766: PetscCall(KSPSetWorkVecs(ksp, 3));
767: if (cheb->emin == 0. || cheb->emax == 0.) { // User did not specify eigenvalues
768: PC pc;
770: PetscCall(KSPGetPC(ksp, &pc));
771: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &flg));
772: if (!flg) { // Provided estimates are only relevant for Jacobi
773: cheb->emax_provided = 0;
774: cheb->emin_provided = 0;
775: }
776: /* We need to estimate eigenvalues */
777: if (!cheb->kspest) PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
778: }
779: if (cheb->kspest) {
780: Mat Pmat, Amat;
781: PetscBool amatsame, pmatsame;
782: MatState amatstate, pmatstate;
784: PetscCall(KSPGetOperators(ksp, &Amat, &Pmat));
785: PetscCall(MatIsSPDKnown(Pmat, &isset, &flg));
786: if (isset && flg) {
787: const char *prefix;
789: PetscCall(KSPGetOptionsPrefix(cheb->kspest, &prefix));
790: PetscCall(PetscOptionsHasName(NULL, prefix, "-ksp_type", &flg));
791: if (!flg) PetscCall(KSPSetType(cheb->kspest, KSPCG));
792: }
793: PetscCall(MatGetState(Amat, &amatstate));
794: PetscCall(MatGetState(Pmat, &pmatstate));
795: PetscCall(MatStateCompare(amatstate, cheb->amatstate, &amatsame));
796: PetscCall(MatStateCompare(pmatstate, cheb->pmatstate, &pmatsame));
797: if (!amatsame || !pmatsame) {
798: PetscReal max = 0.0, min = 0.0;
799: Vec B;
800: KSPConvergedReason reason;
802: PetscCall(KSPSetPC(cheb->kspest, ksp->pc));
803: if (cheb->usenoisy) {
804: B = ksp->work[1];
805: PetscCall(KSPSetNoisy_Private(Amat, B));
806: } else {
807: PetscBool change;
809: PetscCheck(ksp->vec_rhs, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Chebyshev must use a noisy random number generated right-hand side to estimate the eigenvalues when no right-hand side is available");
810: PetscCall(PCPreSolveChangeRHS(ksp->pc, &change));
811: if (change) {
812: B = ksp->work[1];
813: PetscCall(VecCopy(ksp->vec_rhs, B));
814: } else B = ksp->vec_rhs;
815: }
816: if (ksp->setfromoptionscalled && !cheb->kspest->setfromoptionscalled) PetscCall(KSPSetFromOptions(cheb->kspest));
817: PetscCall(KSPSolve(cheb->kspest, B, ksp->work[0]));
818: PetscCall(KSPGetConvergedReason(cheb->kspest, &reason));
819: if (reason == KSP_DIVERGED_ITS) {
820: PetscCall(PetscInfo(ksp, "Eigen estimator ran for prescribed number of iterations\n"));
821: } else if (reason == KSP_DIVERGED_PC_FAILED) {
822: PetscInt its;
823: PCFailedReason pcreason;
825: PetscCall(KSPGetIterationNumber(cheb->kspest, &its));
826: if (ksp->normtype == KSP_NORM_NONE) PetscCall(PCReduceFailedReason(ksp->pc));
827: PetscCall(PCGetFailedReason(ksp->pc, &pcreason));
828: ksp->reason = KSP_DIVERGED_PC_FAILED;
829: PetscCall(PetscInfo(ksp, "Eigen estimator failed: %s %s at iteration %" PetscInt_FMT "\n", KSPConvergedReasons[reason], PCFailedReasons[pcreason], its));
830: PetscFunctionReturn(PETSC_SUCCESS);
831: } else if (reason == KSP_CONVERGED_RTOL || reason == KSP_CONVERGED_ATOL) {
832: PetscCall(PetscInfo(ksp, "Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n"));
833: } else if (reason < 0) {
834: PetscCall(PetscInfo(ksp, "Eigen estimator failed %s, using estimates anyway\n", KSPConvergedReasons[reason]));
835: }
837: PetscCall(KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest, &min, &max));
838: PetscCall(KSPSetPC(cheb->kspest, NULL));
840: cheb->emin_computed = min;
841: cheb->emax_computed = max;
843: cheb->amatstate = amatstate;
844: cheb->pmatstate = pmatstate;
845: }
846: }
847: if (ksp->monitor[0] == (PetscErrorCode (*)(KSP, PetscInt, PetscReal, void *))KSPMonitorResidual && !ksp->normtype) PetscCall(KSPSetNormType(ksp, KSP_NORM_PRECONDITIONED));
848: PetscFunctionReturn(PETSC_SUCCESS);
849: }
851: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
852: {
853: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
855: PetscFunctionBegin;
856: PetscCall(PetscFree(cheb->betas));
857: PetscCall(KSPDestroy(&cheb->kspest));
858: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", NULL));
859: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", NULL));
860: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", NULL));
861: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetKind_C", NULL));
862: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevGetKind_C", NULL));
863: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", NULL));
864: PetscCall(KSPDestroyDefault(ksp));
865: PetscFunctionReturn(PETSC_SUCCESS);
866: }
868: /*MC
869: KSPCHEBYSHEV - The preconditioned Chebyshev iterative method
871: Options Database Keys:
872: + -ksp_chebyshev_eigenvalues emin,emax - set approximations to the smallest and largest eigenvalues
873: of the preconditioned operator. If these are accurate you will get much faster convergence.
874: . -ksp_chebyshev_esteig a,b,c,d - estimate eigenvalues using a Krylov method, then use this
875: transform for Chebyshev eigenvalue bounds (`KSPChebyshevEstEigSet()`)
876: . -ksp_chebyshev_esteig_steps - number of eigenvalue estimation steps
877: - -ksp_chebyshev_esteig_noisy - use a noisy random number generator to create right-hand side for eigenvalue estimator
879: Level: beginner
881: Notes:
882: The Chebyshev method requires both the matrix and preconditioner to be symmetric positive (semi) definite, but it can work
883: as a smoother in other situations
885: Only has support for left preconditioning.
887: Chebyshev is configured as a smoother by default, targeting the "upper" part of the spectrum.
889: By default this uses `KSPGMRES` to estimate the extreme eigenvalues, if the matrix is known to be SPD then it uses `KSPCG` to estimate the eigenvalues.
890: See `MatIsSPDKnown()` for how to indicate a `Mat`, matrix is SPD.
892: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
893: `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigSetUseNoisy()`,
894: `KSPRICHARDSON`, `KSPCG`, `PCMG`
895: M*/
897: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
898: {
899: KSP_Chebyshev *chebyshevP;
901: PetscFunctionBegin;
902: PetscCall(PetscNew(&chebyshevP));
904: ksp->data = (void *)chebyshevP;
905: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
906: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
907: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
908: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
910: chebyshevP->emin = 0.;
911: chebyshevP->emax = 0.;
913: chebyshevP->tform[0] = 0.0;
914: chebyshevP->tform[1] = 0.1;
915: chebyshevP->tform[2] = 0;
916: chebyshevP->tform[3] = 1.1;
917: chebyshevP->eststeps = 10;
918: chebyshevP->usenoisy = PETSC_TRUE;
919: ksp->setupnewmatrix = PETSC_TRUE;
921: ksp->ops->setup = KSPSetUp_Chebyshev;
922: ksp->ops->destroy = KSPDestroy_Chebyshev;
923: ksp->ops->buildsolution = KSPBuildSolutionDefault;
924: ksp->ops->buildresidual = KSPBuildResidualDefault;
925: ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
926: ksp->ops->view = KSPView_Chebyshev;
927: ksp->ops->reset = KSPReset_Chebyshev;
929: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", KSPChebyshevSetEigenvalues_Chebyshev));
930: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", KSPChebyshevEstEigSet_Chebyshev));
931: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", KSPChebyshevEstEigSetUseNoisy_Chebyshev));
932: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetKind_C", KSPChebyshevSetKind_Chebyshev));
933: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevGetKind_C", KSPChebyshevGetKind_Chebyshev));
934: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", KSPChebyshevEstEigGetKSP_Chebyshev));
935: PetscFunctionReturn(PETSC_SUCCESS);
936: }