Actual source code: bicg.c

      1: #define PETSCKSP_DLL

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

      7: PetscErrorCode KSPSetUp_BiCG(KSP ksp)
      8: {

     12:   /* check user parameters and functions */
     13:   if (ksp->pc_side == PC_RIGHT) {
     14:     SETERRQ(PETSC_ERR_SUP,"no right preconditioning for KSPBiCG");
     15:   } else if (ksp->pc_side == PC_SYMMETRIC) {
     16:     SETERRQ(PETSC_ERR_SUP,"no symmetric preconditioning for KSPBiCG");
     17:   }

     19:   /* get work vectors from user code */
     20:   KSPDefaultGetWork(ksp,6);
     21:   return(0);
     22: }

     26: PetscErrorCode  KSPSolve_BiCG(KSP ksp)
     27: {
     29:   PetscInt       i;
     30:   PetscTruth     diagonalscale;
     31:   PetscScalar    dpi,a=1.0,beta,betaold=1.0,b,mone=-1.0,ma;
     32:   PetscReal      dp;
     33:   Vec            X,B,Zl,Zr,Rl,Rr,Pl,Pr;
     34:   Mat            Amat,Pmat;
     35:   MatStructure   pflag;

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

     41:   X       = ksp->vec_sol;
     42:   B       = ksp->vec_rhs;
     43:   Rl      = ksp->work[0];
     44:   Zl      = ksp->work[1];
     45:   Pl      = ksp->work[2];
     46:   Rr      = ksp->work[3];
     47:   Zr      = ksp->work[4];
     48:   Pr      = ksp->work[5];

     50:   PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);

     52:   if (!ksp->guess_zero) {
     53:     KSP_MatMult(ksp,Amat,X,Rr);      /*   r <- b - Ax       */
     54:     VecAYPX(Rr,mone,B);
     55:   } else {
     56:     VecCopy(B,Rr);           /*     r <- b (x is 0) */
     57:   }
     58:   VecCopy(Rr,Rl);
     59:   KSP_PCApply(ksp,Rr,Zr);     /*     z <- Br         */
     60:   VecConjugate(Rl);
     61:   KSP_PCApplyTranspose(ksp,Rl,Zl);
     62:   VecConjugate(Rl);
     63:   VecConjugate(Zl);
     64:   if (ksp->normtype == KSP_PRECONDITIONED_NORM) {
     65:     VecNorm(Zr,NORM_2,&dp);  /*    dp <- z'*z       */
     66:   } else {
     67:     VecNorm(Rr,NORM_2,&dp);  /*    dp <- r'*r       */
     68:   }
     69:   KSPMonitor(ksp,0,dp);
     70:   PetscObjectTakeAccess(ksp);
     71:   ksp->its   = 0;
     72:   ksp->rnorm = dp;
     73:   PetscObjectGrantAccess(ksp);
     74:   KSPLogResidualHistory(ksp,dp);
     75:   (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);
     76:   if (ksp->reason) return(0);

     78:   i = 0;
     79:   do {
     80:      VecDot(Zr,Rl,&beta);       /*     beta <- r'z     */
     81:      if (!i) {
     82:        if (beta == 0.0) {
     83:          ksp->reason = KSP_DIVERGED_BREAKDOWN_BICG;
     84:          return(0);
     85:        }
     86:        VecCopy(Zr,Pr);       /*     p <- z          */
     87:        VecCopy(Zl,Pl);
     88:      } else {
     89:        b = beta/betaold;
     90:        VecAYPX(Pr,b,Zr);  /*     p <- z + b* p   */
     91:        b = PetscConj(b);
     92:        VecAYPX(Pl,b,Zl);
     93:      }
     94:      betaold = beta;
     95:      KSP_MatMult(ksp,Amat,Pr,Zr);    /*     z <- Kp         */
     96:      VecConjugate(Pl);
     97:      KSP_MatMultTranspose(ksp,Amat,Pl,Zl);
     98:      VecConjugate(Pl);
     99:      VecConjugate(Zl);
    100:      VecDot(Zr,Pl,&dpi);               /*     dpi <- z'p      */
    101:      a = beta/dpi;                                 /*     a = beta/p'z    */
    102:      VecAXPY(X,a,Pr);       /*     x <- x + ap     */
    103:      ma = -a;
    104:      VecAXPY(Rr,ma,Zr);CHKERRQ(ierr)
    105:      ma = PetscConj(ma);
    106:      VecAXPY(Rl,ma,Zl);
    107:      if (ksp->normtype == KSP_PRECONDITIONED_NORM) {
    108:        KSP_PCApply(ksp,Rr,Zr);  /*     z <- Br         */
    109:        VecConjugate(Rl);
    110:        KSP_PCApplyTranspose(ksp,Rl,Zl);
    111:        VecConjugate(Rl);
    112:        VecConjugate(Zl);
    113:        VecNorm(Zr,NORM_2,&dp);  /*    dp <- z'*z       */
    114:      } else {
    115:        VecNorm(Rr,NORM_2,&dp);  /*    dp <- r'*r       */
    116:      }
    117:      PetscObjectTakeAccess(ksp);
    118:      ksp->its   = i+1;
    119:      ksp->rnorm = dp;
    120:      PetscObjectGrantAccess(ksp);
    121:      KSPLogResidualHistory(ksp,dp);
    122:      KSPMonitor(ksp,i+1,dp);
    123:      (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);
    124:      if (ksp->reason) break;
    125:      if (ksp->normtype == KSP_UNPRECONDITIONED_NORM) {
    126:        KSP_PCApply(ksp,Rr,Zr);  /* z <- Br  */
    127:        VecConjugate(Rl);
    128:        KSP_PCApplyTranspose(ksp,Rl,Zl);
    129:        VecConjugate(Rl);
    130:        VecConjugate(Zl);
    131:      }
    132:      i++;
    133:   } while (i<ksp->max_it);
    134:   if (i >= ksp->max_it) {
    135:     ksp->reason = KSP_DIVERGED_ITS;
    136:   }
    137:   return(0);
    138: }

    142: PetscErrorCode KSPDestroy_BiCG(KSP ksp)
    143: {

    147:   KSPDefaultFreeWork(ksp);
    148:   return(0);
    149: }

    151: /*MC
    152:      KSPBICG - Implements the Biconjugate gradient method (essentially running the conjugate
    153:          gradient on the normal equations).

    155:    Options Database Keys:
    156: .   see KSPSolve()

    158:    Level: beginner

    160: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPBCGS

    162: M*/
    166: PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_BiCG(KSP ksp)
    167: {
    169:   ksp->data                      = (void*)0;
    170:   ksp->pc_side                   = PC_LEFT;
    171:   ksp->ops->setup                = KSPSetUp_BiCG;
    172:   ksp->ops->solve                = KSPSolve_BiCG;
    173:   ksp->ops->destroy              = KSPDestroy_BiCG;
    174:   ksp->ops->view                 = 0;
    175:   ksp->ops->setfromoptions       = 0;
    176:   ksp->ops->buildsolution        = KSPDefaultBuildSolution;
    177:   ksp->ops->buildresidual        = KSPDefaultBuildResidual;

    179:   return(0);
    180: }