Actual source code: cheby.c

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    115: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp, PetscBool use)
    116: {
    117:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;

    119:   PetscFunctionBegin;
    120:   cheb->usenoisy = use;
    121:   PetscFunctionReturn(PETSC_SUCCESS);
    122: }

    124: static PetscErrorCode KSPChebyshevSetKind_Chebyshev(KSP ksp, KSPChebyshevKind kind)
    125: {
    126:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;

    128:   PetscFunctionBegin;
    129:   cheb->chebykind = kind;
    130:   PetscFunctionReturn(PETSC_SUCCESS);
    131: }

    133: static PetscErrorCode KSPChebyshevGetKind_Chebyshev(KSP ksp, KSPChebyshevKind *kind)
    134: {
    135:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;

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

    144:   Logically Collective

    146:   Input Parameters:
    147: + ksp  - the Krylov space context
    148: . emax - the eigenvalue maximum estimate
    149: - emin - the eigenvalue minimum estimate

    151:   Options Database Key:
    152: . -ksp_chebyshev_eigenvalues emin,emax - extreme eigenvalues

    154:   Level: intermediate

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

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

    165: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`,
    166: @*/
    167: PetscErrorCode KSPChebyshevSetEigenvalues(KSP ksp, PetscReal emax, PetscReal emin)
    168: {
    169:   PetscFunctionBegin;
    173:   PetscTryMethod(ksp, "KSPChebyshevSetEigenvalues_C", (KSP, PetscReal, PetscReal), (ksp, emax, emin));
    174:   PetscFunctionReturn(PETSC_SUCCESS);
    175: }

    177: /*@
    178:   KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev

    180:   Logically Collective

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

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

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

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

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

    205:   The eigenvalues are estimated using the Lanczos (`KSPCG`) or Arnoldi (`KSPGMRES`) process depending on if the operator is
    206:   symmetric definite or not.

    208:   Level: intermediate

    210: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSetUseNoisy()`, `KSPChebyshevEstEigGetKSP()`
    211: @*/
    212: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp, PetscReal a, PetscReal b, PetscReal c, PetscReal d)
    213: {
    214:   PetscFunctionBegin;
    220:   PetscTryMethod(ksp, "KSPChebyshevEstEigSet_C", (KSP, PetscReal, PetscReal, PetscReal, PetscReal), (ksp, a, b, c, d));
    221:   PetscFunctionReturn(PETSC_SUCCESS);
    222: }

    224: /*@
    225:   KSPChebyshevEstEigSetUseNoisy - use a noisy random number generated right-hand side to estimate the extreme eigenvalues instead of the given right-hand side

    227:   Logically Collective

    229:   Input Parameters:
    230: + ksp - linear solver context
    231: - use - `PETSC_TRUE` to use noisy

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

    236:   Level: intermediate

    238:   Note:
    239:   This allegedly works better for multigrid smoothers

    241: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigGetKSP()`
    242: @*/
    243: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp, PetscBool use)
    244: {
    245:   PetscFunctionBegin;
    248:   PetscTryMethod(ksp, "KSPChebyshevEstEigSetUseNoisy_C", (KSP, PetscBool), (ksp, use));
    249:   PetscFunctionReturn(PETSC_SUCCESS);
    250: }

    252: /*@
    253:   KSPChebyshevEstEigGetKSP - Get the Krylov method context used to estimate the eigenvalues for the Chebyshev method.

    255:   Input Parameter:
    256: . ksp - the Krylov space context

    258:   Output Parameter:
    259: . kspest - the eigenvalue estimation Krylov space context

    261:   Level: advanced

    263:   Notes:
    264:   If a Krylov method is not being used for this purpose, `NULL` is returned.  The reference count of the returned `KSP` is
    265:   not incremented: it should not be destroyed by the user.

    267: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`
    268: @*/
    269: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
    270: {
    271:   PetscFunctionBegin;
    273:   PetscAssertPointer(kspest, 2);
    274:   *kspest = NULL;
    275:   PetscTryMethod(ksp, "KSPChebyshevEstEigGetKSP_C", (KSP, KSP *), (ksp, kspest));
    276:   PetscFunctionReturn(PETSC_SUCCESS);
    277: }

    279: /*@
    280:   KSPChebyshevSetKind - set the kind of Chebyshev polynomial to use

    282:   Logically Collective

    284:   Input Parameters:
    285: + ksp  - Linear solver context
    286: - kind - The kind of Chebyshev polynomial to use, see `KSPChebyshevKind`, one of `KSP_CHEBYSHEV_FIRST`, `KSP_CHEBYSHEV_FOURTH`, or `KSP_CHEBYSHEV_OPT_FOURTH`

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

    291:   Level: intermediate

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

    301: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevKind`, `KSPChebyshevGetKind()`, `KSP_CHEBYSHEV_FIRST`, `KSP_CHEBYSHEV_FOURTH`, `KSP_CHEBYSHEV_OPT_FOURTH`
    302: @*/
    303: PetscErrorCode KSPChebyshevSetKind(KSP ksp, KSPChebyshevKind kind)
    304: {
    305:   PetscFunctionBegin;
    308:   PetscUseMethod(ksp, "KSPChebyshevSetKind_C", (KSP, KSPChebyshevKind), (ksp, kind));
    309:   PetscFunctionReturn(PETSC_SUCCESS);
    310: }

    312: /*@
    313:   KSPChebyshevGetKind - get the kind of Chebyshev polynomial to use

    315:   Logically Collective

    317:   Input Parameters:
    318: + ksp  - Linear solver context
    319: - kind - The kind of Chebyshev polynomial used

    321:   Level: intermediate

    323: .seealso: [](ch_ksp), `KSPCHEBYSHEV`, `KSPChebyshevKind`, `KSPChebyshevSetKind()`, `KSP_CHEBYSHEV_FIRST`, `KSP_CHEBYSHEV_FOURTH`, `KSP_CHEBYSHEV_OPT_FOURTH`
    324: @*/
    325: PetscErrorCode KSPChebyshevGetKind(KSP ksp, KSPChebyshevKind *kind)
    326: {
    327:   PetscFunctionBegin;
    329:   PetscUseMethod(ksp, "KSPChebyshevGetKind_C", (KSP, KSPChebyshevKind *), (ksp, kind));
    330:   PetscFunctionReturn(PETSC_SUCCESS);
    331: }

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

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

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

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

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

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

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

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

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

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

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

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

    428:   c[km1] = 1.0;
    429:   c[k]   = mu;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    711: static PetscErrorCode KSPView_Chebyshev(KSP ksp, PetscViewer viewer)
    712: {
    713:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
    714:   PetscBool      isascii;

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

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

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

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

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

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

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

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

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

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

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

    834:         PetscCall(KSPGetIterationNumber(cheb->kspest, &its));
    835:         if (ksp->normtype == KSP_NORM_NONE) PetscCall(PCReduceFailedReason(ksp->pc));
    836:         PetscCall(PCGetFailedReason(ksp->pc, &pcreason));
    837:         ksp->reason = KSP_DIVERGED_PC_FAILED;
    838:         PetscCall(PetscInfo(ksp, "Eigen estimator failed: %s %s at iteration %" PetscInt_FMT "\n", KSPConvergedReasons[reason], PCFailedReasons[pcreason], its));
    839:         PetscFunctionReturn(PETSC_SUCCESS);
    840:       } else if (reason == KSP_CONVERGED_RTOL || reason == KSP_CONVERGED_ATOL) {
    841:         PetscCall(PetscInfo(ksp, "Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n"));
    842:       } else if (reason < 0) {
    843:         PetscCall(PetscInfo(ksp, "Eigen estimator failed %s, using estimates anyway\n", KSPConvergedReasons[reason]));
    844:       }

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

    849:       cheb->emin_computed = min;
    850:       cheb->emax_computed = max;

    852:       cheb->amatid    = amatid;
    853:       cheb->pmatid    = pmatid;
    854:       cheb->amatstate = amatstate;
    855:       cheb->pmatstate = pmatstate;
    856:     }
    857:   }
    858:   if (ksp->monitor[0] == (PetscErrorCode (*)(KSP, PetscInt, PetscReal, void *))KSPMonitorResidual && !ksp->normtype) PetscCall(KSPSetNormType(ksp, KSP_NORM_PRECONDITIONED));
    859:   PetscFunctionReturn(PETSC_SUCCESS);
    860: }

    862: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
    863: {
    864:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;

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

    879: /*MC
    880:      KSPCHEBYSHEV - The preconditioned Chebyshev iterative method

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

    890:    Level: beginner

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

    896:    Only has support for left preconditioning.

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

    900:    By default this uses `KSPGMRES` to estimate the extreme eigenvalues, if the matrix is known to be SPD then it uses `KSPCG` to estimate the eigenvalues.
    901:    See `MatIsSPDKnown()` for how to indicate a `Mat`, matrix is SPD.

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

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

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

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

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

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

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

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