Actual source code: fbcgs.c
petsc-3.4.5 2014-06-29
2: /*
3: This file implements flexible BiCGStab (FBiCGStab).
4: Only allow right preconditioning.
5: */
6: #include <../src/ksp/ksp/impls/bcgs/bcgsimpl.h> /*I "petscksp.h" I*/
10: static PetscErrorCode KSPSetUp_FBCGS(KSP ksp)
11: {
15: KSPSetWorkVecs(ksp,8);
16: return(0);
17: }
19: /* Only need a few hacks from KSPSolve_BCGS */
20: #include <petsc-private/pcimpl.h> /*I "petscksp.h" I*/
23: static PetscErrorCode KSPSolve_FBCGS(KSP ksp)
24: {
26: PetscInt i;
27: PetscScalar rho,rhoold,alpha,beta,omega,omegaold,d1;
28: Vec X,B,V,P,R,RP,T,S,P2,S2;
29: PetscReal dp = 0.0,d2;
30: KSP_BCGS *bcgs = (KSP_BCGS*)ksp->data;
31: PC pc;
34: X = ksp->vec_sol;
35: B = ksp->vec_rhs;
36: R = ksp->work[0];
37: RP = ksp->work[1];
38: V = ksp->work[2];
39: T = ksp->work[3];
40: S = ksp->work[4];
41: P = ksp->work[5];
42: S2 = ksp->work[6];
43: P2 = ksp->work[7];
45: /* Only supports right preconditioning */
46: if (ksp->pc_side != PC_RIGHT) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"KSP fbcgs does not support %s",PCSides[ksp->pc_side]);
47: if (!ksp->guess_zero) {
48: if (!bcgs->guess) {
49: VecDuplicate(X,&bcgs->guess);
50: }
51: VecCopy(X,bcgs->guess);
52: } else {
53: VecSet(X,0.0);
54: }
56: /* Compute initial residual */
57: KSPGetPC(ksp,&pc);
58: PCSetUp(pc);
59: if (!ksp->guess_zero) {
60: MatMult(pc->mat,X,S2);
61: VecCopy(B,R);
62: VecAXPY(R,-1.0,S2);
63: } else {
64: VecCopy(B,R);
65: }
67: /* Test for nothing to do */
68: if (ksp->normtype != KSP_NORM_NONE) {
69: VecNorm(R,NORM_2,&dp);
70: }
71: PetscObjectAMSTakeAccess((PetscObject)ksp);
72: ksp->its = 0;
73: ksp->rnorm = dp;
74: PetscObjectAMSGrantAccess((PetscObject)ksp);
75: KSPLogResidualHistory(ksp,dp);
76: KSPMonitor(ksp,0,dp);
77: (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);
78: if (ksp->reason) return(0);
80: /* Make the initial Rp == R */
81: VecCopy(R,RP);
83: rhoold = 1.0;
84: alpha = 1.0;
85: omegaold = 1.0;
86: VecSet(P,0.0);
87: VecSet(V,0.0);
89: i=0;
90: do {
91: VecDot(R,RP,&rho); /* rho <- (r,rp) */
92: beta = (rho/rhoold) * (alpha/omegaold);
93: VecAXPBYPCZ(P,1.0,-omegaold*beta,beta,R,V); /* p <- r - omega * beta* v + beta * p */
95: PCApply(pc,P,P2); /* p2 <- K p */
96: MatMult(pc->mat,P2,V); /* v <- A p2 */
98: VecDot(V,RP,&d1);
99: if (d1 == 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Divide by zero");
100: alpha = rho / d1; /* alpha <- rho / (v,rp) */
101: VecWAXPY(S,-alpha,V,R); /* s <- r - alpha v */
103: PCApply(pc,S,S2); /* s2 <- K s */
104: MatMult(pc->mat,S2,T); /* t <- A s2 */
106: VecDotNorm2(S,T,&d1,&d2);
107: if (d2 == 0.0) {
108: /* t is 0. if s is 0, then alpha v == r, and hence alpha p may be our solution. Give it a try? */
109: VecDot(S,S,&d1);
110: if (d1 != 0.0) {
111: ksp->reason = KSP_DIVERGED_BREAKDOWN;
112: break;
113: }
114: VecAXPY(X,alpha,P2); /* x <- x + alpha p2 */
115: PetscObjectAMSTakeAccess((PetscObject)ksp);
116: ksp->its++;
117: ksp->rnorm = 0.0;
118: ksp->reason = KSP_CONVERGED_RTOL;
119: PetscObjectAMSGrantAccess((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: PetscObjectAMSTakeAccess((PetscObject)ksp);
136: ksp->its++;
137: ksp->rnorm = dp;
138: PetscObjectAMSGrantAccess((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) {
144: ksp->reason = KSP_DIVERGED_BREAKDOWN;
145: break;
146: }
147: i++;
148: } while (i<ksp->max_it);
150: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
151: return(0);
152: }
154: /*MC
155: KSPFBCGS - Implements flexible BiCGStab method.
157: Options Database Keys:
158: . see KSPSolve()
160: Level: beginner
162: Notes: Only allow right preconditioning
164: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPBICG, KSPFBCGSL, KSPSetPCSide()
165: M*/
168: PETSC_EXTERN PetscErrorCode KSPCreate_FBCGS(KSP ksp)
169: {
171: KSP_BCGS *bcgs;
174: PetscNewLog(ksp,KSP_BCGS,&bcgs);
176: ksp->data = bcgs;
177: ksp->ops->setup = KSPSetUp_FBCGS;
178: ksp->ops->solve = KSPSolve_FBCGS;
179: ksp->ops->destroy = KSPDestroy_BCGS;
180: ksp->ops->reset = KSPReset_BCGS;
181: ksp->ops->buildresidual = KSPBuildResidualDefault;
182: ksp->ops->setfromoptions = KSPSetFromOptions_BCGS;
183: ksp->pc_side = PC_RIGHT; /* set default PC side */
185: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
186: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,1);
187: return(0);
188: }