Actual source code: rich.c

      1: #define PETSCKSP_DLL

      3: /*
      4:             This implements Richardson Iteration.
      5: */
     6:  #include src/ksp/ksp/kspimpl.h
     7:  #include src/ksp/ksp/impls/rich/richctx.h

     11: PetscErrorCode KSPSetUp_Richardson(KSP ksp)
     12: {

     16:   if (ksp->pc_side == PC_RIGHT) {SETERRQ(PETSC_ERR_SUP,"no right preconditioning for KSPRICHARDSON");}
     17:   else if (ksp->pc_side == PC_SYMMETRIC) {SETERRQ(PETSC_ERR_SUP,"no symmetric preconditioning for KSPRICHARDSON");}
     18:   KSPDefaultGetWork(ksp,2);
     19:   return(0);
     20: }

     24: PetscErrorCode  KSPSolve_Richardson(KSP ksp)
     25: {
     27:   PetscInt       i,maxit;
     28:   MatStructure   pflag;
     29:   PetscReal      rnorm = 0.0;
     30:   PetscScalar    scale,mone = -1.0,dt;
     31:   Vec            x,b,r,z;
     32:   Mat            Amat,Pmat;
     33:   KSP_Richardson *richardsonP = (KSP_Richardson*)ksp->data;
     34:   PetscTruth     exists,diagonalscale;

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

     40:   ksp->its = 0;

     42:   PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);
     43:   x       = ksp->vec_sol;
     44:   b       = ksp->vec_rhs;
     45:   r       = ksp->work[0];
     46:   z       = ksp->work[1];
     47:   maxit   = ksp->max_it;

     49:   /* if user has provided fast Richardson code use that */
     50:   PCApplyRichardsonExists(ksp->pc,&exists);
     51:   if (exists && !ksp->numbermonitors && !ksp->transpose_solve) {
     52:     ksp->normtype = KSP_NO_NORM;
     53:     PCApplyRichardson(ksp->pc,b,x,r,ksp->rtol,ksp->abstol,ksp->divtol,maxit);
     54:     if (ksp->normtype != KSP_NO_NORM) {
     55:       KSP_MatMult(ksp,Amat,x,r);
     56:       VecAYPX(r,mone,b);
     57:       if (ksp->normtype == KSP_UNPRECONDITIONED_NORM || ksp->pc_side == PC_RIGHT) {
     58:         VecNorm(r,NORM_2,&rnorm); /*   rnorm <- r'*r     */
     59:       } else {
     60:         PetscExceptionTry1((KSP_PCApply(ksp,r,z)),PETSC_ERR_SUP);
     61:         if (PetscExceptionCaught(ierr,PETSC_ERR_SUP)) {
     62:           ksp->reason = KSP_CONVERGED_ITS;
     63:           return(0);
     64:         }
     65:         if (ksp->normtype == KSP_PRECONDITIONED_NORM) {
     66:           VecNorm(z,NORM_2,&rnorm); /*   rnorm <- r'B'*Br     */
     67:         } else {
     68:           VecDot(r,z,&dt);
     69:           rnorm = PetscAbsScalar(dt);                    /*   rnorm <- z'*r  = r'Br = e'*A'*B*A*e  */
     70:         }
     71:       }
     72:       (*ksp->converged)(ksp,0,rnorm,&ksp->reason,ksp->cnvP);
     73:     }
     74:     if (!ksp->reason) ksp->reason = KSP_CONVERGED_ITS;
     75:     return(0);
     76:   }

     78:   scale   = richardsonP->scale;

     80:   if (!ksp->guess_zero) {                          /*   r <- b - A x     */
     81:     KSP_MatMult(ksp,Amat,x,r);
     82:     VecAYPX(r,mone,b);
     83:   } else {
     84:     VecCopy(b,r);
     85:   }

     87:   for (i=0; i<maxit; i++) {
     88:     PetscObjectTakeAccess(ksp);
     89:     ksp->its++;
     90:     PetscObjectGrantAccess(ksp);

     92:     if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) {
     93:       VecNorm(r,NORM_2,&rnorm); /*   rnorm <- r'*r     */
     94:       KSPMonitor(ksp,i,rnorm);
     95:     }

     97:     KSP_PCApply(ksp,r,z);    /*   z <- B r          */

     99:     if (ksp->normtype == KSP_PRECONDITIONED_NORM) {
    100:       VecNorm(z,NORM_2,&rnorm); /*   rnorm <- z'*z     */
    101:       KSPMonitor(ksp,i,rnorm);
    102:     }

    104:     VecAXPY(x,scale,z);    /*   x  <- x + scale z */
    105:     if (ksp->normtype != KSP_NO_NORM) {
    106:       PetscObjectTakeAccess(ksp);
    107:       ksp->rnorm = rnorm;
    108:       PetscObjectGrantAccess(ksp);
    109:       KSPLogResidualHistory(ksp,rnorm);

    111:       (*ksp->converged)(ksp,i,rnorm,&ksp->reason,ksp->cnvP);
    112:       if (ksp->reason) break;
    113:     }
    114:
    115:     KSP_MatMult(ksp,Amat,x,r);      /*   r  <- b - Ax      */
    116:     VecAYPX(r,mone,b);
    117:   }
    118:   if (!ksp->reason) {
    119:     ksp->reason = KSP_DIVERGED_ITS;
    120:     if (ksp->normtype != KSP_NO_NORM) {
    121:       if (ksp->normtype == KSP_UNPRECONDITIONED_NORM){
    122:         VecNorm(r,NORM_2,&rnorm);     /*   rnorm <- r'*r     */
    123:       } else {
    124:         KSP_PCApply(ksp,r,z);   /*   z <- B r          */
    125:         VecNorm(z,NORM_2,&rnorm);     /*   rnorm <- z'*z     */
    126:       }
    127:     }
    128:     PetscObjectTakeAccess(ksp);
    129:     ksp->rnorm = rnorm;
    130:     PetscObjectGrantAccess(ksp);
    131:     KSPLogResidualHistory(ksp,rnorm);
    132:     KSPMonitor(ksp,i,rnorm);
    133:     i--;
    134:   } else if (!ksp->reason) {
    135:     ksp->reason = KSP_DIVERGED_ITS;
    136:   }

    138:   return(0);
    139: }

    143: PetscErrorCode KSPView_Richardson(KSP ksp,PetscViewer viewer)
    144: {
    145:   KSP_Richardson *richardsonP = (KSP_Richardson*)ksp->data;
    147:   PetscTruth     iascii;

    150:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    151:   if (iascii) {
    152:     PetscViewerASCIIPrintf(viewer,"  Richardson: damping factor=%g\n",richardsonP->scale);
    153:   } else {
    154:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for KSP Richardson",((PetscObject)viewer)->type_name);
    155:   }
    156:   return(0);
    157: }

    161: PetscErrorCode KSPSetFromOptions_Richardson(KSP ksp)
    162: {
    163:   KSP_Richardson *rich = (KSP_Richardson*)ksp->data;
    165:   PetscReal      tmp;
    166:   PetscTruth     flg;

    169:   PetscOptionsHead("KSP Richardson Options");
    170:     PetscOptionsReal("-ksp_richardson_scale","damping factor","KSPRichardsonSetScale",rich->scale,&tmp,&flg);
    171:     if (flg) { KSPRichardsonSetScale(ksp,tmp); }
    172:   PetscOptionsTail();
    173:   return(0);
    174: }

    179: PetscErrorCode PETSCKSP_DLLEXPORT KSPRichardsonSetScale_Richardson(KSP ksp,PetscReal scale)
    180: {
    181:   KSP_Richardson *richardsonP;

    184:   richardsonP = (KSP_Richardson*)ksp->data;
    185:   richardsonP->scale = scale;
    186:   return(0);
    187: }

    190: /*MC
    191:      KSPRICHARDSON - The preconditioned Richardson iterative method

    193:    Options Database Keys:
    194: .   -ksp_richardson_scale - damping factor on the correction (defaults to 1.0)

    196:    Level: beginner

    198:    Notes: x^{n+1} = x^{n} + scale*B(b - A x^{n})

    200:           This method often (usually) will not converge unless scale is very small

    202: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP,
    203:            KSPRichardsonSetScale()

    205: M*/

    210: PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_Richardson(KSP ksp)
    211: {
    213:   KSP_Richardson *richardsonP;

    216:   PetscNew(KSP_Richardson,&richardsonP);
    217:   PetscLogObjectMemory(ksp,sizeof(KSP_Richardson));
    218:   ksp->data                        = (void*)richardsonP;
    219:   richardsonP->scale               = 1.0;
    220:   ksp->ops->setup                  = KSPSetUp_Richardson;
    221:   ksp->ops->solve                  = KSPSolve_Richardson;
    222:   ksp->ops->destroy                = KSPDefaultDestroy;
    223:   ksp->ops->buildsolution          = KSPDefaultBuildSolution;
    224:   ksp->ops->buildresidual          = KSPDefaultBuildResidual;
    225:   ksp->ops->view                   = KSPView_Richardson;
    226:   ksp->ops->setfromoptions         = KSPSetFromOptions_Richardson;

    228:   PetscObjectComposeFunctionDynamic((PetscObject)ksp,"KSPRichardsonSetScale_C",
    229:                                     "KSPRichardsonSetScale_Richardson",
    230:                                     KSPRichardsonSetScale_Richardson);
    231:   return(0);
    232: }