Actual source code: cheby.c

  1: #include "chebyshevimpl.h"
  2: #include <../src/ksp/ksp/impls/cheby/chebyshevimpl.h>

  4: static const char *const KSPChebyshevKinds[] = {"FIRST", "FOURTH", "OPT_FOURTH", "KSPChebyshevKinds", "KSP_CHEBYSHEV_", NULL};

  6: static PetscErrorCode KSPReset_Chebyshev(KSP ksp)
  7: {
  8:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;

 10:   PetscFunctionBegin;
 11:   if (cheb->kspest) PetscCall(KSPReset(cheb->kspest));
 12:   PetscFunctionReturn(PETSC_SUCCESS);
 13: }

 15: /*
 16:  * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
 17:  */
 18: static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest, PetscReal *emin, PetscReal *emax)
 19: {
 20:   PetscInt   n, neig;
 21:   PetscReal *re, *im, min, max;

 23:   PetscFunctionBegin;
 24:   PetscCall(KSPGetIterationNumber(kspest, &n));
 25:   PetscCall(PetscMalloc2(n, &re, n, &im));
 26:   PetscCall(KSPComputeEigenvalues(kspest, n, re, im, &neig));
 27:   min = PETSC_MAX_REAL;
 28:   max = PETSC_MIN_REAL;
 29:   for (n = 0; n < neig; n++) {
 30:     min = PetscMin(min, re[n]);
 31:     max = PetscMax(max, re[n]);
 32:   }
 33:   PetscCall(PetscFree2(re, im));
 34:   *emax = max;
 35:   *emin = min;
 36:   PetscFunctionReturn(PETSC_SUCCESS);
 37: }

 39: static PetscErrorCode KSPChebyshevGetEigenvalues_Chebyshev(KSP ksp, PetscReal *emax, PetscReal *emin)
 40: {
 41:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;

 43:   PetscFunctionBegin;
 44:   *emax = 0;
 45:   *emin = 0;
 46:   if (cheb->emax != 0.) {
 47:     *emax = cheb->emax;
 48:   } else if (cheb->emax_computed != 0.) {
 49:     *emax = cheb->tform[2] * cheb->emin_computed + cheb->tform[3] * cheb->emax_computed;
 50:   } else if (cheb->emax_provided != 0.) {
 51:     *emax = cheb->tform[2] * cheb->emin_provided + cheb->tform[3] * cheb->emax_provided;
 52:   }
 53:   if (cheb->emin != 0.) {
 54:     *emin = cheb->emin;
 55:   } else if (cheb->emin_computed != 0.) {
 56:     *emin = cheb->tform[0] * cheb->emin_computed + cheb->tform[1] * cheb->emax_computed;
 57:   } else if (cheb->emin_provided != 0.) {
 58:     *emin = cheb->tform[0] * cheb->emin_provided + cheb->tform[1] * cheb->emax_provided;
 59:   }
 60:   PetscFunctionReturn(PETSC_SUCCESS);
 61: }

 63: static PetscErrorCode KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp, PetscReal emax, PetscReal emin)
 64: {
 65:   KSP_Chebyshev *chebyshevP = (KSP_Chebyshev *)ksp->data;

 67:   PetscFunctionBegin;
 68:   PetscCheck(emax > emin, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "Maximum eigenvalue must be larger than minimum: max %g min %g", (double)emax, (double)emin);
 69:   PetscCheck(emax * emin > 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "Both eigenvalues must be of the same sign: max %g min %g", (double)emax, (double)emin);
 70:   chebyshevP->emax = emax;
 71:   chebyshevP->emin = emin;

 73:   PetscCall(KSPChebyshevEstEigSet(ksp, 0., 0., 0., 0.)); /* Destroy any estimation setup */
 74:   PetscFunctionReturn(PETSC_SUCCESS);
 75: }

 77: static PetscErrorCode KSPChebyshevEstEigSet_Chebyshev(KSP ksp, PetscReal a, PetscReal b, PetscReal c, PetscReal d)
 78: {
 79:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;

 81:   PetscFunctionBegin;
 82:   if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
 83:     if ((cheb->emin_provided == 0. || cheb->emax_provided == 0.) && !cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
 84:       PetscCall(KSPCreate(PetscObjectComm((PetscObject)ksp), &cheb->kspest));
 85:       PetscCall(KSPSetErrorIfNotConverged(cheb->kspest, ksp->errorifnotconverged));
 86:       PetscCall(PetscObjectIncrementTabLevel((PetscObject)cheb->kspest, (PetscObject)ksp, 1));
 87:       /* use PetscObjectSet/AppendOptionsPrefix() instead of KSPSet/AppendOptionsPrefix() so that the PC prefix is not changed */
 88:       PetscCall(PetscObjectSetOptionsPrefix((PetscObject)cheb->kspest, ((PetscObject)ksp)->prefix));
 89:       PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)cheb->kspest, "esteig_"));
 90:       PetscCall(KSPSetSkipPCSetFromOptions(cheb->kspest, PETSC_TRUE));

 92:       PetscCall(KSPSetComputeEigenvalues(cheb->kspest, PETSC_TRUE));

 94:       /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
 95:       PetscCall(KSPSetTolerances(cheb->kspest, 1.e-12, PETSC_DEFAULT, PETSC_DEFAULT, cheb->eststeps));
 96:     }
 97:     if (a >= 0) cheb->tform[0] = a;
 98:     if (b >= 0) cheb->tform[1] = b;
 99:     if (c >= 0) cheb->tform[2] = c;
100:     if (d >= 0) cheb->tform[3] = d;
101:     cheb->amatid    = 0;
102:     cheb->pmatid    = 0;
103:     cheb->amatstate = -1;
104:     cheb->pmatstate = -1;
105:   } else {
106:     PetscCall(KSPDestroy(&cheb->kspest));
107:   }
108:   PetscFunctionReturn(PETSC_SUCCESS);
109: }

111: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp, PetscBool use)
112: {
113:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;

115:   PetscFunctionBegin;
116:   cheb->usenoisy = use;
117:   PetscFunctionReturn(PETSC_SUCCESS);
118: }

120: static PetscErrorCode KSPChebyshevSetKind_Chebyshev(KSP ksp, KSPChebyshevKind kind)
121: {
122:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;

124:   PetscFunctionBegin;
125:   cheb->chebykind = kind;
126:   PetscFunctionReturn(PETSC_SUCCESS);
127: }

129: static PetscErrorCode KSPChebyshevGetKind_Chebyshev(KSP ksp, KSPChebyshevKind *kind)
130: {
131:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;

133:   PetscFunctionBegin;
134:   *kind = cheb->chebykind;
135:   PetscFunctionReturn(PETSC_SUCCESS);
136: }
137: /*@
138:    KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues
139:    of the preconditioned problem.

141:    Logically Collective

143:    Input Parameters:
144: +  ksp - the Krylov space context
145:    emax - the eigenvalue maximum estimate
146: -  emin - the eigenvalue minimum estimate

148:   Options Database Key:
149: .  -ksp_chebyshev_eigenvalues emin,emax - extreme eigenvalues

151:    Notes:
152:    Call `KSPChebyshevEstEigSet()` or use the option -ksp_chebyshev_esteig a,b,c,d to have the KSP
153:    estimate the eigenvalues and use these estimated values automatically.

155:    When `KSPCHEBYSHEV` is used as a smoother, one often wants to target a portion of the spectrum rather than the entire
156:    spectrum. This function takes the range of target eigenvalues for Chebyshev, which will often slightly over-estimate
157:    the largest eigenvalue of the actual operator (for safety) and greatly overestimate the smallest eigenvalue to
158:    improve the smoothing properties of Chebyshev iteration on the higher frequencies in the spectrum.

160:    Level: intermediate

162: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`,
163: @*/
164: PetscErrorCode KSPChebyshevSetEigenvalues(KSP ksp, PetscReal emax, PetscReal emin)
165: {
166:   PetscFunctionBegin;
170:   PetscTryMethod(ksp, "KSPChebyshevSetEigenvalues_C", (KSP, PetscReal, PetscReal), (ksp, emax, emin));
171:   PetscFunctionReturn(PETSC_SUCCESS);
172: }

174: /*@
175:    KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev

177:    Logically Collective

179:    Input Parameters:
180: +  ksp - the Krylov space context
181: .  a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
182: .  b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
183: .  c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
184: -  d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)

186:   Options Database Key:
187: .  -ksp_chebyshev_esteig a,b,c,d - estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds

189:    Notes:
190:    The Chebyshev bounds are set using
191: .vb
192:    minbound = a*minest + b*maxest
193:    maxbound = c*minest + d*maxest
194: .ve
195:    The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
196:    The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.

198:    If 0.0 is passed for all transform arguments (a,b,c,d), eigenvalue estimation is disabled.

200:    The default transform is (0,0.1; 0,1.1) which targets the "upper" part of the spectrum, as desirable for use with multigrid.

202:    The eigenvalues are estimated using the Lanczo (`KSPCG`) or Arnoldi (`KSPGMRES`) process using a noisy right hand side vector.

204:    Level: intermediate

206: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigSetUseNoisy()`, `KSPChebyshevEstEigGetKSP()`
207: @*/
208: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp, PetscReal a, PetscReal b, PetscReal c, PetscReal d)
209: {
210:   PetscFunctionBegin;
216:   PetscTryMethod(ksp, "KSPChebyshevEstEigSet_C", (KSP, PetscReal, PetscReal, PetscReal, PetscReal), (ksp, a, b, c, d));
217:   PetscFunctionReturn(PETSC_SUCCESS);
218: }

220: /*@
221:    KSPChebyshevEstEigSetUseNoisy - use a noisy right hand side in order to do the estimate instead of the given right hand side

223:    Logically Collective

225:    Input Parameters:
226: +  ksp - linear solver context
227: -  use - `PETSC_TRUE` to use noisy

229:    Options Database Key:
230: .  -ksp_chebyshev_esteig_noisy <true,false> - Use noisy right hand side for estimate

232:    Note:
233:     This allegedly works better for multigrid smoothers

235:   Level: intermediate

237: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigGetKSP()`
238: @*/
239: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp, PetscBool use)
240: {
241:   PetscFunctionBegin;
244:   PetscTryMethod(ksp, "KSPChebyshevEstEigSetUseNoisy_C", (KSP, PetscBool), (ksp, use));
245:   PetscFunctionReturn(PETSC_SUCCESS);
246: }

248: /*@
249:   KSPChebyshevEstEigGetKSP - Get the Krylov method context used to estimate eigenvalues for the Chebyshev method.  If
250:   a Krylov method is not being used for this purpose, NULL is returned.  The reference count of the returned `KSP` is
251:   not incremented: it should not be destroyed by the user.

253:   Input Parameter:
254: . ksp - the Krylov space context

256:   Output Parameter:
257: . kspest - the eigenvalue estimation Krylov space context

259:   Level: advanced

261: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`
262: @*/
263: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
264: {
265:   PetscFunctionBegin;
268:   *kspest = NULL;
269:   PetscTryMethod(ksp, "KSPChebyshevEstEigGetKSP_C", (KSP, KSP *), (ksp, kspest));
270:   PetscFunctionReturn(PETSC_SUCCESS);
271: }

273: /*@
274:   KSPChebyshevSetKind - set the kind of Chebyshev polynomial to use

276:   Logically Collective

278:   Input Parameters:
279: + ksp  - Linear solver context
280: - kind - The kind of Chebyshev polynomial to use

282:   Options Database Key:
283: . -ksp_chebyshev_kind <kind> - which kind of Chebyshev polynomial to use

285:   Level: intermediate

287:   Note:
288:   When using multigrid methods for problems with a poor quality coarse space (e.g., due to anisotropy or aggressive
289:   coarsening), it is necessary for the smoother to handle smaller eigenvalues. With first-kind Chebyshev smoothing, this
290:   requires using higher degree Chebyhev polynomials and reducing the lower end of the target spectrum, at which point
291:   the whole target spectrum experiences about the same damping. Fourth kind Chebyshev polynomials (and the "optimized"
292:   fourth kind) avoid the ad-hoc choice of lower bound and extend smoothing to smaller eigenvalues while preferentially
293:   smoothing higher modes faster as needed to minimize the energy norm of the error.

295:   References:
296: +  * - Malachi Phillips and Paul Fischer, Optimal Chebyshev Smoothers and One-sided V-cycles, https://arxiv.org/abs/2210.03179.
297: -  * - James Lottes, Optimal Polynomial Smoothers for Multigrid V-cycles, https://arxiv.org/abs/2202.08830.

299: .seealso: [](ch_ksp), `KSPCHEBYSHEV` `KSPChebyshevKind`, `KSPChebyshevGetKind()`
300: @*/
301: PetscErrorCode KSPChebyshevSetKind(KSP ksp, KSPChebyshevKind kind)
302: {
303:   PetscFunctionBegin;
306:   PetscUseMethod(ksp, "KSPChebyshevSetKind_C", (KSP, KSPChebyshevKind), (ksp, kind));
307:   PetscFunctionReturn(PETSC_SUCCESS);
308: }

310: /*@
311:   KSPChebyshevGetKind - get the kind of Chebyshev polynomial to use

313:   Logically Collective

315:   Input Parameters:
316: + ksp  - Linear solver context
317: - kind - The kind of Chebyshev polynomial to use

319:   Level: intermediate

321: .seealso: [](ch_ksp), `KSPCHEBYSHEV` `KSPChebyshevKind`, `KSPChebyshevSetKind()`
322: @*/
323: PetscErrorCode KSPChebyshevGetKind(KSP ksp, KSPChebyshevKind *kind)
324: {
325:   PetscFunctionBegin;
327:   PetscUseMethod(ksp, "KSPChebyshevGetKind_C", (KSP, KSPChebyshevKind *), (ksp, kind));
328:   PetscFunctionReturn(PETSC_SUCCESS);
329: }

331: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
332: {
333:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;

335:   PetscFunctionBegin;
336:   *kspest = cheb->kspest;
337:   PetscFunctionReturn(PETSC_SUCCESS);
338: }

340: static PetscErrorCode KSPSetFromOptions_Chebyshev(KSP ksp, PetscOptionItems *PetscOptionsObject)
341: {
342:   KSP_Chebyshev *cheb    = (KSP_Chebyshev *)ksp->data;
343:   PetscInt       neigarg = 2, nestarg = 4;
344:   PetscReal      eminmax[2] = {0., 0.};
345:   PetscReal      tform[4]   = {PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE};
346:   PetscBool      flgeig, flgest;

348:   PetscFunctionBegin;
349:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP Chebyshev Options");
350:   PetscCall(PetscOptionsInt("-ksp_chebyshev_esteig_steps", "Number of est steps in Chebyshev", "", cheb->eststeps, &cheb->eststeps, NULL));
351:   PetscCall(PetscOptionsRealArray("-ksp_chebyshev_eigenvalues", "extreme eigenvalues", "KSPChebyshevSetEigenvalues", eminmax, &neigarg, &flgeig));
352:   if (flgeig) {
353:     PetscCheck(neigarg == 2, PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
354:     PetscCall(KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]));
355:   }
356:   PetscCall(PetscOptionsRealArray("-ksp_chebyshev_esteig", "estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds", "KSPChebyshevEstEigSet", tform, &nestarg, &flgest));
357:   if (flgest) {
358:     switch (nestarg) {
359:     case 0:
360:       PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
361:       break;
362:     case 2: /* Base everything on the max eigenvalues */
363:       PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, tform[0], PETSC_DECIDE, tform[1]));
364:       break;
365:     case 4: /* Use the full 2x2 linear transformation */
366:       PetscCall(KSPChebyshevEstEigSet(ksp, tform[0], tform[1], tform[2], tform[3]));
367:       break;
368:     default:
369:       SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_ARG_INCOMP, "Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
370:     }
371:   }

373:   cheb->chebykind = KSP_CHEBYSHEV_FIRST; /* Default to 1st-kind Chebyshev polynomial */
374:   PetscCall(PetscOptionsEnum("-ksp_chebyshev_kind", "Type of Chebyshev polynomial", "KSPChebyshevKind", KSPChebyshevKinds, (PetscEnum)cheb->chebykind, (PetscEnum *)&cheb->chebykind, NULL));

376:   /* We need to estimate eigenvalues; need to set this here so that KSPSetFromOptions() is called on the estimator */
377:   if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));

379:   if (cheb->kspest) {
380:     PetscCall(PetscOptionsBool("-ksp_chebyshev_esteig_noisy", "Use noisy right hand side for estimate", "KSPChebyshevEstEigSetUseNoisy", cheb->usenoisy, &cheb->usenoisy, NULL));
381:     PetscCall(KSPSetFromOptions(cheb->kspest));
382:   }
383:   PetscOptionsHeadEnd();
384:   PetscFunctionReturn(PETSC_SUCCESS);
385: }

387: static PetscErrorCode KSPSolve_Chebyshev_FirstKind(KSP ksp)
388: {
389:   PetscInt    k, kp1, km1, ktmp, i;
390:   PetscScalar alpha, omegaprod, mu, omega, Gamma, c[3], scale;
391:   PetscReal   rnorm = 0.0, emax, emin;
392:   Vec         sol_orig, b, p[3], r;
393:   Mat         Amat, Pmat;
394:   PetscBool   diagonalscale;

396:   PetscFunctionBegin;
397:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
398:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);

400:   PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
401:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
402:   ksp->its = 0;
403:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
404:   /* These three point to the three active solutions, we
405:      rotate these three at each solution update */
406:   km1      = 0;
407:   k        = 1;
408:   kp1      = 2;
409:   sol_orig = ksp->vec_sol; /* ksp->vec_sol will be assigned to rotating vector p[k], thus save its address */
410:   b        = ksp->vec_rhs;
411:   p[km1]   = sol_orig;
412:   p[k]     = ksp->work[0];
413:   p[kp1]   = ksp->work[1];
414:   r        = ksp->work[2];

416:   PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
417:   /* use scale*B as our preconditioner */
418:   scale = 2.0 / (emax + emin);

420:   /*   -alpha <=  scale*lambda(B^{-1}A) <= alpha   */
421:   alpha     = 1.0 - scale * emin;
422:   Gamma     = 1.0;
423:   mu        = 1.0 / alpha;
424:   omegaprod = 2.0 / alpha;

426:   c[km1] = 1.0;
427:   c[k]   = mu;

429:   if (!ksp->guess_zero) {
430:     PetscCall(KSP_MatMult(ksp, Amat, sol_orig, r)); /*  r = b - A*p[km1] */
431:     PetscCall(VecAYPX(r, -1.0, b));
432:   } else {
433:     PetscCall(VecCopy(b, r));
434:   }

436:   /* calculate residual norm if requested, we have done one iteration */
437:   if (ksp->normtype) {
438:     switch (ksp->normtype) {
439:     case KSP_NORM_PRECONDITIONED:
440:       PetscCall(KSP_PCApply(ksp, r, p[k])); /* p[k] = B^{-1}r */
441:       PetscCall(VecNorm(p[k], NORM_2, &rnorm));
442:       break;
443:     case KSP_NORM_UNPRECONDITIONED:
444:     case KSP_NORM_NATURAL:
445:       PetscCall(VecNorm(r, NORM_2, &rnorm));
446:       break;
447:     default:
448:       SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
449:     }
450:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
451:     ksp->rnorm = rnorm;
452:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
453:     PetscCall(KSPLogResidualHistory(ksp, rnorm));
454:     PetscCall(KSPLogErrorHistory(ksp));
455:     PetscCall(KSPMonitor(ksp, 0, rnorm));
456:     PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
457:   } else ksp->reason = KSP_CONVERGED_ITERATING;
458:   if (ksp->reason || ksp->max_it == 0) {
459:     if (ksp->max_it == 0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
460:     PetscFunctionReturn(PETSC_SUCCESS);
461:   }
462:   if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, p[k])); /* p[k] = B^{-1}r */ }
463:   PetscCall(VecAYPX(p[k], scale, p[km1])); /* p[k] = scale B^{-1}r + p[km1] */
464:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
465:   ksp->its = 1;
466:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));

468:   for (i = 1; i < ksp->max_it; i++) {
469:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
470:     ksp->its++;
471:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));

473:     PetscCall(KSP_MatMult(ksp, Amat, p[k], r)); /*  r = b - Ap[k]    */
474:     PetscCall(VecAYPX(r, -1.0, b));
475:     /* calculate residual norm if requested */
476:     if (ksp->normtype) {
477:       switch (ksp->normtype) {
478:       case KSP_NORM_PRECONDITIONED:
479:         PetscCall(KSP_PCApply(ksp, r, p[kp1])); /*  p[kp1] = B^{-1}r  */
480:         PetscCall(VecNorm(p[kp1], NORM_2, &rnorm));
481:         break;
482:       case KSP_NORM_UNPRECONDITIONED:
483:       case KSP_NORM_NATURAL:
484:         PetscCall(VecNorm(r, NORM_2, &rnorm));
485:         break;
486:       default:
487:         rnorm = 0.0;
488:         break;
489:       }
490:       KSPCheckNorm(ksp, rnorm);
491:       PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
492:       ksp->rnorm = rnorm;
493:       PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
494:       PetscCall(KSPLogResidualHistory(ksp, rnorm));
495:       PetscCall(KSPMonitor(ksp, i, rnorm));
496:       PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
497:       if (ksp->reason) break;
498:       if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, p[kp1])); /*  p[kp1] = B^{-1}r  */ }
499:     } else {
500:       PetscCall(KSP_PCApply(ksp, r, p[kp1])); /*  p[kp1] = B^{-1}r  */
501:     }
502:     ksp->vec_sol = p[k];
503:     PetscCall(KSPLogErrorHistory(ksp));

505:     c[kp1] = 2.0 * mu * c[k] - c[km1];
506:     omega  = omegaprod * c[k] / c[kp1];

508:     /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
509:     PetscCall(VecAXPBYPCZ(p[kp1], 1.0 - omega, omega, omega * Gamma * scale, p[km1], p[k]));

511:     ktmp = km1;
512:     km1  = k;
513:     k    = kp1;
514:     kp1  = ktmp;
515:   }
516:   if (!ksp->reason) {
517:     if (ksp->normtype) {
518:       PetscCall(KSP_MatMult(ksp, Amat, p[k], r)); /*  r = b - Ap[k]    */
519:       PetscCall(VecAYPX(r, -1.0, b));
520:       switch (ksp->normtype) {
521:       case KSP_NORM_PRECONDITIONED:
522:         PetscCall(KSP_PCApply(ksp, r, p[kp1])); /* p[kp1] = B^{-1}r */
523:         PetscCall(VecNorm(p[kp1], NORM_2, &rnorm));
524:         break;
525:       case KSP_NORM_UNPRECONDITIONED:
526:       case KSP_NORM_NATURAL:
527:         PetscCall(VecNorm(r, NORM_2, &rnorm));
528:         break;
529:       default:
530:         rnorm = 0.0;
531:         break;
532:       }
533:       KSPCheckNorm(ksp, rnorm);
534:       PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
535:       ksp->rnorm = rnorm;
536:       PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
537:       PetscCall(KSPLogResidualHistory(ksp, rnorm));
538:       PetscCall(KSPMonitor(ksp, i, rnorm));
539:     }
540:     if (ksp->its >= ksp->max_it) {
541:       if (ksp->normtype != KSP_NORM_NONE) {
542:         PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
543:         if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
544:       } else ksp->reason = KSP_CONVERGED_ITS;
545:     }
546:   }

548:   /* make sure solution is in vector x */
549:   ksp->vec_sol = sol_orig;
550:   if (k) PetscCall(VecCopy(p[k], sol_orig));
551:   if (ksp->reason == KSP_CONVERGED_ITS) PetscCall(KSPLogErrorHistory(ksp));
552:   PetscFunctionReturn(PETSC_SUCCESS);
553: }

555: static PetscErrorCode KSPSolve_Chebyshev_FourthKind(KSP ksp)
556: {
557:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
558:   PetscInt       i;
559:   PetscScalar    scale, rScale, dScale;
560:   PetscReal      rnorm = 0.0, emax, emin;
561:   Vec            x, b, d, r, Br;
562:   Mat            Amat, Pmat;
563:   PetscBool      diagonalscale;
564:   PetscReal     *betas = cheb->betas;

566:   PetscFunctionBegin;
567:   PetscCall(PCGetDiagonalScale(ksp->pc, &diagonalscale));
568:   PetscCheck(!diagonalscale, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Krylov method %s does not support diagonal scaling", ((PetscObject)ksp)->type_name);

570:   PetscCall(PCGetOperators(ksp->pc, &Amat, &Pmat));
571:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
572:   ksp->its = 0;
573:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));

575:   x  = ksp->vec_sol;
576:   b  = ksp->vec_rhs;
577:   r  = ksp->work[0];
578:   d  = ksp->work[1];
579:   Br = ksp->work[2];

581:   PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
582:   /* use scale*B as our preconditioner */
583:   scale = 1.0 / emax;

585:   if (!ksp->guess_zero) {
586:     PetscCall(KSP_MatMult(ksp, Amat, x, r)); /*  r = b - A*x */
587:     PetscCall(VecAYPX(r, -1.0, b));
588:   } else {
589:     PetscCall(VecCopy(b, r));
590:   }

592:   /* calculate residual norm if requested, we have done one iteration */
593:   if (ksp->normtype) {
594:     switch (ksp->normtype) {
595:     case KSP_NORM_PRECONDITIONED:
596:       PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */
597:       PetscCall(VecNorm(Br, NORM_2, &rnorm));
598:       break;
599:     case KSP_NORM_UNPRECONDITIONED:
600:     case KSP_NORM_NATURAL:
601:       PetscCall(VecNorm(r, NORM_2, &rnorm));
602:       break;
603:     default:
604:       SETERRQ(PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "%s", KSPNormTypes[ksp->normtype]);
605:     }
606:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
607:     ksp->rnorm = rnorm;
608:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
609:     PetscCall(KSPLogResidualHistory(ksp, rnorm));
610:     PetscCall(KSPLogErrorHistory(ksp));
611:     PetscCall(KSPMonitor(ksp, 0, rnorm));
612:     PetscCall((*ksp->converged)(ksp, 0, rnorm, &ksp->reason, ksp->cnvP));
613:   } else ksp->reason = KSP_CONVERGED_ITERATING;
614:   if (ksp->reason || ksp->max_it == 0) {
615:     if (ksp->max_it == 0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
616:     PetscFunctionReturn(PETSC_SUCCESS);
617:   }
618:   if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, Br)); /* Br = B^{-1}r */ }
619:   PetscCall(VecAXPBY(d, 4.0 / 3.0 * scale, 0.0, Br)); /* d = 4/3 * scale B^{-1}r */
620:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
621:   ksp->its = 1;
622:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));

624:   for (i = 1; i < ksp->max_it; i++) {
625:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
626:     ksp->its++;
627:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));

629:     PetscCall(VecAXPBY(x, betas[i - 1], 1.0, d)); /* x = x + \beta_k d */

631:     PetscCall(KSP_MatMult(ksp, Amat, d, Br)); /*  r = r - Ad */
632:     PetscCall(VecAXPBY(r, -1.0, 1.0, Br));

634:     /* calculate residual norm if requested */
635:     if (ksp->normtype) {
636:       switch (ksp->normtype) {
637:       case KSP_NORM_PRECONDITIONED:
638:         PetscCall(KSP_PCApply(ksp, r, Br)); /*  Br = B^{-1}r  */
639:         PetscCall(VecNorm(Br, NORM_2, &rnorm));
640:         break;
641:       case KSP_NORM_UNPRECONDITIONED:
642:       case KSP_NORM_NATURAL:
643:         PetscCall(VecNorm(r, NORM_2, &rnorm));
644:         break;
645:       default:
646:         rnorm = 0.0;
647:         break;
648:       }
649:       KSPCheckNorm(ksp, rnorm);
650:       PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
651:       ksp->rnorm = rnorm;
652:       PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
653:       PetscCall(KSPLogResidualHistory(ksp, rnorm));
654:       PetscCall(KSPMonitor(ksp, i, rnorm));
655:       PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
656:       if (ksp->reason) break;
657:       if (ksp->normtype != KSP_NORM_PRECONDITIONED) { PetscCall(KSP_PCApply(ksp, r, Br)); /*  Br = B^{-1}r  */ }
658:     } else {
659:       PetscCall(KSP_PCApply(ksp, r, Br)); /*  Br = B^{-1}r  */
660:     }
661:     PetscCall(KSPLogErrorHistory(ksp));

663:     rScale = scale * (8.0 * i + 4.0) / (2.0 * i + 3.0);
664:     dScale = (2.0 * i - 1.0) / (2.0 * i + 3.0);

666:     /* d_k+1 = \dfrac{2k-1}{2k+3} d_k + \dfrac{8k+4}{2k+3} \dfrac{1}{\rho(SA)} Br */
667:     PetscCall(VecAXPBY(d, rScale, dScale, Br));
668:   }

670:   /* on last pass, update solution vector */
671:   PetscCall(VecAXPBY(x, betas[ksp->max_it - 1], 1.0, d)); /* x = x + d */

673:   if (!ksp->reason) {
674:     if (ksp->normtype) {
675:       PetscCall(KSP_MatMult(ksp, Amat, x, r)); /*  r = b - Ax    */
676:       PetscCall(VecAYPX(r, -1.0, b));
677:       switch (ksp->normtype) {
678:       case KSP_NORM_PRECONDITIONED:
679:         PetscCall(KSP_PCApply(ksp, r, Br)); /* Br= B^{-1}r */
680:         PetscCall(VecNorm(Br, NORM_2, &rnorm));
681:         break;
682:       case KSP_NORM_UNPRECONDITIONED:
683:       case KSP_NORM_NATURAL:
684:         PetscCall(VecNorm(r, NORM_2, &rnorm));
685:         break;
686:       default:
687:         rnorm = 0.0;
688:         break;
689:       }
690:       KSPCheckNorm(ksp, rnorm);
691:       PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
692:       ksp->rnorm = rnorm;
693:       PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
694:       PetscCall(KSPLogResidualHistory(ksp, rnorm));
695:       PetscCall(KSPMonitor(ksp, i, rnorm));
696:     }
697:     if (ksp->its >= ksp->max_it) {
698:       if (ksp->normtype != KSP_NORM_NONE) {
699:         PetscCall((*ksp->converged)(ksp, i, rnorm, &ksp->reason, ksp->cnvP));
700:         if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
701:       } else ksp->reason = KSP_CONVERGED_ITS;
702:     }
703:   }

705:   if (ksp->reason == KSP_CONVERGED_ITS) PetscCall(KSPLogErrorHistory(ksp));
706:   PetscFunctionReturn(PETSC_SUCCESS);
707: }

709: static PetscErrorCode KSPView_Chebyshev(KSP ksp, PetscViewer viewer)
710: {
711:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
712:   PetscBool      iascii;

714:   PetscFunctionBegin;
715:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
716:   if (iascii) {
717:     switch (cheb->chebykind) {
718:     case KSP_CHEBYSHEV_FIRST:
719:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Chebyshev polynomial of first kind\n"));
720:       break;
721:     case KSP_CHEBYSHEV_FOURTH:
722:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Chebyshev polynomial of fourth kind\n"));
723:       break;
724:     case KSP_CHEBYSHEV_OPT_FOURTH:
725:       PetscCall(PetscViewerASCIIPrintf(viewer, "  Chebyshev polynomial of opt. fourth kind\n"));
726:       break;
727:     }
728:     PetscReal emax, emin;
729:     PetscCall(KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin));
730:     PetscCall(PetscViewerASCIIPrintf(viewer, "  eigenvalue targets used: min %g, max %g\n", (double)emin, (double)emax));
731:     if (cheb->kspest) {
732:       PetscCall(PetscViewerASCIIPrintf(viewer, "  eigenvalues estimated via %s: min %g, max %g\n", ((PetscObject)(cheb->kspest))->type_name, (double)cheb->emin_computed, (double)cheb->emax_computed));
733:       PetscCall(PetscViewerASCIIPrintf(viewer, "  eigenvalues estimated using %s with transform: [%g %g; %g %g]\n", ((PetscObject)cheb->kspest)->type_name, (double)cheb->tform[0], (double)cheb->tform[1], (double)cheb->tform[2], (double)cheb->tform[3]));
734:       PetscCall(PetscViewerASCIIPushTab(viewer));
735:       PetscCall(KSPView(cheb->kspest, viewer));
736:       PetscCall(PetscViewerASCIIPopTab(viewer));
737:       if (cheb->usenoisy) PetscCall(PetscViewerASCIIPrintf(viewer, "  estimating eigenvalues using noisy right hand side\n"));
738:     } else if (cheb->emax_provided != 0.) {
739:       PetscCall(PetscViewerASCIIPrintf(viewer, "  eigenvalues provided (min %g, max %g) with transform: [%g %g; %g %g]\n", (double)cheb->emin_provided, (double)cheb->emax_provided, (double)cheb->tform[0], (double)cheb->tform[1], (double)cheb->tform[2],
740:                                        (double)cheb->tform[3]));
741:     }
742:   }
743:   PetscFunctionReturn(PETSC_SUCCESS);
744: }

746: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
747: {
748:   KSP_Chebyshev   *cheb = (KSP_Chebyshev *)ksp->data;
749:   PetscBool        isset, flg;
750:   Mat              Pmat, Amat;
751:   PetscObjectId    amatid, pmatid;
752:   PetscObjectState amatstate, pmatstate;

754:   PetscFunctionBegin;
755:   switch (cheb->chebykind) {
756:   case KSP_CHEBYSHEV_FIRST:
757:     ksp->ops->solve = KSPSolve_Chebyshev_FirstKind;
758:     break;
759:   case KSP_CHEBYSHEV_FOURTH:
760:   case KSP_CHEBYSHEV_OPT_FOURTH:
761:     ksp->ops->solve = KSPSolve_Chebyshev_FourthKind;
762:     break;
763:   }

765:   if (ksp->max_it > cheb->num_betas_alloc) {
766:     PetscCall(PetscFree(cheb->betas));
767:     PetscCall(PetscMalloc1(ksp->max_it, &cheb->betas));
768:     cheb->num_betas_alloc = ksp->max_it;
769:   }

771:   // coefficients for 4th-kind Chebyshev
772:   for (PetscInt i = 0; i < ksp->max_it; i++) cheb->betas[i] = 1.0;

774:   // coefficients for optimized 4th-kind Chebyshev
775:   if (cheb->chebykind == KSP_CHEBYSHEV_OPT_FOURTH) PetscCall(KSPChebyshevGetBetas_Private(ksp));

777:   PetscCall(KSPSetWorkVecs(ksp, 3));
778:   if (cheb->emin == 0. || cheb->emax == 0.) { // User did not specify eigenvalues
779:     PC pc;

781:     PetscCall(KSPGetPC(ksp, &pc));
782:     PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &flg));
783:     if (!flg) { // Provided estimates are only relevant for Jacobi
784:       cheb->emax_provided = 0;
785:       cheb->emin_provided = 0;
786:     }
787:     /* We need to estimate eigenvalues */
788:     if (!cheb->kspest) PetscCall(KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE));
789:   }
790:   if (cheb->kspest) {
791:     PetscCall(KSPGetOperators(ksp, &Amat, &Pmat));
792:     PetscCall(MatIsSPDKnown(Pmat, &isset, &flg));
793:     if (isset && flg) {
794:       const char *prefix;

796:       PetscCall(KSPGetOptionsPrefix(cheb->kspest, &prefix));
797:       PetscCall(PetscOptionsHasName(NULL, prefix, "-ksp_type", &flg));
798:       if (!flg) PetscCall(KSPSetType(cheb->kspest, KSPCG));
799:     }
800:     PetscCall(PetscObjectGetId((PetscObject)Amat, &amatid));
801:     PetscCall(PetscObjectGetId((PetscObject)Pmat, &pmatid));
802:     PetscCall(PetscObjectStateGet((PetscObject)Amat, &amatstate));
803:     PetscCall(PetscObjectStateGet((PetscObject)Pmat, &pmatstate));
804:     if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
805:       PetscReal          max = 0.0, min = 0.0;
806:       Vec                B;
807:       KSPConvergedReason reason;

809:       PetscCall(KSPSetPC(cheb->kspest, ksp->pc));
810:       if (cheb->usenoisy) {
811:         B = ksp->work[1];
812:         PetscCall(KSPSetNoisy_Private(B));
813:       } else {
814:         PetscBool change;

816:         PetscCheck(ksp->vec_rhs, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Chebyshev must use a noisy right hand side to estimate the eigenvalues when no right hand side is available");
817:         PetscCall(PCPreSolveChangeRHS(ksp->pc, &change));
818:         if (change) {
819:           B = ksp->work[1];
820:           PetscCall(VecCopy(ksp->vec_rhs, B));
821:         } else B = ksp->vec_rhs;
822:       }
823:       PetscCall(KSPSolve(cheb->kspest, B, ksp->work[0]));
824:       PetscCall(KSPGetConvergedReason(cheb->kspest, &reason));
825:       if (reason == KSP_DIVERGED_ITS) {
826:         PetscCall(PetscInfo(ksp, "Eigen estimator ran for prescribed number of iterations\n"));
827:       } else if (reason == KSP_DIVERGED_PC_FAILED) {
828:         PetscInt       its;
829:         PCFailedReason pcreason;

831:         PetscCall(KSPGetIterationNumber(cheb->kspest, &its));
832:         if (ksp->normtype == KSP_NORM_NONE) {
833:           PetscInt sendbuf, recvbuf;

835:           PetscCall(PCGetFailedReasonRank(ksp->pc, &pcreason));
836:           sendbuf = (PetscInt)pcreason;
837:           PetscCall(MPIU_Allreduce(&sendbuf, &recvbuf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)ksp)));
838:           PetscCall(PCSetFailedReason(ksp->pc, (PCFailedReason)recvbuf));
839:         }
840:         PetscCall(PCGetFailedReason(ksp->pc, &pcreason));
841:         ksp->reason = KSP_DIVERGED_PC_FAILED;
842:         PetscCall(PetscInfo(ksp, "Eigen estimator failed: %s %s at iteration %" PetscInt_FMT, KSPConvergedReasons[reason], PCFailedReasons[pcreason], its));
843:         PetscFunctionReturn(PETSC_SUCCESS);
844:       } else if (reason == KSP_CONVERGED_RTOL || reason == KSP_CONVERGED_ATOL) {
845:         PetscCall(PetscInfo(ksp, "Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n"));
846:       } else if (reason < 0) {
847:         PetscCall(PetscInfo(ksp, "Eigen estimator failed %s, using estimates anyway\n", KSPConvergedReasons[reason]));
848:       }

850:       PetscCall(KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest, &min, &max));
851:       PetscCall(KSPSetPC(cheb->kspest, NULL));

853:       cheb->emin_computed = min;
854:       cheb->emax_computed = max;

856:       cheb->amatid    = amatid;
857:       cheb->pmatid    = pmatid;
858:       cheb->amatstate = amatstate;
859:       cheb->pmatstate = pmatstate;
860:     }
861:   }
862:   PetscFunctionReturn(PETSC_SUCCESS);
863: }

865: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
866: {
867:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;

869:   PetscFunctionBegin;
870:   PetscCall(PetscFree(cheb->betas));
871:   PetscCall(KSPDestroy(&cheb->kspest));
872:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", NULL));
873:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", NULL));
874:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", NULL));
875:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetKind_C", NULL));
876:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevGetKind_C", NULL));
877:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", NULL));
878:   PetscCall(KSPDestroyDefault(ksp));
879:   PetscFunctionReturn(PETSC_SUCCESS);
880: }

882: /*MC
883:      KSPCHEBYSHEV - The preconditioned Chebyshev iterative method

885:    Options Database Keys:
886: +   -ksp_chebyshev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
887:                   of the preconditioned operator. If these are accurate you will get much faster convergence.
888: .   -ksp_chebyshev_esteig <a,b,c,d> - estimate eigenvalues using a Krylov method, then use this
889:                          transform for Chebyshev eigenvalue bounds (`KSPChebyshevEstEigSet()`)
890: .   -ksp_chebyshev_esteig_steps - number of estimation steps
891: -   -ksp_chebyshev_esteig_noisy - use noisy number generator to create right hand side for eigenvalue estimator

893:    Level: beginner

895:    Notes:
896:    The Chebyshev method requires both the matrix and preconditioner to be symmetric positive (semi) definite, but it can work as a smoother in other situations

898:    Only support for left preconditioning.

900:    Chebyshev is configured as a smoother by default, targeting the "upper" part of the spectrum.

902:    The user should call `KSPChebyshevSetEigenvalues()` to get eigenvalue estimates.

904: .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
905:           `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigSetUseNoisy()`
906:           `KSPRICHARDSON`, `KSPCG`, `PCMG`
907: M*/

909: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
910: {
911:   KSP_Chebyshev *chebyshevP;

913:   PetscFunctionBegin;
914:   PetscCall(PetscNew(&chebyshevP));

916:   ksp->data = (void *)chebyshevP;
917:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
918:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2));
919:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
920:   PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));

922:   chebyshevP->emin = 0.;
923:   chebyshevP->emax = 0.;

925:   chebyshevP->tform[0] = 0.0;
926:   chebyshevP->tform[1] = 0.1;
927:   chebyshevP->tform[2] = 0;
928:   chebyshevP->tform[3] = 1.1;
929:   chebyshevP->eststeps = 10;
930:   chebyshevP->usenoisy = PETSC_TRUE;
931:   ksp->setupnewmatrix  = PETSC_TRUE;

933:   ksp->ops->setup          = KSPSetUp_Chebyshev;
934:   ksp->ops->destroy        = KSPDestroy_Chebyshev;
935:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
936:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
937:   ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
938:   ksp->ops->view           = KSPView_Chebyshev;
939:   ksp->ops->reset          = KSPReset_Chebyshev;

941:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", KSPChebyshevSetEigenvalues_Chebyshev));
942:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", KSPChebyshevEstEigSet_Chebyshev));
943:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", KSPChebyshevEstEigSetUseNoisy_Chebyshev));
944:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetKind_C", KSPChebyshevSetKind_Chebyshev));
945:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevGetKind_C", KSPChebyshevGetKind_Chebyshev));
946:   PetscCall(PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", KSPChebyshevEstEigGetKSP_Chebyshev));
947:   PetscFunctionReturn(PETSC_SUCCESS);
948: }