Actual source code: cheby.c
2: #include <private/kspimpl.h> /*I "petscksp.h" I*/
3: #include <../src/ksp/ksp/impls/cheby/chebychevimpl.h>
7: PetscErrorCode KSPSetUp_Chebychev(KSP ksp)
8: {
9: KSP_Chebychev *cheb = (KSP_Chebychev*)ksp->data;
13: KSPDefaultGetWork(ksp,3);
14: cheb->estimate_current = PETSC_FALSE;
15: return(0);
16: }
18: EXTERN_C_BEGIN
21: PetscErrorCode KSPChebychevSetEigenvalues_Chebychev(KSP ksp,PetscReal emax,PetscReal emin)
22: {
23: KSP_Chebychev *chebychevP = (KSP_Chebychev*)ksp->data;
26: if (emax <= emin) SETERRQ2(((PetscObject)ksp)->comm,PETSC_ERR_ARG_INCOMP,"Maximum eigenvalue must be larger than minimum: max %g min %G",emax,emin);
27: if (emax*emin <= 0.0) SETERRQ2(((PetscObject)ksp)->comm,PETSC_ERR_ARG_INCOMP,"Both eigenvalues must be of the same sign: max %G min %G",emax,emin);
28: chebychevP->emax = emax;
29: chebychevP->emin = emin;
30: return(0);
31: }
32: EXTERN_C_END
34: EXTERN_C_BEGIN
37: PetscErrorCode KSPChebychevSetEstimateEigenvalues_Chebychev(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
38: {
39: KSP_Chebychev *cheb = (KSP_Chebychev*)ksp->data;
43: if (a != 0.0 || b != 0.0 || c != 0.0 || d != 0.0) {
44: if (!cheb->kspest) {
45: PetscBool nonzero;
47: KSPCreate(((PetscObject)ksp)->comm,&cheb->kspest);
49: KSPGetPC(cheb->kspest,&cheb->pcnone);
50: PetscObjectReference((PetscObject)cheb->pcnone);
51: PCSetType(cheb->pcnone,PCNONE);
52: KSPSetPC(cheb->kspest,ksp->pc);
54: KSPGetInitialGuessNonzero(ksp,&nonzero);
55: KSPSetInitialGuessNonzero(cheb->kspest,nonzero);
56: KSPSetComputeSingularValues(cheb->kspest,PETSC_TRUE);
58: /* Estimate with a fixed number of iterations */
59: KSPSetConvergenceTest(cheb->kspest,KSPSkipConverged,0,0);
60: KSPSetNormType(cheb->kspest,KSP_NORM_NONE);
61: KSPSetTolerances(cheb->kspest,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,10);
62: }
63: if (a >= 0) cheb->tform[0] = a;
64: if (b >= 0) cheb->tform[1] = b;
65: if (c >= 0) cheb->tform[2] = c;
66: if (d >= 0) cheb->tform[3] = d;
67: } else {
68: KSPDestroy(&cheb->kspest);
69: PCDestroy(&cheb->pcnone);
70: }
71: return(0);
72: }
73: EXTERN_C_END
77: /*@
78: KSPChebychevSetEigenvalues - Sets estimates for the extreme eigenvalues
79: of the preconditioned problem.
81: Logically Collective on KSP
83: Input Parameters:
84: + ksp - the Krylov space context
85: - emax, emin - the eigenvalue estimates
87: Options Database:
88: . -ksp_chebychev_eigenvalues emin,emax
90: Note: If you run with the Krylov method of KSPCG with the option -ksp_monitor_singular_value it will
91: for that given matrix and preconditioner an estimate of the extreme eigenvalues.
93: Level: intermediate
95: .keywords: KSP, Chebyshev, set, eigenvalues
96: @*/
97: PetscErrorCode KSPChebychevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
98: {
105: PetscTryMethod(ksp,"KSPChebychevSetEigenvalues_C",(KSP,PetscReal,PetscReal),(ksp,emax,emin));
106: return(0);
107: }
111: /*@
112: KSPChebychevSetEstimateEigenvalues - Automatically estimate the eigenvalues to use for Chebychev
114: Logically Collective on KSP
116: Input Parameters:
117: + ksp - the Krylov space context
118: . a - multiple of min eigenvalue estimate to use for min Chebychev bound (or PETSC_DECIDE)
119: . b - multiple of max eigenvalue estimate to use for min Chebychev bound (or PETSC_DECIDE)
120: . c - multiple of min eigenvalue estimate to use for max Chebychev bound (or PETSC_DECIDE)
121: - d - multiple of max eigenvalue estimate to use for max Chebychev bound (or PETSC_DECIDE)
123: Options Database:
124: . -ksp_chebychev_estimate_eigenvalues a,b,c,d
126: Notes:
127: The Chebychev bounds are estimated using
128: .vb
129: minbound = a*minest + b*maxest
130: maxbound = c*minest + d*maxest
131: .ve
132: The default configuration targets the upper part of the spectrum for use as a multigrid smoother, so only the maximum eigenvalue estimate is used.
133: The minimum eigenvalue estimate obtained by Krylov iteration is typically not accurate until the method has converged.
135: If 0.0 is passed for all transform arguments (a,b,c,d), eigenvalue estimation is disabled.
137: Level: intermediate
139: .keywords: KSP, Chebyshev, set, eigenvalues
140: @*/
141: PetscErrorCode KSPChebychevSetEstimateEigenvalues(KSP ksp,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
142: {
151: PetscTryMethod(ksp,"KSPChebychevSetEstimateEigenvalues_C",(KSP,PetscReal,PetscReal,PetscReal,PetscReal),(ksp,a,b,c,d));
152: return(0);
153: }
157: PetscErrorCode KSPSetFromOptions_Chebychev(KSP ksp)
158: {
159: KSP_Chebychev *cheb = (KSP_Chebychev*)ksp->data;
161: PetscInt two = 2,four = 4;
162: PetscReal tform[4] = {PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE};
163: PetscBool flg;
166: PetscOptionsHead("KSP Chebychev Options");
167: PetscOptionsRealArray("-ksp_chebychev_eigenvalues","extreme eigenvalues","KSPChebychevSetEigenvalues",&cheb->emin,&two,0);
168: PetscOptionsRealArray("-ksp_chebychev_estimate_eigenvalues","estimate eigenvalues using a Krylov method, then use this transform for Chebychev eigenvalue bounds","KSPChebychevSetEstimateEigenvalues",tform,&four,&flg);
169: if (flg) {KSPChebychevSetEstimateEigenvalues(ksp,tform[0],tform[1],tform[2],tform[3]);}
170: if (cheb->kspest) {
171: /* Mask the PC so that PCSetFromOptions does not do anything */
172: KSPSetPC(cheb->kspest,cheb->pcnone);
173: KSPSetOptionsPrefix(cheb->kspest,((PetscObject)ksp)->prefix);
174: KSPAppendOptionsPrefix(cheb->kspest,"est_");
175: if (!((PetscObject)cheb->kspest)->type_name) {
176: KSPSetType(cheb->kspest,KSPGMRES);
177: }
178: KSPSetFromOptions(cheb->kspest);
179: KSPSetPC(cheb->kspest,ksp->pc);
180: }
181: PetscOptionsTail();
182: return(0);
183: }
187: PetscErrorCode KSPSolve_Chebychev(KSP ksp)
188: {
189: KSP_Chebychev *cheb = (KSP_Chebychev*)ksp->data;
191: PetscInt k,kp1,km1,maxit,ktmp,i;
192: PetscScalar alpha,omegaprod,mu,omega,Gamma,c[3],scale;
193: PetscReal rnorm = 0.0;
194: Vec x,b,p[3],r;
195: Mat Amat,Pmat;
196: MatStructure pflag;
197: PetscBool diagonalscale;
200: PCGetDiagonalScale(ksp->pc,&diagonalscale);
201: if (diagonalscale) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
203: if (cheb->kspest && !cheb->estimate_current) {
204: PetscReal max,min;
205: PetscBool nonzero;
206: Vec X = ksp->vec_sol;
207: KSPGetInitialGuessNonzero(ksp,&nonzero);
208: if (nonzero) {VecDuplicate(ksp->vec_sol,&X);}
209: KSPSolve(cheb->kspest,ksp->vec_rhs,X);
210: if (nonzero) {VecDestroy(&X);}
211: else {VecZeroEntries(X);}
212: KSPComputeExtremeSingularValues(cheb->kspest,&max,&min);
213: cheb->emin = cheb->tform[0]*min + cheb->tform[1]*max;
214: cheb->emax = cheb->tform[2]*min + cheb->tform[3]*max;
215: cheb->estimate_current = PETSC_TRUE;
216: }
218: ksp->its = 0;
219: PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);
220: maxit = ksp->max_it;
222: /* These three point to the three active solutions, we
223: rotate these three at each solution update */
224: km1 = 0; k = 1; kp1 = 2;
225: x = ksp->vec_sol;
226: b = ksp->vec_rhs;
227: p[km1] = x;
228: p[k] = ksp->work[0];
229: p[kp1] = ksp->work[1];
230: r = ksp->work[2];
232: /* use scale*B as our preconditioner */
233: scale = 2.0/(cheb->emax + cheb->emin);
235: /* -alpha <= scale*lambda(B^{-1}A) <= alpha */
236: alpha = 1.0 - scale*(cheb->emin); ;
237: Gamma = 1.0;
238: mu = 1.0/alpha;
239: omegaprod = 2.0/alpha;
241: c[km1] = 1.0;
242: c[k] = mu;
244: if (!ksp->guess_zero) {
245: KSP_MatMult(ksp,Amat,x,r); /* r = b - Ax */
246: VecAYPX(r,-1.0,b);
247: } else {
248: VecCopy(b,r);
249: }
250:
251: KSP_PCApply(ksp,r,p[k]); /* p[k] = scale B^{-1}r + x */
252: VecAYPX(p[k],scale,x);
254: for (i=0; i<maxit; i++) {
255: PetscObjectTakeAccess(ksp);
256: ksp->its++;
257: PetscObjectGrantAccess(ksp);
258: c[kp1] = 2.0*mu*c[k] - c[km1];
259: omega = omegaprod*c[k]/c[kp1];
261: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
262: VecAYPX(r,-1.0,b);
263: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}z */
265: /* calculate residual norm if requested */
266: if (ksp->normtype != KSP_NORM_NONE) {
267: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {VecNorm(r,NORM_2,&rnorm);}
268: else {VecNorm(p[kp1],NORM_2,&rnorm);}
269: PetscObjectTakeAccess(ksp);
270: ksp->rnorm = rnorm;
271: PetscObjectGrantAccess(ksp);
272: ksp->vec_sol = p[k];
273: KSPLogResidualHistory(ksp,rnorm);
274: KSPMonitor(ksp,i,rnorm);
275: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
276: if (ksp->reason) break;
277: }
279: /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
280: VecScale(p[kp1],omega*Gamma*scale);
281: VecAXPY(p[kp1],1.0-omega,p[km1]);
282: VecAXPY(p[kp1],omega,p[k]);
284: ktmp = km1;
285: km1 = k;
286: k = kp1;
287: kp1 = ktmp;
288: }
289: if (!ksp->reason) {
290: if (ksp->normtype != KSP_NORM_NONE) {
291: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
292: VecAYPX(r,-1.0,b);
293: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED) {
294: VecNorm(r,NORM_2,&rnorm);
295: } else {
296: KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}z */
297: VecNorm(p[kp1],NORM_2,&rnorm);
298: }
299: PetscObjectTakeAccess(ksp);
300: ksp->rnorm = rnorm;
301: PetscObjectGrantAccess(ksp);
302: ksp->vec_sol = p[k];
303: KSPLogResidualHistory(ksp,rnorm);
304: KSPMonitor(ksp,i,rnorm);
305: }
306: if (ksp->its >= ksp->max_it) {
307: if (ksp->normtype != KSP_NORM_NONE) {
308: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
309: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
310: } else {
311: ksp->reason = KSP_CONVERGED_ITS;
312: }
313: }
314: }
316: /* make sure solution is in vector x */
317: ksp->vec_sol = x;
318: if (k) {
319: VecCopy(p[k],x);
320: }
321: return(0);
322: }
326: PetscErrorCode KSPView_Chebychev(KSP ksp,PetscViewer viewer)
327: {
328: KSP_Chebychev *cheb = (KSP_Chebychev*)ksp->data;
330: PetscBool iascii;
333: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
334: if (iascii) {
335: PetscViewerASCIIPrintf(viewer," Chebychev: eigenvalue estimates: min = %G, max = %G\n",cheb->emin,cheb->emax);
336: if (cheb->kspest) {
337: PetscViewerASCIIPushTab(viewer);
338: KSPView(cheb->kspest,viewer);
339: PetscViewerASCIIPopTab(viewer);
340: }
341: } else {
342: SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for KSP Chebychev",((PetscObject)viewer)->type_name);
343: }
344: return(0);
345: }
349: PetscErrorCode KSPDestroy_Chebychev(KSP ksp)
350: {
351: KSP_Chebychev *cheb = (KSP_Chebychev*)ksp->data;
355: KSPDestroy(&cheb->kspest);
356: PCDestroy(&cheb->pcnone);
357: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPChebychevSetEigenvalues_C","",PETSC_NULL);
358: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPChebychevSetEstimateEigenvalues_C","",PETSC_NULL);
359: KSPDefaultDestroy(ksp);
360: return(0);
361: }
363: /*MC
364: KSPCHEBYCHEV - The preconditioned Chebychev iterative method
366: Options Database Keys:
367: + -ksp_chebychev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
368: of the preconditioned operator. If these are accurate you will get much faster convergence.
369: - -ksp_chebychev_estimate_eigenvalues <a,b,c,d> - estimate eigenvalues using a Krylov method, then use this
370: transform for Chebychev eigenvalue bounds (KSPChebychevSetEstimateEigenvalues)
373: Level: beginner
375: Notes: The Chebychev method requires both the matrix and preconditioner to
376: be symmetric positive (semi) definite.
377: Only support for left preconditioning.
379: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
380: KSPChebychevSetEigenvalues(), KSPRICHARDSON, KSPCG
382: M*/
384: EXTERN_C_BEGIN
387: PetscErrorCode KSPCreate_Chebychev(KSP ksp)
388: {
390: KSP_Chebychev *chebychevP;
393: PetscNewLog(ksp,KSP_Chebychev,&chebychevP);
395: ksp->data = (void*)chebychevP;
396: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
397: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_LEFT,1);
399: chebychevP->emin = 1.e-2;
400: chebychevP->emax = 1.e+2;
402: chebychevP->tform[0] = 0.0;
403: chebychevP->tform[1] = 0.02;
404: chebychevP->tform[1] = 0;
405: chebychevP->tform[2] = 1.1;
407: ksp->ops->setup = KSPSetUp_Chebychev;
408: ksp->ops->solve = KSPSolve_Chebychev;
409: ksp->ops->destroy = KSPDestroy_Chebychev;
410: ksp->ops->buildsolution = KSPDefaultBuildSolution;
411: ksp->ops->buildresidual = KSPDefaultBuildResidual;
412: ksp->ops->setfromoptions = KSPSetFromOptions_Chebychev;
413: ksp->ops->view = KSPView_Chebychev;
415: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPChebychevSetEigenvalues_C",
416: "KSPChebychevSetEigenvalues_Chebychev",
417: KSPChebychevSetEigenvalues_Chebychev);
418: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPChebychevSetEstimateEigenvalues_C",
419: "KSPChebychevSetEstimateEigenvalues_Chebychev",
420: KSPChebychevSetEstimateEigenvalues_Chebychev);
421: return(0);
422: }
423: EXTERN_C_END