Actual source code: cheby.c
petsc-3.14.3 2021-01-09
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: PETSC_STATIC_INLINE PetscScalar chebyhash(PetscInt xx)
17: {
18: unsigned int x = xx;
19: x = ((x >> 16) ^ x) * 0x45d9f3b;
20: x = ((x >> 16) ^ x) * 0x45d9f3b;
21: x = ((x >> 16) ^ x);
22: return (PetscScalar)((PetscInt64)x-2147483648)*5.e-10; /* center around zero, scaled about -1. to 1.*/
23: }
25: /*
26: * Must be passed a KSP solver that has "converged", with KSPSetComputeEigenvalues() called before the solve
27: */
28: static PetscErrorCode KSPChebyshevComputeExtremeEigenvalues_Private(KSP kspest,PetscReal *emin,PetscReal *emax)
29: {
31: PetscInt n,neig;
32: PetscReal *re,*im,min,max;
35: KSPGetIterationNumber(kspest,&n);
36: PetscMalloc2(n,&re,n,&im);
37: KSPComputeEigenvalues(kspest,n,re,im,&neig);
38: min = PETSC_MAX_REAL;
39: max = PETSC_MIN_REAL;
40: for (n=0; n<neig; n++) {
41: min = PetscMin(min,re[n]);
42: max = PetscMax(max,re[n]);
43: }
44: PetscFree2(re,im);
45: *emax = max;
46: *emin = min;
47: return(0);
48: }
50: static PetscErrorCode KSPSetUp_Chebyshev(KSP ksp)
51: {
52: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
53: PetscErrorCode ierr;
54: PetscBool flg;
55: Mat Pmat,Amat;
56: PetscObjectId amatid, pmatid;
57: PetscObjectState amatstate, pmatstate;
59: KSPSetWorkVecs(ksp,3);
60: if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) { /* We need to estimate eigenvalues */
61: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
62: }
63: if (cheb->kspest) {
64: KSPGetOperators(ksp,&Amat,&Pmat);
65: MatGetOption(Pmat, MAT_SPD, &flg);
66: if (flg) {
67: const char *prefix;
68: KSPGetOptionsPrefix(cheb->kspest,&prefix);
69: PetscOptionsHasName(NULL,prefix,"-ksp_type",&flg);
70: if (!flg) {
71: KSPSetType(cheb->kspest, KSPCG);
72: }
73: }
74: PetscObjectGetId((PetscObject)Amat,&amatid);
75: PetscObjectGetId((PetscObject)Pmat,&pmatid);
76: PetscObjectStateGet((PetscObject)Amat,&amatstate);
77: PetscObjectStateGet((PetscObject)Pmat,&pmatstate);
78: if (amatid != cheb->amatid || pmatid != cheb->pmatid || amatstate != cheb->amatstate || pmatstate != cheb->pmatstate) {
79: PetscReal max=0.0,min=0.0;
80: Vec B;
81: KSPConvergedReason reason;
82: KSPSetPC(cheb->kspest,ksp->pc);
83: if (cheb->usenoisy) {
84: PetscInt n,i,istart;
85: PetscScalar *xx;
87: B = ksp->work[1];
88: VecGetOwnershipRange(B,&istart,NULL);
89: VecGetLocalSize(B,&n);
90: VecGetArrayWrite(B,&xx);
91: for (i=0; i<n; i++) xx[i] = chebyhash(i+istart);
92: VecRestoreArrayWrite(B,&xx);
93: } else {
94: PetscBool change;
96: 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");
97: PCPreSolveChangeRHS(ksp->pc,&change);
98: if (change) {
99: B = ksp->work[1];
100: VecCopy(ksp->vec_rhs,B);
101: } else B = ksp->vec_rhs;
102: }
103: KSPSolve(cheb->kspest,B,ksp->work[0]);
104: KSPGetConvergedReason(cheb->kspest,&reason);
105: if (reason == KSP_DIVERGED_ITS) {
106: PetscInfo(ksp,"Eigen estimator ran for prescribed number of iterations\n");
107: } else if (reason == KSP_DIVERGED_PC_FAILED) {
108: PetscInt its;
109: PCFailedReason pcreason;
111: KSPGetIterationNumber(cheb->kspest,&its);
112: if (ksp->normtype == KSP_NORM_NONE) {
113: PetscInt sendbuf,recvbuf;
114: PCGetFailedReasonRank(ksp->pc,&pcreason);
115: sendbuf = (PetscInt)pcreason;
116: MPI_Allreduce(&sendbuf,&recvbuf,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)ksp));
117: PCSetFailedReason(ksp->pc,(PCFailedReason) recvbuf);
118: }
119: PCGetFailedReason(ksp->pc,&pcreason);
120: ksp->reason = KSP_DIVERGED_PC_FAILED;
121: PetscInfo3(ksp,"Eigen estimator failed: %s %s at iteration %D",KSPConvergedReasons[reason],PCFailedReasons[pcreason],its);
122: return(0);
123: } else if (reason == KSP_CONVERGED_RTOL || reason == KSP_CONVERGED_ATOL) {
124: PetscInfo(ksp,"Eigen estimator converged prematurely. Should not happen except for small or low rank problem\n");
125: } else if (reason < 0) {
126: PetscInfo1(ksp,"Eigen estimator failed %s, using estimates anyway\n",KSPConvergedReasons[reason]);
127: }
129: KSPChebyshevComputeExtremeEigenvalues_Private(cheb->kspest,&min,&max);
130: KSPSetPC(cheb->kspest,NULL);
132: cheb->emin_computed = min;
133: cheb->emax_computed = max;
134: cheb->emin = cheb->tform[0]*min + cheb->tform[1]*max;
135: cheb->emax = cheb->tform[2]*min + cheb->tform[3]*max;
137: cheb->amatid = amatid;
138: cheb->pmatid = pmatid;
139: cheb->amatstate = amatstate;
140: cheb->pmatstate = pmatstate;
141: }
142: }
143: return(0);
144: }
146: static PetscErrorCode KSPChebyshevSetEigenvalues_Chebyshev(KSP ksp,PetscReal emax,PetscReal emin)
147: {
148: KSP_Chebyshev *chebyshevP = (KSP_Chebyshev*)ksp->data;
152: 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);
153: 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);
154: chebyshevP->emax = emax;
155: chebyshevP->emin = emin;
157: KSPChebyshevEstEigSet(ksp,0.,0.,0.,0.); /* Destroy any estimation setup */
158: return(0);
159: }
161: static PetscErrorCode KSPChebyshevEstEigSet_Chebyshev(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
162: {
163: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
167: if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
168: if (!cheb->kspest) { /* should this block of code be moved to KSPSetUp_Chebyshev()? */
169: KSPCreate(PetscObjectComm((PetscObject)ksp),&cheb->kspest);
170: KSPSetErrorIfNotConverged(cheb->kspest,ksp->errorifnotconverged);
171: PetscObjectIncrementTabLevel((PetscObject)cheb->kspest,(PetscObject)ksp,1);
172: /* use PetscObjectSet/AppendOptionsPrefix() instead of KSPSet/AppendOptionsPrefix() so that the PC prefix is not changed */
173: PetscObjectSetOptionsPrefix((PetscObject)cheb->kspest,((PetscObject)ksp)->prefix);
174: PetscObjectAppendOptionsPrefix((PetscObject)cheb->kspest,"esteig_");
175: KSPSetSkipPCSetFromOptions(cheb->kspest,PETSC_TRUE);
177: KSPSetComputeEigenvalues(cheb->kspest,PETSC_TRUE);
179: /* We cannot turn off convergence testing because GMRES will break down if you attempt to keep iterating after a zero norm is obtained */
180: KSPSetTolerances(cheb->kspest,1.e-12,PETSC_DEFAULT,PETSC_DEFAULT,cheb->eststeps);
181: }
182: if (a >= 0) cheb->tform[0] = a;
183: if (b >= 0) cheb->tform[1] = b;
184: if (c >= 0) cheb->tform[2] = c;
185: if (d >= 0) cheb->tform[3] = d;
186: cheb->amatid = 0;
187: cheb->pmatid = 0;
188: cheb->amatstate = -1;
189: cheb->pmatstate = -1;
190: } else {
191: KSPDestroy(&cheb->kspest);
192: }
193: return(0);
194: }
196: static PetscErrorCode KSPChebyshevEstEigSetUseNoisy_Chebyshev(KSP ksp,PetscBool use)
197: {
198: 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 on ksp
211: Input Parameters:
212: + ksp - the Krylov space context
213: - emax, emin - the eigenvalue estimates
215: Options Database:
216: . -ksp_chebyshev_eigenvalues emin,emax
218: Note: Call KSPChebyshevEstEigSet() or use the option -ksp_chebyshev_esteig a,b,c,d to have the KSP
219: estimate the eigenvalues and use these estimated values automatically
221: Level: intermediate
223: @*/
224: PetscErrorCode KSPChebyshevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
225: {
232: PetscTryMethod(ksp,"KSPChebyshevSetEigenvalues_C",(KSP,PetscReal,PetscReal),(ksp,emax,emin));
233: return(0);
234: }
236: /*@
237: KSPChebyshevEstEigSet - Automatically estimate the eigenvalues to use for Chebyshev
239: Logically Collective on ksp
241: Input Parameters:
242: + ksp - the Krylov space context
243: . a - multiple of min eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
244: . b - multiple of max eigenvalue estimate to use for min Chebyshev bound (or PETSC_DECIDE)
245: . c - multiple of min eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
246: - d - multiple of max eigenvalue estimate to use for max Chebyshev bound (or PETSC_DECIDE)
248: Options Database:
249: . -ksp_chebyshev_esteig a,b,c,d
251: Notes:
252: The Chebyshev bounds are set using
253: .vb
254: minbound = a*minest + b*maxest
255: maxbound = c*minest + d*maxest
256: .ve
257: The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
258: The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.
260: If 0.0 is passed for all transform arguments (a,b,c,d), eigenvalue estimation is disabled.
262: The default transform is (0,0.1; 0,1.1) which targets the "upper" part of the spectrum, as desirable for use with multigrid.
264: The eigenvalues are estimated using the Lanczo (KSPCG) or Arnoldi (KSPGMRES) process using a noisy right hand side vector.
266: Level: intermediate
268: @*/
269: PetscErrorCode KSPChebyshevEstEigSet(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
270: {
279: PetscTryMethod(ksp,"KSPChebyshevEstEigSet_C",(KSP,PetscReal,PetscReal,PetscReal,PetscReal),(ksp,a,b,c,d));
280: return(0);
281: }
283: /*@
284: KSPChebyshevEstEigSetUseNoisy - use a noisy right hand side in order to do the estimate instead of the given right hand side
286: Logically Collective
288: Input Arguments:
289: + ksp - linear solver context
290: - use - PETSC_TRUE to use noisy
292: Options Database:
293: . -ksp_chebyshev_esteig_noisy <true,false>
295: Notes:
296: This alledgely works better for multigrid smoothers
298: Level: intermediate
300: .seealso: KSPChebyshevEstEigSet()
301: @*/
302: PetscErrorCode KSPChebyshevEstEigSetUseNoisy(KSP ksp,PetscBool use)
303: {
308: PetscTryMethod(ksp,"KSPChebyshevEstEigSetUseNoisy_C",(KSP,PetscBool),(ksp,use));
309: return(0);
310: }
312: /*@
313: KSPChebyshevEstEigGetKSP - Get the Krylov method context used to estimate eigenvalues for the Chebyshev method. If
314: a Krylov method is not being used for this purpose, NULL is returned. The reference count of the returned KSP is
315: not incremented: it should not be destroyed by the user.
317: Input Parameters:
318: . ksp - the Krylov space context
320: Output Parameters:
321: . kspest the eigenvalue estimation Krylov space context
323: Level: intermediate
325: .seealso: KSPChebyshevEstEigSet()
326: @*/
327: PetscErrorCode KSPChebyshevEstEigGetKSP(KSP ksp, KSP *kspest)
328: {
334: *kspest = NULL;
335: PetscTryMethod(ksp,"KSPChebyshevEstEigGetKSP_C",(KSP,KSP*),(ksp,kspest));
336: return(0);
337: }
339: static PetscErrorCode KSPChebyshevEstEigGetKSP_Chebyshev(KSP ksp, KSP *kspest)
340: {
341: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
344: *kspest = cheb->kspest;
345: return(0);
346: }
348: static PetscErrorCode KSPSetFromOptions_Chebyshev(PetscOptionItems *PetscOptionsObject,KSP ksp)
349: {
350: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
352: PetscInt neigarg = 2, nestarg = 4;
353: PetscReal eminmax[2] = {0., 0.};
354: PetscReal tform[4] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
355: PetscBool flgeig, flgest;
358: PetscOptionsHead(PetscOptionsObject,"KSP Chebyshev Options");
359: PetscOptionsInt("-ksp_chebyshev_esteig_steps","Number of est steps in Chebyshev","",cheb->eststeps,&cheb->eststeps,NULL);
360: PetscOptionsRealArray("-ksp_chebyshev_eigenvalues","extreme eigenvalues","KSPChebyshevSetEigenvalues",eminmax,&neigarg,&flgeig);
361: if (flgeig) {
362: if (neigarg != 2) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"-ksp_chebyshev_eigenvalues: must specify 2 parameters, min and max eigenvalues");
363: KSPChebyshevSetEigenvalues(ksp, eminmax[1], eminmax[0]);
364: }
365: PetscOptionsRealArray("-ksp_chebyshev_esteig","estimate eigenvalues using a Krylov method, then use this transform for Chebyshev eigenvalue bounds","KSPChebyshevEstEigSet",tform,&nestarg,&flgest);
366: if (flgest) {
367: switch (nestarg) {
368: case 0:
369: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
370: break;
371: case 2: /* Base everything on the max eigenvalues */
372: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,tform[0],PETSC_DECIDE,tform[1]);
373: break;
374: case 4: /* Use the full 2x2 linear transformation */
375: KSPChebyshevEstEigSet(ksp,tform[0],tform[1],tform[2],tform[3]);
376: break;
377: default: SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"Must specify either 0, 2, or 4 parameters for eigenvalue estimation");
378: }
379: }
381: /* We need to estimate eigenvalues; need to set this here so that KSPSetFromOptions() is called on the estimator */
382: if ((cheb->emin == 0. || cheb->emax == 0.) && !cheb->kspest) {
383: KSPChebyshevEstEigSet(ksp,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE);
384: }
386: if (cheb->kspest) {
387: PetscOptionsBool("-ksp_chebyshev_esteig_noisy","Use noisy right hand side for estimate","KSPChebyshevEstEigSetUseNoisy",cheb->usenoisy,&cheb->usenoisy,NULL);
388: KSPSetFromOptions(cheb->kspest);
389: }
390: PetscOptionsTail();
391: return(0);
392: }
394: static PetscErrorCode KSPSolve_Chebyshev(KSP ksp)
395: {
396: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
398: PetscInt k,kp1,km1,ktmp,i;
399: PetscScalar alpha,omegaprod,mu,omega,Gamma,c[3],scale;
400: PetscReal rnorm = 0.0;
401: Vec sol_orig,b,p[3],r;
402: Mat Amat,Pmat;
403: PetscBool diagonalscale;
406: PCGetDiagonalScale(ksp->pc,&diagonalscale);
407: if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
409: PCGetOperators(ksp->pc,&Amat,&Pmat);
410: PetscObjectSAWsTakeAccess((PetscObject)ksp);
411: ksp->its = 0;
412: PetscObjectSAWsGrantAccess((PetscObject)ksp);
413: /* These three point to the three active solutions, we
414: rotate these three at each solution update */
415: km1 = 0; k = 1; kp1 = 2;
416: sol_orig = ksp->vec_sol; /* ksp->vec_sol will be asigned to rotating vector p[k], thus save its address */
417: b = ksp->vec_rhs;
418: p[km1] = sol_orig;
419: p[k] = ksp->work[0];
420: p[kp1] = ksp->work[1];
421: r = ksp->work[2];
423: /* use scale*B as our preconditioner */
424: scale = 2.0/(cheb->emax + cheb->emin);
426: /* -alpha <= scale*lambda(B^{-1}A) <= alpha */
427: alpha = 1.0 - scale*(cheb->emin);
428: Gamma = 1.0;
429: mu = 1.0/alpha;
430: omegaprod = 2.0/alpha;
432: c[km1] = 1.0;
433: c[k] = mu;
435: if (!ksp->guess_zero) {
436: KSP_MatMult(ksp,Amat,sol_orig,r); /* r = b - A*p[km1] */
437: VecAYPX(r,-1.0,b);
438: } else {
439: VecCopy(b,r);
440: }
442: /* calculate residual norm if requested, we have done one iteration */
443: if (ksp->normtype) {
444: switch (ksp->normtype) {
445: case KSP_NORM_PRECONDITIONED:
446: KSP_PCApply(ksp,r,p[k]); /* p[k] = B^{-1}r */
447: VecNorm(p[k],NORM_2,&rnorm);
448: break;
449: case KSP_NORM_UNPRECONDITIONED:
450: case KSP_NORM_NATURAL:
451: VecNorm(r,NORM_2,&rnorm);
452: break;
453: default: SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"%s",KSPNormTypes[ksp->normtype]);
454: }
455: PetscObjectSAWsTakeAccess((PetscObject)ksp);
456: ksp->rnorm = rnorm;
457: PetscObjectSAWsGrantAccess((PetscObject)ksp);
458: KSPLogResidualHistory(ksp,rnorm);
459: KSPMonitor(ksp,0,rnorm);
460: (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
461: } else ksp->reason = KSP_CONVERGED_ITERATING;
462: if (ksp->reason || ksp->max_it==0) {
463: if (ksp->max_it==0) ksp->reason = KSP_DIVERGED_ITS; /* This for a V(0,x) cycle */
464: return(0);
465: }
466: if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
467: KSP_PCApply(ksp,r,p[k]); /* p[k] = B^{-1}r */
468: }
469: VecAYPX(p[k],scale,p[km1]); /* p[k] = scale B^{-1}r + p[km1] */
470: PetscObjectSAWsTakeAccess((PetscObject)ksp);
471: ksp->its = 1;
472: PetscObjectSAWsGrantAccess((PetscObject)ksp);
474: for (i=1; i<ksp->max_it; i++) {
475: PetscObjectSAWsTakeAccess((PetscObject)ksp);
476: ksp->its++;
477: PetscObjectSAWsGrantAccess((PetscObject)ksp);
479: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
480: VecAYPX(r,-1.0,b);
481: /* calculate residual norm if requested */
482: if (ksp->normtype) {
483: switch (ksp->normtype) {
484: case KSP_NORM_PRECONDITIONED:
485: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
486: VecNorm(p[kp1],NORM_2,&rnorm);
487: break;
488: case KSP_NORM_UNPRECONDITIONED:
489: case KSP_NORM_NATURAL:
490: VecNorm(r,NORM_2,&rnorm);
491: break;
492: default:
493: rnorm = 0.0;
494: break;
495: }
496: KSPCheckNorm(ksp,rnorm);
497: PetscObjectSAWsTakeAccess((PetscObject)ksp);
498: ksp->rnorm = rnorm;
499: PetscObjectSAWsGrantAccess((PetscObject)ksp);
500: KSPLogResidualHistory(ksp,rnorm);
501: KSPMonitor(ksp,i,rnorm);
502: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
503: if (ksp->reason) break;
504: if (ksp->normtype != KSP_NORM_PRECONDITIONED) {
505: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
506: }
507: } else {
508: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
509: }
510: ksp->vec_sol = p[k];
512: c[kp1] = 2.0*mu*c[k] - c[km1];
513: omega = omegaprod*c[k]/c[kp1];
515: /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
516: VecAXPBYPCZ(p[kp1],1.0-omega,omega,omega*Gamma*scale,p[km1],p[k]);
518: ktmp = km1;
519: km1 = k;
520: k = kp1;
521: kp1 = ktmp;
522: }
523: if (!ksp->reason) {
524: if (ksp->normtype) {
525: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
526: VecAYPX(r,-1.0,b);
527: switch (ksp->normtype) {
528: case KSP_NORM_PRECONDITIONED:
529: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}r */
530: VecNorm(p[kp1],NORM_2,&rnorm);
531: break;
532: case KSP_NORM_UNPRECONDITIONED:
533: case KSP_NORM_NATURAL:
534: VecNorm(r,NORM_2,&rnorm);
535: break;
536: default:
537: rnorm = 0.0;
538: break;
539: }
540: KSPCheckNorm(ksp,rnorm);
541: PetscObjectSAWsTakeAccess((PetscObject)ksp);
542: ksp->rnorm = rnorm;
543: PetscObjectSAWsGrantAccess((PetscObject)ksp);
544: KSPLogResidualHistory(ksp,rnorm);
545: KSPMonitor(ksp,i,rnorm);
546: }
547: if (ksp->its >= ksp->max_it) {
548: if (ksp->normtype != KSP_NORM_NONE) {
549: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
550: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
551: } else ksp->reason = KSP_CONVERGED_ITS;
552: }
553: }
555: /* make sure solution is in vector x */
556: ksp->vec_sol = sol_orig;
557: if (k) {
558: VecCopy(p[k],sol_orig);
559: }
560: return(0);
561: }
563: static PetscErrorCode KSPView_Chebyshev(KSP ksp,PetscViewer viewer)
564: {
565: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
567: PetscBool iascii;
570: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
571: if (iascii) {
572: PetscViewerASCIIPrintf(viewer," eigenvalue estimates used: min = %g, max = %g\n",(double)cheb->emin,(double)cheb->emax);
573: if (cheb->kspest) {
574: PetscViewerASCIIPrintf(viewer," eigenvalues estimate via %s min %g, max %g\n",((PetscObject)(cheb->kspest))->type_name,(double)cheb->emin_computed,(double)cheb->emax_computed);
575: 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]);
576: PetscViewerASCIIPushTab(viewer);
577: KSPView(cheb->kspest,viewer);
578: PetscViewerASCIIPopTab(viewer);
579: if (cheb->usenoisy) {
580: PetscViewerASCIIPrintf(viewer," estimating eigenvalues using noisy right hand side\n");
581: }
582: }
583: }
584: return(0);
585: }
587: static PetscErrorCode KSPDestroy_Chebyshev(KSP ksp)
588: {
589: KSP_Chebyshev *cheb = (KSP_Chebyshev*)ksp->data;
593: KSPDestroy(&cheb->kspest);
594: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",NULL);
595: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",NULL);
596: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",NULL);
597: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",NULL);
598: KSPDestroyDefault(ksp);
599: return(0);
600: }
602: /*MC
603: KSPCHEBYSHEV - The preconditioned Chebyshev iterative method
605: Options Database Keys:
606: + -ksp_chebyshev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
607: of the preconditioned operator. If these are accurate you will get much faster convergence.
608: . -ksp_chebyshev_esteig <a,b,c,d> - estimate eigenvalues using a Krylov method, then use this
609: transform for Chebyshev eigenvalue bounds (KSPChebyshevEstEigSet())
610: . -ksp_chebyshev_esteig_steps - number of estimation steps
611: - -ksp_chebyshev_esteig_noisy - use noisy number generator to create right hand side for eigenvalue estimator
613: Level: beginner
615: Notes:
616: The Chebyshev method requires both the matrix and preconditioner to
617: be symmetric positive (semi) definite.
618: Only support for left preconditioning.
620: Chebyshev is configured as a smoother by default, targetting the "upper" part of the spectrum.
621: The user should call KSPChebyshevSetEigenvalues() if they have eigenvalue estimates.
623: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
624: KSPChebyshevSetEigenvalues(), KSPChebyshevEstEigSet(), KSPChebyshevEstEigSetUseNoisy()
625: KSPRICHARDSON, KSPCG, PCMG
627: M*/
629: PETSC_EXTERN PetscErrorCode KSPCreate_Chebyshev(KSP ksp)
630: {
632: KSP_Chebyshev *chebyshevP;
635: PetscNewLog(ksp,&chebyshevP);
637: ksp->data = (void*)chebyshevP;
638: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
639: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,2);
640: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
641: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
643: chebyshevP->emin = 0.;
644: chebyshevP->emax = 0.;
646: chebyshevP->tform[0] = 0.0;
647: chebyshevP->tform[1] = 0.1;
648: chebyshevP->tform[2] = 0;
649: chebyshevP->tform[3] = 1.1;
650: chebyshevP->eststeps = 10;
651: chebyshevP->usenoisy = PETSC_TRUE;
652: ksp->setupnewmatrix = PETSC_TRUE;
654: ksp->ops->setup = KSPSetUp_Chebyshev;
655: ksp->ops->solve = KSPSolve_Chebyshev;
656: ksp->ops->destroy = KSPDestroy_Chebyshev;
657: ksp->ops->buildsolution = KSPBuildSolutionDefault;
658: ksp->ops->buildresidual = KSPBuildResidualDefault;
659: ksp->ops->setfromoptions = KSPSetFromOptions_Chebyshev;
660: ksp->ops->view = KSPView_Chebyshev;
661: ksp->ops->reset = KSPReset_Chebyshev;
663: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevSetEigenvalues_C",KSPChebyshevSetEigenvalues_Chebyshev);
664: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSet_C",KSPChebyshevEstEigSet_Chebyshev);
665: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigSetUseNoisy_C",KSPChebyshevEstEigSetUseNoisy_Chebyshev);
666: PetscObjectComposeFunction((PetscObject)ksp,"KSPChebyshevEstEigGetKSP_C",KSPChebyshevEstEigGetKSP_Chebyshev);
667: return(0);
668: }