Actual source code: fbcgs.c


  2: /*
  3:     This file implements flexible BiCGStab (FBiCGStab).
  4:     Only allow right preconditioning.
  5: */
  6: #include <../src/ksp/ksp/impls/bcgs/bcgsimpl.h>

  8: static PetscErrorCode KSPSetUp_FBCGS(KSP ksp)
  9: {
 10:   KSPSetWorkVecs(ksp,8);
 11:   return 0;
 12: }

 14: /* Only need a few hacks from KSPSolve_BCGS */

 16: static PetscErrorCode  KSPSolve_FBCGS(KSP ksp)
 17: {
 18:   PetscInt       i;
 19:   PetscScalar    rho,rhoold,alpha,beta,omega,omegaold,d1;
 20:   Vec            X,B,V,P,R,RP,T,S,P2,S2;
 21:   PetscReal      dp    = 0.0,d2;
 22:   KSP_BCGS       *bcgs = (KSP_BCGS*)ksp->data;
 23:   PC             pc;
 24:   Mat            mat;

 26:   X  = ksp->vec_sol;
 27:   B  = ksp->vec_rhs;
 28:   R  = ksp->work[0];
 29:   RP = ksp->work[1];
 30:   V  = ksp->work[2];
 31:   T  = ksp->work[3];
 32:   S  = ksp->work[4];
 33:   P  = ksp->work[5];
 34:   S2 = ksp->work[6];
 35:   P2 = ksp->work[7];

 37:   /* Only supports right preconditioning */
 39:   if (!ksp->guess_zero) {
 40:     if (!bcgs->guess) {
 41:       VecDuplicate(X,&bcgs->guess);
 42:     }
 43:     VecCopy(X,bcgs->guess);
 44:   } else {
 45:     VecSet(X,0.0);
 46:   }

 48:   /* Compute initial residual */
 49:   KSPGetPC(ksp,&pc);
 50:   PCSetUp(pc);
 51:   PCGetOperators(pc,&mat,NULL);
 52:   if (!ksp->guess_zero) {
 53:     KSP_MatMult(ksp,mat,X,S2);
 54:     VecCopy(B,R);
 55:     VecAXPY(R,-1.0,S2);
 56:   } else {
 57:     VecCopy(B,R);
 58:   }

 60:   /* Test for nothing to do */
 61:   if (ksp->normtype != KSP_NORM_NONE) {
 62:     VecNorm(R,NORM_2,&dp);
 63:   }
 64:   PetscObjectSAWsTakeAccess((PetscObject)ksp);
 65:   ksp->its   = 0;
 66:   ksp->rnorm = dp;
 67:   PetscObjectSAWsGrantAccess((PetscObject)ksp);
 68:   KSPLogResidualHistory(ksp,dp);
 69:   KSPMonitor(ksp,0,dp);
 70:   (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);
 71:   if (ksp->reason) return 0;

 73:   /* Make the initial Rp == R */
 74:   VecCopy(R,RP);

 76:   rhoold   = 1.0;
 77:   alpha    = 1.0;
 78:   omegaold = 1.0;
 79:   VecSet(P,0.0);
 80:   VecSet(V,0.0);

 82:   i=0;
 83:   do {
 84:     VecDot(R,RP,&rho); /* rho <- (r,rp) */
 85:     beta = (rho/rhoold) * (alpha/omegaold);
 86:     VecAXPBYPCZ(P,1.0,-omegaold*beta,beta,R,V); /* p <- r - omega * beta* v + beta * p */

 88:     KSP_PCApply(ksp,P,P2); /* p2 <- K p */
 89:     KSP_MatMult(ksp,mat,P2,V); /* v <- A p2 */

 91:     VecDot(V,RP,&d1);
 92:     if (d1 == 0.0) {
 94:       else ksp->reason = KSP_DIVERGED_BREAKDOWN;
 95:       PetscInfo(ksp,"Breakdown due to zero inner product\n");
 96:       break;
 97:     }
 98:     alpha = rho / d1; /* alpha <- rho / (v,rp) */
 99:     VecWAXPY(S,-alpha,V,R); /* s <- r - alpha v */

101:     KSP_PCApply(ksp,S,S2); /* s2 <- K s */
102:     KSP_MatMult(ksp,mat,S2,T); /* t <- A s2 */

104:     VecDotNorm2(S,T,&d1,&d2);
105:     if (d2 == 0.0) {
106:       /* t is 0. if s is 0, then alpha v == r, and hence alpha p may be our solution. Give it a try? */
107:       VecDot(S,S,&d1);
108:       if (d1 != 0.0) {
110:         else ksp->reason = KSP_DIVERGED_BREAKDOWN;
111:         PetscInfo(ksp,"Failed due to singular preconditioned operator\n");
112:         break;
113:       }
114:       VecAXPY(X,alpha,P2);   /* x <- x + alpha p2 */
115:       PetscObjectSAWsTakeAccess((PetscObject)ksp);
116:       ksp->its++;
117:       ksp->rnorm  = 0.0;
118:       ksp->reason = KSP_CONVERGED_RTOL;
119:       PetscObjectSAWsGrantAccess((PetscObject)ksp);
120:       KSPLogResidualHistory(ksp,dp);
121:       KSPMonitor(ksp,i+1,0.0);
122:       break;
123:     }
124:     omega = d1 / d2; /* omega <- (t's) / (t't) */
125:     VecAXPBYPCZ(X,alpha,omega,1.0,P2,S2); /* x <- alpha * p2 + omega * s2 + x */

127:     VecWAXPY(R,-omega,T,S);  /* r <- s - omega t */
128:     if (ksp->normtype != KSP_NORM_NONE && ksp->chknorm < i+2) {
129:       VecNorm(R,NORM_2,&dp);
130:     }

132:     rhoold   = rho;
133:     omegaold = omega;

135:     PetscObjectSAWsTakeAccess((PetscObject)ksp);
136:     ksp->its++;
137:     ksp->rnorm = dp;
138:     PetscObjectSAWsGrantAccess((PetscObject)ksp);
139:     KSPLogResidualHistory(ksp,dp);
140:     KSPMonitor(ksp,i+1,dp);
141:     (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);
142:     if (ksp->reason) break;
143:     if (rho == 0.0) {
145:       else ksp->reason = KSP_DIVERGED_BREAKDOWN;
146:       PetscInfo(ksp,"Breakdown due to zero rho inner product\n");
147:       break;
148:     }
149:     i++;
150:   } while (i<ksp->max_it);

152:   if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
153:   return 0;
154: }

156: /*MC
157:      KSPFBCGS - Implements flexible BiCGStab method.

159:    Options Database Keys:
160:     see KSPSolve()

162:    Level: beginner

164:    Notes:
165:     Only allow right preconditioning

167: .seealso:  KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPBICG, KSPFBCGSL, KSPSetPCSide()
168: M*/
169: PETSC_EXTERN PetscErrorCode KSPCreate_FBCGS(KSP ksp)
170: {
171:   KSP_BCGS       *bcgs;

173:   PetscNewLog(ksp,&bcgs);

175:   ksp->data                = bcgs;
176:   ksp->ops->setup          = KSPSetUp_FBCGS;
177:   ksp->ops->solve          = KSPSolve_FBCGS;
178:   ksp->ops->destroy        = KSPDestroy_BCGS;
179:   ksp->ops->reset          = KSPReset_BCGS;
180:   ksp->ops->buildresidual  = KSPBuildResidualDefault;
181:   ksp->ops->setfromoptions = KSPSetFromOptions_BCGS;
182:   ksp->pc_side             = PC_RIGHT;  /* set default PC side */

184:   KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
185:   KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
186:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
187:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
188:   return 0;
189: }