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, 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, 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));
92: PetscCall(KSPSetComputeEigenvalues(cheb->kspest, PETSC_TRUE));
94: /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
95: PetscCall(KSPSetTolerances(cheb->kspest, 1.e-12, PETSC_DEFAULT, PETSC_DEFAULT, cheb->eststeps));
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
139: of the preconditioned problem.
141: Logically Collective
143: Input Parameters:
144: + ksp - the Krylov space context
145: emax - the eigenvalue maximum estimate
146: - emin - the eigenvalue minimum estimate
148: Options Database Key:
149: . -ksp_chebyshev_eigenvalues emin,emax - extreme eigenvalues
151: Notes:
152: Call `KSPChebyshevEstEigSet()` or use the option -ksp_chebyshev_esteig a,b,c,d to have the KSP
153: estimate the eigenvalues and use these estimated values automatically.
155: When `KSPCHEBYSHEV` is used as a smoother, one often wants to target a portion of the spectrum rather than the entire
156: spectrum. This function takes the range of target eigenvalues for Chebyshev, which will often slightly over-estimate
157: the largest eigenvalue of the actual operator (for safety) and greatly overestimate the smallest eigenvalue to
158: improve the smoothing properties of Chebyshev iteration on the higher frequencies in the spectrum.
160: Level: intermediate
162: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`,
163: @*/
164: PetscErrorCode KSPChebyshevSetEigenvalues(KSP ksp, PetscReal emax, PetscReal emin)
165: {
166: PetscFunctionBegin;
170: PetscTryMethod(ksp, "KSPChebyshevSetEigenvalues_C", (KSP, PetscReal, PetscReal), (ksp, emax, emin));
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: /*@
175: KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev
177: Logically Collective
179: Input Parameters:
180: + ksp - the Krylov space context
181: . a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
182: . b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
183: . c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
184: - d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
186: Options Database Key:
187: . -ksp_chebyshev_esteig a,b,c,d - estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds
189: Notes:
190: The Chebyshev bounds are set using
191: .vb
192: minbound = a*minest + b*maxest
193: maxbound = c*minest + d*maxest
194: .ve
195: The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
196: The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.
198: If 0.0 is passed for all transform arguments (a,b,c,d), eigenvalue estimation is disabled.
200: The default transform is (0,0.1; 0,1.1) which targets the "upper" part of the spectrum, as desirable for use with multigrid.
202: The eigenvalues are estimated using the Lanczo (`KSPCG`) or Arnoldi (`KSPGMRES`) process using a noisy right hand side vector.
204: Level: intermediate
206: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigSetUseNoisy()`, `KSPChebyshevEstEigGetKSP()`
207: @*/
208: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp, PetscReal a, PetscReal b, PetscReal c, PetscReal d)
209: {
210: PetscFunctionBegin;
216: PetscTryMethod(ksp, "KSPChebyshevEstEigSet_C", (KSP, PetscReal, PetscReal, PetscReal, PetscReal), (ksp, a, b, c, d));
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: /*@
221: KSPChebyshevEstEigSetUseNoisy - use a noisy right hand side in order to do the estimate instead of the given right hand side
223: Logically Collective
225: Input Parameters:
226: + ksp - linear solver context
227: - use - `PETSC_TRUE` to use noisy
229: Options Database Key:
230: . -ksp_chebyshev_esteig_noisy <true,false> - Use noisy right hand side for estimate
232: Note:
233: This allegedly works better for multigrid smoothers
235: Level: intermediate
237: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigGetKSP()`
238: @*/
239: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp, PetscBool use)
240: {
241: PetscFunctionBegin;
244: PetscTryMethod(ksp, "KSPChebyshevEstEigSetUseNoisy_C", (KSP, PetscBool), (ksp, use));
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: /*@
249: KSPChebyshevEstEigGetKSP - Get the Krylov method context used to estimate eigenvalues for the Chebyshev method. If
250: a Krylov method is not being used for this purpose, NULL is returned. The reference count of the returned `KSP` is
251: not incremented: it should not be destroyed by the user.
253: Input Parameter:
254: . ksp - the Krylov space context
256: Output Parameter:
257: . kspest - the eigenvalue estimation Krylov space context
259: Level: advanced
261: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`
262: @*/
263: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
264: {
265: PetscFunctionBegin;
268: *kspest = NULL;
269: PetscTryMethod(ksp, "KSPChebyshevEstEigGetKSP_C", (KSP, KSP *), (ksp, kspest));
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: /*@
274: KSPChebyshevSetKind - set the kind of Chebyshev polynomial to use
276: Logically Collective
278: Input Parameters:
279: + ksp - Linear solver context
280: - kind - The kind of Chebyshev polynomial to use
282: Options Database Key:
283: . -ksp_chebyshev_kind <kind> - which kind of Chebyshev polynomial to use
285: Level: intermediate
287: Note:
288: When using multigrid methods for problems with a poor quality coarse space (e.g., due to anisotropy or aggressive
289: coarsening), it is necessary for the smoother to handle smaller eigenvalues. With first-kind Chebyshev smoothing, this
290: requires using higher degree Chebyhev polynomials and reducing the lower end of the target spectrum, at which point
291: the whole target spectrum experiences about the same damping. Fourth kind Chebyshev polynomials (and the "optimized"
292: fourth kind) avoid the ad-hoc choice of lower bound and extend smoothing to smaller eigenvalues while preferentially
293: smoothing higher modes faster as needed to minimize the energy norm of the error.
295: References:
296: + * - Malachi Phillips and Paul Fischer, Optimal Chebyshev Smoothers and One-sided V-cycles, https://arxiv.org/abs/2210.03179.
297: - * - James Lottes, Optimal Polynomial Smoothers for Multigrid V-cycles, https://arxiv.org/abs/2202.08830.
299: .seealso: [](ch_ksp), `KSPCHEBYSHEV` `KSPChebyshevKind`, `KSPChebyshevGetKind()`
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 to use
319: Level: intermediate
321: .seealso: [](ch_ksp), `KSPCHEBYSHEV` `KSPChebyshevKind`, `KSPChebyshevSetKind()`
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 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;
394: PetscBool diagonalscale;
396: PetscFunctionBegin;
397: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
398: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
400: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
401: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
402: ksp->its = 0;
403: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
404: /* These three point to the three active solutions, we
405: rotate these three at each solution update */
406: km1 = 0;
407: k = 1;
408: kp1 = 2;
409: sol_orig = ksp->vec_sol; /* ksp->vec_sol will be assigned to rotating vector p[k], thus save its address */
410: b = ksp->vec_rhs;
411: p[km1] = sol_orig;
412: p[k] = ksp->work[0];
413: p[kp1] = ksp->work[1];
414: r = ksp->work[2];
416: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
417: /* use scale*B as our preconditioner */
418: scale = 2.0 / (emax + emin);
420: /* -alpha <= scale*lambda(B^{-1}A) <= alpha */
421: alpha = 1.0 - scale * emin;
422: Gamma = 1.0;
423: mu = 1.0 / alpha;
424: omegaprod = 2.0 / alpha;
426: c[km1] = 1.0;
427: c[k] = mu;
429: if (!ksp->guess_zero) {
430: PetscCall(KSP_MatMult(ksp, Amat, sol_orig, r)); /* r = b - A*p[km1] */
431: PetscCall(VecAYPX(r, -1.0, b));
432: } else {
433: PetscCall(VecCopy(b, r));
434: }
436: /* calculate residual norm if requested, we have done one iteration */
437: if (ksp->normtype) {
438: switch (ksp->normtype) {
439: case KSP_NORM_PRECONDITIONED:
440: PetscCall(KSP_PCApply(ksp, r, p[k])); /* p[k] = B^{-1}r */
441: PetscCall(VecNorm(p[k], NORM_2, &rnorm));
442: break;
443: case KSP_NORM_UNPRECONDITIONED:
444: case KSP_NORM_NATURAL:
445: PetscCall(VecNorm(r, NORM_2, &rnorm));
446: break;
447: default:
448: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
449: }
450: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
451: ksp->rnorm = rnorm;
452: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
453: PetscCall(KSPLogResidualHistory(ksp, rnorm));
454: PetscCall(KSPLogErrorHistory(ksp));
455: PetscCall(KSPMonitor(ksp, 0, rnorm));
456: PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
457: } else ksp->reason = KSP_CONVERGED_ITERATING;
458: if (ksp->reason || ksp->max_it == 0) {
459: if (ksp->max_it == 0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }
462: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, p[k])); /* p[k] = B^{-1}r */ }
463: PetscCall(VecAYPX(p[k], scale, p[km1])); /* p[k] = scale B^{-1}r + p[km1] */
464: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
465: ksp->its = 1;
466: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
468: for (i = 1; i < ksp->max_it; i++) {
469: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
470: ksp->its++;
471: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
473: PetscCall(KSP_MatMult(ksp, Amat, p[k], r)); /* r = b - Ap[k] */
474: PetscCall(VecAYPX(r, -1.0, b));
475: /* calculate residual norm if requested */
476: if (ksp->normtype) {
477: switch (ksp->normtype) {
478: case KSP_NORM_PRECONDITIONED:
479: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
480: PetscCall(VecNorm(p[kp1], NORM_2, &rnorm));
481: break;
482: case KSP_NORM_UNPRECONDITIONED:
483: case KSP_NORM_NATURAL:
484: PetscCall(VecNorm(r, NORM_2, &rnorm));
485: break;
486: default:
487: rnorm = 0.0;
488: break;
489: }
490: KSPCheckNorm(ksp, rnorm);
491: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
492: ksp->rnorm = rnorm;
493: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
494: PetscCall(KSPLogResidualHistory(ksp, rnorm));
495: PetscCall(KSPMonitor(ksp, i, rnorm));
496: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
497: if (ksp->reason) break;
498: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */ }
499: } else {
500: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
501: }
502: ksp->vec_sol = p[k];
503: PetscCall(KSPLogErrorHistory(ksp));
505: c[kp1] = 2.0 * mu * c[k] - c[km1];
506: omega = omegaprod * c[k] / c[kp1];
508: /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
509: PetscCall(VecAXPBYPCZ(p[kp1], 1.0 - omega, omega, omega * Gamma * scale, p[km1], p[k]));
511: ktmp = km1;
512: km1 = k;
513: k = kp1;
514: kp1 = ktmp;
515: }
516: if (!ksp->reason) {
517: if (ksp->normtype) {
518: PetscCall(KSP_MatMult(ksp, Amat, p[k], r)); /* r = b - Ap[k] */
519: PetscCall(VecAYPX(r, -1.0, b));
520: switch (ksp->normtype) {
521: case KSP_NORM_PRECONDITIONED:
522: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
523: PetscCall(VecNorm(p[kp1], NORM_2, &rnorm));
524: break;
525: case KSP_NORM_UNPRECONDITIONED:
526: case KSP_NORM_NATURAL:
527: PetscCall(VecNorm(r, NORM_2, &rnorm));
528: break;
529: default:
530: rnorm = 0.0;
531: break;
532: }
533: KSPCheckNorm(ksp, rnorm);
534: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
535: ksp->rnorm = rnorm;
536: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
537: PetscCall(KSPLogResidualHistory(ksp, rnorm));
538: PetscCall(KSPMonitor(ksp, i, rnorm));
539: }
540: if (ksp->its >= ksp->max_it) {
541: if (ksp->normtype != KSP_NORM_NONE) {
542: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
543: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
544: } else ksp->reason = KSP_CONVERGED_ITS;
545: }
546: }
548: /* make sure solution is in vector x */
549: ksp->vec_sol = sol_orig;
550: if (k) PetscCall(VecCopy(p[k], sol_orig));
551: if (ksp->reason == KSP_CONVERGED_ITS) PetscCall(KSPLogErrorHistory(ksp));
552: PetscFunctionReturn(PETSC_SUCCESS);
553: }
555: static PetscErrorCode KSPSolve_Chebyshev_FourthKind(KSP ksp)
556: {
557: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
558: PetscInt i;
559: PetscScalar scale, rScale, dScale;
560: PetscReal rnorm = 0.0, emax, emin;
561: Vec x, b, d, r, Br;
562: Mat Amat, Pmat;
563: PetscBool diagonalscale;
564: PetscReal *betas = cheb->betas;
566: PetscFunctionBegin;
567: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
568: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
570: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
571: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
572: ksp->its = 0;
573: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
575: x = ksp->vec_sol;
576: b = ksp->vec_rhs;
577: r = ksp->work[0];
578: d = ksp->work[1];
579: Br = ksp->work[2];
581: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
582: /* use scale*B as our preconditioner */
583: scale = 1.0 / emax;
585: if (!ksp->guess_zero) {
586: PetscCall(KSP_MatMult(ksp, Amat, x, r)); /* r = b - A*x */
587: PetscCall(VecAYPX(r, -1.0, b));
588: } else {
589: PetscCall(VecCopy(b, r));
590: }
592: /* calculate residual norm if requested, we have done one iteration */
593: if (ksp->normtype) {
594: switch (ksp->normtype) {
595: case KSP_NORM_PRECONDITIONED:
596: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
597: PetscCall(VecNorm(Br, NORM_2, &rnorm));
598: break;
599: case KSP_NORM_UNPRECONDITIONED:
600: case KSP_NORM_NATURAL:
601: PetscCall(VecNorm(r, NORM_2, &rnorm));
602: break;
603: default:
604: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
605: }
606: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
607: ksp->rnorm = rnorm;
608: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
609: PetscCall(KSPLogResidualHistory(ksp, rnorm));
610: PetscCall(KSPLogErrorHistory(ksp));
611: PetscCall(KSPMonitor(ksp, 0, rnorm));
612: PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
613: } else ksp->reason = KSP_CONVERGED_ITERATING;
614: if (ksp->reason || ksp->max_it == 0) {
615: if (ksp->max_it == 0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
618: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */ }
619: PetscCall(VecAXPBY(d, 4.0 / 3.0 * scale, 0.0, Br)); /* d = 4/3 * scale B^{-1}r */
620: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
621: ksp->its = 1;
622: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
624: for (i = 1; i < ksp->max_it; i++) {
625: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
626: ksp->its++;
627: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
629: PetscCall(VecAXPBY(x, betas[i - 1], 1.0, d)); /* x = x + \beta_k d */
631: PetscCall(KSP_MatMult(ksp, Amat, d, Br)); /* r = r - Ad */
632: PetscCall(VecAXPBY(r, -1.0, 1.0, Br));
634: /* calculate residual norm if requested */
635: if (ksp->normtype) {
636: switch (ksp->normtype) {
637: case KSP_NORM_PRECONDITIONED:
638: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
639: PetscCall(VecNorm(Br, NORM_2, &rnorm));
640: break;
641: case KSP_NORM_UNPRECONDITIONED:
642: case KSP_NORM_NATURAL:
643: PetscCall(VecNorm(r, NORM_2, &rnorm));
644: break;
645: default:
646: rnorm = 0.0;
647: break;
648: }
649: KSPCheckNorm(ksp, rnorm);
650: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
651: ksp->rnorm = rnorm;
652: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
653: PetscCall(KSPLogResidualHistory(ksp, rnorm));
654: PetscCall(KSPMonitor(ksp, i, rnorm));
655: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
656: if (ksp->reason) break;
657: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */ }
658: } else {
659: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
660: }
661: PetscCall(KSPLogErrorHistory(ksp));
663: rScale = scale * (8.0 * i + 4.0) / (2.0 * i + 3.0);
664: dScale = (2.0 * i - 1.0) / (2.0 * i + 3.0);
666: /* d_k+1 = \dfrac{2k-1}{2k+3} d_k + \dfrac{8k+4}{2k+3} \dfrac{1}{\rho(SA)} Br */
667: PetscCall(VecAXPBY(d, rScale, dScale, Br));
668: }
670: /* on last pass, update solution vector */
671: PetscCall(VecAXPBY(x, betas[ksp->max_it - 1], 1.0, d)); /* x = x + d */
673: if (!ksp->reason) {
674: if (ksp->normtype) {
675: PetscCall(KSP_MatMult(ksp, Amat, x, r)); /* r = b - Ax */
676: PetscCall(VecAYPX(r, -1.0, b));
677: switch (ksp->normtype) {
678: case KSP_NORM_PRECONDITIONED:
679: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br= B^{-1}r */
680: PetscCall(VecNorm(Br, NORM_2, &rnorm));
681: break;
682: case KSP_NORM_UNPRECONDITIONED:
683: case KSP_NORM_NATURAL:
684: PetscCall(VecNorm(r, NORM_2, &rnorm));
685: break;
686: default:
687: rnorm = 0.0;
688: break;
689: }
690: KSPCheckNorm(ksp, rnorm);
691: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
692: ksp->rnorm = rnorm;
693: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
694: PetscCall(KSPLogResidualHistory(ksp, rnorm));
695: PetscCall(KSPMonitor(ksp, i, rnorm));
696: }
697: if (ksp->its >= ksp->max_it) {
698: if (ksp->normtype != KSP_NORM_NONE) {
699: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
700: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
701: } else ksp->reason = KSP_CONVERGED_ITS;
702: }
703: }
705: if (ksp->reason == KSP_CONVERGED_ITS) PetscCall(KSPLogErrorHistory(ksp));
706: PetscFunctionReturn(PETSC_SUCCESS);
707: }
709: static PetscErrorCode KSPView_Chebyshev(KSP ksp, PetscViewer viewer)
710: {
711: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
712: PetscBool iascii;
714: PetscFunctionBegin;
715: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
716: if (iascii) {
717: switch (cheb->chebykind) {
718: case KSP_CHEBYSHEV_FIRST:
719: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of first kind\n"));
720: break;
721: case KSP_CHEBYSHEV_FOURTH:
722: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of fourth kind\n"));
723: break;
724: case KSP_CHEBYSHEV_OPT_FOURTH:
725: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of opt. fourth kind\n"));
726: break;
727: }
728: PetscReal emax, emin;
729: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
730: PetscCall(PetscViewerASCIIPrintf(viewer, " eigenvalue targets used: min %g, max %g\n", (double)emin, (double)emax));
731: if (cheb->kspest) {
732: 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));
733: 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]));
734: PetscCall(PetscViewerASCIIPushTab(viewer));
735: PetscCall(KSPView(cheb->kspest, viewer));
736: PetscCall(PetscViewerASCIIPopTab(viewer));
737: if (cheb->usenoisy) PetscCall(PetscViewerASCIIPrintf(viewer, " estimating eigenvalues using noisy right hand side\n"));
738: } else if (cheb->emax_provided != 0.) {
739: 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],
740: (double)cheb->tform[3]));
741: }
742: }
743: PetscFunctionReturn(PETSC_SUCCESS);
744: }
746: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
747: {
748: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
749: PetscBool isset, flg;
750: Mat Pmat, Amat;
751: PetscObjectId amatid, pmatid;
752: PetscObjectState amatstate, pmatstate;
754: PetscFunctionBegin;
755: switch (cheb->chebykind) {
756: case KSP_CHEBYSHEV_FIRST:
757: ksp->ops->solve = KSPSolve_Chebyshev_FirstKind;
758: break;
759: case KSP_CHEBYSHEV_FOURTH:
760: case KSP_CHEBYSHEV_OPT_FOURTH:
761: ksp->ops->solve = KSPSolve_Chebyshev_FourthKind;
762: break;
763: }
765: if (ksp->max_it > cheb->num_betas_alloc) {
766: PetscCall(PetscFree(cheb->betas));
767: PetscCall(PetscMalloc1(ksp->max_it, &cheb->betas));
768: cheb->num_betas_alloc = ksp->max_it;
769: }
771: // coefficients for 4th-kind Chebyshev
772: for (PetscInt i = 0; i < ksp->max_it; i++) cheb->betas[i] = 1.0;
774: // coefficients for optimized 4th-kind Chebyshev
775: if (cheb->chebykind == KSP_CHEBYSHEV_OPT_FOURTH) PetscCall(KSPChebyshevGetBetas_Private(ksp));
777: PetscCall(KSPSetWorkVecs(ksp, 3));
778: if (cheb->emin == 0. || cheb->emax == 0.) { // User did not specify eigenvalues
779: PC pc;
781: PetscCall(KSPGetPC(ksp, &pc));
782: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &flg));
783: if (!flg) { // Provided estimates are only relevant for Jacobi
784: cheb->emax_provided = 0;
785: cheb->emin_provided = 0;
786: }
787: /* We need to estimate eigenvalues */
788: if (!cheb->kspest) PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
789: }
790: if (cheb->kspest) {
791: PetscCall(KSPGetOperators(ksp, &Amat, &Pmat));
792: PetscCall(MatIsSPDKnown(Pmat, &isset, &flg));
793: if (isset && flg) {
794: const char *prefix;
796: PetscCall(KSPGetOptionsPrefix(cheb->kspest, &prefix));
797: PetscCall(PetscOptionsHasName(NULL, prefix, "-ksp_type", &flg));
798: if (!flg) PetscCall(KSPSetType(cheb->kspest, KSPCG));
799: }
800: PetscCall(PetscObjectGetId((PetscObject)Amat, &amatid));
801: PetscCall(PetscObjectGetId((PetscObject)Pmat, &pmatid));
802: PetscCall(PetscObjectStateGet((PetscObject)Amat, &amatstate));
803: PetscCall(PetscObjectStateGet((PetscObject)Pmat, &pmatstate));
804: if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
805: PetscReal max = 0.0, min = 0.0;
806: Vec B;
807: KSPConvergedReason reason;
809: PetscCall(KSPSetPC(cheb->kspest, ksp->pc));
810: if (cheb->usenoisy) {
811: B = ksp->work[1];
812: PetscCall(KSPSetNoisy_Private(B));
813: } else {
814: PetscBool change;
816: 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");
817: PetscCall(PCPreSolveChangeRHS(ksp->pc, &change));
818: if (change) {
819: B = ksp->work[1];
820: PetscCall(VecCopy(ksp->vec_rhs, B));
821: } else B = ksp->vec_rhs;
822: }
823: PetscCall(KSPSolve(cheb->kspest, B, ksp->work[0]));
824: PetscCall(KSPGetConvergedReason(cheb->kspest, &reason));
825: if (reason == KSP_DIVERGED_ITS) {
826: PetscCall(PetscInfo(ksp, "Eigen estimator ran for prescribed number of iterations\n"));
827: } else if (reason == KSP_DIVERGED_PC_FAILED) {
828: PetscInt its;
829: PCFailedReason pcreason;
831: PetscCall(KSPGetIterationNumber(cheb->kspest, &its));
832: if (ksp->normtype == KSP_NORM_NONE) {
833: PetscInt sendbuf, recvbuf;
835: PetscCall(PCGetFailedReasonRank(ksp->pc, &pcreason));
836: sendbuf = (PetscInt)pcreason;
837: PetscCall(MPIU_Allreduce(&sendbuf, &recvbuf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)ksp)));
838: PetscCall(PCSetFailedReason(ksp->pc, (PCFailedReason)recvbuf));
839: }
840: PetscCall(PCGetFailedReason(ksp->pc, &pcreason));
841: ksp->reason = KSP_DIVERGED_PC_FAILED;
842: PetscCall(PetscInfo(ksp, "Eigen estimator failed: %s %s at iteration %" PetscInt_FMT, KSPConvergedReasons[reason], PCFailedReasons[pcreason], its));
843: PetscFunctionReturn(PETSC_SUCCESS);
844: } else if (reason == KSP_CONVERGED_RTOL || reason == KSP_CONVERGED_ATOL) {
845: PetscCall(PetscInfo(ksp, "Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n"));
846: } else if (reason < 0) {
847: PetscCall(PetscInfo(ksp, "Eigen estimator failed %s, using estimates anyway\n", KSPConvergedReasons[reason]));
848: }
850: PetscCall(KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest, &min, &max));
851: PetscCall(KSPSetPC(cheb->kspest, NULL));
853: cheb->emin_computed = min;
854: cheb->emax_computed = max;
856: cheb->amatid = amatid;
857: cheb->pmatid = pmatid;
858: cheb->amatstate = amatstate;
859: cheb->pmatstate = pmatstate;
860: }
861: }
862: PetscFunctionReturn(PETSC_SUCCESS);
863: }
865: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
866: {
867: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
869: PetscFunctionBegin;
870: PetscCall(PetscFree(cheb->betas));
871: PetscCall(KSPDestroy(&cheb->kspest));
872: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", NULL));
873: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", NULL));
874: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", NULL));
875: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetKind_C", NULL));
876: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevGetKind_C", NULL));
877: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", NULL));
878: PetscCall(KSPDestroyDefault(ksp));
879: PetscFunctionReturn(PETSC_SUCCESS);
880: }
882: /*MC
883: KSPCHEBYSHEV - The preconditioned Chebyshev iterative method
885: Options Database Keys:
886: + -ksp_chebyshev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
887: of the preconditioned operator. If these are accurate you will get much faster convergence.
888: . -ksp_chebyshev_esteig <a,b,c,d> - estimate eigenvalues using a Krylov method, then use this
889: transform for Chebyshev eigenvalue bounds (`KSPChebyshevEstEigSet()`)
890: . -ksp_chebyshev_esteig_steps - number of estimation steps
891: - -ksp_chebyshev_esteig_noisy - use noisy number generator to create right hand side for eigenvalue estimator
893: Level: beginner
895: Notes:
896: 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
898: Only support for left preconditioning.
900: Chebyshev is configured as a smoother by default, targeting the "upper" part of the spectrum.
902: The user should call `KSPChebyshevSetEigenvalues()` to get eigenvalue estimates.
904: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
905: `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigSetUseNoisy()`
906: `KSPRICHARDSON`, `KSPCG`, `PCMG`
907: M*/
909: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
910: {
911: KSP_Chebyshev *chebyshevP;
913: PetscFunctionBegin;
914: PetscCall(PetscNew(&chebyshevP));
916: ksp->data = (void *)chebyshevP;
917: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
918: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
919: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
920: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
922: chebyshevP->emin = 0.;
923: chebyshevP->emax = 0.;
925: chebyshevP->tform[0] = 0.0;
926: chebyshevP->tform[1] = 0.1;
927: chebyshevP->tform[2] = 0;
928: chebyshevP->tform[3] = 1.1;
929: chebyshevP->eststeps = 10;
930: chebyshevP->usenoisy = PETSC_TRUE;
931: ksp->setupnewmatrix = PETSC_TRUE;
933: ksp->ops->setup = KSPSetUp_Chebyshev;
934: ksp->ops->destroy = KSPDestroy_Chebyshev;
935: ksp->ops->buildsolution = KSPBuildSolutionDefault;
936: ksp->ops->buildresidual = KSPBuildResidualDefault;
937: ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
938: ksp->ops->view = KSPView_Chebyshev;
939: ksp->ops->reset = KSPReset_Chebyshev;
941: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", KSPChebyshevSetEigenvalues_Chebyshev));
942: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", KSPChebyshevEstEigSet_Chebyshev));
943: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", KSPChebyshevEstEigSetUseNoisy_Chebyshev));
944: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetKind_C", KSPChebyshevSetKind_Chebyshev));
945: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevGetKind_C", KSPChebyshevGetKind_Chebyshev));
946: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", KSPChebyshevEstEigGetKSP_Chebyshev));
947: PetscFunctionReturn(PETSC_SUCCESS);
948: }