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;
82: PetscFunctionBegin;
83: if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
84: if ((cheb->emin_provided == 0. || cheb->emax_provided == 0.) && !cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
85: PetscCall(KSPCreate(PetscObjectComm((PetscObject)ksp), &cheb->kspest));
86: PetscCall(KSPSetErrorIfNotConverged(cheb->kspest, ksp->errorifnotconverged));
87: PetscCall(PetscObjectIncrementTabLevel((PetscObject)cheb->kspest, (PetscObject)ksp, 1));
88: /* use PetscObjectSet/AppendOptionsPrefix() instead of KSPSet/AppendOptionsPrefix() so that the PC prefix is not changed */
89: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)cheb->kspest, ((PetscObject)ksp)->prefix));
90: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)cheb->kspest, "esteig_"));
91: 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: PetscCall(PetscInfo(ksp, "Created eigen estimator KSP\n"));
97: }
98: if (a >= 0) cheb->tform[0] = a;
99: if (b >= 0) cheb->tform[1] = b;
100: if (c >= 0) cheb->tform[2] = c;
101: if (d >= 0) cheb->tform[3] = d;
102: cheb->amatid = 0;
103: cheb->pmatid = 0;
104: cheb->amatstate = -1;
105: cheb->pmatstate = -1;
106: } else {
107: PetscCall(KSPDestroy(&cheb->kspest));
108: }
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
112: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp, PetscBool use)
113: {
114: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
116: PetscFunctionBegin;
117: cheb->usenoisy = use;
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
121: static PetscErrorCode KSPChebyshevSetKind_Chebyshev(KSP ksp, KSPChebyshevKind kind)
122: {
123: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
125: PetscFunctionBegin;
126: cheb->chebykind = kind;
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
130: static PetscErrorCode KSPChebyshevGetKind_Chebyshev(KSP ksp, KSPChebyshevKind *kind)
131: {
132: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
134: PetscFunctionBegin;
135: *kind = cheb->chebykind;
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
138: /*@
139: KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues 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: Level: intermediate
153: Notes:
154: Call `KSPChebyshevEstEigSet()` or use the option `-ksp_chebyshev_esteig a,b,c,d` to have the `KSP`
155: estimate the eigenvalues and use these estimated values automatically.
157: When `KSPCHEBYSHEV` is used as a smoother, one often wants to target a portion of the spectrum rather than the entire
158: spectrum. This function takes the range of target eigenvalues for Chebyshev, which will often slightly over-estimate
159: the largest eigenvalue of the actual operator (for safety) and greatly overestimate the smallest eigenvalue to
160: improve the smoothing properties of Chebyshev iteration on the higher frequencies in the spectrum.
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 Lanczos (`KSPCG`) or Arnoldi (`KSPGMRES`) process depending on if the operator is
203: symmetric definite or not.
205: Level: intermediate
207: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSetUseNoisy()`, `KSPChebyshevEstEigGetKSP()`
208: @*/
209: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp, PetscReal a, PetscReal b, PetscReal c, PetscReal d)
210: {
211: PetscFunctionBegin;
217: PetscTryMethod(ksp, "KSPChebyshevEstEigSet_C", (KSP, PetscReal, PetscReal, PetscReal, PetscReal), (ksp, a, b, c, d));
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: /*@
222: KSPChebyshevEstEigSetUseNoisy - use a noisy random number generated right-hand side to estimate the extreme eigenvalues instead of the given right-hand side
224: Logically Collective
226: Input Parameters:
227: + ksp - linear solver context
228: - use - `PETSC_TRUE` to use noisy
230: Options Database Key:
231: . -ksp_chebyshev_esteig_noisy <true,false> - Use noisy right-hand side for estimate
233: Level: intermediate
235: Note:
236: This allegedly works better for multigrid smoothers
238: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigGetKSP()`
239: @*/
240: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp, PetscBool use)
241: {
242: PetscFunctionBegin;
245: PetscTryMethod(ksp, "KSPChebyshevEstEigSetUseNoisy_C", (KSP, PetscBool), (ksp, use));
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
249: /*@
250: KSPChebyshevEstEigGetKSP - Get the Krylov method context used to estimate the eigenvalues for the Chebyshev method.
252: Input Parameter:
253: . ksp - the Krylov space context
255: Output Parameter:
256: . kspest - the eigenvalue estimation Krylov space context
258: Level: advanced
260: Notes:
261: If a Krylov method is not being used for this purpose, `NULL` is returned. The reference count of the returned `KSP` is
262: not incremented: it should not be destroyed by the user.
264: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`
265: @*/
266: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
267: {
268: PetscFunctionBegin;
270: PetscAssertPointer(kspest, 2);
271: *kspest = NULL;
272: PetscTryMethod(ksp, "KSPChebyshevEstEigGetKSP_C", (KSP, KSP *), (ksp, kspest));
273: PetscFunctionReturn(PETSC_SUCCESS);
274: }
276: /*@
277: KSPChebyshevSetKind - set the kind of Chebyshev polynomial to use
279: Logically Collective
281: Input Parameters:
282: + ksp - Linear solver context
283: - kind - The kind of Chebyshev polynomial to use, see `KSPChebyshevKind`, one of `KSP_CHEBYSHEV_FIRST`, `KSP_CHEBYSHEV_FOURTH`, or `KSP_CHEBYSHEV_OPT_FOURTH`
285: Options Database Key:
286: . -ksp_chebyshev_kind <kind> - which kind of Chebyshev polynomial to use
288: Level: intermediate
290: Note:
291: When using multigrid methods for problems with a poor quality coarse space (e.g., due to anisotropy or aggressive
292: coarsening), it is necessary for the smoother to handle smaller eigenvalues. With first-kind Chebyshev smoothing, this
293: requires using higher degree Chebyhev polynomials and reducing the lower end of the target spectrum, at which point
294: the whole target spectrum experiences about the same damping. Fourth kind Chebyshev polynomials (and the "optimized"
295: fourth kind) avoid the ad-hoc choice of lower bound and extend smoothing to smaller eigenvalues while preferentially
296: smoothing higher modes faster as needed to minimize the energy norm of the error. {cite}`phillips2022optimal`, {cite}`lottes2023optimal`
298: .seealso: [](ch_ksp), `KSPCHEBYSHEV` `KSPChebyshevKind`, `KSPChebyshevGetKind()`, `KSP_CHEBYSHEV_FIRST`, `KSP_CHEBYSHEV_FOURTH`, `KSP_CHEBYSHEV_OPT_FOURTH`
299: @*/
300: PetscErrorCode KSPChebyshevSetKind(KSP ksp, KSPChebyshevKind kind)
301: {
302: PetscFunctionBegin;
305: PetscUseMethod(ksp, "KSPChebyshevSetKind_C", (KSP, KSPChebyshevKind), (ksp, kind));
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: /*@
310: KSPChebyshevGetKind - get the kind of Chebyshev polynomial to use
312: Logically Collective
314: Input Parameters:
315: + ksp - Linear solver context
316: - kind - The kind of Chebyshev polynomial used
318: Level: intermediate
320: .seealso: [](ch_ksp), `KSPCHEBYSHEV` `KSPChebyshevKind`, `KSPChebyshevSetKind()`, `KSP_CHEBYSHEV_FIRST`, `KSP_CHEBYSHEV_FOURTH`, `KSP_CHEBYSHEV_OPT_FOURTH`
321: @*/
322: PetscErrorCode KSPChebyshevGetKind(KSP ksp, KSPChebyshevKind *kind)
323: {
324: PetscFunctionBegin;
326: PetscUseMethod(ksp, "KSPChebyshevGetKind_C", (KSP, KSPChebyshevKind *), (ksp, kind));
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }
330: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
331: {
332: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
334: PetscFunctionBegin;
335: *kspest = cheb->kspest;
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: static PetscErrorCode KSPSetFromOptions_Chebyshev(KSP ksp, PetscOptionItems *PetscOptionsObject)
340: {
341: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
342: PetscInt neigarg = 2, nestarg = 4;
343: PetscReal eminmax[2] = {0., 0.};
344: PetscReal tform[4] = {PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE};
345: PetscBool flgeig, flgest;
347: PetscFunctionBegin;
348: PetscOptionsHeadBegin(PetscOptionsObject, "KSP Chebyshev Options");
349: PetscCall(PetscOptionsInt("-ksp_chebyshev_esteig_steps", "Number of est steps in Chebyshev", "", cheb->eststeps, &cheb->eststeps, NULL));
350: PetscCall(PetscOptionsRealArray("-ksp_chebyshev_eigenvalues", "extreme eigenvalues", "KSPChebyshevSetEigenvalues", eminmax, &neigarg, &flgeig));
351: if (flgeig) {
352: PetscCheck(neigarg == 2, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
353: PetscCall(KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]));
354: }
355: PetscCall(PetscOptionsRealArray("-ksp_chebyshev_esteig", "estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds", "KSPChebyshevEstEigSet", tform, &nestarg, &flgest));
356: if (flgest) {
357: switch (nestarg) {
358: case 0:
359: PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
360: break;
361: case 2: /* Base everything on the max eigenvalues */
362: PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, tform[0], PETSC_DECIDE, tform[1]));
363: break;
364: case 4: /* Use the full 2x2 linear transformation */
365: PetscCall(KSPChebyshevEstEigSet(ksp, tform[0], tform[1], tform[2], tform[3]));
366: break;
367: default:
368: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
369: }
370: }
372: cheb->chebykind = KSP_CHEBYSHEV_FIRST; /* Default to 1st-kind Chebyshev polynomial */
373: PetscCall(PetscOptionsEnum("-ksp_chebyshev_kind", "Type of Chebyshev polynomial", "KSPChebyshevKind", KSPChebyshevKinds, (PetscEnum)cheb->chebykind, (PetscEnum *)&cheb->chebykind, NULL));
375: /* We need to estimate eigenvalues; need to set this here so that KSPSetFromOptions() is called on the estimator */
376: if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
378: if (cheb->kspest) {
379: PetscCall(PetscOptionsBool("-ksp_chebyshev_esteig_noisy", "Use noisy random number generated right-hand side for estimate", "KSPChebyshevEstEigSetUseNoisy", cheb->usenoisy, &cheb->usenoisy, NULL));
380: PetscCall(KSPSetFromOptions(cheb->kspest));
381: }
382: PetscOptionsHeadEnd();
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: static PetscErrorCode KSPSolve_Chebyshev_FirstKind(KSP ksp)
387: {
388: PetscInt k, kp1, km1, ktmp, i;
389: PetscScalar alpha, omegaprod, mu, omega, Gamma, c[3], scale;
390: PetscReal rnorm = 0.0, emax, emin;
391: Vec sol_orig, b, p[3], r;
392: Mat Amat, Pmat;
393: PetscBool diagonalscale;
395: PetscFunctionBegin;
396: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
397: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
399: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
400: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
401: ksp->its = 0;
402: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
403: /* These three point to the three active solutions, we
404: rotate these three at each solution update */
405: km1 = 0;
406: k = 1;
407: kp1 = 2;
408: sol_orig = ksp->vec_sol; /* ksp->vec_sol will be assigned to rotating vector p[k], thus save its address */
409: b = ksp->vec_rhs;
410: p[km1] = sol_orig;
411: p[k] = ksp->work[0];
412: p[kp1] = ksp->work[1];
413: r = ksp->work[2];
415: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
416: /* use scale*B as our preconditioner */
417: scale = 2.0 / (emax + emin);
419: /* -alpha <= scale*lambda(B^{-1}A) <= alpha */
420: alpha = 1.0 - scale * emin;
421: Gamma = 1.0;
422: mu = 1.0 / alpha;
423: omegaprod = 2.0 / alpha;
425: c[km1] = 1.0;
426: c[k] = mu;
428: if (!ksp->guess_zero) {
429: PetscCall(KSP_MatMult(ksp, Amat, sol_orig, r)); /* r = b - A*p[km1] */
430: PetscCall(VecAYPX(r, -1.0, b));
431: } else {
432: PetscCall(VecCopy(b, r));
433: }
435: /* calculate residual norm if requested, we have done one iteration */
436: if (ksp->normtype) {
437: switch (ksp->normtype) {
438: case KSP_NORM_PRECONDITIONED:
439: PetscCall(KSP_PCApply(ksp, r, p[k])); /* p[k] = B^{-1}r */
440: PetscCall(VecNorm(p[k], NORM_2, &rnorm));
441: break;
442: case KSP_NORM_UNPRECONDITIONED:
443: case KSP_NORM_NATURAL:
444: PetscCall(VecNorm(r, NORM_2, &rnorm));
445: break;
446: default:
447: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
448: }
449: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
450: ksp->rnorm = rnorm;
451: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
452: PetscCall(KSPLogResidualHistory(ksp, rnorm));
453: PetscCall(KSPLogErrorHistory(ksp));
454: PetscCall(KSPMonitor(ksp, 0, rnorm));
455: PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
456: } else ksp->reason = KSP_CONVERGED_ITERATING;
457: if (ksp->reason || ksp->max_it == 0) {
458: if (ksp->max_it == 0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }
461: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, p[k])); /* p[k] = B^{-1}r */ }
462: PetscCall(VecAYPX(p[k], scale, p[km1])); /* p[k] = scale B^{-1}r + p[km1] */
463: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
464: ksp->its = 1;
465: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
467: for (i = 1; i < ksp->max_it; i++) {
468: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
469: ksp->its++;
470: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
472: PetscCall(KSP_MatMult(ksp, Amat, p[k], r)); /* r = b - Ap[k] */
473: PetscCall(VecAYPX(r, -1.0, b));
474: /* calculate residual norm if requested */
475: if (ksp->normtype) {
476: switch (ksp->normtype) {
477: case KSP_NORM_PRECONDITIONED:
478: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
479: PetscCall(VecNorm(p[kp1], NORM_2, &rnorm));
480: break;
481: case KSP_NORM_UNPRECONDITIONED:
482: case KSP_NORM_NATURAL:
483: PetscCall(VecNorm(r, NORM_2, &rnorm));
484: break;
485: default:
486: rnorm = 0.0;
487: break;
488: }
489: KSPCheckNorm(ksp, rnorm);
490: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
491: ksp->rnorm = rnorm;
492: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
493: PetscCall(KSPLogResidualHistory(ksp, rnorm));
494: PetscCall(KSPMonitor(ksp, i, rnorm));
495: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
496: if (ksp->reason) break;
497: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */ }
498: } else {
499: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
500: }
501: ksp->vec_sol = p[k];
502: PetscCall(KSPLogErrorHistory(ksp));
504: c[kp1] = 2.0 * mu * c[k] - c[km1];
505: omega = omegaprod * c[k] / c[kp1];
507: /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
508: PetscCall(VecAXPBYPCZ(p[kp1], 1.0 - omega, omega, omega * Gamma * scale, p[km1], p[k]));
510: ktmp = km1;
511: km1 = k;
512: k = kp1;
513: kp1 = ktmp;
514: }
515: if (!ksp->reason) {
516: if (ksp->normtype) {
517: PetscCall(KSP_MatMult(ksp, Amat, p[k], r)); /* r = b - Ap[k] */
518: PetscCall(VecAYPX(r, -1.0, b));
519: switch (ksp->normtype) {
520: case KSP_NORM_PRECONDITIONED:
521: PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
522: PetscCall(VecNorm(p[kp1], NORM_2, &rnorm));
523: break;
524: case KSP_NORM_UNPRECONDITIONED:
525: case KSP_NORM_NATURAL:
526: PetscCall(VecNorm(r, NORM_2, &rnorm));
527: break;
528: default:
529: rnorm = 0.0;
530: break;
531: }
532: KSPCheckNorm(ksp, rnorm);
533: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
534: ksp->rnorm = rnorm;
535: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
536: PetscCall(KSPLogResidualHistory(ksp, rnorm));
537: PetscCall(KSPMonitor(ksp, i, rnorm));
538: }
539: if (ksp->its >= ksp->max_it) {
540: if (ksp->normtype != KSP_NORM_NONE) {
541: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
542: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
543: } else ksp->reason = KSP_CONVERGED_ITS;
544: }
545: }
547: /* make sure solution is in vector x */
548: ksp->vec_sol = sol_orig;
549: if (k) PetscCall(VecCopy(p[k], sol_orig));
550: if (ksp->reason == KSP_CONVERGED_ITS) PetscCall(KSPLogErrorHistory(ksp));
551: PetscFunctionReturn(PETSC_SUCCESS);
552: }
554: static PetscErrorCode KSPSolve_Chebyshev_FourthKind(KSP ksp)
555: {
556: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
557: PetscInt i;
558: PetscScalar scale, rScale, dScale;
559: PetscReal rnorm = 0.0, emax, emin;
560: Vec x, b, d, r, Br;
561: Mat Amat, Pmat;
562: PetscBool diagonalscale;
563: PetscReal *betas = cheb->betas;
565: PetscFunctionBegin;
566: PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
567: PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);
569: PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
570: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
571: ksp->its = 0;
572: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
574: x = ksp->vec_sol;
575: b = ksp->vec_rhs;
576: r = ksp->work[0];
577: d = ksp->work[1];
578: Br = ksp->work[2];
580: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
581: /* use scale*B as our preconditioner */
582: scale = 1.0 / emax;
584: if (!ksp->guess_zero) {
585: PetscCall(KSP_MatMult(ksp, Amat, x, r)); /* r = b - A*x */
586: PetscCall(VecAYPX(r, -1.0, b));
587: } else {
588: PetscCall(VecCopy(b, r));
589: }
591: /* calculate residual norm if requested, we have done one iteration */
592: if (ksp->normtype) {
593: switch (ksp->normtype) {
594: case KSP_NORM_PRECONDITIONED:
595: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
596: PetscCall(VecNorm(Br, NORM_2, &rnorm));
597: break;
598: case KSP_NORM_UNPRECONDITIONED:
599: case KSP_NORM_NATURAL:
600: PetscCall(VecNorm(r, NORM_2, &rnorm));
601: break;
602: default:
603: SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
604: }
605: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
606: ksp->rnorm = rnorm;
607: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
608: PetscCall(KSPLogResidualHistory(ksp, rnorm));
609: PetscCall(KSPLogErrorHistory(ksp));
610: PetscCall(KSPMonitor(ksp, 0, rnorm));
611: PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
612: } else ksp->reason = KSP_CONVERGED_ITERATING;
613: if (ksp->reason || ksp->max_it == 0) {
614: if (ksp->max_it == 0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
615: PetscFunctionReturn(PETSC_SUCCESS);
616: }
617: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */ }
618: PetscCall(VecAXPBY(d, 4.0 / 3.0 * scale, 0.0, Br)); /* d = 4/3 * scale B^{-1}r */
619: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
620: ksp->its = 1;
621: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
623: for (i = 1; i < ksp->max_it; i++) {
624: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
625: ksp->its++;
626: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
628: PetscCall(VecAXPBY(x, betas[i - 1], 1.0, d)); /* x = x + \beta_k d */
630: PetscCall(KSP_MatMult(ksp, Amat, d, Br)); /* r = r - Ad */
631: PetscCall(VecAXPBY(r, -1.0, 1.0, Br));
633: /* calculate residual norm if requested */
634: if (ksp->normtype) {
635: switch (ksp->normtype) {
636: case KSP_NORM_PRECONDITIONED:
637: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
638: PetscCall(VecNorm(Br, NORM_2, &rnorm));
639: break;
640: case KSP_NORM_UNPRECONDITIONED:
641: case KSP_NORM_NATURAL:
642: PetscCall(VecNorm(r, NORM_2, &rnorm));
643: break;
644: default:
645: rnorm = 0.0;
646: break;
647: }
648: KSPCheckNorm(ksp, rnorm);
649: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
650: ksp->rnorm = rnorm;
651: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
652: PetscCall(KSPLogResidualHistory(ksp, rnorm));
653: PetscCall(KSPMonitor(ksp, i, rnorm));
654: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
655: if (ksp->reason) break;
656: if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */ }
657: } else {
658: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
659: }
660: PetscCall(KSPLogErrorHistory(ksp));
662: rScale = scale * (8.0 * i + 4.0) / (2.0 * i + 3.0);
663: dScale = (2.0 * i - 1.0) / (2.0 * i + 3.0);
665: /* d_k+1 = \dfrac{2k-1}{2k+3} d_k + \dfrac{8k+4}{2k+3} \dfrac{1}{\rho(SA)} Br */
666: PetscCall(VecAXPBY(d, rScale, dScale, Br));
667: }
669: /* on last pass, update solution vector */
670: PetscCall(VecAXPBY(x, betas[ksp->max_it - 1], 1.0, d)); /* x = x + d */
672: if (!ksp->reason) {
673: if (ksp->normtype) {
674: PetscCall(KSP_MatMult(ksp, Amat, x, r)); /* r = b - Ax */
675: PetscCall(VecAYPX(r, -1.0, b));
676: switch (ksp->normtype) {
677: case KSP_NORM_PRECONDITIONED:
678: PetscCall(KSP_PCApply(ksp, r, Br)); /* Br= B^{-1}r */
679: PetscCall(VecNorm(Br, NORM_2, &rnorm));
680: break;
681: case KSP_NORM_UNPRECONDITIONED:
682: case KSP_NORM_NATURAL:
683: PetscCall(VecNorm(r, NORM_2, &rnorm));
684: break;
685: default:
686: rnorm = 0.0;
687: break;
688: }
689: KSPCheckNorm(ksp, rnorm);
690: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
691: ksp->rnorm = rnorm;
692: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
693: PetscCall(KSPLogResidualHistory(ksp, rnorm));
694: PetscCall(KSPMonitor(ksp, i, rnorm));
695: }
696: if (ksp->its >= ksp->max_it) {
697: if (ksp->normtype != KSP_NORM_NONE) {
698: PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
699: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
700: } else ksp->reason = KSP_CONVERGED_ITS;
701: }
702: }
704: if (ksp->reason == KSP_CONVERGED_ITS) PetscCall(KSPLogErrorHistory(ksp));
705: PetscFunctionReturn(PETSC_SUCCESS);
706: }
708: static PetscErrorCode KSPView_Chebyshev(KSP ksp, PetscViewer viewer)
709: {
710: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
711: PetscBool iascii;
713: PetscFunctionBegin;
714: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
715: if (iascii) {
716: switch (cheb->chebykind) {
717: case KSP_CHEBYSHEV_FIRST:
718: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of first kind\n"));
719: break;
720: case KSP_CHEBYSHEV_FOURTH:
721: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of fourth kind\n"));
722: break;
723: case KSP_CHEBYSHEV_OPT_FOURTH:
724: PetscCall(PetscViewerASCIIPrintf(viewer, " Chebyshev polynomial of opt. fourth kind\n"));
725: break;
726: }
727: PetscReal emax, emin;
728: PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
729: PetscCall(PetscViewerASCIIPrintf(viewer, " eigenvalue targets used: min %g, max %g\n", (double)emin, (double)emax));
730: if (cheb->kspest) {
731: 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));
732: 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]));
733: PetscCall(PetscViewerASCIIPushTab(viewer));
734: PetscCall(KSPView(cheb->kspest, viewer));
735: PetscCall(PetscViewerASCIIPopTab(viewer));
736: if (cheb->usenoisy) PetscCall(PetscViewerASCIIPrintf(viewer, " estimating eigenvalues using a noisy random number generated right-hand side\n"));
737: } else if (cheb->emax_provided != 0.) {
738: 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],
739: (double)cheb->tform[3]));
740: }
741: }
742: PetscFunctionReturn(PETSC_SUCCESS);
743: }
745: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
746: {
747: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
748: PetscBool isset, flg;
749: Mat Pmat, Amat;
750: PetscObjectId amatid, pmatid;
751: PetscObjectState amatstate, pmatstate;
753: PetscFunctionBegin;
754: switch (cheb->chebykind) {
755: case KSP_CHEBYSHEV_FIRST:
756: ksp->ops->solve = KSPSolve_Chebyshev_FirstKind;
757: break;
758: case KSP_CHEBYSHEV_FOURTH:
759: case KSP_CHEBYSHEV_OPT_FOURTH:
760: ksp->ops->solve = KSPSolve_Chebyshev_FourthKind;
761: break;
762: }
764: if (ksp->max_it > cheb->num_betas_alloc) {
765: PetscCall(PetscFree(cheb->betas));
766: PetscCall(PetscMalloc1(ksp->max_it, &cheb->betas));
767: cheb->num_betas_alloc = ksp->max_it;
768: }
770: // coefficients for 4th-kind Chebyshev
771: for (PetscInt i = 0; i < ksp->max_it; i++) cheb->betas[i] = 1.0;
773: // coefficients for optimized 4th-kind Chebyshev
774: if (cheb->chebykind == KSP_CHEBYSHEV_OPT_FOURTH) PetscCall(KSPChebyshevGetBetas_Private(ksp));
776: PetscCall(KSPSetWorkVecs(ksp, 3));
777: if (cheb->emin == 0. || cheb->emax == 0.) { // User did not specify eigenvalues
778: PC pc;
780: PetscCall(KSPGetPC(ksp, &pc));
781: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &flg));
782: if (!flg) { // Provided estimates are only relevant for Jacobi
783: cheb->emax_provided = 0;
784: cheb->emin_provided = 0;
785: }
786: /* We need to estimate eigenvalues */
787: if (!cheb->kspest) PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
788: }
789: if (cheb->kspest) {
790: PetscCall(KSPGetOperators(ksp, &Amat, &Pmat));
791: PetscCall(MatIsSPDKnown(Pmat, &isset, &flg));
792: if (isset && flg) {
793: const char *prefix;
795: PetscCall(KSPGetOptionsPrefix(cheb->kspest, &prefix));
796: PetscCall(PetscOptionsHasName(NULL, prefix, "-ksp_type", &flg));
797: if (!flg) PetscCall(KSPSetType(cheb->kspest, KSPCG));
798: }
799: PetscCall(PetscObjectGetId((PetscObject)Amat, &amatid));
800: PetscCall(PetscObjectGetId((PetscObject)Pmat, &pmatid));
801: PetscCall(PetscObjectStateGet((PetscObject)Amat, &amatstate));
802: PetscCall(PetscObjectStateGet((PetscObject)Pmat, &pmatstate));
803: if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
804: PetscReal max = 0.0, min = 0.0;
805: Vec B;
806: KSPConvergedReason reason;
808: PetscCall(KSPSetPC(cheb->kspest, ksp->pc));
809: if (cheb->usenoisy) {
810: B = ksp->work[1];
811: PetscCall(KSPSetNoisy_Private(B));
812: } else {
813: PetscBool change;
815: 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");
816: PetscCall(PCPreSolveChangeRHS(ksp->pc, &change));
817: if (change) {
818: B = ksp->work[1];
819: PetscCall(VecCopy(ksp->vec_rhs, B));
820: } else B = ksp->vec_rhs;
821: }
822: if (ksp->setfromoptionscalled && !cheb->kspest->setfromoptionscalled) PetscCall(KSPSetFromOptions(cheb->kspest));
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) PetscCall(PCReduceFailedReason(ksp->pc));
833: PetscCall(PCGetFailedReason(ksp->pc, &pcreason));
834: ksp->reason = KSP_DIVERGED_PC_FAILED;
835: PetscCall(PetscInfo(ksp, "Eigen estimator failed: %s %s at iteration %" PetscInt_FMT "\n", KSPConvergedReasons[reason], PCFailedReasons[pcreason], its));
836: PetscFunctionReturn(PETSC_SUCCESS);
837: } else if (reason == KSP_CONVERGED_RTOL || reason == KSP_CONVERGED_ATOL) {
838: PetscCall(PetscInfo(ksp, "Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n"));
839: } else if (reason < 0) {
840: PetscCall(PetscInfo(ksp, "Eigen estimator failed %s, using estimates anyway\n", KSPConvergedReasons[reason]));
841: }
843: PetscCall(KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest, &min, &max));
844: PetscCall(KSPSetPC(cheb->kspest, NULL));
846: cheb->emin_computed = min;
847: cheb->emax_computed = max;
849: cheb->amatid = amatid;
850: cheb->pmatid = pmatid;
851: cheb->amatstate = amatstate;
852: cheb->pmatstate = pmatstate;
853: }
854: }
855: PetscFunctionReturn(PETSC_SUCCESS);
856: }
858: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
859: {
860: KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
862: PetscFunctionBegin;
863: PetscCall(PetscFree(cheb->betas));
864: PetscCall(KSPDestroy(&cheb->kspest));
865: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", NULL));
866: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", NULL));
867: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", NULL));
868: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetKind_C", NULL));
869: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevGetKind_C", NULL));
870: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", NULL));
871: PetscCall(KSPDestroyDefault(ksp));
872: PetscFunctionReturn(PETSC_SUCCESS);
873: }
875: /*MC
876: KSPCHEBYSHEV - The preconditioned Chebyshev iterative method
878: Options Database Keys:
879: + -ksp_chebyshev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
880: of the preconditioned operator. If these are accurate you will get much faster convergence.
881: . -ksp_chebyshev_esteig <a,b,c,d> - estimate eigenvalues using a Krylov method, then use this
882: transform for Chebyshev eigenvalue bounds (`KSPChebyshevEstEigSet()`)
883: . -ksp_chebyshev_esteig_steps - number of estimation steps
884: - -ksp_chebyshev_esteig_noisy - use a noisy random number generator to create right-hand side for eigenvalue estimator
886: Level: beginner
888: Notes:
889: The Chebyshev method requires both the matrix and preconditioner to be symmetric positive (semi) definite, but it can work
890: as a smoother in other situations
892: Only support for left preconditioning.
894: Chebyshev is configured as a smoother by default, targeting the "upper" part of the spectrum.
896: The user should call `KSPChebyshevSetEigenvalues()` to get eigenvalue estimates.
898: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
899: `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigSetUseNoisy()`
900: `KSPRICHARDSON`, `KSPCG`, `PCMG`
901: M*/
903: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
904: {
905: KSP_Chebyshev *chebyshevP;
907: PetscFunctionBegin;
908: PetscCall(PetscNew(&chebyshevP));
910: ksp->data = (void *)chebyshevP;
911: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
912: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
913: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
914: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
916: chebyshevP->emin = 0.;
917: chebyshevP->emax = 0.;
919: chebyshevP->tform[0] = 0.0;
920: chebyshevP->tform[1] = 0.1;
921: chebyshevP->tform[2] = 0;
922: chebyshevP->tform[3] = 1.1;
923: chebyshevP->eststeps = 10;
924: chebyshevP->usenoisy = PETSC_TRUE;
925: ksp->setupnewmatrix = PETSC_TRUE;
927: ksp->ops->setup = KSPSetUp_Chebyshev;
928: ksp->ops->destroy = KSPDestroy_Chebyshev;
929: ksp->ops->buildsolution = KSPBuildSolutionDefault;
930: ksp->ops->buildresidual = KSPBuildResidualDefault;
931: ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
932: ksp->ops->view = KSPView_Chebyshev;
933: ksp->ops->reset = KSPReset_Chebyshev;
935: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", KSPChebyshevSetEigenvalues_Chebyshev));
936: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", KSPChebyshevEstEigSet_Chebyshev));
937: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", KSPChebyshevEstEigSetUseNoisy_Chebyshev));
938: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetKind_C", KSPChebyshevSetKind_Chebyshev));
939: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevGetKind_C", KSPChebyshevGetKind_Chebyshev));
940: PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", KSPChebyshevEstEigGetKSP_Chebyshev));
941: PetscFunctionReturn(PETSC_SUCCESS);
942: }