Actual source code: fbcgs.c
1: /*
2: This file implements flexible BiCGStab (FBiCGStab).
3: */
4: #include <../src/ksp/ksp/impls/bcgs/bcgsimpl.h>
6: static PetscErrorCode KSPSetUp_FBCGS(KSP ksp)
7: {
8: PetscFunctionBegin;
9: PetscCall(KSPSetWorkVecs(ksp, 8));
10: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
26: X = ksp->vec_sol;
27: B = ksp->vec_rhs;
28: R = ksp->work[0];
29: RP = ksp->work[1];
30: V = ksp->work[2];
31: T = ksp->work[3];
32: S = ksp->work[4];
33: P = ksp->work[5];
34: S2 = ksp->work[6];
35: P2 = ksp->work[7];
37: if (!ksp->guess_zero) {
38: if (!bcgs->guess) PetscCall(VecDuplicate(X, &bcgs->guess));
39: PetscCall(VecCopy(X, bcgs->guess));
40: } else {
41: PetscCall(VecSet(X, 0.0));
42: }
44: /* Compute initial residual */
45: PetscCall(KSPGetPC(ksp, &pc));
46: PetscCall(PCGetOperators(pc, &mat, NULL));
47: if (!ksp->guess_zero) {
48: PetscCall(KSP_MatMult(ksp, mat, X, S2));
49: PetscCall(VecCopy(B, R));
50: PetscCall(VecAXPY(R, -1.0, S2));
51: } else {
52: PetscCall(VecCopy(B, R));
53: }
55: /* Test for nothing to do */
56: if (ksp->normtype != KSP_NORM_NONE) PetscCall(VecNorm(R, NORM_2, &dp));
57: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
58: ksp->its = 0;
59: ksp->rnorm = dp;
60: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
61: PetscCall(KSPLogResidualHistory(ksp, dp));
62: PetscCall(KSPMonitor(ksp, 0, dp));
63: PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP));
64: if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
66: /* Make the initial Rp == R */
67: PetscCall(VecCopy(R, RP));
69: rhoold = 1.0;
70: alpha = 1.0;
71: omegaold = 1.0;
72: PetscCall(VecSet(P, 0.0));
73: PetscCall(VecSet(V, 0.0));
75: i = 0;
76: do {
77: PetscCall(VecDot(R, RP, &rho)); /* rho <- (r,rp) */
78: beta = (rho / rhoold) * (alpha / omegaold);
79: PetscCall(VecAXPBYPCZ(P, 1.0, -omegaold * beta, beta, R, V)); /* p <- r - omega * beta* v + beta * p */
81: PetscCall(KSP_PCApply(ksp, P, P2)); /* p2 <- K p */
82: PetscCall(KSP_MatMult(ksp, mat, P2, V)); /* v <- A p2 */
84: PetscCall(VecDot(V, RP, &d1));
85: if (d1 == 0.0) {
86: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve breakdown due to zero inner product");
87: ksp->reason = KSP_DIVERGED_BREAKDOWN;
88: PetscCall(PetscInfo(ksp, "Breakdown due to zero inner product\n"));
89: break;
90: }
91: alpha = rho / d1; /* alpha <- rho / (v,rp) */
92: PetscCall(VecWAXPY(S, -alpha, V, R)); /* s <- r - alpha v */
94: PetscCall(KSP_PCApply(ksp, S, S2)); /* s2 <- K s */
95: PetscCall(KSP_MatMult(ksp, mat, S2, T)); /* t <- A s2 */
97: PetscCall(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: PetscCall(VecDot(S, S, &d1));
101: if (d1 != 0.0) {
102: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has failed due to singular preconditioned operator");
103: ksp->reason = KSP_DIVERGED_BREAKDOWN;
104: PetscCall(PetscInfo(ksp, "Failed due to singular preconditioned operator\n"));
105: break;
106: }
107: PetscCall(VecAXPY(X, alpha, P2)); /* x <- x + alpha p2 */
108: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
109: ksp->its++;
110: ksp->rnorm = 0.0;
111: ksp->reason = KSP_CONVERGED_RTOL;
112: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
113: PetscCall(KSPLogResidualHistory(ksp, dp));
114: PetscCall(KSPMonitor(ksp, i + 1, 0.0));
115: break;
116: }
117: omega = d1 / d2; /* omega <- (t's) / (t't) */
118: PetscCall(VecAXPBYPCZ(X, alpha, omega, 1.0, P2, S2)); /* x <- alpha * p2 + omega * s2 + x */
120: PetscCall(VecWAXPY(R, -omega, T, S)); /* r <- s - omega t */
121: if (ksp->normtype != KSP_NORM_NONE && ksp->chknorm < i + 2) PetscCall(VecNorm(R, NORM_2, &dp));
123: rhoold = rho;
124: omegaold = omega;
126: PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
127: ksp->its++;
128: ksp->rnorm = dp;
129: PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
130: PetscCall(KSPLogResidualHistory(ksp, dp));
131: PetscCall(KSPMonitor(ksp, i + 1, dp));
132: PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
133: if (ksp->reason) break;
134: if (rho == 0.0) {
135: PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve breakdown due to zero rho inner product");
136: ksp->reason = KSP_DIVERGED_BREAKDOWN;
137: PetscCall(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: PetscFunctionReturn(PETSC_SUCCESS);
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: [](ch_ksp), [](sec_flexibleksp), `KSPFBCGSR`, `KSPPIPEBCGS`, `KSPBCGSL`, `KSPFBCGS`, `KSPCreate()`, `KSPFGMRES`, `KSPSetType()`, `KSPType`, `KSP`, `KSPBICG`, `KSPSetPCSide()`
163: M*/
164: PETSC_EXTERN PetscErrorCode KSPCreate_FBCGS(KSP ksp)
165: {
166: KSP_BCGS *bcgs;
168: PetscFunctionBegin;
169: PetscCall(PetscNew(&bcgs));
171: ksp->data = bcgs;
172: ksp->ops->setup = KSPSetUp_FBCGS;
173: ksp->ops->solve = KSPSolve_FBCGS;
174: ksp->ops->destroy = KSPDestroy_BCGS;
175: ksp->ops->reset = KSPReset_BCGS;
176: ksp->ops->buildresidual = KSPBuildResidualDefault;
177: ksp->ops->setfromoptions = KSPSetFromOptions_BCGS;
178: ksp->pc_side = PC_RIGHT;
180: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
181: PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }