Actual source code: cheby.c
1: /*$Id: cheby.c,v 1.94 2001/08/07 03:03:53 balay Exp $*/
2: /*
3: This is a first attempt at a Chebychev routine, it is not
4: necessarily well optimized.
5: */
6: #include src/sles/ksp/kspimpl.h
7: #include src/sles/ksp/impls/cheby/chebctx.h
9: #undef __FUNCT__
11: int KSPSetUp_Chebychev(KSP ksp)
12: {
16: if (ksp->pc_side == PC_SYMMETRIC) SETERRQ(2,"no symmetric preconditioning for KSPCHEBYCHEV");
17: KSPDefaultGetWork(ksp,3);
18: return(0);
19: }
21: EXTERN_C_BEGIN
22: #undef __FUNCT__
24: int KSPChebychevSetEigenvalues_Chebychev(KSP ksp,PetscReal emax,PetscReal emin)
25: {
26: KSP_Chebychev *chebychevP = (KSP_Chebychev*)ksp->data;
29: chebychevP->emax = emax;
30: chebychevP->emin = emin;
31: return(0);
32: }
33: EXTERN_C_END
35: #undef __FUNCT__
37: /*@
38: KSPChebychevSetEigenvalues - Sets estimates for the extreme eigenvalues
39: of the preconditioned problem.
41: Collective on KSP
43: Input Parameters:
44: + ksp - the Krylov space context
45: - emax, emin - the eigenvalue estimates
47: Level: intermediate
49: .keywords: KSP, Chebyshev, set, eigenvalues
50: @*/
51: int KSPChebychevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
52: {
53: int ierr,(*f)(KSP,PetscReal,PetscReal);
57: PetscObjectQueryFunction((PetscObject)ksp,"KSPChebychevSetEigenvalues_C",(void (**)(void))&f);
58: if (f) {
59: (*f)(ksp,emax,emin);
60: }
61: return(0);
62: }
64: #undef __FUNCT__
66: int KSPSolve_Chebychev(KSP ksp,int *its)
67: {
68: int k,kp1,km1,maxit,ktmp,i,ierr;
69: PetscScalar alpha,omegaprod,mu,omega,Gamma,c[3],scale,mone = -1.0,tmp;
70: PetscReal rnorm;
71: Vec x,b,p[3],r;
72: KSP_Chebychev *chebychevP = (KSP_Chebychev*)ksp->data;
73: Mat Amat,Pmat;
74: MatStructure pflag;
75: PetscTruth diagonalscale;
78: ierr = PCDiagonalScale(ksp->B,&diagonalscale);
79: if (diagonalscale) SETERRQ1(1,"Krylov method %s does not support diagonal scaling",ksp->type_name);
81: ksp->its = 0;
82: ierr = PCGetOperators(ksp->B,&Amat,&Pmat,&pflag);
83: maxit = ksp->max_it;
85: /* These three point to the three active solutions, we
86: rotate these three at each solution update */
87: km1 = 0; k = 1; kp1 = 2;
88: x = ksp->vec_sol;
89: b = ksp->vec_rhs;
90: p[km1] = x;
91: p[k] = ksp->work[0];
92: p[kp1] = ksp->work[1];
93: r = ksp->work[2];
95: /* use scale*B as our preconditioner */
96: scale = 2.0/(chebychevP->emax + chebychevP->emin);
98: /* -alpha <= scale*lambda(B^{-1}A) <= alpha */
99: alpha = 1.0 - scale*(chebychevP->emin); ;
100: Gamma = 1.0;
101: mu = 1.0/alpha;
102: omegaprod = 2.0/alpha;
104: c[km1] = 1.0;
105: c[k] = mu;
107: if (!ksp->guess_zero) {
108: KSP_MatMult(ksp,Amat,x,r); /* r = b - Ax */
109: VecAYPX(&mone,b,r);
110: } else {
111: VecCopy(b,r);
112: }
113:
114: KSP_PCApply(ksp,ksp->B,r,p[k]); /* p[k] = scale B^{-1}r + x */
115: VecAYPX(&scale,x,p[k]);
117: for (i=0; i<maxit; i++) {
118: PetscObjectTakeAccess(ksp);
119: ksp->its++;
120: PetscObjectGrantAccess(ksp);
121: c[kp1] = 2.0*mu*c[k] - c[km1];
122: omega = omegaprod*c[k]/c[kp1];
124: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
125: VecAYPX(&mone,b,r);
126: KSP_PCApply(ksp,ksp->B,r,p[kp1]); /* p[kp1] = B^{-1}z */
128: /* calculate residual norm if requested */
129: if (ksp->normtype != KSP_NO_NORM) {
130: if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) {VecNorm(r,NORM_2,&rnorm);}
131: else {VecNorm(p[kp1],NORM_2,&rnorm);}
132: PetscObjectTakeAccess(ksp);
133: ksp->rnorm = rnorm;
134: PetscObjectGrantAccess(ksp);
135: ksp->vec_sol = p[k];
136: KSPLogResidualHistory(ksp,rnorm);
137: KSPMonitor(ksp,i,rnorm);
138: (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
139: if (ksp->reason) break;
140: }
142: /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
143: tmp = omega*Gamma*scale;
144: VecScale(&tmp,p[kp1]);
145: tmp = 1.0-omega; VecAXPY(&tmp,p[km1],p[kp1]);
146: VecAXPY(&omega,p[k],p[kp1]);
148: ktmp = km1;
149: km1 = k;
150: k = kp1;
151: kp1 = ktmp;
152: }
153: if (!ksp->reason && ksp->normtype != KSP_NO_NORM) {
154: ksp->reason = KSP_DIVERGED_ITS;
155: KSP_MatMult(ksp,Amat,p[k],r); /* r = b - Ap[k] */
156: VecAYPX(&mone,b,r);
157: if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) {VecNorm(r,NORM_2,&rnorm);}
158: else {
159: KSP_PCApply(ksp,ksp->B,r,p[kp1]); /* p[kp1] = B^{-1}z */
160: VecNorm(p[kp1],NORM_2,&rnorm);
161: }
162: PetscObjectTakeAccess(ksp);
163: ksp->rnorm = rnorm;
164: PetscObjectGrantAccess(ksp);
165: ksp->vec_sol = p[k];
166: KSPLogResidualHistory(ksp,rnorm);
167: KSPMonitor(ksp,i,rnorm);
168: }
170: /* make sure solution is in vector x */
171: ksp->vec_sol = x;
172: if (k) {
173: VecCopy(p[k],x);
174: }
175: *its = ksp->its;
176: return(0);
177: }
179: #undef __FUNCT__
181: int KSPView_Chebychev(KSP ksp,PetscViewer viewer)
182: {
183: KSP_Chebychev *cheb = (KSP_Chebychev*)ksp->data;
184: int ierr;
185: PetscTruth isascii;
188: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
189: if (isascii) {
190: PetscViewerASCIIPrintf(viewer," Chebychev: eigenvalue estimates: min = %g, max = %gn",cheb->emin,cheb->emax);
191: } else {
192: SETERRQ1(1,"Viewer type %s not supported for KSP Chebychev",((PetscObject)viewer)->type_name);
193: }
194: return(0);
195: }
197: EXTERN_C_BEGIN
198: #undef __FUNCT__
200: int KSPCreate_Chebychev(KSP ksp)
201: {
202: int ierr;
203: KSP_Chebychev *chebychevP;
206: PetscNew(KSP_Chebychev,&chebychevP);
207: PetscLogObjectMemory(ksp,sizeof(KSP_Chebychev));
209: ksp->data = (void*)chebychevP;
210: ksp->pc_side = PC_LEFT;
212: chebychevP->emin = 1.e-2;
213: chebychevP->emax = 1.e+2;
215: ksp->ops->setup = KSPSetUp_Chebychev;
216: ksp->ops->solve = KSPSolve_Chebychev;
217: ksp->ops->destroy = KSPDefaultDestroy;
218: ksp->ops->buildsolution = KSPDefaultBuildSolution;
219: ksp->ops->buildresidual = KSPDefaultBuildResidual;
220: ksp->ops->setfromoptions = 0;
221: ksp->ops->view = KSPView_Chebychev;
223: PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPChebychevSetEigenvalues_C",
224: "KSPChebychevSetEigenvalues_Chebychev",
225: KSPChebychevSetEigenvalues_Chebychev);
226: return(0);
227: }
228: EXTERN_C_END