Actual source code: bcgs.c
petsc-3.3-p7 2013-05-11
2: #include <../src/ksp/ksp/impls/bcgs/bcgsimpl.h> /*I "petscksp.h" I*/
6: PetscErrorCode KSPSetFromOptions_BCGS(KSP ksp)
7: {
9:
11: PetscOptionsHead("KSP BCGS Options");
12: PetscOptionsTail();
13: return(0);
14: }
18: PetscErrorCode KSPView_BCGS(KSP ksp,PetscViewer viewer)
19: {
21: PetscBool iascii;
24: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
25: if (!iascii){
26: SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for KSP cg",((PetscObject)viewer)->type_name);
27: }
28: return(0);
29: }
33: PetscErrorCode KSPSetUp_BCGS(KSP ksp)
34: {
38: KSPDefaultGetWork(ksp,6);
39: return(0);
40: }
45: PetscErrorCode KSPSolve_BCGS(KSP ksp)
46: {
48: PetscInt i;
49: PetscScalar rho,rhoold,alpha,beta,omega,omegaold,d1,d2;
50: Vec X,B,V,P,R,RP,T,S;
51: PetscReal dp = 0.0;
52: KSP_BCGS *bcgs = (KSP_BCGS*)ksp->data;
55: X = ksp->vec_sol;
56: B = ksp->vec_rhs;
57: R = ksp->work[0];
58: RP = ksp->work[1];
59: V = ksp->work[2];
60: T = ksp->work[3];
61: S = ksp->work[4];
62: P = ksp->work[5];
64: /* Compute initial preconditioned residual */
65: KSPInitialResidual(ksp,X,V,T,R,B);
67: /* with right preconditioning need to save initial guess to add to final solution */
68: if (ksp->pc_side == PC_RIGHT && !ksp->guess_zero) {
69: if (!bcgs->guess) {
70: VecDuplicate(X,&bcgs->guess);
71: }
72: VecCopy(X,bcgs->guess);
73: VecSet(X,0.0);
74: }
76: /* Test for nothing to do */
77: if (ksp->normtype != KSP_NORM_NONE) {
78: VecNorm(R,NORM_2,&dp);
79: }
80: PetscObjectTakeAccess(ksp);
81: ksp->its = 0;
82: ksp->rnorm = dp;
83: PetscObjectGrantAccess(ksp);
84: KSPLogResidualHistory(ksp,dp);
85: KSPMonitor(ksp,0,dp);
86: (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);
87: if (ksp->reason) return(0);
89: /* Make the initial Rp == R */
90: VecCopy(R,RP);
92: rhoold = 1.0;
93: alpha = 1.0;
94: omegaold = 1.0;
95: VecSet(P,0.0);
96: VecSet(V,0.0);
98: i=0;
99: do {
100: VecDot(R,RP,&rho); /* rho <- (r,rp) */
101: beta = (rho/rhoold) * (alpha/omegaold);
102: VecAXPBYPCZ(P,1.0,-omegaold*beta,beta,R,V); /* p <- r - omega * beta* v + beta * p */
103: KSP_PCApplyBAorAB(ksp,P,V,T); /* v <- K p */
104: VecDot(V,RP,&d1);
105: if (d1 == 0.0) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_PLIB,"Divide by zero");
106: alpha = rho / d1; /* a <- rho / (v,rp) */
107: VecWAXPY(S,-alpha,V,R); /* s <- r - a v */
108: KSP_PCApplyBAorAB(ksp,S,T,R);/* t <- K s */
109: VecDotNorm2(S,T,&d1,&d2);
110: if (d2 == 0.0) {
111: /* t is 0. if s is 0, then alpha v == r, and hence alpha p
112: may be our solution. Give it a try? */
113: VecDot(S,S,&d1);
114: if (d1 != 0.0) {
115: ksp->reason = KSP_DIVERGED_BREAKDOWN;
116: break;
117: }
118: VecAXPY(X,alpha,P); /* x <- x + a p */
119: PetscObjectTakeAccess(ksp);
120: ksp->its++;
121: ksp->rnorm = 0.0;
122: ksp->reason = KSP_CONVERGED_RTOL;
123: PetscObjectGrantAccess(ksp);
124: KSPLogResidualHistory(ksp,dp);
125: KSPMonitor(ksp,i+1,0.0);
126: break;
127: }
128: omega = d1 / d2; /* w <- (t's) / (t't) */
129: VecAXPBYPCZ(X,alpha,omega,1.0,P,S); /* x <- alpha * p + omega * s + x */
130: VecWAXPY(R,-omega,T,S); /* r <- s - w t */
131: if (ksp->normtype != KSP_NORM_NONE && ksp->chknorm < i+2) {
132: VecNorm(R,NORM_2,&dp);
133: }
135: rhoold = rho;
136: omegaold = omega;
138: PetscObjectTakeAccess(ksp);
139: ksp->its++;
140: ksp->rnorm = dp;
141: PetscObjectGrantAccess(ksp);
142: KSPLogResidualHistory(ksp,dp);
143: KSPMonitor(ksp,i+1,dp);
144: (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);
145: if (ksp->reason) break;
146: if (rho == 0.0) {
147: ksp->reason = KSP_DIVERGED_BREAKDOWN;
148: break;
149: }
150: i++;
151: } while (i<ksp->max_it);
153: if (i >= ksp->max_it) {
154: ksp->reason = KSP_DIVERGED_ITS;
155: }
157: KSPUnwindPreconditioner(ksp,X,T);
158: if (bcgs->guess) {
159: VecAXPY(X,1.0,bcgs->guess);
160: }
161: return(0);
162: }
166: PetscErrorCode KSPBuildSolution_BCGS(KSP ksp,Vec v,Vec *V)
167: {
169: KSP_BCGS *bcgs = (KSP_BCGS*)ksp->data;
172: if (ksp->pc_side == PC_RIGHT) {
173: if (v) {
174: KSP_PCApply(ksp,ksp->vec_sol,v);
175: if (bcgs->guess) {
176: VecAXPY(v,1.0,bcgs->guess);
177: }
178: *V = v;
179: } else SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Not working with right preconditioner");
180: } else {
181: if (v) {VecCopy(ksp->vec_sol,v); *V = v;}
182: else { *V = ksp->vec_sol; }
183: }
184: return(0);
185: }
189: PetscErrorCode KSPReset_BCGS(KSP ksp)
190: {
191: KSP_BCGS *cg = (KSP_BCGS*)ksp->data;
195: VecDestroy(&cg->guess);
196: return(0);
197: }
201: PetscErrorCode KSPDestroy_BCGS(KSP ksp)
202: {
206: KSPReset_BCGS(ksp);
207: KSPDefaultDestroy(ksp);
208: return(0);
209: }
211: /*MC
212: KSPBCGS - Implements the BiCGStab (Stabilized version of BiConjugate Gradient Squared) method.
214: Options Database Keys:
215: . see KSPSolve()
217: Level: beginner
219: Notes: See KSPBCGSL for additional stabilization
220: Supports left and right preconditioning but not symmetric
222: References: van der Vorst, SIAM J. Sci. Stat. Comput., 1992.
224: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPBICG, KSPBCGSL, KSPFBICG, KSPSetPCSide()
225: M*/
226: EXTERN_C_BEGIN
229: PetscErrorCode KSPCreate_BCGS(KSP ksp)
230: {
232: KSP_BCGS *bcgs;
235: PetscNewLog(ksp,KSP_BCGS,&bcgs);
236: ksp->data = bcgs;
237: ksp->ops->setup = KSPSetUp_BCGS;
238: ksp->ops->solve = KSPSolve_BCGS;
239: ksp->ops->destroy = KSPDestroy_BCGS;
240: ksp->ops->reset = KSPReset_BCGS;
241: ksp->ops->buildsolution = KSPBuildSolution_BCGS;
242: ksp->ops->buildresidual = KSPDefaultBuildResidual;
243: ksp->ops->setfromoptions = KSPSetFromOptions_BCGS;
244: ksp->ops->view = KSPView_BCGS;
245:
246: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);
247: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,1);
248: return(0);
249: }
250: EXTERN_C_END