Actual source code: fbcgs.c
petsc-3.10.5 2019-03-28
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: {
13: KSPSetWorkVecs(ksp,8);
14: return(0);
15: }
17: /* Only need a few hacks from KSPSolve_BCGS */
19: static PetscErrorCode KSPSolve_FBCGS(KSP ksp)
20: {
22: PetscInt i;
23: PetscScalar rho,rhoold,alpha,beta,omega,omegaold,d1;
24: Vec X,B,V,P,R,RP,T,S,P2,S2;
25: PetscReal dp = 0.0,d2;
26: KSP_BCGS *bcgs = (KSP_BCGS*)ksp->data;
27: PC pc;
28: Mat mat;
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];
39: S2 = ksp->work[6];
40: P2 = ksp->work[7];
42: /* Only supports right preconditioning */
43: if (ksp->pc_side != PC_RIGHT) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"KSP fbcgs does not support %s",PCSides[ksp->pc_side]);
44: if (!ksp->guess_zero) {
45: if (!bcgs->guess) {
46: VecDuplicate(X,&bcgs->guess);
47: }
48: VecCopy(X,bcgs->guess);
49: } else {
50: VecSet(X,0.0);
51: }
53: /* Compute initial residual */
54: KSPGetPC(ksp,&pc);
55: PCSetUp(pc);
56: PCGetOperators(pc,&mat,NULL);
57: if (!ksp->guess_zero) {
58: KSP_MatMult(ksp,mat,X,S2);
59: VecCopy(B,R);
60: VecAXPY(R,-1.0,S2);
61: } else {
62: VecCopy(B,R);
63: }
65: /* Test for nothing to do */
66: if (ksp->normtype != KSP_NORM_NONE) {
67: VecNorm(R,NORM_2,&dp);
68: }
69: PetscObjectSAWsTakeAccess((PetscObject)ksp);
70: ksp->its = 0;
71: ksp->rnorm = dp;
72: PetscObjectSAWsGrantAccess((PetscObject)ksp);
73: KSPLogResidualHistory(ksp,dp);
74: KSPMonitor(ksp,0,dp);
75: (*ksp->converged)(ksp,0,dp,&ksp->reason,ksp->cnvP);
76: if (ksp->reason) return(0);
78: /* Make the initial Rp == R */
79: VecCopy(R,RP);
81: rhoold = 1.0;
82: alpha = 1.0;
83: omegaold = 1.0;
84: VecSet(P,0.0);
85: VecSet(V,0.0);
87: i=0;
88: do {
89: VecDot(R,RP,&rho); /* rho <- (r,rp) */
90: beta = (rho/rhoold) * (alpha/omegaold);
91: VecAXPBYPCZ(P,1.0,-omegaold*beta,beta,R,V); /* p <- r - omega * beta* v + beta * p */
93: KSP_PCApply(ksp,P,P2); /* p2 <- K p */
94: KSP_MatMult(ksp,mat,P2,V); /* v <- A p2 */
96: VecDot(V,RP,&d1);
97: if (d1 == 0.0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Divide by zero");
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) {
109: ksp->reason = KSP_DIVERGED_BREAKDOWN;
110: break;
111: }
112: VecAXPY(X,alpha,P2); /* x <- x + alpha p2 */
113: PetscObjectSAWsTakeAccess((PetscObject)ksp);
114: ksp->its++;
115: ksp->rnorm = 0.0;
116: ksp->reason = KSP_CONVERGED_RTOL;
117: PetscObjectSAWsGrantAccess((PetscObject)ksp);
118: KSPLogResidualHistory(ksp,dp);
119: KSPMonitor(ksp,i+1,0.0);
120: break;
121: }
122: omega = d1 / d2; /* omega <- (t's) / (t't) */
123: VecAXPBYPCZ(X,alpha,omega,1.0,P2,S2); /* x <- alpha * p2 + omega * s2 + x */
125: VecWAXPY(R,-omega,T,S); /* r <- s - omega t */
126: if (ksp->normtype != KSP_NORM_NONE && ksp->chknorm < i+2) {
127: VecNorm(R,NORM_2,&dp);
128: }
130: rhoold = rho;
131: omegaold = omega;
133: PetscObjectSAWsTakeAccess((PetscObject)ksp);
134: ksp->its++;
135: ksp->rnorm = dp;
136: PetscObjectSAWsGrantAccess((PetscObject)ksp);
137: KSPLogResidualHistory(ksp,dp);
138: KSPMonitor(ksp,i+1,dp);
139: (*ksp->converged)(ksp,i+1,dp,&ksp->reason,ksp->cnvP);
140: if (ksp->reason) break;
141: if (rho == 0.0) {
142: ksp->reason = KSP_DIVERGED_BREAKDOWN;
143: break;
144: }
145: i++;
146: } while (i<ksp->max_it);
148: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
149: return(0);
150: }
152: /*MC
153: KSPFBCGS - Implements flexible BiCGStab method.
155: Options Database Keys:
156: . see KSPSolve()
158: Level: beginner
160: Notes:
161: Only allow right preconditioning
163: .seealso: KSPCreate(), KSPSetType(), KSPType (for list of available types), KSP, KSPBICG, KSPFBCGSL, KSPSetPCSide()
164: M*/
165: PETSC_EXTERN PetscErrorCode KSPCreate_FBCGS(KSP ksp)
166: {
168: KSP_BCGS *bcgs;
171: PetscNewLog(ksp,&bcgs);
173: ksp->data = bcgs;
174: ksp->ops->setup = KSPSetUp_FBCGS;
175: ksp->ops->solve = KSPSolve_FBCGS;
176: ksp->ops->destroy = KSPDestroy_BCGS;
177: ksp->ops->reset = KSPReset_BCGS;
178: ksp->ops->buildresidual = KSPBuildResidualDefault;
179: ksp->ops->setfromoptions = KSPSetFromOptions_BCGS;
180: ksp->pc_side = PC_RIGHT; /* set default PC side */
182: KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,3);
183: KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,2);
184: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
185: KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
186: return(0);
187: }