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 depending on if the operator is
202: symmetric definite or not.
204: Level: intermediate
206: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `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 of the extreme eigenvalues 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: Level: intermediate
234: Note:
235: This allegedly works better for multigrid smoothers
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 the eigenvalues for the Chebyshev method.
251: Input Parameter:
252: . ksp - the Krylov space context
254: Output Parameter:
255: . kspest - the eigenvalue estimation Krylov space context
257: Level: advanced
259: Notes:
260: If a Krylov method is not being used for this purpose, `NULL` is returned. The reference count of the returned `KSP` is
261: not incremented: it should not be destroyed by the user.
263: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`
264: @*/
265: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
266: {
267: PetscFunctionBegin;
269: PetscAssertPointer(kspest, 2);
270: *kspest = NULL;
271: PetscTryMethod(ksp, "KSPChebyshevEstEigGetKSP_C", (KSP, KSP *), (ksp, kspest));
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: /*@
276: KSPChebyshevSetKind - set the kind of Chebyshev polynomial to use
278: Logically Collective
280: Input Parameters:
281: + ksp - Linear solver context
282: - kind - The kind of Chebyshev polynomial to use, see `KSPChebyshevKind`, one of `KSP_CHEBYSHEV_FIRST`, `KSP_CHEBYSHEV_FOURTH`, or `KSP_CHEBYSHEV_OPT_FOURTH`
284: Options Database Key:
285: . -ksp_chebyshev_kind <kind> - which kind of Chebyshev polynomial to use
287: Level: intermediate
289: Note:
290: When using multigrid methods for problems with a poor quality coarse space (e.g., due to anisotropy or aggressive
291: coarsening), it is necessary for the smoother to handle smaller eigenvalues. With first-kind Chebyshev smoothing, this
292: requires using higher degree Chebyhev polynomials and reducing the lower end of the target spectrum, at which point
293: the whole target spectrum experiences about the same damping. Fourth kind Chebyshev polynomials (and the "optimized"
294: fourth kind) avoid the ad-hoc choice of lower bound and extend smoothing to smaller eigenvalues while preferentially
295: smoothing higher modes faster as needed to minimize the energy norm of the error. {cite}`phillips2022optimal`, {cite}`lottes2023optimal`
297: .seealso: [](ch_ksp), `KSPCHEBYSHEV` `KSPChebyshevKind`, `KSPChebyshevGetKind()`, `KSP_CHEBYSHEV_FIRST`, `KSP_CHEBYSHEV_FOURTH`, `KSP_CHEBYSHEV_OPT_FOURTH`
298: @*/
299: PetscErrorCode KSPChebyshevSetKind(KSP ksp, KSPChebyshevKind kind)
300: {
301: PetscFunctionBegin;
304: PetscUseMethod(ksp, "KSPChebyshevSetKind_C", (KSP, KSPChebyshevKind), (ksp, kind));
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }
308: /*@
309: KSPChebyshevGetKind - get the kind of Chebyshev polynomial to use
311: Logically Collective
313: Input Parameters:
314: + ksp - Linear solver context
315: - kind - The kind of Chebyshev polynomial used
317: Level: intermediate
319: .seealso: [](ch_ksp), `KSPCHEBYSHEV` `KSPChebyshevKind`, `KSPChebyshevSetKind()`, `KSP_CHEBYSHEV_FIRST`, `KSP_CHEBYSHEV_FOURTH`, `KSP_CHEBYSHEV_OPT_FOURTH`
320: @*/
321: PetscErrorCode KSPChebyshevGetKind(KSP ksp, KSPChebyshevKind *kind)
322: {
323: PetscFunctionBegin;
325: PetscUseMethod(ksp, "KSPChebyshevGetKind_C", (KSP, KSPChebyshevKind *), (ksp, kind));
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }
329: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
330: {
331: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
333: PetscFunctionBegin;
334: *kspest = cheb->kspest;
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: static PetscErrorCode KSPSetFromOptions_Chebyshev(KSP ksp, PetscOptionItems *PetscOptionsObject)
339: {
340: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
341: PetscInt neigarg = 2, nestarg = 4;
342: PetscReal eminmax[2] = {0., 0.};
343: PetscReal tform[4] = {PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE};
344: PetscBool flgeig, flgest;
346: PetscFunctionBegin;
347: PetscOptionsHeadBegin(PetscOptionsObject, "KSP Chebyshev Options");
348: PetscCall(PetscOptionsInt("-ksp_chebyshev_esteig_steps", "Number of est steps in Chebyshev", "", cheb->eststeps, &cheb->eststeps, NULL));
349: PetscCall(PetscOptionsRealArray("-ksp_chebyshev_eigenvalues", "extreme eigenvalues", "KSPChebyshevSetEigenvalues", eminmax, &neigarg, &flgeig));
350: if (flgeig) {
351: PetscCheck(neigarg == 2, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
352: PetscCall(KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]));
353: }
354: PetscCall(PetscOptionsRealArray("-ksp_chebyshev_esteig", "estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds", "KSPChebyshevEstEigSet", tform, &nestarg, &flgest));
355: if (flgest) {
356: switch (nestarg) {
357: case 0:
358: PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
359: break;
360: case 2: /* Base everything on the max eigenvalues */
361: PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, tform[0], PETSC_DECIDE, tform[1]));
362: break;
363: case 4: /* Use the full 2x2 linear transformation */
364: PetscCall(KSPChebyshevEstEigSet(ksp, tform[0], tform[1], tform[2], tform[3]));
365: break;
366: default:
367: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
368: }
369: }
371: cheb->chebykind = KSP_CHEBYSHEV_FIRST; /* Default to 1st-kind Chebyshev polynomial */
372: PetscCall(PetscOptionsEnum("-ksp_chebyshev_kind", "Type of Chebyshev polynomial", "KSPChebyshevKind", KSPChebyshevKinds, (PetscEnum)cheb->chebykind, (PetscEnum *)&cheb->chebykind, NULL));
374: /* We need to estimate eigenvalues; need to set this here so that KSPSetFromOptions() is called on the estimator */
375: if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
377: if (cheb->kspest) {
378: PetscCall(PetscOptionsBool("-ksp_chebyshev_esteig_noisy", "Use noisy right hand side for estimate", "KSPChebyshevEstEigSetUseNoisy", cheb->usenoisy, &cheb->usenoisy, NULL));
379: PetscCall(KSPSetFromOptions(cheb->kspest));
380: }
381: PetscOptionsHeadEnd();
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: static PetscErrorCode KSPSolve_Chebyshev_FirstKind(KSP ksp)
386: {
387: PetscInt k, kp1, km1, ktmp, i;
388: PetscScalar alpha, omegaprod, mu, omega, Gamma, c[3], scale;
389: PetscReal rnorm = 0.0, emax, emin;
390: Vec sol_orig, b, p[3], r;
391: Mat Amat, Pmat;
392: PetscBool diagonalscale;
394: PetscFunctionBegin;
395: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
396: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
398: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
399: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
400: ksp->its = 0;
401: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
402: /* These three point to the three active solutions, we
403: rotate these three at each solution update */
404: km1 = 0;
405: k = 1;
406: kp1 = 2;
407: sol_orig = ksp->vec_sol; /* ksp->vec_sol will be assigned to rotating vector p[k], thus save its address */
408: b = ksp->vec_rhs;
409: p[km1] = sol_orig;
410: p[k] = ksp->work[0];
411: p[kp1] = ksp->work[1];
412: r = ksp->work[2];
414: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
415: /* use scale*B as our preconditioner */
416: scale = 2.0 / (emax + emin);
418: /* -alpha <= scale*lambda(B^{-1}A) <= alpha */
419: alpha = 1.0 - scale * emin;
420: Gamma = 1.0;
421: mu = 1.0 / alpha;
422: omegaprod = 2.0 / alpha;
424: c[km1] = 1.0;
425: c[k] = mu;
427: if (!ksp->guess_zero) {
428: PetscCall(KSP_MatMult(ksp, Amat, sol_orig, r)); /* r = b - A*p[km1] */
429: PetscCall(VecAYPX(r, -1.0, b));
430: } else {
431: PetscCall(VecCopy(b, r));
432: }
434: /* calculate residual norm if requested, we have done one iteration */
435: if (ksp->normtype) {
436: switch (ksp->normtype) {
437: case KSP_NORM_PRECONDITIONED:
438: PetscCall(KSP_PCApply(ksp, r, p[k])); /* p[k] = B^{-1}r */
439: PetscCall(VecNorm(p[k], NORM_2, &rnorm));
440: break;
441: case KSP_NORM_UNPRECONDITIONED:
442: case KSP_NORM_NATURAL:
443: PetscCall(VecNorm(r, NORM_2, &rnorm));
444: break;
445: default:
446: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
447: }
448: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
449: ksp->rnorm = rnorm;
450: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
451: PetscCall(KSPLogResidualHistory(ksp, rnorm));
452: PetscCall(KSPLogErrorHistory(ksp));
453: PetscCall(KSPMonitor(ksp, 0, rnorm));
454: PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
455: } else ksp->reason = KSP_CONVERGED_ITERATING;
456: if (ksp->reason || ksp->max_it == 0) {
457: if (ksp->max_it == 0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
460: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, p[k])); /* p[k] = B^{-1}r */ }
461: PetscCall(VecAYPX(p[k], scale, p[km1])); /* p[k] = scale B^{-1}r + p[km1] */
462: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
463: ksp->its = 1;
464: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
466: for (i = 1; i < ksp->max_it; i++) {
467: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
468: ksp->its++;
469: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
471: PetscCall(KSP_MatMult(ksp, Amat, p[k], r)); /* r = b - Ap[k] */
472: PetscCall(VecAYPX(r, -1.0, b));
473: /* calculate residual norm if requested */
474: if (ksp->normtype) {
475: switch (ksp->normtype) {
476: case KSP_NORM_PRECONDITIONED:
477: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
478: PetscCall(VecNorm(p[kp1], NORM_2, &rnorm));
479: break;
480: case KSP_NORM_UNPRECONDITIONED:
481: case KSP_NORM_NATURAL:
482: PetscCall(VecNorm(r, NORM_2, &rnorm));
483: break;
484: default:
485: rnorm = 0.0;
486: break;
487: }
488: KSPCheckNorm(ksp, rnorm);
489: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
490: ksp->rnorm = rnorm;
491: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
492: PetscCall(KSPLogResidualHistory(ksp, rnorm));
493: PetscCall(KSPMonitor(ksp, i, rnorm));
494: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
495: if (ksp->reason) break;
496: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */ }
497: } else {
498: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
499: }
500: ksp->vec_sol = p[k];
501: PetscCall(KSPLogErrorHistory(ksp));
503: c[kp1] = 2.0 * mu * c[k] - c[km1];
504: omega = omegaprod * c[k] / c[kp1];
506: /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
507: PetscCall(VecAXPBYPCZ(p[kp1], 1.0 - omega, omega, omega * Gamma * scale, p[km1], p[k]));
509: ktmp = km1;
510: km1 = k;
511: k = kp1;
512: kp1 = ktmp;
513: }
514: if (!ksp->reason) {
515: if (ksp->normtype) {
516: PetscCall(KSP_MatMult(ksp, Amat, p[k], r)); /* r = b - Ap[k] */
517: PetscCall(VecAYPX(r, -1.0, b));
518: switch (ksp->normtype) {
519: case KSP_NORM_PRECONDITIONED:
520: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
521: PetscCall(VecNorm(p[kp1], NORM_2, &rnorm));
522: break;
523: case KSP_NORM_UNPRECONDITIONED:
524: case KSP_NORM_NATURAL:
525: PetscCall(VecNorm(r, NORM_2, &rnorm));
526: break;
527: default:
528: rnorm = 0.0;
529: break;
530: }
531: KSPCheckNorm(ksp, rnorm);
532: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
533: ksp->rnorm = rnorm;
534: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
535: PetscCall(KSPLogResidualHistory(ksp, rnorm));
536: PetscCall(KSPMonitor(ksp, i, rnorm));
537: }
538: if (ksp->its >= ksp->max_it) {
539: if (ksp->normtype != KSP_NORM_NONE) {
540: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
541: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
542: } else ksp->reason = KSP_CONVERGED_ITS;
543: }
544: }
546: /* make sure solution is in vector x */
547: ksp->vec_sol = sol_orig;
548: if (k) PetscCall(VecCopy(p[k], sol_orig));
549: if (ksp->reason == KSP_CONVERGED_ITS) PetscCall(KSPLogErrorHistory(ksp));
550: PetscFunctionReturn(PETSC_SUCCESS);
551: }
553: static PetscErrorCode KSPSolve_Chebyshev_FourthKind(KSP ksp)
554: {
555: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
556: PetscInt i;
557: PetscScalar scale, rScale, dScale;
558: PetscReal rnorm = 0.0, emax, emin;
559: Vec x, b, d, r, Br;
560: Mat Amat, Pmat;
561: PetscBool diagonalscale;
562: PetscReal *betas = cheb->betas;
564: PetscFunctionBegin;
565: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
566: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
568: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
569: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
570: ksp->its = 0;
571: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
573: x = ksp->vec_sol;
574: b = ksp->vec_rhs;
575: r = ksp->work[0];
576: d = ksp->work[1];
577: Br = ksp->work[2];
579: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
580: /* use scale*B as our preconditioner */
581: scale = 1.0 / emax;
583: if (!ksp->guess_zero) {
584: PetscCall(KSP_MatMult(ksp, Amat, x, r)); /* r = b - A*x */
585: PetscCall(VecAYPX(r, -1.0, b));
586: } else {
587: PetscCall(VecCopy(b, r));
588: }
590: /* calculate residual norm if requested, we have done one iteration */
591: if (ksp->normtype) {
592: switch (ksp->normtype) {
593: case KSP_NORM_PRECONDITIONED:
594: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
595: PetscCall(VecNorm(Br, NORM_2, &rnorm));
596: break;
597: case KSP_NORM_UNPRECONDITIONED:
598: case KSP_NORM_NATURAL:
599: PetscCall(VecNorm(r, NORM_2, &rnorm));
600: break;
601: default:
602: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
603: }
604: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
605: ksp->rnorm = rnorm;
606: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
607: PetscCall(KSPLogResidualHistory(ksp, rnorm));
608: PetscCall(KSPLogErrorHistory(ksp));
609: PetscCall(KSPMonitor(ksp, 0, rnorm));
610: PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
611: } else ksp->reason = KSP_CONVERGED_ITERATING;
612: if (ksp->reason || ksp->max_it == 0) {
613: if (ksp->max_it == 0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
616: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */ }
617: PetscCall(VecAXPBY(d, 4.0 / 3.0 * scale, 0.0, Br)); /* d = 4/3 * scale B^{-1}r */
618: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
619: ksp->its = 1;
620: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
622: for (i = 1; i < ksp->max_it; i++) {
623: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
624: ksp->its++;
625: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
627: PetscCall(VecAXPBY(x, betas[i - 1], 1.0, d)); /* x = x + \beta_k d */
629: PetscCall(KSP_MatMult(ksp, Amat, d, Br)); /* r = r - Ad */
630: PetscCall(VecAXPBY(r, -1.0, 1.0, Br));
632: /* calculate residual norm if requested */
633: if (ksp->normtype) {
634: switch (ksp->normtype) {
635: case KSP_NORM_PRECONDITIONED:
636: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
637: PetscCall(VecNorm(Br, NORM_2, &rnorm));
638: break;
639: case KSP_NORM_UNPRECONDITIONED:
640: case KSP_NORM_NATURAL:
641: PetscCall(VecNorm(r, NORM_2, &rnorm));
642: break;
643: default:
644: rnorm = 0.0;
645: break;
646: }
647: KSPCheckNorm(ksp, rnorm);
648: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
649: ksp->rnorm = rnorm;
650: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
651: PetscCall(KSPLogResidualHistory(ksp, rnorm));
652: PetscCall(KSPMonitor(ksp, i, rnorm));
653: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
654: if (ksp->reason) break;
655: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */ }
656: } else {
657: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
658: }
659: PetscCall(KSPLogErrorHistory(ksp));
661: rScale = scale * (8.0 * i + 4.0) / (2.0 * i + 3.0);
662: dScale = (2.0 * i - 1.0) / (2.0 * i + 3.0);
664: /* d_k+1 = \dfrac{2k-1}{2k+3} d_k + \dfrac{8k+4}{2k+3} \dfrac{1}{\rho(SA)} Br */
665: PetscCall(VecAXPBY(d, rScale, dScale, Br));
666: }
668: /* on last pass, update solution vector */
669: PetscCall(VecAXPBY(x, betas[ksp->max_it - 1], 1.0, d)); /* x = x + d */
671: if (!ksp->reason) {
672: if (ksp->normtype) {
673: PetscCall(KSP_MatMult(ksp, Amat, x, r)); /* r = b - Ax */
674: PetscCall(VecAYPX(r, -1.0, b));
675: switch (ksp->normtype) {
676: case KSP_NORM_PRECONDITIONED:
677: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br= B^{-1}r */
678: PetscCall(VecNorm(Br, NORM_2, &rnorm));
679: break;
680: case KSP_NORM_UNPRECONDITIONED:
681: case KSP_NORM_NATURAL:
682: PetscCall(VecNorm(r, NORM_2, &rnorm));
683: break;
684: default:
685: rnorm = 0.0;
686: break;
687: }
688: KSPCheckNorm(ksp, rnorm);
689: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
690: ksp->rnorm = rnorm;
691: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
692: PetscCall(KSPLogResidualHistory(ksp, rnorm));
693: PetscCall(KSPMonitor(ksp, i, rnorm));
694: }
695: if (ksp->its >= ksp->max_it) {
696: if (ksp->normtype != KSP_NORM_NONE) {
697: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
698: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
699: } else ksp->reason = KSP_CONVERGED_ITS;
700: }
701: }
703: if (ksp->reason == KSP_CONVERGED_ITS) PetscCall(KSPLogErrorHistory(ksp));
704: PetscFunctionReturn(PETSC_SUCCESS);
705: }
707: static PetscErrorCode KSPView_Chebyshev(KSP ksp, PetscViewer viewer)
708: {
709: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
710: PetscBool iascii;
712: PetscFunctionBegin;
713: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
714: if (iascii) {
715: switch (cheb->chebykind) {
716: case KSP_CHEBYSHEV_FIRST:
717: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of first kind\n"));
718: break;
719: case KSP_CHEBYSHEV_FOURTH:
720: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of fourth kind\n"));
721: break;
722: case KSP_CHEBYSHEV_OPT_FOURTH:
723: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of opt. fourth kind\n"));
724: break;
725: }
726: PetscReal emax, emin;
727: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
728: PetscCall(PetscViewerASCIIPrintf(viewer, " eigenvalue targets used: min %g, max %g\n", (double)emin, (double)emax));
729: if (cheb->kspest) {
730: 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));
731: 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]));
732: PetscCall(PetscViewerASCIIPushTab(viewer));
733: PetscCall(KSPView(cheb->kspest, viewer));
734: PetscCall(PetscViewerASCIIPopTab(viewer));
735: if (cheb->usenoisy) PetscCall(PetscViewerASCIIPrintf(viewer, " estimating eigenvalues using noisy right hand side\n"));
736: } else if (cheb->emax_provided != 0.) {
737: 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],
738: (double)cheb->tform[3]));
739: }
740: }
741: PetscFunctionReturn(PETSC_SUCCESS);
742: }
744: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
745: {
746: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
747: PetscBool isset, flg;
748: Mat Pmat, Amat;
749: PetscObjectId amatid, pmatid;
750: PetscObjectState amatstate, pmatstate;
752: PetscFunctionBegin;
753: switch (cheb->chebykind) {
754: case KSP_CHEBYSHEV_FIRST:
755: ksp->ops->solve = KSPSolve_Chebyshev_FirstKind;
756: break;
757: case KSP_CHEBYSHEV_FOURTH:
758: case KSP_CHEBYSHEV_OPT_FOURTH:
759: ksp->ops->solve = KSPSolve_Chebyshev_FourthKind;
760: break;
761: }
763: if (ksp->max_it > cheb->num_betas_alloc) {
764: PetscCall(PetscFree(cheb->betas));
765: PetscCall(PetscMalloc1(ksp->max_it, &cheb->betas));
766: cheb->num_betas_alloc = ksp->max_it;
767: }
769: // coefficients for 4th-kind Chebyshev
770: for (PetscInt i = 0; i < ksp->max_it; i++) cheb->betas[i] = 1.0;
772: // coefficients for optimized 4th-kind Chebyshev
773: if (cheb->chebykind == KSP_CHEBYSHEV_OPT_FOURTH) PetscCall(KSPChebyshevGetBetas_Private(ksp));
775: PetscCall(KSPSetWorkVecs(ksp, 3));
776: if (cheb->emin == 0. || cheb->emax == 0.) { // User did not specify eigenvalues
777: PC pc;
779: PetscCall(KSPGetPC(ksp, &pc));
780: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &flg));
781: if (!flg) { // Provided estimates are only relevant for Jacobi
782: cheb->emax_provided = 0;
783: cheb->emin_provided = 0;
784: }
785: /* We need to estimate eigenvalues */
786: if (!cheb->kspest) PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
787: }
788: if (cheb->kspest) {
789: PetscCall(KSPGetOperators(ksp, &Amat, &Pmat));
790: PetscCall(MatIsSPDKnown(Pmat, &isset, &flg));
791: if (isset && flg) {
792: const char *prefix;
794: PetscCall(KSPGetOptionsPrefix(cheb->kspest, &prefix));
795: PetscCall(PetscOptionsHasName(NULL, prefix, "-ksp_type", &flg));
796: if (!flg) PetscCall(KSPSetType(cheb->kspest, KSPCG));
797: }
798: PetscCall(PetscObjectGetId((PetscObject)Amat, &amatid));
799: PetscCall(PetscObjectGetId((PetscObject)Pmat, &pmatid));
800: PetscCall(PetscObjectStateGet((PetscObject)Amat, &amatstate));
801: PetscCall(PetscObjectStateGet((PetscObject)Pmat, &pmatstate));
802: if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
803: PetscReal max = 0.0, min = 0.0;
804: Vec B;
805: KSPConvergedReason reason;
807: PetscCall(KSPSetPC(cheb->kspest, ksp->pc));
808: if (cheb->usenoisy) {
809: B = ksp->work[1];
810: PetscCall(KSPSetNoisy_Private(B));
811: } else {
812: PetscBool change;
814: 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");
815: PetscCall(PCPreSolveChangeRHS(ksp->pc, &change));
816: if (change) {
817: B = ksp->work[1];
818: PetscCall(VecCopy(ksp->vec_rhs, B));
819: } else B = ksp->vec_rhs;
820: }
821: if (ksp->setfromoptionscalled && !cheb->kspest->setfromoptionscalled) PetscCall(KSPSetFromOptions(cheb->kspest));
822: PetscCall(KSPSolve(cheb->kspest, B, ksp->work[0]));
823: PetscCall(KSPGetConvergedReason(cheb->kspest, &reason));
824: if (reason == KSP_DIVERGED_ITS) {
825: PetscCall(PetscInfo(ksp, "Eigen estimator ran for prescribed number of iterations\n"));
826: } else if (reason == KSP_DIVERGED_PC_FAILED) {
827: PetscInt its;
828: PCFailedReason pcreason;
830: PetscCall(KSPGetIterationNumber(cheb->kspest, &its));
831: if (ksp->normtype == KSP_NORM_NONE) PetscCall(PCReduceFailedReason(ksp->pc));
832: PetscCall(PCGetFailedReason(ksp->pc, &pcreason));
833: ksp->reason = KSP_DIVERGED_PC_FAILED;
834: PetscCall(PetscInfo(ksp, "Eigen estimator failed: %s %s at iteration %" PetscInt_FMT "\n", KSPConvergedReasons[reason], PCFailedReasons[pcreason], its));
835: PetscFunctionReturn(PETSC_SUCCESS);
836: } else if (reason == KSP_CONVERGED_RTOL || reason == KSP_CONVERGED_ATOL) {
837: PetscCall(PetscInfo(ksp, "Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n"));
838: } else if (reason < 0) {
839: PetscCall(PetscInfo(ksp, "Eigen estimator failed %s, using estimates anyway\n", KSPConvergedReasons[reason]));
840: }
842: PetscCall(KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest, &min, &max));
843: PetscCall(KSPSetPC(cheb->kspest, NULL));
845: cheb->emin_computed = min;
846: cheb->emax_computed = max;
848: cheb->amatid = amatid;
849: cheb->pmatid = pmatid;
850: cheb->amatstate = amatstate;
851: cheb->pmatstate = pmatstate;
852: }
853: }
854: PetscFunctionReturn(PETSC_SUCCESS);
855: }
857: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
858: {
859: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
861: PetscFunctionBegin;
862: PetscCall(PetscFree(cheb->betas));
863: PetscCall(KSPDestroy(&cheb->kspest));
864: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", NULL));
865: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", NULL));
866: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", NULL));
867: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetKind_C", NULL));
868: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevGetKind_C", NULL));
869: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", NULL));
870: PetscCall(KSPDestroyDefault(ksp));
871: PetscFunctionReturn(PETSC_SUCCESS);
872: }
874: /*MC
875: KSPCHEBYSHEV - The preconditioned Chebyshev iterative method
877: Options Database Keys:
878: + -ksp_chebyshev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
879: of the preconditioned operator. If these are accurate you will get much faster convergence.
880: . -ksp_chebyshev_esteig <a,b,c,d> - estimate eigenvalues using a Krylov method, then use this
881: transform for Chebyshev eigenvalue bounds (`KSPChebyshevEstEigSet()`)
882: . -ksp_chebyshev_esteig_steps - number of estimation steps
883: - -ksp_chebyshev_esteig_noisy - use noisy number generator to create right hand side for eigenvalue estimator
885: Level: beginner
887: Notes:
888: The Chebyshev method requires both the matrix and preconditioner to be symmetric positive (semi) definite, but it can work
889: as a smoother in other situations
891: Only support for left preconditioning.
893: Chebyshev is configured as a smoother by default, targeting the "upper" part of the spectrum.
895: The user should call `KSPChebyshevSetEigenvalues()` to get eigenvalue estimates.
897: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
898: `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigSetUseNoisy()`
899: `KSPRICHARDSON`, `KSPCG`, `PCMG`
900: M*/
902: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
903: {
904: KSP_Chebyshev *chebyshevP;
906: PetscFunctionBegin;
907: PetscCall(PetscNew(&chebyshevP));
909: ksp->data = (void *)chebyshevP;
910: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
911: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
912: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
913: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
915: chebyshevP->emin = 0.;
916: chebyshevP->emax = 0.;
918: chebyshevP->tform[0] = 0.0;
919: chebyshevP->tform[1] = 0.1;
920: chebyshevP->tform[2] = 0;
921: chebyshevP->tform[3] = 1.1;
922: chebyshevP->eststeps = 10;
923: chebyshevP->usenoisy = PETSC_TRUE;
924: ksp->setupnewmatrix = PETSC_TRUE;
926: ksp->ops->setup = KSPSetUp_Chebyshev;
927: ksp->ops->destroy = KSPDestroy_Chebyshev;
928: ksp->ops->buildsolution = KSPBuildSolutionDefault;
929: ksp->ops->buildresidual = KSPBuildResidualDefault;
930: ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
931: ksp->ops->view = KSPView_Chebyshev;
932: ksp->ops->reset = KSPReset_Chebyshev;
934: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", KSPChebyshevSetEigenvalues_Chebyshev));
935: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", KSPChebyshevEstEigSet_Chebyshev));
936: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", KSPChebyshevEstEigSetUseNoisy_Chebyshev));
937: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetKind_C", KSPChebyshevSetKind_Chebyshev));
938: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevGetKind_C", KSPChebyshevGetKind_Chebyshev));
939: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", KSPChebyshevEstEigGetKSP_Chebyshev));
940: PetscFunctionReturn(PETSC_SUCCESS);
941: }