Actual source code: cheby.c

      1: #define PETSCKSP_DLL

      3: /*
      4:     This is a first attempt at a Chebychev routine, it is not
      5:     necessarily well optimized.
      6: */
     7:  #include src/ksp/ksp/kspimpl.h
     8:  #include src/ksp/ksp/impls/cheby/chebctx.h

     12: PetscErrorCode KSPSetUp_Chebychev(KSP ksp)
     13: {

     17:   if (ksp->pc_side == PC_SYMMETRIC) SETERRQ(PETSC_ERR_SUP,"no symmetric preconditioning for KSPCHEBYCHEV");
     18:   KSPDefaultGetWork(ksp,3);
     19:   return(0);
     20: }

     25: PetscErrorCode PETSCKSP_DLLEXPORT KSPChebychevSetEigenvalues_Chebychev(KSP ksp,PetscReal emax,PetscReal emin)
     26: {
     27:   KSP_Chebychev *chebychevP = (KSP_Chebychev*)ksp->data;

     30:   if (emax <= emin) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Maximum eigenvalue must be larger than minimum: max %g min %g",emax,emin);
     31:   if (emax*emin <= 0.0) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Both eigenvalues must be of the same sign: max %g min %g",emax,emin);
     32:   chebychevP->emax = emax;
     33:   chebychevP->emin = emin;
     34:   return(0);
     35: }

     40: /*@
     41:    KSPChebychevSetEigenvalues - Sets estimates for the extreme eigenvalues
     42:    of the preconditioned problem.

     44:    Collective on KSP

     46:    Input Parameters:
     47: +  ksp - the Krylov space context
     48: -  emax, emin - the eigenvalue estimates

     50:   Options Database:
     51: .  -ksp_chebychev_eigenvalues emin,emax

     53:    Level: intermediate

     55: .keywords: KSP, Chebyshev, set, eigenvalues
     56: @*/
     57: PetscErrorCode PETSCKSP_DLLEXPORT KSPChebychevSetEigenvalues(KSP ksp,PetscReal emax,PetscReal emin)
     58: {
     59:   PetscErrorCode ierr,(*f)(KSP,PetscReal,PetscReal);

     63:   PetscObjectQueryFunction((PetscObject)ksp,"KSPChebychevSetEigenvalues_C",(void (**)(void))&f);
     64:   if (f) {
     65:     (*f)(ksp,emax,emin);
     66:   }
     67:   return(0);
     68: }

     72: PetscErrorCode KSPSetFromOptions_Chebychev(KSP ksp)
     73: {
     74:   KSP_Chebychev *cheb = (KSP_Chebychev*)ksp->data;
     76:   PetscInt       two = 2;

     79:   PetscOptionsHead("KSP Chebychev Options");
     80:     PetscOptionsRealArray("-ksp_chebychev_eigenvalues","extreme eigenvalues","KSPChebychevSetEigenvalues",&cheb->emin,&two,0);
     81:   PetscOptionsTail();
     82:   return(0);
     83: }

     87: PetscErrorCode KSPSolve_Chebychev(KSP ksp)
     88: {
     90:   PetscInt       k,kp1,km1,maxit,ktmp,i;
     91:   PetscScalar    alpha,omegaprod,mu,omega,Gamma,c[3],scale,mone = -1.0,tmp;
     92:   PetscReal      rnorm;
     93:   Vec            x,b,p[3],r;
     94:   KSP_Chebychev  *chebychevP = (KSP_Chebychev*)ksp->data;
     95:   Mat            Amat,Pmat;
     96:   MatStructure   pflag;
     97:   PetscTruth     diagonalscale;

    100:   PCDiagonalScale(ksp->pc,&diagonalscale);
    101:   if (diagonalscale) SETERRQ1(PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",ksp->type_name);

    103:   ksp->its = 0;
    104:   PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);
    105:   maxit    = ksp->max_it;

    107:   /* These three point to the three active solutions, we
    108:      rotate these three at each solution update */
    109:   km1    = 0; k = 1; kp1 = 2;
    110:   x      = ksp->vec_sol;
    111:   b      = ksp->vec_rhs;
    112:   p[km1] = x;
    113:   p[k]   = ksp->work[0];
    114:   p[kp1] = ksp->work[1];
    115:   r      = ksp->work[2];

    117:   /* use scale*B as our preconditioner */
    118:   scale  = 2.0/(chebychevP->emax + chebychevP->emin);

    120:   /*   -alpha <=  scale*lambda(B^{-1}A) <= alpha   */
    121:   alpha  = 1.0 - scale*(chebychevP->emin); ;
    122:   Gamma  = 1.0;
    123:   mu     = 1.0/alpha;
    124:   omegaprod = 2.0/alpha;

    126:   c[km1] = 1.0;
    127:   c[k]   = mu;

    129:   if (!ksp->guess_zero) {
    130:     KSP_MatMult(ksp,Amat,x,r);     /*  r = b - Ax     */
    131:     VecAYPX(r,mone,b);
    132:   } else {
    133:     VecCopy(b,r);
    134:   }
    135:
    136:   KSP_PCApply(ksp,r,p[k]);  /* p[k] = scale B^{-1}r + x */
    137:   VecAYPX(p[k],scale,x);

    139:   for (i=0; i<maxit; i++) {
    140:     PetscObjectTakeAccess(ksp);
    141:     ksp->its++;
    142:     PetscObjectGrantAccess(ksp);
    143:     c[kp1] = 2.0*mu*c[k] - c[km1];
    144:     omega = omegaprod*c[k]/c[kp1];

    146:     KSP_MatMult(ksp,Amat,p[k],r);                 /*  r = b - Ap[k]    */
    147:     VecAYPX(r,mone,b);
    148:     KSP_PCApply(ksp,r,p[kp1]);             /*  p[kp1] = B^{-1}z  */

    150:     /* calculate residual norm if requested */
    151:     if (ksp->normtype != KSP_NO_NORM) {
    152:       if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) {VecNorm(r,NORM_2,&rnorm);}
    153:       else {VecNorm(p[kp1],NORM_2,&rnorm);}
    154:       PetscObjectTakeAccess(ksp);
    155:       ksp->rnorm                              = rnorm;
    156:       PetscObjectGrantAccess(ksp);
    157:       ksp->vec_sol = p[k];
    158:       KSPLogResidualHistory(ksp,rnorm);
    159:       KSPMonitor(ksp,i,rnorm);
    160:       (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
    161:       if (ksp->reason) break;
    162:     }

    164:     /* y^{k+1} = omega(y^{k} - y^{k-1} + Gamma*r^{k}) + y^{k-1} */
    165:     tmp  = omega*Gamma*scale;
    166:     VecScale(p[kp1],tmp);
    167:     tmp  = 1.0-omega; VecAXPY(p[kp1],tmp,p[km1]);
    168:     VecAXPY(p[kp1],omega,p[k]);

    170:     ktmp = km1;
    171:     km1  = k;
    172:     k    = kp1;
    173:     kp1  = ktmp;
    174:   }
    175:   if (!ksp->reason && ksp->normtype != KSP_NO_NORM) {
    176:     ksp->reason = KSP_DIVERGED_ITS;
    177:     KSP_MatMult(ksp,Amat,p[k],r);       /*  r = b - Ap[k]    */
    178:     VecAYPX(r,mone,b);
    179:     if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) {VecNorm(r,NORM_2,&rnorm);}
    180:     else {
    181:       KSP_PCApply(ksp,r,p[kp1]); /* p[kp1] = B^{-1}z */
    182:       VecNorm(p[kp1],NORM_2,&rnorm);
    183:     }
    184:     PetscObjectTakeAccess(ksp);
    185:     ksp->rnorm                              = rnorm;
    186:     PetscObjectGrantAccess(ksp);
    187:     ksp->vec_sol = p[k];
    188:     KSPLogResidualHistory(ksp,rnorm);
    189:     KSPMonitor(ksp,i,rnorm);
    190:   }

    192:   /* make sure solution is in vector x */
    193:   ksp->vec_sol = x;
    194:   if (k) {
    195:     VecCopy(p[k],x);
    196:   }
    197:   return(0);
    198: }

    202: PetscErrorCode KSPView_Chebychev(KSP ksp,PetscViewer viewer)
    203: {
    204:   KSP_Chebychev  *cheb = (KSP_Chebychev*)ksp->data;
    206:   PetscTruth     iascii;

    209:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    210:   if (iascii) {
    211:     PetscViewerASCIIPrintf(viewer,"  Chebychev: eigenvalue estimates:  min = %g, max = %g\n",cheb->emin,cheb->emax);
    212:   } else {
    213:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for KSP Chebychev",((PetscObject)viewer)->type_name);
    214:   }
    215:   return(0);
    216: }

    218: /*MC
    219:      KSPCHEBYCHEV - The preconditioned Chebychev iterative method

    221:    Options Database Keys:
    222: .   -ksp_chebychev_eigenvalues <emin,emax> - set approximations to the smallest and largest eigenvalues
    223:                   of the preconditioned operator. If these are accurate you will get much faster convergence.

    225:    Level: beginner

    227:    Notes: The Chebychev method requires both the matrix and preconditioner to
    228:           be symmetric positive (semi) definite

    230: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
    231:            KSPChebychevSetEigenvalues()

    233: M*/

    238: PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_Chebychev(KSP ksp)
    239: {
    241:   KSP_Chebychev  *chebychevP;

    244:   PetscNew(KSP_Chebychev,&chebychevP);
    245:   PetscLogObjectMemory(ksp,sizeof(KSP_Chebychev));

    247:   ksp->data                      = (void*)chebychevP;
    248:   ksp->pc_side                   = PC_LEFT;

    250:   chebychevP->emin               = 1.e-2;
    251:   chebychevP->emax               = 1.e+2;

    253:   ksp->ops->setup                = KSPSetUp_Chebychev;
    254:   ksp->ops->solve                = KSPSolve_Chebychev;
    255:   ksp->ops->destroy              = KSPDefaultDestroy;
    256:   ksp->ops->buildsolution        = KSPDefaultBuildSolution;
    257:   ksp->ops->buildresidual        = KSPDefaultBuildResidual;
    258:   ksp->ops->setfromoptions       = KSPSetFromOptions_Chebychev;
    259:   ksp->ops->view                 = KSPView_Chebychev;

    261:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPChebychevSetEigenvalues_C",
    262:                                     "KSPChebychevSetEigenvalues_Chebychev",
    263:                                     KSPChebychevSetEigenvalues_Chebychev);
    264:   return(0);
    265: }