Actual source code: bcgs.c
2: #include <../src/ksp/ksp/impls/bcgs/bcgsimpl.h>
4: PetscErrorCode KSPSetFromOptions_BCGS(PetscOptionItems *PetscOptionsObject,KSP ksp)
5: {
6: PetscOptionsHead(PetscOptionsObject,"KSP BCGS Options");
7: PetscOptionsTail();
8: return 0;
9: }
11: PetscErrorCode KSPSetUp_BCGS(KSP ksp)
12: {
13: KSPSetWorkVecs(ksp,6);
14: return 0;
15: }
17: PetscErrorCode KSPSolve_BCGS(KSP ksp)
18: {
19: PetscInt i;
20: PetscScalar rho,rhoold,alpha,beta,omega,omegaold,d1;
21: Vec X,B,V,P,R,RP,T,S;
22: PetscReal dp = 0.0,d2;
23: KSP_BCGS *bcgs = (KSP_BCGS*)ksp->data;
25: X = ksp->vec_sol;
26: B = ksp->vec_rhs;
27: R = ksp->work[0];
28: RP = ksp->work[1];
29: V = ksp->work[2];
30: T = ksp->work[3];
31: S = ksp->work[4];
32: P = ksp->work[5];
34: /* Compute initial preconditioned residual */
35: KSPInitialResidual(ksp,X,V,T,R,B);
37: /* with right preconditioning need to save initial guess to add to final solution */
38: if (ksp->pc_side == PC_RIGHT && !ksp->guess_zero) {
39: if (!bcgs->guess) {
40: VecDuplicate(X,&bcgs->guess);
41: }
42: VecCopy(X,bcgs->guess);
43: VecSet(X,0.0);
44: }
46: /* Test for nothing to do */
47: if (ksp->normtype != KSP_NORM_NONE) {
48: VecNorm(R,NORM_2,&dp);
49: KSPCheckNorm(ksp,dp);
50: }
51: PetscObjectSAWsTakeAccess((PetscObject)ksp);
52: ksp->its = 0;
53: ksp->rnorm = dp;
54: PetscObjectSAWsGrantAccess((PetscObject)ksp);
55: KSPLogResidualHistory(ksp,dp);
56: KSPMonitor(ksp,0,dp);
57: (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);
58: if (ksp->reason) {
59: if (bcgs->guess) {
60: VecAXPY(X,1.0,bcgs->guess);
61: }
62: return 0;
63: }
65: /* Make the initial Rp == R */
66: VecCopy(R,RP);
68: rhoold = 1.0;
69: alpha = 1.0;
70: omegaold = 1.0;
71: VecSet(P,0.0);
72: VecSet(V,0.0);
74: i = 0;
75: do {
76: VecDot(R,RP,&rho); /* rho <- (r,rp) */
77: beta = (rho/rhoold) * (alpha/omegaold);
78: VecAXPBYPCZ(P,1.0,-omegaold*beta,beta,R,V); /* p <- r - omega * beta* v + beta * p */
79: KSP_PCApplyBAorAB(ksp,P,V,T); /* v <- K p */
80: VecDot(V,RP,&d1);
81: KSPCheckDot(ksp,d1);
82: if (d1 == 0.0) {
84: else ksp->reason = KSP_DIVERGED_BREAKDOWN;
85: PetscInfo(ksp,"Breakdown due to zero inner product\n");
86: break;
87: }
88: alpha = rho / d1; /* a <- rho / (v,rp) */
89: VecWAXPY(S,-alpha,V,R); /* s <- r - a v */
90: KSP_PCApplyBAorAB(ksp,S,T,R); /* t <- K s */
91: VecDotNorm2(S,T,&d1,&d2);
92: if (d2 == 0.0) {
93: /* t is 0. if s is 0, then alpha v == r, and hence alpha p
94: may be our solution. Give it a try? */
95: VecDot(S,S,&d1);
96: if (d1 != 0.0) {
98: else ksp->reason = KSP_DIVERGED_BREAKDOWN;
99: PetscInfo(ksp,"Failed due to singular preconditioned operator\n");
100: break;
101: }
102: VecAXPY(X,alpha,P); /* x <- x + a p */
103: PetscObjectSAWsTakeAccess((PetscObject)ksp);
104: ksp->its++;
105: ksp->rnorm = 0.0;
106: ksp->reason = KSP_CONVERGED_RTOL;
107: PetscObjectSAWsGrantAccess((PetscObject)ksp);
108: KSPLogResidualHistory(ksp,dp);
109: KSPMonitor(ksp,i+1,0.0);
110: break;
111: }
112: omega = d1 / d2; /* w <- (t's) / (t't) */
113: VecAXPBYPCZ(X,alpha,omega,1.0,P,S); /* x <- alpha * p + omega * s + x */
114: VecWAXPY(R,-omega,T,S); /* r <- s - w t */
115: if (ksp->normtype != KSP_NORM_NONE && ksp->chknorm < i+2) {
116: VecNorm(R,NORM_2,&dp);
117: KSPCheckNorm(ksp,dp);
118: }
120: rhoold = rho;
121: omegaold = omega;
123: PetscObjectSAWsTakeAccess((PetscObject)ksp);
124: ksp->its++;
125: ksp->rnorm = dp;
126: PetscObjectSAWsGrantAccess((PetscObject)ksp);
127: KSPLogResidualHistory(ksp,dp);
128: KSPMonitor(ksp,i+1,dp);
129: (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);
130: if (ksp->reason) break;
131: if (rho == 0.0) {
133: else ksp->reason = KSP_DIVERGED_BREAKDOWN;
134: PetscInfo(ksp,"Breakdown due to zero rho inner product\n");
135: break;
136: }
137: i++;
138: } while (i<ksp->max_it);
140: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
142: KSPUnwindPreconditioner(ksp,X,T);
143: if (bcgs->guess) {
144: VecAXPY(X,1.0,bcgs->guess);
145: }
146: return 0;
147: }
149: PetscErrorCode KSPBuildSolution_BCGS(KSP ksp,Vec v,Vec *V)
150: {
151: KSP_BCGS *bcgs = (KSP_BCGS*)ksp->data;
153: if (ksp->pc_side == PC_RIGHT) {
154: if (v) {
155: KSP_PCApply(ksp,ksp->vec_sol,v);
156: if (bcgs->guess) {
157: VecAXPY(v,1.0,bcgs->guess);
158: }
159: *V = v;
160: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with right preconditioner");
161: } else {
162: if (v) {
163: VecCopy(ksp->vec_sol,v); *V = v;
164: } else *V = ksp->vec_sol;
165: }
166: return 0;
167: }
169: PetscErrorCode KSPReset_BCGS(KSP ksp)
170: {
171: KSP_BCGS *cg = (KSP_BCGS*)ksp->data;
173: VecDestroy(&cg->guess);
174: return 0;
175: }
177: PetscErrorCode KSPDestroy_BCGS(KSP ksp)
178: {
179: KSPReset_BCGS(ksp);
180: KSPDestroyDefault(ksp);
181: return 0;
182: }
184: /*MC
185: KSPBCGS - Implements the BiCGStab (Stabilized version of BiConjugate Gradient) method.
187: Options Database Keys:
188: see KSPSolve()
190: Level: beginner
192: Notes:
193: See KSPBCGSL for additional stabilization
194: Supports left and right preconditioning but not symmetric
196: References:
197: . * - van der Vorst, SIAM J. Sci. Stat. Comput., 1992.
199: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPBICG, KSPBCGSL, KSPFBICG, KSPQMRCGS, KSPSetPCSide()
200: M*/
201: PETSC_EXTERN PetscErrorCode KSPCreate_BCGS(KSP ksp)
202: {
203: KSP_BCGS *bcgs;
205: PetscNewLog(ksp,&bcgs);
207: ksp->data = bcgs;
208: ksp->ops->setup = KSPSetUp_BCGS;
209: ksp->ops->solve = KSPSolve_BCGS;
210: ksp->ops->destroy = KSPDestroy_BCGS;
211: ksp->ops->reset = KSPReset_BCGS;
212: ksp->ops->buildsolution = KSPBuildSolution_BCGS;
213: ksp->ops->buildresidual = KSPBuildResidualDefault;
214: ksp->ops->setfromoptions = KSPSetFromOptions_BCGS;
216: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
217: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
218: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
219: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
220: return 0;
221: }