1 /*
2 This file implements flexible BiCGStab (FBiCGStab).
3 */
4 #include <../src/ksp/ksp/impls/bcgs/bcgsimpl.h> /*I "petscksp.h" I*/
5
KSPSetUp_FBCGS(KSP ksp)6 static PetscErrorCode KSPSetUp_FBCGS(KSP ksp)
7 {
8 PetscFunctionBegin;
9 PetscCall(KSPSetWorkVecs(ksp, 8));
10 PetscFunctionReturn(PETSC_SUCCESS);
11 }
12
13 /* Only need a few hacks from KSPSolve_BCGS */
14
KSPSolve_FBCGS(KSP ksp)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;
24
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];
36
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 }
43
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 }
54
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);
65
66 /* Make the initial Rp == R */
67 PetscCall(VecCopy(R, RP));
68
69 rhoold = 1.0;
70 alpha = 1.0;
71 omegaold = 1.0;
72 PetscCall(VecSet(P, 0.0));
73 PetscCall(VecSet(V, 0.0));
74
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 */
80
81 PetscCall(KSP_PCApply(ksp, P, P2)); /* p2 <- K p */
82 PetscCall(KSP_MatMult(ksp, mat, P2, V)); /* v <- A p2 */
83
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 */
93
94 PetscCall(KSP_PCApply(ksp, S, S2)); /* s2 <- K s */
95 PetscCall(KSP_MatMult(ksp, mat, S2, T)); /* t <- A s2 */
96
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 */
119
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));
122
123 rhoold = rho;
124 omegaold = omega;
125
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);
142
143 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
144 PetscFunctionReturn(PETSC_SUCCESS);
145 }
146
147 /*MC
148 KSPFBCGS - Implements the flexible BiCGStab method. [](sec_flexibleksp)
149
150 Level: beginner
151
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.
155
156 `KSPFBCGSR` provides another variant of this algorithm that requires fewer `MPI_Allreduce()` calls and my converge faster
157
158 See `KSPPIPEBCGS` for a pipelined version of the algorithm
159
160 Only supportst right preconditioning
161
162 .seealso: [](ch_ksp), [](sec_flexibleksp), `KSPFBCGSR`, `KSPPIPEBCGS`, `KSPBCGSL`, `KSPFBCGS`, `KSPCreate()`, `KSPFGMRES`, `KSPSetType()`, `KSPType`, `KSP`, `KSPBICG`, `KSPSetPCSide()`
163 M*/
KSPCreate_FBCGS(KSP ksp)164 PETSC_EXTERN PetscErrorCode KSPCreate_FBCGS(KSP ksp)
165 {
166 KSP_BCGS *bcgs;
167
168 PetscFunctionBegin;
169 PetscCall(PetscNew(&bcgs));
170
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;
179
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 }
184