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