Actual source code: cr.c

      1: #define PETSCKSP_DLL

     3:  #include src/ksp/ksp/kspimpl.h

      7: static PetscErrorCode KSPSetUp_CR(KSP ksp)
      8: {

     12:   if (ksp->pc_side == PC_RIGHT) {SETERRQ(PETSC_ERR_SUP,"no right preconditioning for KSPCR");}
     13:   else if (ksp->pc_side == PC_SYMMETRIC) {SETERRQ(PETSC_ERR_SUP,"no symmetric preconditioning for KSPCR");}
     14:   KSPDefaultGetWork(ksp,6);
     15:   return(0);
     16: }

     20: static PetscErrorCode  KSPSolve_CR(KSP ksp)
     21: {
     23:   PetscInt       i = 0;
     24:   MatStructure   pflag;
     25:   PetscReal      dp;
     26:   PetscScalar    ai, bi;
     27:   PetscScalar    apq,btop, bbot, mone = -1.0;
     28:   Vec            X,B,R,RT,P,AP,ART,Q;
     29:   Mat            Amat, Pmat;

     32:   X       = ksp->vec_sol;
     33:   B       = ksp->vec_rhs;
     34:   R       = ksp->work[0];
     35:   RT      = ksp->work[1];
     36:   P       = ksp->work[2];
     37:   AP      = ksp->work[3];
     38:   ART     = ksp->work[4];
     39:   Q       = ksp->work[5];

     41:   /* R is the true residual norm, RT is the preconditioned residual norm */
     42:   PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);
     43:   if (!ksp->guess_zero) {
     44:     KSP_MatMult(ksp,Amat,X,R);     /*   R <- A*X           */
     45:     VecAYPX(R,mone,B);            /*   R <- B-R == B-A*X  */
     46:   } else {
     47:     VecCopy(B,R);                  /*   R <- B (X is 0)    */
     48:   }
     49:   KSP_PCApply(ksp,R,P);     /*   P   <- B*R         */
     50:   KSP_MatMult(ksp,Amat,P,AP);      /*   AP  <- A*P         */
     51:   VecCopy(P,RT);                   /*   RT  <- P           */
     52:   VecCopy(AP,ART);                 /*   ART <- AP          */
     53:   VecDot(RT,ART,&btop);          /*   (RT,ART)           */
     54:   if (PetscAbsScalar(btop) < 0.0) {
     55:     ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
     56:     PetscLogInfo((ksp,"KSPSolve_CR:diverging due to indefinite or negative definite matrix\n"));
     57:     return(0);
     58:   }
     59:
     60:   if (ksp->normtype == KSP_PRECONDITIONED_NORM) {
     61:     VecNorm(RT,NORM_2,&dp);        /*   dp <- RT'*RT       */
     62:   } else if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) {
     63:     VecNorm(R,NORM_2,&dp);         /*   dp <- R'*R         */
     64:   } else if (ksp->normtype == KSP_NATURAL_NORM) {
     65:     dp = sqrt(PetscAbsScalar(btop));                    /* dp = sqrt(R,AR)      */
     66:   }
     67:   ksp->its = 0;
     68:   KSPMonitor(ksp,0,dp);
     69:   PetscObjectTakeAccess(ksp);
     70:   ksp->rnorm              = dp;
     71:   PetscObjectGrantAccess(ksp);
     72:   KSPLogResidualHistory(ksp,dp);
     73:   (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);
     74:   if (ksp->reason) return(0);

     76:   i = 0;
     77:   do {
     78:     KSP_PCApply(ksp,AP,Q);/*   Q <- B* AP          */

     80:     VecDot(AP,Q,&apq);
     81:     if (PetscAbsScalar(apq) <= 0.0) {
     82:       ksp->reason = KSP_DIVERGED_INDEFINITE_PC;
     83:       PetscLogInfo((ksp,"KSPSolve_CR:diverging due to indefinite or negative definite PC\n"));
     84:       break;
     85:     }
     86:     ai = btop/apq;                                      /* ai = (RT,ART)/(AP,Q)  */

     88:     VecAXPY(X,ai,P);            /*   X   <- X + ai*P     */
     89:     ai     = -ai;
     90:     VecAXPY(RT,ai,Q);           /*   RT  <- RT - ai*Q    */
     91:     KSP_MatMult(ksp,Amat,RT,ART);/*   ART <-   A*RT       */
     92:     bbot = btop;
     93:     VecDot(RT,ART,&btop);
     94:     if (PetscAbsScalar(btop) < 0.0) {
     95:       ksp->reason = KSP_DIVERGED_INDEFINITE_MAT;
     96:       PetscLogInfo((ksp,"KSPSolve_CR:diverging due to indefinite or negative definite matrix\n"));
     97:       break;
     98:     }

    100:     if (ksp->normtype == KSP_PRECONDITIONED_NORM) {
    101:       VecNorm(RT,NORM_2,&dp);      /*   dp <- || RT ||      */
    102:     } else if (ksp->normtype == KSP_NATURAL_NORM) {
    103:       dp = sqrt(PetscAbsScalar(btop));                  /* dp = sqrt(R,AR)       */
    104:     } else if (ksp->normtype == KSP_NO_NORM) {
    105:       dp = 0.0;
    106:     } else if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) {
    107:       VecAXPY(R,ai,AP);           /*   R   <- R - ai*AP    */
    108:       VecNorm(R,NORM_2,&dp);       /*   dp <- R'*R          */
    109:     } else {
    110:       SETERRQ1(PETSC_ERR_SUP,"KSPNormType of %d not supported",(int)ksp->normtype);
    111:     }

    113:     PetscObjectTakeAccess(ksp);
    114:     ksp->its++;
    115:     ksp->rnorm = dp;
    116:     PetscObjectGrantAccess(ksp);

    118:     KSPLogResidualHistory(ksp,dp);
    119:     KSPMonitor(ksp,i+1,dp);
    120:     (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);
    121:     if (ksp->reason) break;

    123:     bi = btop/bbot;
    124:     VecAYPX(P,bi,RT);              /*   P <- RT + Bi P     */
    125:     VecAYPX(AP,bi,ART);            /*   AP <- ART + Bi AP  */
    126:     i++;
    127:   } while (i<ksp->max_it);
    128:   if (i >= ksp->max_it) {
    129:     ksp->reason =  KSP_DIVERGED_ITS;
    130:   }
    131:   return(0);
    132: }

    135: /*MC
    136:      KSPCR - This code implements the (preconditioned) conjugate residuals method

    138:    Options Database Keys:
    139: .   see KSPSolve()

    141:    Level: beginner

    143:    Notes: The operator and the preconditioner must be symmetric for this method. The
    144:           preconditioner must be POSITIVE-DEFINITE and the operator POSITIVE-SEMIDEFINITE

    146: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPCG
    147: M*/
    151: PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_CR(KSP ksp)
    152: {
    154:   ksp->pc_side                   = PC_LEFT;
    155:   ksp->ops->setup                = KSPSetUp_CR;
    156:   ksp->ops->solve                = KSPSolve_CR;
    157:   ksp->ops->destroy              = KSPDefaultDestroy;
    158:   ksp->ops->buildsolution        = KSPDefaultBuildSolution;
    159:   ksp->ops->buildresidual        = KSPDefaultBuildResidual;
    160:   ksp->ops->setfromoptions       = 0;
    161:   ksp->ops->view                 = 0;
    162:   return(0);
    163: }