Actual source code: fbcgs.c
2: /*
3: This file implements flexible BiCGStab (FBiCGStab).
4: */
5: #include <../src/ksp/ksp/impls/bcgs/bcgsimpl.h>
7: static PetscErrorCode KSPSetUp_FBCGS(KSP ksp)
8: {
9: KSPSetWorkVecs(ksp, 8);
10: return 0;
11: }
13: /* Only need a few hacks from KSPSolve_BCGS */
15: static PetscErrorCode KSPSolve_FBCGS(KSP ksp)
16: {
17: PetscInt i;
18: PetscScalar rho, rhoold, alpha, beta, omega, omegaold, d1;
19: Vec X, B, V, P, R, RP, T, S, P2, S2;
20: PetscReal dp = 0.0, d2;
21: KSP_BCGS *bcgs = (KSP_BCGS *)ksp->data;
22: PC pc;
23: Mat mat;
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];
33: S2 = ksp->work[6];
34: P2 = ksp->work[7];
36: if (!ksp->guess_zero) {
37: if (!bcgs->guess) VecDuplicate(X, &bcgs->guess);
38: VecCopy(X, bcgs->guess);
39: } else {
40: VecSet(X, 0.0);
41: }
43: /* Compute initial residual */
44: KSPGetPC(ksp, &pc);
45: PCSetUp(pc);
46: PCGetOperators(pc, &mat, NULL);
47: if (!ksp->guess_zero) {
48: KSP_MatMult(ksp, mat, X, S2);
49: VecCopy(B, R);
50: VecAXPY(R, -1.0, S2);
51: } else {
52: VecCopy(B, R);
53: }
55: /* Test for nothing to do */
56: if (ksp->normtype != KSP_NORM_NONE) VecNorm(R, NORM_2, &dp);
57: PetscObjectSAWsTakeAccess((PetscObject)ksp);
58: ksp->its = 0;
59: ksp->rnorm = dp;
60: PetscObjectSAWsGrantAccess((PetscObject)ksp);
61: KSPLogResidualHistory(ksp, dp);
62: KSPMonitor(ksp, 0, dp);
63: (*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP);
64: if (ksp->reason) return 0;
66: /* Make the initial Rp == R */
67: VecCopy(R, RP);
69: rhoold = 1.0;
70: alpha = 1.0;
71: omegaold = 1.0;
72: VecSet(P, 0.0);
73: VecSet(V, 0.0);
75: i = 0;
76: do {
77: VecDot(R, RP, &rho); /* rho <- (r,rp) */
78: beta = (rho / rhoold) * (alpha / omegaold);
79: VecAXPBYPCZ(P, 1.0, -omegaold * beta, beta, R, V); /* p <- r - omega * beta* v + beta * p */
81: KSP_PCApply(ksp, P, P2); /* p2 <- K p */
82: KSP_MatMult(ksp, mat, P2, V); /* v <- A p2 */
84: VecDot(V, RP, &d1);
85: if (d1 == 0.0) {
87: ksp->reason = KSP_DIVERGED_BREAKDOWN;
88: PetscInfo(ksp, "Breakdown due to zero inner product\n");
89: break;
90: }
91: alpha = rho / d1; /* alpha <- rho / (v,rp) */
92: VecWAXPY(S, -alpha, V, R); /* s <- r - alpha v */
94: KSP_PCApply(ksp, S, S2); /* s2 <- K s */
95: KSP_MatMult(ksp, mat, S2, T); /* t <- A s2 */
97: VecDotNorm2(S, T, &d1, &d2);
98: if (d2 == 0.0) {
99: /* t is 0. if s is 0, then alpha v == r, and hence alpha p may be our solution. Give it a try? */
100: VecDot(S, S, &d1);
101: if (d1 != 0.0) {
103: ksp->reason = KSP_DIVERGED_BREAKDOWN;
104: PetscInfo(ksp, "Failed due to singular preconditioned operator\n");
105: break;
106: }
107: VecAXPY(X, alpha, P2); /* x <- x + alpha p2 */
108: PetscObjectSAWsTakeAccess((PetscObject)ksp);
109: ksp->its++;
110: ksp->rnorm = 0.0;
111: ksp->reason = KSP_CONVERGED_RTOL;
112: PetscObjectSAWsGrantAccess((PetscObject)ksp);
113: KSPLogResidualHistory(ksp, dp);
114: KSPMonitor(ksp, i + 1, 0.0);
115: break;
116: }
117: omega = d1 / d2; /* omega <- (t's) / (t't) */
118: VecAXPBYPCZ(X, alpha, omega, 1.0, P2, S2); /* x <- alpha * p2 + omega * s2 + x */
120: VecWAXPY(R, -omega, T, S); /* r <- s - omega t */
121: if (ksp->normtype != KSP_NORM_NONE && ksp->chknorm < i + 2) VecNorm(R, NORM_2, &dp);
123: rhoold = rho;
124: omegaold = omega;
126: PetscObjectSAWsTakeAccess((PetscObject)ksp);
127: ksp->its++;
128: ksp->rnorm = dp;
129: PetscObjectSAWsGrantAccess((PetscObject)ksp);
130: KSPLogResidualHistory(ksp, dp);
131: KSPMonitor(ksp, i + 1, dp);
132: (*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP);
133: if (ksp->reason) break;
134: if (rho == 0.0) {
136: ksp->reason = KSP_DIVERGED_BREAKDOWN;
137: PetscInfo(ksp, "Breakdown due to zero rho inner product\n");
138: break;
139: }
140: i++;
141: } while (i < ksp->max_it);
143: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
144: return 0;
145: }
147: /*MC
148: KSPFBCGS - Implements the flexible BiCGStab method. [](sec_flexibleksp)
150: Level: beginner
152: Notes:
153: Flexible BiCGStab, unlike most Krylov methods allows the preconditioner to be nonlinear, that is the action of the preconditioner to a vector need not be linear
154: in the vector entries.
156: `KSPFBCGSR` provides another variant of this algorithm that requires fewer `MPI_Allreduce()` calls and my converge faster
158: See `KSPPIPEBCGS` for a pipelined version of the algorithm
160: Only supportst right preconditioning
162: .seealso: [](chapter_ksp), [](sec_flexibleksp), `KSPFBCGSR`, `KSPPIPEBCGS`, `KSPBCGSL`, `KSPFBCGS`, `KSPCreate()`, `KSPFGMRES`, `KSPSetType()`, `KSPType`, `KSP`, `KSPBICG`, `KSPFBCGSL`, `KSPSetPCSide()`
163: M*/
164: PETSC_EXTERN PetscErrorCode KSPCreate_FBCGS(KSP ksp)
165: {
166: KSP_BCGS *bcgs;
168: PetscNew(&bcgs);
170: ksp->data = bcgs;
171: ksp->ops->setup = KSPSetUp_FBCGS;
172: ksp->ops->solve = KSPSolve_FBCGS;
173: ksp->ops->destroy = KSPDestroy_BCGS;
174: ksp->ops->reset = KSPReset_BCGS;
175: ksp->ops->buildresidual = KSPBuildResidualDefault;
176: ksp->ops->setfromoptions = KSPSetFromOptions_BCGS;
177: ksp->pc_side = PC_RIGHT;
179: KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2);
180: KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1);
181: return 0;
182: }