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: PetscFunctionBegin;
10: PetscCall(KSPSetWorkVecs(ksp, 8));
11: PetscFunctionReturn(PETSC_SUCCESS);
12: }
14: /* Only need a few hacks from KSPSolve_BCGS */
16: static PetscErrorCode KSPSolve_FBCGS(KSP ksp)
17: {
18: PetscInt i;
19: PetscScalar rho, rhoold, alpha, beta, omega, omegaold, d1;
20: Vec X, B, V, P, R, RP, T, S, P2, S2;
21: PetscReal dp = 0.0, d2;
22: KSP_BCGS *bcgs = (KSP_BCGS *)ksp->data;
23: PC pc;
24: Mat mat;
26: PetscFunctionBegin;
27: X = ksp->vec_sol;
28: B = ksp->vec_rhs;
29: R = ksp->work[0];
30: RP = ksp->work[1];
31: V = ksp->work[2];
32: T = ksp->work[3];
33: S = ksp->work[4];
34: P = ksp->work[5];
35: S2 = ksp->work[6];
36: P2 = ksp->work[7];
38: if (!ksp->guess_zero) {
39: if (!bcgs->guess) PetscCall(VecDuplicate(X, &bcgs->guess));
40: PetscCall(VecCopy(X, bcgs->guess));
41: } else {
42: PetscCall(VecSet(X, 0.0));
43: }
45: /* Compute initial residual */
46: PetscCall(KSPGetPC(ksp, &pc));
47: PetscCall(PCSetUp(pc));
48: PetscCall(PCGetOperators(pc, &mat, NULL));
49: if (!ksp->guess_zero) {
50: PetscCall(KSP_MatMult(ksp, mat, X, S2));
51: PetscCall(VecCopy(B, R));
52: PetscCall(VecAXPY(R, -1.0, S2));
53: } else {
54: PetscCall(VecCopy(B, R));
55: }
57: /* Test for nothing to do */
58: if (ksp->normtype != KSP_NORM_NONE) PetscCall(VecNorm(R, NORM_2, &dp));
59: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
60: ksp->its = 0;
61: ksp->rnorm = dp;
62: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
63: PetscCall(KSPLogResidualHistory(ksp, dp));
64: PetscCall(KSPMonitor(ksp, 0, dp));
65: PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP));
66: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
68: /* Make the initial Rp == R */
69: PetscCall(VecCopy(R, RP));
71: rhoold = 1.0;
72: alpha = 1.0;
73: omegaold = 1.0;
74: PetscCall(VecSet(P, 0.0));
75: PetscCall(VecSet(V, 0.0));
77: i = 0;
78: do {
79: PetscCall(VecDot(R, RP, &rho)); /* rho <- (r,rp) */
80: beta = (rho / rhoold) * (alpha / omegaold);
81: PetscCall(VecAXPBYPCZ(P, 1.0, -omegaold * beta, beta, R, V)); /* p <- r - omega * beta* v + beta * p */
83: PetscCall(KSP_PCApply(ksp, P, P2)); /* p2 <- K p */
84: PetscCall(KSP_MatMult(ksp, mat, P2, V)); /* v <- A p2 */
86: PetscCall(VecDot(V, RP, &d1));
87: if (d1 == 0.0) {
88: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve breakdown due to zero inner product");
89: ksp->reason = KSP_DIVERGED_BREAKDOWN;
90: PetscCall(PetscInfo(ksp, "Breakdown due to zero inner product\n"));
91: break;
92: }
93: alpha = rho / d1; /* alpha <- rho / (v,rp) */
94: PetscCall(VecWAXPY(S, -alpha, V, R)); /* s <- r - alpha v */
96: PetscCall(KSP_PCApply(ksp, S, S2)); /* s2 <- K s */
97: PetscCall(KSP_MatMult(ksp, mat, S2, T)); /* t <- A s2 */
99: PetscCall(VecDotNorm2(S, T, &d1, &d2));
100: if (d2 == 0.0) {
101: /* t is 0. if s is 0, then alpha v == r, and hence alpha p may be our solution. Give it a try? */
102: PetscCall(VecDot(S, S, &d1));
103: if (d1 != 0.0) {
104: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has failed due to singular preconditioned operator");
105: ksp->reason = KSP_DIVERGED_BREAKDOWN;
106: PetscCall(PetscInfo(ksp, "Failed due to singular preconditioned operator\n"));
107: break;
108: }
109: PetscCall(VecAXPY(X, alpha, P2)); /* x <- x + alpha p2 */
110: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
111: ksp->its++;
112: ksp->rnorm = 0.0;
113: ksp->reason = KSP_CONVERGED_RTOL;
114: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
115: PetscCall(KSPLogResidualHistory(ksp, dp));
116: PetscCall(KSPMonitor(ksp, i + 1, 0.0));
117: break;
118: }
119: omega = d1 / d2; /* omega <- (t's) / (t't) */
120: PetscCall(VecAXPBYPCZ(X, alpha, omega, 1.0, P2, S2)); /* x <- alpha * p2 + omega * s2 + x */
122: PetscCall(VecWAXPY(R, -omega, T, S)); /* r <- s - omega t */
123: if (ksp->normtype != KSP_NORM_NONE && ksp->chknorm < i + 2) PetscCall(VecNorm(R, NORM_2, &dp));
125: rhoold = rho;
126: omegaold = omega;
128: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
129: ksp->its++;
130: ksp->rnorm = dp;
131: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
132: PetscCall(KSPLogResidualHistory(ksp, dp));
133: PetscCall(KSPMonitor(ksp, i + 1, dp));
134: PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
135: if (ksp->reason) break;
136: if (rho == 0.0) {
137: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve breakdown due to zero rho inner product");
138: ksp->reason = KSP_DIVERGED_BREAKDOWN;
139: PetscCall(PetscInfo(ksp, "Breakdown due to zero rho inner product\n"));
140: break;
141: }
142: i++;
143: } while (i < ksp->max_it);
145: if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*MC
150: KSPFBCGS - Implements the flexible BiCGStab method. [](sec_flexibleksp)
152: Level: beginner
154: Notes:
155: 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
156: in the vector entries.
158: `KSPFBCGSR` provides another variant of this algorithm that requires fewer `MPI_Allreduce()` calls and my converge faster
160: See `KSPPIPEBCGS` for a pipelined version of the algorithm
162: Only supportst right preconditioning
164: .seealso: [](ch_ksp), [](sec_flexibleksp), `KSPFBCGSR`, `KSPPIPEBCGS`, `KSPBCGSL`, `KSPFBCGS`, `KSPCreate()`, `KSPFGMRES`, `KSPSetType()`, `KSPType`, `KSP`, `KSPBICG`, `KSPSetPCSide()`
165: M*/
166: PETSC_EXTERN PetscErrorCode KSPCreate_FBCGS(KSP ksp)
167: {
168: KSP_BCGS *bcgs;
170: PetscFunctionBegin;
171: PetscCall(PetscNew(&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;
182: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
183: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }