Actual source code: bcgs.c

      1: #define PETSCKSP_DLL

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

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

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

     21: static PetscErrorCode  KSPSolve_BCGS(KSP ksp)
     22: {
     24:   PetscInt       i;
     25:   PetscScalar    rho,rhoold,alpha,beta,omega,omegaold,d1,d2,zero = 0.0,tmp;
     26:   Vec            X,B,V,P,R,RP,T,S;
     27:   PetscReal      dp = 0.0;

     31:   X       = ksp->vec_sol;
     32:   B       = ksp->vec_rhs;
     33:   R       = ksp->work[0];
     34:   RP      = ksp->work[1];
     35:   V       = ksp->work[2];
     36:   T       = ksp->work[3];
     37:   S       = ksp->work[4];
     38:   P       = ksp->work[5];

     40:   /* Compute initial preconditioned residual */
     41:   KSPInitialResidual(ksp,X,V,T,R,B);

     43:   /* Test for nothing to do */
     44:   if (ksp->normtype != KSP_NO_NORM) {
     45:     VecNorm(R,NORM_2,&dp);
     46:   }
     47:   PetscObjectTakeAccess(ksp);
     48:   ksp->its   = 0;
     49:   ksp->rnorm = dp;
     50:   PetscObjectGrantAccess(ksp);
     51:   KSPLogResidualHistory(ksp,dp);
     52:   KSPMonitor(ksp,0,dp);
     53:   (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);
     54:   if (ksp->reason) return(0);

     56:   /* Make the initial Rp == R */
     57:   VecCopy(R,RP);

     59:   rhoold   = 1.0;
     60:   alpha    = 1.0;
     61:   omegaold = 1.0;
     62:   VecSet(P,zero);
     63:   VecSet(V,zero);

     65:   i=0;
     66:   do {
     67:     VecDot(R,RP,&rho);       /*   rho <- (r,rp)      */
     68:     if (rho == 0.0) {
     69:       ksp->reason = KSP_DIVERGED_BREAKDOWN;
     70:       break;
     71:     }
     72:     beta = (rho/rhoold) * (alpha/omegaold);
     73:     tmp = -omegaold; VecAXPY(P,tmp,V);            /*   p <- p - w v       */
     74:     VecAYPX(P,beta,R);      /*   p <- r + p beta    */
     75:     KSP_PCApplyBAorAB(ksp,P,V,T);  /*   v <- K p           */
     76:     VecDot(V,RP,&d1);
     77:     alpha = rho / d1; tmp = -alpha;                /*   a <- rho / (v,rp)  */
     78:     VecWAXPY(S,tmp,V,R);      /*   s <- r - a v       */
     79:     KSP_PCApplyBAorAB(ksp,S,T,R);/*   t <- K s    */
     80:     VecDot(S,T,&d1);
     81:     VecDot(T,T,&d2);
     82:     if (d2 == 0.0) {
     83:       /* t is 0.  if s is 0, then alpha v == r, and hence alpha p
     84:          may be our solution.  Give it a try? */
     85:       VecDot(S,S,&d1);
     86:       if (d1 == 0.0) {
     87:         ksp->reason = KSP_DIVERGED_BREAKDOWN;
     88:         break;
     89:       }
     90:       VecAXPY(X,alpha,P);   /*   x <- x + a p       */
     91:       PetscObjectTakeAccess(ksp);
     92:       ksp->its++;
     93:       ksp->rnorm  = 0.0;
     94:       ksp->reason = KSP_CONVERGED_RTOL;
     95:       PetscObjectGrantAccess(ksp);
     96:       KSPLogResidualHistory(ksp,dp);
     97:       KSPMonitor(ksp,i+1,0.0);
     98:       break;
     99:     }
    100:     omega = d1 / d2;                               /*   w <- (t's) / (t't) */
    101:     VecAXPY(X,alpha,P);     /*   x <- x + a p       */
    102:     VecAXPY(X,omega,S);     /*   x <- x + w s       */
    103:     tmp   = -omega;
    104:     VecWAXPY(R,tmp,T,S);     /*   r <- s - w t       */
    105:     if (ksp->normtype != KSP_NO_NORM) {
    106:       VecNorm(R,NORM_2,&dp);
    107:     }

    109:     rhoold   = rho;
    110:     omegaold = omega;

    112:     PetscObjectTakeAccess(ksp);
    113:     ksp->its++;
    114:     ksp->rnorm = dp;
    115:     PetscObjectGrantAccess(ksp);
    116:     KSPLogResidualHistory(ksp,dp);
    117:     KSPMonitor(ksp,i+1,dp);
    118:     (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);
    119:     if (ksp->reason) break;
    120:     i++;
    121:   } while (i<ksp->max_it);

    123:   if (i >= ksp->max_it) {
    124:     ksp->reason = KSP_DIVERGED_ITS;
    125:   }

    127:   KSPUnwindPreconditioner(ksp,X,T);
    128:   return(0);
    129: }

    131: /*MC
    132:      KSPBCGS - Implements the BiCGStab (Stabilized version of BiConjugate Gradient Squared) method.

    134:    Options Database Keys:
    135: .   see KSPSolve()

    137:    Level: beginner

    139:    Notes: Reference: van der Vorst, SIAM J. Sci. Stat. Comput., 1992.
    140:           See KSPBCGSL for additional stabilization

    142: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPBICG, KSPBCGSL
    143: M*/
    147: PetscErrorCode PETSCKSP_DLLEXPORT KSPCreate_BCGS(KSP ksp)
    148: {
    150:   ksp->data                 = (void*)0;
    151:   ksp->pc_side              = PC_LEFT;
    152:   ksp->ops->setup           = KSPSetUp_BCGS;
    153:   ksp->ops->solve           = KSPSolve_BCGS;
    154:   ksp->ops->destroy         = KSPDefaultDestroy;
    155:   ksp->ops->buildsolution   = KSPDefaultBuildSolution;
    156:   ksp->ops->buildresidual   = KSPDefaultBuildResidual;
    157:   ksp->ops->setfromoptions  = 0;
    158:   ksp->ops->view            = 0;
    159:   return(0);
    160: }