Actual source code: cheby.c
petsc-3.15.0 2021-04-05
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;
10: if (cheb->kspest) {
11: KSPReset(cheb->kspest);
12: }
13: return(0);
14: }
16: /*
17: * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
18: */
19: static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest,PetscReal *emin,PetscReal *emax)
20: {
22: PetscInt n,neig;
23: PetscReal *re,*im,min,max;
26: KSPGetIterationNumber(kspest,&n);
27: PetscMalloc2(n,&re,n,&im);
28: KSPComputeEigenvalues(kspest,n,re,im,&neig);
29: min = PETSC_MAX_REAL;
30: max = PETSC_MIN_REAL;
31: for (n=0; n<neig; n++) {
32: min = PetscMin(min,re[n]);
33: max = PetscMax(max,re[n]);
34: }
35: PetscFree2(re,im);
36: *emax = max;
37: *emin = min;
38: return(0);
39: }
41: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
42: {
43: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
44: PetscErrorCode ierr;
45: PetscBool flg;
46: Mat Pmat,Amat;
47: PetscObjectId amatid, pmatid;
48: PetscObjectState amatstate, pmatstate;
50: KSPSetWorkVecs(ksp,3);
51: if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) { /* We need to estimate eigenvalues */
52: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
53: }
54: if (cheb->kspest) {
55: KSPGetOperators(ksp,&Amat,&Pmat);
56: MatGetOption(Pmat, MAT_SPD, &flg);
57: if (flg) {
58: const char *prefix;
59: KSPGetOptionsPrefix(cheb->kspest,&prefix);
60: PetscOptionsHasName(NULL,prefix,"-ksp_type",&flg);
61: if (!flg) {
62: KSPSetType(cheb->kspest, KSPCG);
63: }
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;
80: if (!ksp->vec_rhs) SETERRQ(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");
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: PetscInfo3(ksp,"Eigen estimator failed: %s %s at iteration %D",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: PetscInfo1(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;
118: cheb->emin = cheb->tform[0]*min + cheb->tform[1]*max;
119: cheb->emax = cheb->tform[2]*min + cheb->tform[3]*max;
121: cheb->amatid = amatid;
122: cheb->pmatid = pmatid;
123: cheb->amatstate = amatstate;
124: cheb->pmatstate = pmatstate;
125: }
126: }
127: return(0);
128: }
130: static PetscErrorCode KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp,PetscReal emax,PetscReal emin)
131: {
132: KSP_Chebyshev *chebyshevP = (KSP_Chebyshev*)ksp->data;
136: if (emax <= emin) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Maximum eigenvalue must be larger than minimum: max %g min %g",(double)emax,(double)emin);
137: if (emax*emin <= 0.0) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Both eigenvalues must be of the same sign: max %g min %g",(double)emax,(double)emin);
138: chebyshevP->emax = emax;
139: chebyshevP->emin = emin;
141: KSPChebyshevEstEigSet(ksp,0.,0.,0.,0.); /* Destroy any estimation setup */
142: return(0);
143: }
145: static PetscErrorCode KSPChebyshevEstEigSet_Chebyshev(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
146: {
147: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
151: if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
152: if (!cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
153: KSPCreate(PetscObjectComm((PetscObject)ksp),&cheb->kspest);
154: KSPSetErrorIfNotConverged(cheb->kspest,ksp->errorifnotconverged);
155: PetscObjectIncrementTabLevel((PetscObject)cheb->kspest,(PetscObject)ksp,1);
156: /* use PetscObjectSet/AppendOptionsPrefix() instead of KSPSet/AppendOptionsPrefix() so that the PC prefix is not changed */
157: PetscObjectSetOptionsPrefix((PetscObject)cheb->kspest,((PetscObject)ksp)->prefix);
158: PetscObjectAppendOptionsPrefix((PetscObject)cheb->kspest,"esteig_");
159: KSPSetSkipPCSetFromOptions(cheb->kspest,PETSC_TRUE);
161: KSPSetComputeEigenvalues(cheb->kspest,PETSC_TRUE);
163: /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
164: KSPSetTolerances(cheb->kspest,1.e-12,PETSC_DEFAULT,PETSC_DEFAULT,cheb->eststeps);
165: }
166: if (a >= 0) cheb->tform[0] = a;
167: if (b >= 0) cheb->tform[1] = b;
168: if (c >= 0) cheb->tform[2] = c;
169: if (d >= 0) cheb->tform[3] = d;
170: cheb->amatid = 0;
171: cheb->pmatid = 0;
172: cheb->amatstate = -1;
173: cheb->pmatstate = -1;
174: } else {
175: KSPDestroy(&cheb->kspest);
176: }
177: return(0);
178: }
180: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp,PetscBool use)
181: {
182: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
185: cheb->usenoisy = use;
186: return(0);
187: }
189: /*@
190: KSPChebyshevSetEigenvalues - Sets estimates for the extreme eigenvalues
191: of the preconditioned problem.
193: Logically Collective on ksp
195: Input Parameters:
196: + ksp - the Krylov space context
197: - emax, emin - the eigenvalue estimates
199: Options Database:
200: . -ksp_chebyshev_eigenvalues emin,emax
202: Note: Call KSPChebyshevEstEigSet() or use the option -ksp_chebyshev_esteig a,b,c,d to have the KSP
203: estimate the eigenvalues and use these estimated values automatically
205: Level: intermediate
207: @*/
208: PetscErrorCode KSPChebyshevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
209: {
216: PetscTryMethod(ksp,"KSPChebyshevSetEigenvalues_C",(KSP,PetscReal,PetscReal),(ksp,emax,emin));
217: return(0);
218: }
220: /*@
221: KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev
223: Logically Collective on ksp
225: Input Parameters:
226: + ksp - the Krylov space context
227: . a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
228: . b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
229: . c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
230: - d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
232: Options Database:
233: . -ksp_chebyshev_esteig a,b,c,d
235: Notes:
236: The Chebyshev bounds are set using
237: .vb
238: minbound = a*minest + b*maxest
239: maxbound = c*minest + d*maxest
240: .ve
241: The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
242: The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.
244: If 0.0 is passed for all transform arguments (a,b,c,d), eigenvalue estimation is disabled.
246: The default transform is (0,0.1; 0,1.1) which targets the "upper" part of the spectrum, as desirable for use with multigrid.
248: The eigenvalues are estimated using the Lanczo (KSPCG) or Arnoldi (KSPGMRES) process using a noisy right hand side vector.
250: Level: intermediate
252: @*/
253: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
254: {
263: PetscTryMethod(ksp,"KSPChebyshevEstEigSet_C",(KSP,PetscReal,PetscReal,PetscReal,PetscReal),(ksp,a,b,c,d));
264: return(0);
265: }
267: /*@
268: KSPChebyshevEstEigSetUseNoisy - use a noisy right hand side in order to do the estimate instead of the given right hand side
270: Logically Collective
272: Input Arguments:
273: + ksp - linear solver context
274: - use - PETSC_TRUE to use noisy
276: Options Database:
277: . -ksp_chebyshev_esteig_noisy <true,false>
279: Notes:
280: This alledgely works better for multigrid smoothers
282: Level: intermediate
284: .seealso: KSPChebyshevEstEigSet()
285: @*/
286: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp,PetscBool use)
287: {
292: PetscTryMethod(ksp,"KSPChebyshevEstEigSetUseNoisy_C",(KSP,PetscBool),(ksp,use));
293: return(0);
294: }
296: /*@
297: KSPChebyshevEstEigGetKSP - Get the Krylov method context used to estimate eigenvalues for the Chebyshev method. If
298: a Krylov method is not being used for this purpose, NULL is returned. The reference count of the returned KSP is
299: not incremented: it should not be destroyed by the user.
301: Input Parameters:
302: . ksp - the Krylov space context
304: Output Parameters:
305: . kspest the eigenvalue estimation Krylov space context
307: Level: intermediate
309: .seealso: KSPChebyshevEstEigSet()
310: @*/
311: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
312: {
318: *kspest = NULL;
319: PetscTryMethod(ksp,"KSPChebyshevEstEigGetKSP_C",(KSP,KSP*),(ksp,kspest));
320: return(0);
321: }
323: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
324: {
325: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
328: *kspest = cheb->kspest;
329: return(0);
330: }
332: static PetscErrorCode KSPSetFromOptions_Chebyshev(PetscOptionItems *PetscOptionsObject,KSP ksp)
333: {
334: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
336: PetscInt neigarg = 2, nestarg = 4;
337: PetscReal eminmax[2] = {0., 0.};
338: PetscReal tform[4] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
339: PetscBool flgeig, flgest;
342: PetscOptionsHead(PetscOptionsObject,"KSP Chebyshev Options");
343: PetscOptionsInt("-ksp_chebyshev_esteig_steps","Number of est steps in Chebyshev","",cheb->eststeps,&cheb->eststeps,NULL);
344: PetscOptionsRealArray("-ksp_chebyshev_eigenvalues","extreme eigenvalues","KSPChebyshevSetEigenvalues",eminmax,&neigarg,&flgeig);
345: if (flgeig) {
346: if (neigarg != 2) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
347: KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]);
348: }
349: PetscOptionsRealArray("-ksp_chebyshev_esteig","estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds","KSPChebyshevEstEigSet",tform,&nestarg,&flgest);
350: if (flgest) {
351: switch (nestarg) {
352: case 0:
353: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
354: break;
355: case 2: /* Base everything on the max eigenvalues */
356: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,tform[0],PETSC_DECIDE,tform[1]);
357: break;
358: case 4: /* Use the full 2x2 linear transformation */
359: KSPChebyshevEstEigSet(ksp,tform[0],tform[1],tform[2],tform[3]);
360: break;
361: default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
362: }
363: }
365: /* We need to estimate eigenvalues; need to set this here so that KSPSetFromOptions() is called on the estimator */
366: if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) {
367: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
368: }
370: if (cheb->kspest) {
371: PetscOptionsBool("-ksp_chebyshev_esteig_noisy","Use noisy right hand side for estimate","KSPChebyshevEstEigSetUseNoisy",cheb->usenoisy,&cheb->usenoisy,NULL);
372: KSPSetFromOptions(cheb->kspest);
373: }
374: PetscOptionsTail();
375: return(0);
376: }
378: static PetscErrorCode KSPSolve_Chebyshev(KSP ksp)
379: {
380: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
382: PetscInt k,kp1,km1,ktmp,i;
383: PetscScalar alpha,omegaprod,mu,omega,Gamma,c[3],scale;
384: PetscReal rnorm = 0.0;
385: Vec sol_orig,b,p[3],r;
386: Mat Amat,Pmat;
387: PetscBool diagonalscale;
390: PCGetDiagonalScale(ksp->pc,&diagonalscale);
391: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
393: PCGetOperators(ksp->pc,&Amat,&Pmat);
394: PetscObjectSAWsTakeAccess((PetscObject)ksp);
395: ksp->its = 0;
396: PetscObjectSAWsGrantAccess((PetscObject)ksp);
397: /* These three point to the three active solutions, we
398: rotate these three at each solution update */
399: km1 = 0; k = 1; kp1 = 2;
400: sol_orig = ksp->vec_sol; /* ksp->vec_sol will be assigned to rotating vector p[k], thus save its address */
401: b = ksp->vec_rhs;
402: p[km1] = sol_orig;
403: p[k] = ksp->work[0];
404: p[kp1] = ksp->work[1];
405: r = ksp->work[2];
407: /* use scale*B as our preconditioner */
408: scale = 2.0/(cheb->emax + cheb->emin);
410: /* -alpha <= scale*lambda(B^{-1}A) <= alpha */
411: alpha = 1.0 - scale*(cheb->emin);
412: Gamma = 1.0;
413: mu = 1.0/alpha;
414: omegaprod = 2.0/alpha;
416: c[km1] = 1.0;
417: c[k] = mu;
419: if (!ksp->guess_zero) {
420: KSP_MatMult(ksp,Amat,sol_orig,r); /* r = b - A*p[km1] */
421: VecAYPX(r,-1.0,b);
422: } else {
423: VecCopy(b,r);
424: }
426: /* calculate residual norm if requested, we have done one iteration */
427: if (ksp->normtype) {
428: switch (ksp->normtype) {
429: case KSP_NORM_PRECONDITIONED:
430: KSP_PCApply(ksp,r,p[k]); /* p[k] = B^{-1}r */
431: VecNorm(p[k],NORM_2,&rnorm);
432: break;
433: case KSP_NORM_UNPRECONDITIONED:
434: case KSP_NORM_NATURAL:
435: VecNorm(r,NORM_2,&rnorm);
436: break;
437: default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
438: }
439: PetscObjectSAWsTakeAccess((PetscObject)ksp);
440: ksp->rnorm = rnorm;
441: PetscObjectSAWsGrantAccess((PetscObject)ksp);
442: KSPLogResidualHistory(ksp,rnorm);
443: KSPLogErrorHistory(ksp);
444: KSPMonitor(ksp,0,rnorm);
445: (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
446: } else ksp->reason = KSP_CONVERGED_ITERATING;
447: if (ksp->reason || ksp->max_it==0) {
448: if (ksp->max_it==0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
449: return(0);
450: }
451: if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
452: KSP_PCApply(ksp,r,p[k]); /* p[k] = B^{-1}r */
453: }
454: VecAYPX(p[k],scale,p[km1]); /* p[k] = scale B^{-1}r + p[km1] */
455: PetscObjectSAWsTakeAccess((PetscObject)ksp);
456: ksp->its = 1;
457: PetscObjectSAWsGrantAccess((PetscObject)ksp);
459: for (i=1; i<ksp->max_it; i++) {
460: PetscObjectSAWsTakeAccess((PetscObject)ksp);
461: ksp->its++;
462: PetscObjectSAWsGrantAccess((PetscObject)ksp);
464: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
465: VecAYPX(r,-1.0,b);
466: /* calculate residual norm if requested */
467: if (ksp->normtype) {
468: switch (ksp->normtype) {
469: case KSP_NORM_PRECONDITIONED:
470: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
471: VecNorm(p[kp1],NORM_2,&rnorm);
472: break;
473: case KSP_NORM_UNPRECONDITIONED:
474: case KSP_NORM_NATURAL:
475: VecNorm(r,NORM_2,&rnorm);
476: break;
477: default:
478: rnorm = 0.0;
479: break;
480: }
481: KSPCheckNorm(ksp,rnorm);
482: PetscObjectSAWsTakeAccess((PetscObject)ksp);
483: ksp->rnorm = rnorm;
484: PetscObjectSAWsGrantAccess((PetscObject)ksp);
485: KSPLogResidualHistory(ksp,rnorm);
486: KSPMonitor(ksp,i,rnorm);
487: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
488: if (ksp->reason) break;
489: if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
490: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
491: }
492: } else {
493: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
494: }
495: ksp->vec_sol = p[k];
496: KSPLogErrorHistory(ksp);
498: c[kp1] = 2.0*mu*c[k] - c[km1];
499: omega = omegaprod*c[k]/c[kp1];
501: /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
502: VecAXPBYPCZ(p[kp1],1.0-omega,omega,omega*Gamma*scale,p[km1],p[k]);
504: ktmp = km1;
505: km1 = k;
506: k = kp1;
507: kp1 = ktmp;
508: }
509: if (!ksp->reason) {
510: if (ksp->normtype) {
511: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
512: VecAYPX(r,-1.0,b);
513: switch (ksp->normtype) {
514: case KSP_NORM_PRECONDITIONED:
515: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
516: VecNorm(p[kp1],NORM_2,&rnorm);
517: break;
518: case KSP_NORM_UNPRECONDITIONED:
519: case KSP_NORM_NATURAL:
520: VecNorm(r,NORM_2,&rnorm);
521: break;
522: default:
523: rnorm = 0.0;
524: break;
525: }
526: KSPCheckNorm(ksp,rnorm);
527: PetscObjectSAWsTakeAccess((PetscObject)ksp);
528: ksp->rnorm = rnorm;
529: PetscObjectSAWsGrantAccess((PetscObject)ksp);
530: KSPLogResidualHistory(ksp,rnorm);
531: KSPMonitor(ksp,i,rnorm);
532: }
533: if (ksp->its >= ksp->max_it) {
534: if (ksp->normtype != KSP_NORM_NONE) {
535: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
536: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
537: } else ksp->reason = KSP_CONVERGED_ITS;
538: }
539: }
541: /* make sure solution is in vector x */
542: ksp->vec_sol = sol_orig;
543: if (k) {
544: VecCopy(p[k],sol_orig);
545: }
546: if (ksp->reason == KSP_CONVERGED_ITS) {
547: KSPLogErrorHistory(ksp);
548: }
549: return(0);
550: }
552: static PetscErrorCode KSPView_Chebyshev(KSP ksp,PetscViewer viewer)
553: {
554: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
556: PetscBool iascii;
559: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
560: if (iascii) {
561: PetscViewerASCIIPrintf(viewer," eigenvalue estimates used: min = %g, max = %g\n",(double)cheb->emin,(double)cheb->emax);
562: if (cheb->kspest) {
563: PetscViewerASCIIPrintf(viewer," eigenvalues estimate via %s min %g, max %g\n",((PetscObject)(cheb->kspest))->type_name,(double)cheb->emin_computed,(double)cheb->emax_computed);
564: PetscViewerASCIIPrintf(viewer," eigenvalues estimated using %s with translations [%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]);
565: PetscViewerASCIIPushTab(viewer);
566: KSPView(cheb->kspest,viewer);
567: PetscViewerASCIIPopTab(viewer);
568: if (cheb->usenoisy) {
569: PetscViewerASCIIPrintf(viewer," estimating eigenvalues using noisy right hand side\n");
570: }
571: }
572: }
573: return(0);
574: }
576: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
577: {
578: 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
606: be symmetric positive (semi) definite.
607: Only support for left preconditioning.
609: Chebyshev is configured as a smoother by default, targetting the "upper" part of the spectrum.
610: The user should call KSPChebyshevSetEigenvalues() if they have eigenvalue estimates.
612: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
613: KSPChebyshevSetEigenvalues(), KSPChebyshevEstEigSet(), KSPChebyshevEstEigSetUseNoisy()
614: KSPRICHARDSON, KSPCG, PCMG
616: M*/
618: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
619: {
621: KSP_Chebyshev *chebyshevP;
624: PetscNewLog(ksp,&chebyshevP);
626: ksp->data = (void*)chebyshevP;
627: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
628: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
629: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
630: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
632: chebyshevP->emin = 0.;
633: chebyshevP->emax = 0.;
635: chebyshevP->tform[0] = 0.0;
636: chebyshevP->tform[1] = 0.1;
637: chebyshevP->tform[2] = 0;
638: chebyshevP->tform[3] = 1.1;
639: chebyshevP->eststeps = 10;
640: chebyshevP->usenoisy = PETSC_TRUE;
641: ksp->setupnewmatrix = PETSC_TRUE;
643: ksp->ops->setup = KSPSetUp_Chebyshev;
644: ksp->ops->solve = KSPSolve_Chebyshev;
645: ksp->ops->destroy = KSPDestroy_Chebyshev;
646: ksp->ops->buildsolution = KSPBuildSolutionDefault;
647: ksp->ops->buildresidual = KSPBuildResidualDefault;
648: ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
649: ksp->ops->view = KSPView_Chebyshev;
650: ksp->ops->reset = KSPReset_Chebyshev;
652: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",KSPChebyshevSetEigenvalues_Chebyshev);
653: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",KSPChebyshevEstEigSet_Chebyshev);
654: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",KSPChebyshevEstEigSetUseNoisy_Chebyshev);
655: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",KSPChebyshevEstEigGetKSP_Chebyshev);
656: return(0);
657: }