Actual source code: cheby.c


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

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

  8:   if (cheb->kspest) KSPReset(cheb->kspest);
  9:   return 0;
 10: }

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

 20:   KSPGetIterationNumber(kspest, &n);
 21:   PetscMalloc2(n, &re, n, &im);
 22:   KSPComputeEigenvalues(kspest, n, re, im, &neig);
 23:   min = PETSC_MAX_REAL;
 24:   max = PETSC_MIN_REAL;
 25:   for (n = 0; n < neig; n++) {
 26:     min = PetscMin(min, re[n]);
 27:     max = PetscMax(max, re[n]);
 28:   }
 29:   PetscFree2(re, im);
 30:   *emax = max;
 31:   *emin = min;
 32:   return 0;
 33: }

 35: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
 36: {
 37:   KSP_Chebyshev   *cheb = (KSP_Chebyshev *)ksp->data;
 38:   PetscBool        isset, flg;
 39:   Mat              Pmat, Amat;
 40:   PetscObjectId    amatid, pmatid;
 41:   PetscObjectState amatstate, pmatstate;

 43:   KSPSetWorkVecs(ksp, 3);
 44:   if (cheb->emin == 0. || cheb->emax == 0.) { // User did not specify eigenvalues
 45:     PC pc;
 46:     KSPGetPC(ksp, &pc);
 47:     PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &flg);
 48:     if (!flg) { // Provided estimates are only relevant for Jacobi
 49:       cheb->emax_provided = 0;
 50:       cheb->emin_provided = 0;
 51:     }
 52:     if (!cheb->kspest) { /* We need to estimate eigenvalues */
 53:       KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);
 54:     }
 55:   }
 56:   if (cheb->kspest) {
 57:     KSPGetOperators(ksp, &Amat, &Pmat);
 58:     MatIsSPDKnown(Pmat, &isset, &flg);
 59:     if (isset && flg) {
 60:       const char *prefix;
 61:       KSPGetOptionsPrefix(cheb->kspest, &prefix);
 62:       PetscOptionsHasName(NULL, prefix, "-ksp_type", &flg);
 63:       if (!flg) KSPSetType(cheb->kspest, KSPCG);
 64:     }
 65:     PetscObjectGetId((PetscObject)Amat, &amatid);
 66:     PetscObjectGetId((PetscObject)Pmat, &pmatid);
 67:     PetscObjectStateGet((PetscObject)Amat, &amatstate);
 68:     PetscObjectStateGet((PetscObject)Pmat, &pmatstate);
 69:     if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
 70:       PetscReal          max = 0.0, min = 0.0;
 71:       Vec                B;
 72:       KSPConvergedReason reason;
 73:       KSPSetPC(cheb->kspest, ksp->pc);
 74:       if (cheb->usenoisy) {
 75:         B = ksp->work[1];
 76:         KSPSetNoisy_Private(B);
 77:       } else {
 78:         PetscBool change;

 81:         PCPreSolveChangeRHS(ksp->pc, &change);
 82:         if (change) {
 83:           B = ksp->work[1];
 84:           VecCopy(ksp->vec_rhs, B);
 85:         } else B = ksp->vec_rhs;
 86:       }
 87:       KSPSolve(cheb->kspest, B, ksp->work[0]);
 88:       KSPGetConvergedReason(cheb->kspest, &reason);
 89:       if (reason == KSP_DIVERGED_ITS) {
 90:         PetscInfo(ksp, "Eigen estimator ran for prescribed number of iterations\n");
 91:       } else if (reason == KSP_DIVERGED_PC_FAILED) {
 92:         PetscInt       its;
 93:         PCFailedReason pcreason;

 95:         KSPGetIterationNumber(cheb->kspest, &its);
 96:         if (ksp->normtype == KSP_NORM_NONE) {
 97:           PetscInt sendbuf, recvbuf;
 98:           PCGetFailedReasonRank(ksp->pc, &pcreason);
 99:           sendbuf = (PetscInt)pcreason;
100:           MPI_Allreduce(&sendbuf, &recvbuf, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)ksp));
101:           PCSetFailedReason(ksp->pc, (PCFailedReason)recvbuf);
102:         }
103:         PCGetFailedReason(ksp->pc, &pcreason);
104:         ksp->reason = KSP_DIVERGED_PC_FAILED;
105:         PetscInfo(ksp, "Eigen estimator failed: %s %s at iteration %" PetscInt_FMT, KSPConvergedReasons[reason], PCFailedReasons[pcreason], its);
106:         return 0;
107:       } else if (reason == KSP_CONVERGED_RTOL || reason == KSP_CONVERGED_ATOL) {
108:         PetscInfo(ksp, "Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n");
109:       } else if (reason < 0) {
110:         PetscInfo(ksp, "Eigen estimator failed %s, using estimates anyway\n", KSPConvergedReasons[reason]);
111:       }

113:       KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest, &min, &max);
114:       KSPSetPC(cheb->kspest, NULL);

116:       cheb->emin_computed = min;
117:       cheb->emax_computed = max;

119:       cheb->amatid    = amatid;
120:       cheb->pmatid    = pmatid;
121:       cheb->amatstate = amatstate;
122:       cheb->pmatstate = pmatstate;
123:     }
124:   }
125:   return 0;
126: }

128: static PetscErrorCode KSPChebyshevGetEigenvalues_Chebyshev(KSP ksp, PetscReal *emax, PetscReal *emin)
129: {
130:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;

132:   *emax = 0;
133:   *emin = 0;
134:   if (cheb->emax != 0.) {
135:     *emax = cheb->emax;
136:   } else if (cheb->emax_computed != 0.) {
137:     *emax = cheb->tform[2] * cheb->emin_computed + cheb->tform[3] * cheb->emax_computed;
138:   } else if (cheb->emax_provided != 0.) {
139:     *emax = cheb->tform[2] * cheb->emin_provided + cheb->tform[3] * cheb->emax_provided;
140:   }
141:   if (cheb->emin != 0.) {
142:     *emin = cheb->emin;
143:   } else if (cheb->emin_computed != 0.) {
144:     *emin = cheb->tform[0] * cheb->emin_computed + cheb->tform[1] * cheb->emax_computed;
145:   } else if (cheb->emin_provided != 0.) {
146:     *emin = cheb->tform[0] * cheb->emin_provided + cheb->tform[1] * cheb->emax_provided;
147:   }
148:   return 0;
149: }

151: static PetscErrorCode KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp, PetscReal emax, PetscReal emin)
152: {
153:   KSP_Chebyshev *chebyshevP = (KSP_Chebyshev *)ksp->data;

157:   chebyshevP->emax = emax;
158:   chebyshevP->emin = emin;

160:   KSPChebyshevEstEigSet(ksp, 0., 0., 0., 0.); /* Destroy any estimation setup */
161:   return 0;
162: }

164: static PetscErrorCode KSPChebyshevEstEigSet_Chebyshev(KSP ksp, PetscReal a, PetscReal b, PetscReal c, PetscReal d)
165: {
166:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;

168:   if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
169:     if ((cheb->emin_provided == 0. || cheb->emax_provided == 0.) && !cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
170:       KSPCreate(PetscObjectComm((PetscObject)ksp), &cheb->kspest);
171:       KSPSetErrorIfNotConverged(cheb->kspest, ksp->errorifnotconverged);
172:       PetscObjectIncrementTabLevel((PetscObject)cheb->kspest, (PetscObject)ksp, 1);
173:       /* use PetscObjectSet/AppendOptionsPrefix() instead of KSPSet/AppendOptionsPrefix() so that the PC prefix is not changed */
174:       PetscObjectSetOptionsPrefix((PetscObject)cheb->kspest, ((PetscObject)ksp)->prefix);
175:       PetscObjectAppendOptionsPrefix((PetscObject)cheb->kspest, "esteig_");
176:       KSPSetSkipPCSetFromOptions(cheb->kspest, PETSC_TRUE);

178:       KSPSetComputeEigenvalues(cheb->kspest, PETSC_TRUE);

180:       /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
181:       KSPSetTolerances(cheb->kspest, 1.e-12, PETSC_DEFAULT, PETSC_DEFAULT, cheb->eststeps);
182:     }
183:     if (a >= 0) cheb->tform[0] = a;
184:     if (b >= 0) cheb->tform[1] = b;
185:     if (c >= 0) cheb->tform[2] = c;
186:     if (d >= 0) cheb->tform[3] = d;
187:     cheb->amatid    = 0;
188:     cheb->pmatid    = 0;
189:     cheb->amatstate = -1;
190:     cheb->pmatstate = -1;
191:   } else {
192:     KSPDestroy(&cheb->kspest);
193:   }
194:   return 0;
195: }

197: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp, PetscBool use)
198: {
199:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;

201:   cheb->usenoisy = use;
202:   return 0;
203: }

205: /*@
206:    KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues
207:    of the preconditioned problem.

209:    Logically Collective

211:    Input Parameters:
212: +  ksp - the Krylov space context
213: -  emax, emin - the eigenvalue estimates

215:   Options Database Key:
216: .  -ksp_chebyshev_eigenvalues emin,emax - extreme eigenvalues

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

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

227:    Level: intermediate

229: .seealso: [](chapter_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`,
230: @*/
231: PetscErrorCode KSPChebyshevSetEigenvalues(KSP ksp, PetscReal emax, PetscReal emin)
232: {
236:   PetscTryMethod(ksp, "KSPChebyshevSetEigenvalues_C", (KSP, PetscReal, PetscReal), (ksp, emax, emin));
237:   return 0;
238: }

240: /*@
241:    KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev

243:    Logically Collective

245:    Input Parameters:
246: +  ksp - the Krylov space context
247: .  a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
248: .  b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
249: .  c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
250: -  d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)

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

255:    Notes:
256:    The Chebyshev bounds are set using
257: .vb
258:    minbound = a*minest + b*maxest
259:    maxbound = c*minest + d*maxest
260: .ve
261:    The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
262:    The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.

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

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

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

270:    Level: intermediate

272: .seealso: [](chapter_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigSetUseNoisy()`, `KSPChebyshevEstEigGetKSP()`
273: @*/
274: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp, PetscReal a, PetscReal b, PetscReal c, PetscReal d)
275: {
281:   PetscTryMethod(ksp, "KSPChebyshevEstEigSet_C", (KSP, PetscReal, PetscReal, PetscReal, PetscReal), (ksp, a, b, c, d));
282:   return 0;
283: }

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

288:    Logically Collective

290:    Input Parameters:
291: +  ksp - linear solver context
292: -  use - `PETSC_TRUE` to use noisy

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

297:    Note:
298:     This allegedly works better for multigrid smoothers

300:   Level: intermediate

302: .seealso: [](chapter_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigGetKSP()`
303: @*/
304: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp, PetscBool use)
305: {
307:   PetscTryMethod(ksp, "KSPChebyshevEstEigSetUseNoisy_C", (KSP, PetscBool), (ksp, use));
308:   return 0;
309: }

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

316:   Input Parameters:
317: . ksp - the Krylov space context

319:   Output Parameters:
320: . kspest - the eigenvalue estimation Krylov space context

322:   Level: advanced

324: .seealso: [](chapter_ksp), `KSPCHEBYSHEV`, `KSPChebyshevEstEigSet()`
325: @*/
326: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
327: {
330:   *kspest = NULL;
331:   PetscTryMethod(ksp, "KSPChebyshevEstEigGetKSP_C", (KSP, KSP *), (ksp, kspest));
332:   return 0;
333: }

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

339:   *kspest = cheb->kspest;
340:   return 0;
341: }

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

351:   PetscOptionsHeadBegin(PetscOptionsObject, "KSP Chebyshev Options");
352:   PetscOptionsInt("-ksp_chebyshev_esteig_steps", "Number of est steps in Chebyshev", "", cheb->eststeps, &cheb->eststeps, NULL);
353:   PetscOptionsRealArray("-ksp_chebyshev_eigenvalues", "extreme eigenvalues", "KSPChebyshevSetEigenvalues", eminmax, &neigarg, &flgeig);
354:   if (flgeig) {
356:     KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]);
357:   }
358:   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:       KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);
363:       break;
364:     case 2: /* Base everything on the max eigenvalues */
365:       KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, tform[0], PETSC_DECIDE, tform[1]);
366:       break;
367:     case 4: /* Use the full 2x2 linear transformation */
368:       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:   /* 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) KSPChebyshevEstEigSet(ksp, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE);

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

386: static PetscErrorCode KSPSolve_Chebyshev(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:   PCGetDiagonalScale(ksp->pc, &diagonalscale);

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

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

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

424:   c[km1] = 1.0;
425:   c[k]   = mu;

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

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

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

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

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

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

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

546:   /* make sure solution is in vector x */
547:   ksp->vec_sol = sol_orig;
548:   if (k) VecCopy(p[k], sol_orig);
549:   if (ksp->reason == KSP_CONVERGED_ITS) KSPLogErrorHistory(ksp);
550:   return 0;
551: }

553: static PetscErrorCode KSPView_Chebyshev(KSP ksp, PetscViewer viewer)
554: {
555:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;
556:   PetscBool      iascii;

558:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
559:   if (iascii) {
560:     PetscReal emax, emin;
561:     KSPChebyshevGetEigenvalues_Chebyshev(ksp, &emax, &emin);
562:     PetscViewerASCIIPrintf(viewer, "  eigenvalue targets used: min %g, max %g\n", (double)emin, (double)emax);
563:     if (cheb->kspest) {
564:       PetscViewerASCIIPrintf(viewer, "  eigenvalues estimated via %s: min %g, max %g\n", ((PetscObject)(cheb->kspest))->type_name, (double)cheb->emin_computed, (double)cheb->emax_computed);
565:       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]);
566:       PetscViewerASCIIPushTab(viewer);
567:       KSPView(cheb->kspest, viewer);
568:       PetscViewerASCIIPopTab(viewer);
569:       if (cheb->usenoisy) PetscViewerASCIIPrintf(viewer, "  estimating eigenvalues using noisy right hand side\n");
570:     } else if (cheb->emax_provided != 0.) {
571:       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],
572:                                        (double)cheb->tform[3]));
573:     }
574:   }
575:   return 0;
576: }

578: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
579: {
580:   KSP_Chebyshev *cheb = (KSP_Chebyshev *)ksp->data;

582:   KSPDestroy(&cheb->kspest);
583:   PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", NULL);
584:   PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", NULL);
585:   PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", NULL);
586:   PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", NULL);
587:   KSPDestroyDefault(ksp);
588:   return 0;
589: }

591: /*MC
592:      KSPCHEBYSHEV - The preconditioned Chebyshev iterative method

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

602:    Level: beginner

604:    Notes:
605:    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

607:    Only support for left preconditioning.

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

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

613: .seealso: [](chapter_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`,
614:           `KSPChebyshevSetEigenvalues()`, `KSPChebyshevEstEigSet()`, `KSPChebyshevEstEigSetUseNoisy()`
615:           `KSPRICHARDSON`, `KSPCG`, `PCMG`
616: M*/

618: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
619: {
620:   KSP_Chebyshev *chebyshevP;

622:   PetscNew(&chebyshevP);

624:   ksp->data = (void *)chebyshevP;
625:   KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3);
626:   KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_LEFT, 2);
627:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1);
628:   KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1);

630:   chebyshevP->emin = 0.;
631:   chebyshevP->emax = 0.;

633:   chebyshevP->tform[0] = 0.0;
634:   chebyshevP->tform[1] = 0.1;
635:   chebyshevP->tform[2] = 0;
636:   chebyshevP->tform[3] = 1.1;
637:   chebyshevP->eststeps = 10;
638:   chebyshevP->usenoisy = PETSC_TRUE;
639:   ksp->setupnewmatrix  = PETSC_TRUE;

641:   ksp->ops->setup          = KSPSetUp_Chebyshev;
642:   ksp->ops->solve          = KSPSolve_Chebyshev;
643:   ksp->ops->destroy        = KSPDestroy_Chebyshev;
644:   ksp->ops->buildsolution  = KSPBuildSolutionDefault;
645:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
646:   ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
647:   ksp->ops->view           = KSPView_Chebyshev;
648:   ksp->ops->reset          = KSPReset_Chebyshev;

650:   PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevSetEigenvalues_C", KSPChebyshevSetEigenvalues_Chebyshev);
651:   PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSet_C", KSPChebyshevEstEigSet_Chebyshev);
652:   PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigSetUseNoisy_C", KSPChebyshevEstEigSetUseNoisy_Chebyshev);
653:   PetscObjectComposeFunction((PetscObject)ksp, "KSPChebyshevEstEigGetKSP_C", KSPChebyshevEstEigGetKSP_Chebyshev);
654:   return 0;
655: }