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: cheb->amatid = 0;
106: cheb->pmatid = 0;
107: cheb->amatstate = -1;
108: cheb->pmatstate = -1;
109: } else {
110: PetscCall(KSPDestroy(&cheb->kspest));
111: }
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp, PetscBool use)
116: {
117: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
119: PetscFunctionBegin;
120: cheb->usenoisy = use;
121: PetscFunctionReturn(PETSC_SUCCESS);
122: }
124: static PetscErrorCode KSPChebyshevSetKind_Chebyshev(KSP ksp, KSPChebyshevKind kind)
125: {
126: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
128: PetscFunctionBegin;
129: cheb->chebykind = kind;
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: static PetscErrorCode KSPChebyshevGetKind_Chebyshev(KSP ksp, KSPChebyshevKind *kind)
134: {
135: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
137: PetscFunctionBegin;
138: *kind = cheb->chebykind;
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
141: /*@
142: KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues of the preconditioned problem.
144: Logically Collective
146: Input Parameters:
147: + ksp - the Krylov space context
148: . emax - the eigenvalue maximum estimate
149: - emin - the eigenvalue minimum estimate
151: Options Database Key:
152: . -ksp_chebyshev_eigenvalues emin,emax - extreme eigenvalues
154: Level: intermediate
156: Notes:
157: Call `KSPChebyshevEstEigSet()` or use the option `-ksp_chebyshev_esteig a,b,c,d` to have the `KSP`
158: estimate the eigenvalues and use these estimated values automatically.
160: When `KSPCHEBYSHEV` is used as a smoother, one often wants to target a portion of the spectrum rather than the entire
161: spectrum. This function takes the range of target eigenvalues for Chebyshev, which will often slightly over-estimate
162: the largest eigenvalue of the actual operator (for safety) and greatly overestimate the smallest eigenvalue to
163: improve the smoothing properties of Chebyshev iteration on the higher frequencies in the spectrum.
165: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`,
166: @*/
167: PetscErrorCode KSPChebyshevSetEigenvalues(KSP ksp, PetscReal emax, PetscReal emin)
168: {
169: PetscFunctionBegin;
173: PetscTryMethod(ksp, "KSPChebyshevSetEigenvalues_C", (KSP, PetscReal, PetscReal), (ksp, emax, emin));
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: /*@
178: KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev
180: Logically Collective
182: Input Parameters:
183: + ksp - the Krylov space context
184: . a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or `PETSC_DECIDE`)
185: . b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or `PETSC_DECIDE`)
186: . c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or `PETSC_DECIDE`)
187: - d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or `PETSC_DECIDE`)
189: Options Database Key:
190: . -ksp_chebyshev_esteig a,b,c,d - estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds
192: Notes:
193: The Chebyshev bounds are set using
194: .vb
195: minbound = a*minest + b*maxest
196: maxbound = c*minest + d*maxest
197: .ve
198: The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
199: The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.
201: If 0.0 is passed for all transform arguments (a,b,c,d), eigenvalue estimation is disabled.
203: The default transform is (0,0.1; 0,1.1) which targets the "upper" part of the spectrum, as desirable for use with multigrid.
205: The eigenvalues are estimated using the Lanczos (`KSPCG`) or Arnoldi (`KSPGMRES`) process depending on if the operator is
206: symmetric definite or not.
208: Level: intermediate
210: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSetUseNoisy()`, `KSPChebyshevEstEigGetKSP()`
211: @*/
212: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp, PetscReal a, PetscReal b, PetscReal c, PetscReal d)
213: {
214: PetscFunctionBegin;
220: PetscTryMethod(ksp, "KSPChebyshevEstEigSet_C", (KSP, PetscReal, PetscReal, PetscReal, PetscReal), (ksp, a, b, c, d));
221: PetscFunctionReturn(PETSC_SUCCESS);
222: }
224: /*@
225: KSPChebyshevEstEigSetUseNoisy - use a noisy random number generated right-hand side to estimate the extreme eigenvalues instead of the given right-hand side
227: Logically Collective
229: Input Parameters:
230: + ksp - linear solver context
231: - use - `PETSC_TRUE` to use noisy
233: Options Database Key:
234: . -ksp_chebyshev_esteig_noisy <true,false> - Use noisy right-hand side for estimate
236: Level: intermediate
238: Note:
239: This allegedly works better for multigrid smoothers
241: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigGetKSP()`
242: @*/
243: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp, PetscBool use)
244: {
245: PetscFunctionBegin;
248: PetscTryMethod(ksp, "KSPChebyshevEstEigSetUseNoisy_C", (KSP, PetscBool), (ksp, use));
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: /*@
253: KSPChebyshevEstEigGetKSP - Get the Krylov method context used to estimate the eigenvalues for the Chebyshev method.
255: Input Parameter:
256: . ksp - the Krylov space context
258: Output Parameter:
259: . kspest - the eigenvalue estimation Krylov space context
261: Level: advanced
263: Notes:
264: If a Krylov method is not being used for this purpose, `NULL` is returned. The reference count of the returned `KSP` is
265: not incremented: it should not be destroyed by the user.
267: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`
268: @*/
269: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
270: {
271: PetscFunctionBegin;
273: PetscAssertPointer(kspest, 2);
274: *kspest = NULL;
275: PetscTryMethod(ksp, "KSPChebyshevEstEigGetKSP_C", (KSP, KSP *), (ksp, kspest));
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: /*@
280: KSPChebyshevSetKind - set the kind of Chebyshev polynomial to use
282: Logically Collective
284: Input Parameters:
285: + ksp - Linear solver context
286: - kind - The kind of Chebyshev polynomial to use, see `KSPChebyshevKind`, one of `KSP_CHEBYSHEV_FIRST`, `KSP_CHEBYSHEV_FOURTH`, or `KSP_CHEBYSHEV_OPT_FOURTH`
288: Options Database Key:
289: . -ksp_chebyshev_kind <kind> - which kind of Chebyshev polynomial to use
291: Level: intermediate
293: Note:
294: When using multigrid methods for problems with a poor quality coarse space (e.g., due to anisotropy or aggressive
295: coarsening), it is necessary for the smoother to handle smaller eigenvalues. With first-kind Chebyshev smoothing, this
296: requires using higher degree Chebyhev polynomials and reducing the lower end of the target spectrum, at which point
297: the whole target spectrum experiences about the same damping. Fourth kind Chebyshev polynomials (and the "optimized"
298: fourth kind) avoid the ad-hoc choice of lower bound and extend smoothing to smaller eigenvalues while preferentially
299: smoothing higher modes faster as needed to minimize the energy norm of the error. {cite}`phillips2022optimal`, {cite}`lottes2023optimal`
301: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevKind`, `KSPChebyshevGetKind()`, `KSP_CHEBYSHEV_FIRST`, `KSP_CHEBYSHEV_FOURTH`, `KSP_CHEBYSHEV_OPT_FOURTH`
302: @*/
303: PetscErrorCode KSPChebyshevSetKind(KSP ksp, KSPChebyshevKind kind)
304: {
305: PetscFunctionBegin;
308: PetscUseMethod(ksp, "KSPChebyshevSetKind_C", (KSP, KSPChebyshevKind), (ksp, kind));
309: PetscFunctionReturn(PETSC_SUCCESS);
310: }
312: /*@
313: KSPChebyshevGetKind - get the kind of Chebyshev polynomial to use
315: Logically Collective
317: Input Parameters:
318: + ksp - Linear solver context
319: - kind - The kind of Chebyshev polynomial used
321: Level: intermediate
323: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevKind`, `KSPChebyshevSetKind()`, `KSP_CHEBYSHEV_FIRST`, `KSP_CHEBYSHEV_FOURTH`, `KSP_CHEBYSHEV_OPT_FOURTH`
324: @*/
325: PetscErrorCode KSPChebyshevGetKind(KSP ksp, KSPChebyshevKind *kind)
326: {
327: PetscFunctionBegin;
329: PetscUseMethod(ksp, "KSPChebyshevGetKind_C", (KSP, KSPChebyshevKind *), (ksp, kind));
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
334: {
335: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
337: PetscFunctionBegin;
338: *kspest = cheb->kspest;
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: static PetscErrorCode KSPSetFromOptions_Chebyshev(KSP ksp, PetscOptionItems PetscOptionsObject)
343: {
344: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
345: PetscInt neigarg = 2, nestarg = 4;
346: PetscReal eminmax[2] = {0., 0.};
347: PetscReal tform[4] = {PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE};
348: PetscBool flgeig, flgest;
350: PetscFunctionBegin;
351: PetscOptionsHeadBegin(PetscOptionsObject, "KSP Chebyshev Options");
352: PetscCall(PetscOptionsInt("-ksp_chebyshev_esteig_steps", "Number of est steps in Chebyshev", "", cheb->eststeps, &cheb->eststeps, NULL));
353: PetscCall(PetscOptionsRealArray("-ksp_chebyshev_eigenvalues", "extreme eigenvalues", "KSPChebyshevSetEigenvalues", eminmax, &neigarg, &flgeig));
354: if (flgeig) {
355: PetscCheck(neigarg == 2, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
356: PetscCall(KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]));
357: }
358: PetscCall(PetscOptionsRealArray("-ksp_chebyshev_esteig", "estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds", "KSPChebyshevEstEigSet", tform, &nestarg, &flgest));
359: if (flgest) {
360: switch (nestarg) {
361: case 0:
362: PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
363: break;
364: case 2: /* Base everything on the max eigenvalues */
365: PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, tform[0], PETSC_DECIDE, tform[1]));
366: break;
367: case 4: /* Use the full 2x2 linear transformation */
368: PetscCall(KSPChebyshevEstEigSet(ksp, tform[0], tform[1], tform[2], tform[3]));
369: break;
370: default:
371: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
372: }
373: }
375: cheb->chebykind = KSP_CHEBYSHEV_FIRST; /* Default to 1st-kind Chebyshev polynomial */
376: PetscCall(PetscOptionsEnum("-ksp_chebyshev_kind", "Type of Chebyshev polynomial", "KSPChebyshevKind", KSPChebyshevKinds, (PetscEnum)cheb->chebykind, (PetscEnum *)&cheb->chebykind, NULL));
378: /* We need to estimate eigenvalues; need to set this here so that KSPSetFromOptions() is called on the estimator */
379: if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
381: if (cheb->kspest) {
382: PetscCall(PetscOptionsBool("-ksp_chebyshev_esteig_noisy", "Use noisy random number generated right-hand side for estimate", "KSPChebyshevEstEigSetUseNoisy", cheb->usenoisy, &cheb->usenoisy, NULL));
383: PetscCall(KSPSetFromOptions(cheb->kspest));
384: }
385: PetscOptionsHeadEnd();
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: static PetscErrorCode KSPSolve_Chebyshev_FirstKind(KSP ksp)
390: {
391: PetscInt k, kp1, km1, ktmp, i;
392: PetscScalar alpha, omegaprod, mu, omega, Gamma, c[3], scale;
393: PetscReal rnorm = 0.0, emax, emin;
394: Vec sol_orig, b, p[3], r;
395: Mat Amat, Pmat;
396: PetscBool diagonalscale;
398: PetscFunctionBegin;
399: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
400: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
402: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
403: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
404: ksp->its = 0;
405: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
406: /* These three point to the three active solutions, we
407: rotate these three at each solution update */
408: km1 = 0;
409: k = 1;
410: kp1 = 2;
411: sol_orig = ksp->vec_sol; /* ksp->vec_sol will be assigned to rotating vector p[k], thus save its address */
412: b = ksp->vec_rhs;
413: p[km1] = sol_orig;
414: p[k] = ksp->work[0];
415: p[kp1] = ksp->work[1];
416: r = ksp->work[2];
418: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
419: /* use scale*B as our preconditioner */
420: scale = 2.0 / (emax + emin);
422: /* -alpha <= scale*lambda(B^{-1}A) <= alpha */
423: alpha = 1.0 - scale * emin;
424: Gamma = 1.0;
425: mu = 1.0 / alpha;
426: omegaprod = 2.0 / alpha;
428: c[km1] = 1.0;
429: c[k] = mu;
431: if (!ksp->guess_zero) {
432: PetscCall(KSP_MatMult(ksp, Amat, sol_orig, r)); /* r = b - A*p[km1] */
433: PetscCall(VecAYPX(r, -1.0, b));
434: } else {
435: PetscCall(VecCopy(b, r));
436: }
438: /* calculate residual norm if requested, we have done one iteration */
439: if (ksp->normtype) {
440: switch (ksp->normtype) {
441: case KSP_NORM_PRECONDITIONED:
442: PetscCall(KSP_PCApply(ksp, r, p[k])); /* p[k] = B^{-1}r */
443: PetscCall(VecNorm(p[k], NORM_2, &rnorm));
444: break;
445: case KSP_NORM_UNPRECONDITIONED:
446: case KSP_NORM_NATURAL:
447: PetscCall(VecNorm(r, NORM_2, &rnorm));
448: break;
449: default:
450: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
451: }
452: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
453: ksp->rnorm = rnorm;
454: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
455: PetscCall(KSPLogResidualHistory(ksp, rnorm));
456: PetscCall(KSPLogErrorHistory(ksp));
457: PetscCall(KSPMonitor(ksp, 0, rnorm));
458: PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
459: } else ksp->reason = KSP_CONVERGED_ITERATING;
460: if (ksp->reason || ksp->max_it == 0) {
461: if (ksp->max_it == 0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
462: PetscFunctionReturn(PETSC_SUCCESS);
463: }
464: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, p[k])); /* p[k] = B^{-1}r */ }
465: PetscCall(VecAYPX(p[k], scale, p[km1])); /* p[k] = scale B^{-1}r + p[km1] */
466: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
467: ksp->its = 1;
468: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
470: for (i = 1; i < ksp->max_it; i++) {
471: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
472: ksp->its++;
473: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
475: PetscCall(KSP_MatMult(ksp, Amat, p[k], r)); /* r = b - Ap[k] */
476: PetscCall(VecAYPX(r, -1.0, b));
477: /* calculate residual norm if requested */
478: if (ksp->normtype) {
479: switch (ksp->normtype) {
480: case KSP_NORM_PRECONDITIONED:
481: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
482: PetscCall(VecNorm(p[kp1], NORM_2, &rnorm));
483: break;
484: case KSP_NORM_UNPRECONDITIONED:
485: case KSP_NORM_NATURAL:
486: PetscCall(VecNorm(r, NORM_2, &rnorm));
487: break;
488: default:
489: rnorm = 0.0;
490: break;
491: }
492: KSPCheckNorm(ksp, rnorm);
493: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
494: ksp->rnorm = rnorm;
495: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
496: PetscCall(KSPLogResidualHistory(ksp, rnorm));
497: PetscCall(KSPMonitor(ksp, i, rnorm));
498: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
499: if (ksp->reason) break;
500: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */ }
501: } else {
502: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
503: }
504: ksp->vec_sol = p[k];
505: PetscCall(KSPLogErrorHistory(ksp));
507: c[kp1] = 2.0 * mu * c[k] - c[km1];
508: omega = omegaprod * c[k] / c[kp1];
510: /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
511: PetscCall(VecAXPBYPCZ(p[kp1], 1.0 - omega, omega, omega * Gamma * scale, p[km1], p[k]));
513: ktmp = km1;
514: km1 = k;
515: k = kp1;
516: kp1 = ktmp;
517: }
518: if (!ksp->reason) {
519: if (ksp->normtype) {
520: PetscCall(KSP_MatMult(ksp, Amat, p[k], r)); /* r = b - Ap[k] */
521: PetscCall(VecAYPX(r, -1.0, b));
522: switch (ksp->normtype) {
523: case KSP_NORM_PRECONDITIONED:
524: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
525: PetscCall(VecNorm(p[kp1], NORM_2, &rnorm));
526: break;
527: case KSP_NORM_UNPRECONDITIONED:
528: case KSP_NORM_NATURAL:
529: PetscCall(VecNorm(r, NORM_2, &rnorm));
530: break;
531: default:
532: rnorm = 0.0;
533: break;
534: }
535: KSPCheckNorm(ksp, rnorm);
536: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
537: ksp->rnorm = rnorm;
538: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
539: PetscCall(KSPLogResidualHistory(ksp, rnorm));
540: PetscCall(KSPMonitor(ksp, i, rnorm));
541: }
542: if (ksp->its >= ksp->max_it) {
543: if (ksp->normtype != KSP_NORM_NONE) {
544: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
545: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
546: } else ksp->reason = KSP_CONVERGED_ITS;
547: }
548: }
550: /* make sure solution is in vector x */
551: ksp->vec_sol = sol_orig;
552: if (k) PetscCall(VecCopy(p[k], sol_orig));
553: if (ksp->reason == KSP_CONVERGED_ITS) PetscCall(KSPLogErrorHistory(ksp));
554: PetscFunctionReturn(PETSC_SUCCESS);
555: }
557: static PetscErrorCode KSPSolve_Chebyshev_FourthKind(KSP ksp)
558: {
559: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
560: PetscInt i;
561: PetscScalar scale, rScale, dScale;
562: PetscReal rnorm = 0.0, emax, emin;
563: Vec x, b, d, r, Br;
564: Mat Amat, Pmat;
565: PetscBool diagonalscale;
566: PetscReal *betas = cheb->betas;
568: PetscFunctionBegin;
569: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
570: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
572: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
573: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
574: ksp->its = 0;
575: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
577: x = ksp->vec_sol;
578: b = ksp->vec_rhs;
579: r = ksp->work[0];
580: d = ksp->work[1];
581: Br = ksp->work[2];
583: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
584: /* use scale*B as our preconditioner */
585: scale = 1.0 / emax;
587: if (!ksp->guess_zero) {
588: PetscCall(KSP_MatMult(ksp, Amat, x, r)); /* r = b - A*x */
589: PetscCall(VecAYPX(r, -1.0, b));
590: } else {
591: PetscCall(VecCopy(b, r));
592: }
594: /* calculate residual norm if requested, we have done one iteration */
595: if (ksp->normtype) {
596: switch (ksp->normtype) {
597: case KSP_NORM_PRECONDITIONED:
598: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
599: PetscCall(VecNorm(Br, NORM_2, &rnorm));
600: break;
601: case KSP_NORM_UNPRECONDITIONED:
602: case KSP_NORM_NATURAL:
603: PetscCall(VecNorm(r, NORM_2, &rnorm));
604: break;
605: default:
606: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
607: }
608: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
609: ksp->rnorm = rnorm;
610: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
611: PetscCall(KSPLogResidualHistory(ksp, rnorm));
612: PetscCall(KSPLogErrorHistory(ksp));
613: PetscCall(KSPMonitor(ksp, 0, rnorm));
614: PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
615: } else ksp->reason = KSP_CONVERGED_ITERATING;
616: if (ksp->reason || ksp->max_it == 0) {
617: if (ksp->max_it == 0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
618: PetscFunctionReturn(PETSC_SUCCESS);
619: }
620: if (ksp->normtype != KSP_NORM_PRECONDITIONED) PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
621: PetscCall(VecAXPBY(d, 4.0 / 3.0 * scale, 0.0, Br)); /* d = 4/3 * scale B^{-1}r */
622: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
623: ksp->its = 1;
624: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
626: for (i = 1; i < ksp->max_it; i++) {
627: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
628: ksp->its++;
629: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
631: PetscCall(VecAXPBY(x, betas[i - 1], 1.0, d)); /* x = x + \beta_k d */
633: PetscCall(KSP_MatMult(ksp, Amat, d, Br)); /* r = r - Ad */
634: PetscCall(VecAXPBY(r, -1.0, 1.0, Br));
636: /* calculate residual norm if requested */
637: if (ksp->normtype) {
638: switch (ksp->normtype) {
639: case KSP_NORM_PRECONDITIONED:
640: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
641: PetscCall(VecNorm(Br, NORM_2, &rnorm));
642: break;
643: case KSP_NORM_UNPRECONDITIONED:
644: case KSP_NORM_NATURAL:
645: PetscCall(VecNorm(r, NORM_2, &rnorm));
646: break;
647: default:
648: rnorm = 0.0;
649: break;
650: }
651: KSPCheckNorm(ksp, rnorm);
652: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
653: ksp->rnorm = rnorm;
654: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
655: PetscCall(KSPLogResidualHistory(ksp, rnorm));
656: PetscCall(KSPMonitor(ksp, i, rnorm));
657: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
658: if (ksp->reason) break;
659: if (ksp->normtype != KSP_NORM_PRECONDITIONED) PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
660: } else {
661: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
662: }
663: PetscCall(KSPLogErrorHistory(ksp));
665: rScale = scale * (8.0 * i + 4.0) / (2.0 * i + 3.0);
666: dScale = (2.0 * i - 1.0) / (2.0 * i + 3.0);
668: /* d_k+1 = \dfrac{2k-1}{2k+3} d_k + \dfrac{8k+4}{2k+3} \dfrac{1}{\rho(SA)} Br */
669: PetscCall(VecAXPBY(d, rScale, dScale, Br));
670: }
672: /* on last pass, update solution vector */
673: PetscCall(VecAXPBY(x, betas[ksp->max_it - 1], 1.0, d)); /* x = x + d */
675: if (!ksp->reason) {
676: if (ksp->normtype) {
677: PetscCall(KSP_MatMult(ksp, Amat, x, r)); /* r = b - Ax */
678: PetscCall(VecAYPX(r, -1.0, b));
679: switch (ksp->normtype) {
680: case KSP_NORM_PRECONDITIONED:
681: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br= B^{-1}r */
682: PetscCall(VecNorm(Br, NORM_2, &rnorm));
683: break;
684: case KSP_NORM_UNPRECONDITIONED:
685: case KSP_NORM_NATURAL:
686: PetscCall(VecNorm(r, NORM_2, &rnorm));
687: break;
688: default:
689: rnorm = 0.0;
690: break;
691: }
692: KSPCheckNorm(ksp, rnorm);
693: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
694: ksp->rnorm = rnorm;
695: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
696: PetscCall(KSPLogResidualHistory(ksp, rnorm));
697: PetscCall(KSPMonitor(ksp, i, rnorm));
698: }
699: if (ksp->its >= ksp->max_it) {
700: if (ksp->normtype != KSP_NORM_NONE) {
701: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
702: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
703: } else ksp->reason = KSP_CONVERGED_ITS;
704: }
705: }
707: if (ksp->reason == KSP_CONVERGED_ITS) PetscCall(KSPLogErrorHistory(ksp));
708: PetscFunctionReturn(PETSC_SUCCESS);
709: }
711: static PetscErrorCode KSPView_Chebyshev(KSP ksp, PetscViewer viewer)
712: {
713: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
714: PetscBool iascii;
716: PetscFunctionBegin;
717: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
718: if (iascii) {
719: switch (cheb->chebykind) {
720: case KSP_CHEBYSHEV_FIRST:
721: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of first kind\n"));
722: break;
723: case KSP_CHEBYSHEV_FOURTH:
724: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of fourth kind\n"));
725: break;
726: case KSP_CHEBYSHEV_OPT_FOURTH:
727: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of opt. fourth kind\n"));
728: break;
729: }
730: PetscReal emax, emin;
731: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
732: PetscCall(PetscViewerASCIIPrintf(viewer, " eigenvalue targets used: min %g, max %g\n", (double)emin, (double)emax));
733: if (cheb->kspest) {
734: 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));
735: 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]));
736: PetscCall(PetscViewerASCIIPushTab(viewer));
737: PetscCall(KSPView(cheb->kspest, viewer));
738: PetscCall(PetscViewerASCIIPopTab(viewer));
739: if (cheb->usenoisy) PetscCall(PetscViewerASCIIPrintf(viewer, " estimating eigenvalues using a noisy random number generated right-hand side\n"));
740: } else if (cheb->emax_provided != 0.) {
741: 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],
742: (double)cheb->tform[3]));
743: }
744: }
745: PetscFunctionReturn(PETSC_SUCCESS);
746: }
748: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
749: {
750: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
751: PetscBool isset, flg;
752: Mat Pmat, Amat;
753: PetscObjectId amatid, pmatid;
754: PetscObjectState amatstate, pmatstate;
756: PetscFunctionBegin;
757: switch (cheb->chebykind) {
758: case KSP_CHEBYSHEV_FIRST:
759: ksp->ops->solve = KSPSolve_Chebyshev_FirstKind;
760: break;
761: case KSP_CHEBYSHEV_FOURTH:
762: case KSP_CHEBYSHEV_OPT_FOURTH:
763: ksp->ops->solve = KSPSolve_Chebyshev_FourthKind;
764: break;
765: }
767: if (ksp->max_it > cheb->num_betas_alloc) {
768: PetscCall(PetscFree(cheb->betas));
769: PetscCall(PetscMalloc1(ksp->max_it, &cheb->betas));
770: cheb->num_betas_alloc = ksp->max_it;
771: }
773: // coefficients for 4th-kind Chebyshev
774: for (PetscInt i = 0; i < ksp->max_it; i++) cheb->betas[i] = 1.0;
776: // coefficients for optimized 4th-kind Chebyshev
777: if (cheb->chebykind == KSP_CHEBYSHEV_OPT_FOURTH) PetscCall(KSPChebyshevGetBetas_Private(ksp));
779: PetscCall(KSPSetWorkVecs(ksp, 3));
780: if (cheb->emin == 0. || cheb->emax == 0.) { // User did not specify eigenvalues
781: PC pc;
783: PetscCall(KSPGetPC(ksp, &pc));
784: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &flg));
785: if (!flg) { // Provided estimates are only relevant for Jacobi
786: cheb->emax_provided = 0;
787: cheb->emin_provided = 0;
788: }
789: /* We need to estimate eigenvalues */
790: if (!cheb->kspest) PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
791: }
792: if (cheb->kspest) {
793: PetscCall(KSPGetOperators(ksp, &Amat, &Pmat));
794: PetscCall(MatIsSPDKnown(Pmat, &isset, &flg));
795: if (isset && flg) {
796: const char *prefix;
798: PetscCall(KSPGetOptionsPrefix(cheb->kspest, &prefix));
799: PetscCall(PetscOptionsHasName(NULL, prefix, "-ksp_type", &flg));
800: if (!flg) PetscCall(KSPSetType(cheb->kspest, KSPCG));
801: }
802: PetscCall(PetscObjectGetId((PetscObject)Amat, &amatid));
803: PetscCall(PetscObjectGetId((PetscObject)Pmat, &pmatid));
804: PetscCall(PetscObjectStateGet((PetscObject)Amat, &amatstate));
805: PetscCall(PetscObjectStateGet((PetscObject)Pmat, &pmatstate));
806: if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
807: PetscReal max = 0.0, min = 0.0;
808: Vec B;
809: KSPConvergedReason reason;
811: PetscCall(KSPSetPC(cheb->kspest, ksp->pc));
812: if (cheb->usenoisy) {
813: B = ksp->work[1];
814: PetscCall(KSPSetNoisy_Private(Amat, B));
815: } else {
816: PetscBool change;
818: 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");
819: PetscCall(PCPreSolveChangeRHS(ksp->pc, &change));
820: if (change) {
821: B = ksp->work[1];
822: PetscCall(VecCopy(ksp->vec_rhs, B));
823: } else B = ksp->vec_rhs;
824: }
825: if (ksp->setfromoptionscalled && !cheb->kspest->setfromoptionscalled) PetscCall(KSPSetFromOptions(cheb->kspest));
826: PetscCall(KSPSolve(cheb->kspest, B, ksp->work[0]));
827: PetscCall(KSPGetConvergedReason(cheb->kspest, &reason));
828: if (reason == KSP_DIVERGED_ITS) {
829: PetscCall(PetscInfo(ksp, "Eigen estimator ran for prescribed number of iterations\n"));
830: } else if (reason == KSP_DIVERGED_PC_FAILED) {
831: PetscInt its;
832: PCFailedReason pcreason;
834: PetscCall(KSPGetIterationNumber(cheb->kspest, &its));
835: if (ksp->normtype == KSP_NORM_NONE) PetscCall(PCReduceFailedReason(ksp->pc));
836: PetscCall(PCGetFailedReason(ksp->pc, &pcreason));
837: ksp->reason = KSP_DIVERGED_PC_FAILED;
838: PetscCall(PetscInfo(ksp, "Eigen estimator failed: %s %s at iteration %" PetscInt_FMT "\n", KSPConvergedReasons[reason], PCFailedReasons[pcreason], its));
839: PetscFunctionReturn(PETSC_SUCCESS);
840: } else if (reason == KSP_CONVERGED_RTOL || reason == KSP_CONVERGED_ATOL) {
841: PetscCall(PetscInfo(ksp, "Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n"));
842: } else if (reason < 0) {
843: PetscCall(PetscInfo(ksp, "Eigen estimator failed %s, using estimates anyway\n", KSPConvergedReasons[reason]));
844: }
846: PetscCall(KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest, &min, &max));
847: PetscCall(KSPSetPC(cheb->kspest, NULL));
849: cheb->emin_computed = min;
850: cheb->emax_computed = max;
852: cheb->amatid = amatid;
853: cheb->pmatid = pmatid;
854: cheb->amatstate = amatstate;
855: cheb->pmatstate = pmatstate;
856: }
857: }
858: if (ksp->monitor[0] == (PetscErrorCode (*)(KSP, PetscInt, PetscReal, void *))KSPMonitorResidual && !ksp->normtype) PetscCall(KSPSetNormType(ksp, KSP_NORM_PRECONDITIONED));
859: PetscFunctionReturn(PETSC_SUCCESS);
860: }
862: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
863: {
864: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
866: PetscFunctionBegin;
867: PetscCall(PetscFree(cheb->betas));
868: PetscCall(KSPDestroy(&cheb->kspest));
869: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", NULL));
870: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", NULL));
871: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", NULL));
872: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetKind_C", NULL));
873: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevGetKind_C", NULL));
874: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", NULL));
875: PetscCall(KSPDestroyDefault(ksp));
876: PetscFunctionReturn(PETSC_SUCCESS);
877: }
879: /*MC
880: KSPCHEBYSHEV - The preconditioned Chebyshev iterative method
882: Options Database Keys:
883: + -ksp_chebyshev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
884: of the preconditioned operator. If these are accurate you will get much faster convergence.
885: . -ksp_chebyshev_esteig <a,b,c,d> - estimate eigenvalues using a Krylov method, then use this
886: transform for Chebyshev eigenvalue bounds (`KSPChebyshevEstEigSet()`)
887: . -ksp_chebyshev_esteig_steps - number of eigenvalue estimation steps
888: - -ksp_chebyshev_esteig_noisy - use a noisy random number generator to create right-hand side for eigenvalue estimator
890: Level: beginner
892: Notes:
893: The Chebyshev method requires both the matrix and preconditioner to be symmetric positive (semi) definite, but it can work
894: as a smoother in other situations
896: Only has support for left preconditioning.
898: Chebyshev is configured as a smoother by default, targeting the "upper" part of the spectrum.
900: 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.
901: See `MatIsSPDKnown()` for how to indicate a `Mat`, matrix is SPD.
903: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
904: `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigSetUseNoisy()`
905: `KSPRICHARDSON`, `KSPCG`, `PCMG`
906: M*/
908: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
909: {
910: KSP_Chebyshev *chebyshevP;
912: PetscFunctionBegin;
913: PetscCall(PetscNew(&chebyshevP));
915: ksp->data = (void *)chebyshevP;
916: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
917: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
918: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
919: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
921: chebyshevP->emin = 0.;
922: chebyshevP->emax = 0.;
924: chebyshevP->tform[0] = 0.0;
925: chebyshevP->tform[1] = 0.1;
926: chebyshevP->tform[2] = 0;
927: chebyshevP->tform[3] = 1.1;
928: chebyshevP->eststeps = 10;
929: chebyshevP->usenoisy = PETSC_TRUE;
930: ksp->setupnewmatrix = PETSC_TRUE;
932: ksp->ops->setup = KSPSetUp_Chebyshev;
933: ksp->ops->destroy = KSPDestroy_Chebyshev;
934: ksp->ops->buildsolution = KSPBuildSolutionDefault;
935: ksp->ops->buildresidual = KSPBuildResidualDefault;
936: ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
937: ksp->ops->view = KSPView_Chebyshev;
938: ksp->ops->reset = KSPReset_Chebyshev;
940: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", KSPChebyshevSetEigenvalues_Chebyshev));
941: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", KSPChebyshevEstEigSet_Chebyshev));
942: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", KSPChebyshevEstEigSetUseNoisy_Chebyshev));
943: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetKind_C", KSPChebyshevSetKind_Chebyshev));
944: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevGetKind_C", KSPChebyshevGetKind_Chebyshev));
945: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", KSPChebyshevEstEigGetKSP_Chebyshev));
946: PetscFunctionReturn(PETSC_SUCCESS);
947: }