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: }