1 #include <../src/ksp/ksp/impls/bcgs/bcgsimpl.h> /*I "petscksp.h" I*/
2
KSPSetFromOptions_BCGS(KSP ksp,PetscOptionItems PetscOptionsObject)3 PetscErrorCode KSPSetFromOptions_BCGS(KSP ksp, PetscOptionItems PetscOptionsObject)
4 {
5 PetscFunctionBegin;
6 PetscOptionsHeadBegin(PetscOptionsObject, "KSP BCGS Options");
7 PetscOptionsHeadEnd();
8 PetscFunctionReturn(PETSC_SUCCESS);
9 }
10
KSPSetUp_BCGS(KSP ksp)11 PetscErrorCode KSPSetUp_BCGS(KSP ksp)
12 {
13 PetscFunctionBegin;
14 PetscCall(KSPSetWorkVecs(ksp, 6));
15 PetscFunctionReturn(PETSC_SUCCESS);
16 }
17
KSPSolve_BCGS(KSP ksp)18 PetscErrorCode KSPSolve_BCGS(KSP ksp)
19 {
20 PetscInt i;
21 PetscScalar rho, rhoold, alpha, beta, omega, omegaold, d1;
22 Vec X, B, V, P, R, RP, T, S;
23 PetscReal dp = 0.0, d2;
24 KSP_BCGS *bcgs = (KSP_BCGS *)ksp->data;
25
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
36 /* Compute initial preconditioned residual */
37 PetscCall(KSPInitialResidual(ksp, X, V, T, R, B));
38
39 /* with right preconditioning need to save initial guess to add to final solution */
40 if (ksp->pc_side == PC_RIGHT && !ksp->guess_zero) {
41 if (!bcgs->guess) PetscCall(VecDuplicate(X, &bcgs->guess));
42 PetscCall(VecCopy(X, bcgs->guess));
43 PetscCall(VecSet(X, 0.0));
44 }
45
46 /* Test for nothing to do */
47 if (ksp->normtype != KSP_NORM_NONE) {
48 PetscCall(VecNorm(R, NORM_2, &dp));
49 KSPCheckNorm(ksp, dp);
50 }
51 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
52 ksp->its = 0;
53 ksp->rnorm = dp;
54 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
55 PetscCall(KSPLogResidualHistory(ksp, dp));
56 PetscCall(KSPMonitor(ksp, 0, dp));
57 PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP));
58 if (ksp->reason) {
59 if (bcgs->guess) PetscCall(VecAXPY(X, 1.0, bcgs->guess));
60 PetscFunctionReturn(PETSC_SUCCESS);
61 }
62
63 /* Make the initial Rp == R */
64 PetscCall(VecCopy(R, RP));
65
66 rhoold = 1.0;
67 alpha = 1.0;
68 omegaold = 1.0;
69 PetscCall(VecSet(P, 0.0));
70 PetscCall(VecSet(V, 0.0));
71
72 i = 0;
73 do {
74 PetscCall(VecDot(R, RP, &rho)); /* rho <- (r,rp) */
75 beta = (rho / rhoold) * (alpha / omegaold);
76 PetscCall(VecAXPBYPCZ(P, 1.0, -omegaold * beta, beta, R, V)); /* p <- r - omega * beta* v + beta * p */
77 PetscCall(KSP_PCApplyBAorAB(ksp, P, V, T)); /* v <- K p */
78 PetscCall(VecDot(V, RP, &d1));
79 KSPCheckDot(ksp, d1);
80 if (d1 == 0.0) {
81 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve breakdown due to zero inner product");
82 ksp->reason = KSP_DIVERGED_BREAKDOWN;
83 PetscCall(PetscInfo(ksp, "Breakdown due to zero inner product\n"));
84 break;
85 }
86 alpha = rho / d1; /* a <- rho / (v,rp) */
87 PetscCall(VecWAXPY(S, -alpha, V, R)); /* s <- r - a v */
88 PetscCall(KSP_PCApplyBAorAB(ksp, S, T, R)); /* t <- K s */
89 PetscCall(VecDotNorm2(S, T, &d1, &d2));
90 if (d2 == 0.0) {
91 /* t is 0. if s is 0, then alpha v == r, and hence alpha p
92 may be our solution. Give it a try? */
93 PetscCall(VecDot(S, S, &d1));
94 if (d1 != 0.0) {
95 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has failed due to singular preconditioned operator");
96 ksp->reason = KSP_DIVERGED_BREAKDOWN;
97 PetscCall(PetscInfo(ksp, "Failed due to singular preconditioned operator\n"));
98 break;
99 }
100 PetscCall(VecAXPY(X, alpha, P)); /* x <- x + a p */
101 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
102 ksp->its++;
103 ksp->rnorm = 0.0;
104 ksp->reason = KSP_CONVERGED_RTOL;
105 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
106 PetscCall(KSPLogResidualHistory(ksp, dp));
107 PetscCall(KSPMonitor(ksp, i + 1, 0.0));
108 break;
109 }
110 omega = d1 / d2; /* w <- (t's) / (t't) */
111 PetscCall(VecAXPBYPCZ(X, alpha, omega, 1.0, P, S)); /* x <- alpha * p + omega * s + x */
112 PetscCall(VecWAXPY(R, -omega, T, S)); /* r <- s - w t */
113 if (ksp->normtype != KSP_NORM_NONE && ksp->chknorm < i + 2) {
114 PetscCall(VecNorm(R, NORM_2, &dp));
115 KSPCheckNorm(ksp, dp);
116 }
117
118 rhoold = rho;
119 omegaold = omega;
120
121 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
122 ksp->its++;
123 ksp->rnorm = dp;
124 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
125 PetscCall(KSPLogResidualHistory(ksp, dp));
126 PetscCall(KSPMonitor(ksp, i + 1, dp));
127 PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
128 if (ksp->reason) break;
129 if (rho == 0.0) {
130 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve breakdown due to zero inner product");
131 ksp->reason = KSP_DIVERGED_BREAKDOWN;
132 PetscCall(PetscInfo(ksp, "Breakdown due to zero rho inner product\n"));
133 break;
134 }
135 i++;
136 } while (i < ksp->max_it);
137
138 if (i >= ksp->max_it) ksp->reason = KSP_DIVERGED_ITS;
139
140 PetscCall(KSPUnwindPreconditioner(ksp, X, T));
141 if (bcgs->guess) PetscCall(VecAXPY(X, 1.0, bcgs->guess));
142 PetscFunctionReturn(PETSC_SUCCESS);
143 }
144
KSPBuildSolution_BCGS(KSP ksp,Vec v,Vec * V)145 PetscErrorCode KSPBuildSolution_BCGS(KSP ksp, Vec v, Vec *V)
146 {
147 KSP_BCGS *bcgs = (KSP_BCGS *)ksp->data;
148
149 PetscFunctionBegin;
150 if (ksp->pc_side == PC_RIGHT) {
151 PetscCheck(v, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "Not working with right preconditioner");
152 PetscCall(KSP_PCApply(ksp, ksp->vec_sol, v));
153 if (bcgs->guess) PetscCall(VecAXPY(v, 1.0, bcgs->guess));
154 *V = v;
155 } else {
156 if (v) {
157 PetscCall(VecCopy(ksp->vec_sol, v));
158 *V = v;
159 } else *V = ksp->vec_sol;
160 }
161 PetscFunctionReturn(PETSC_SUCCESS);
162 }
163
KSPReset_BCGS(KSP ksp)164 PetscErrorCode KSPReset_BCGS(KSP ksp)
165 {
166 KSP_BCGS *cg = (KSP_BCGS *)ksp->data;
167
168 PetscFunctionBegin;
169 PetscCall(VecDestroy(&cg->guess));
170 PetscFunctionReturn(PETSC_SUCCESS);
171 }
172
KSPDestroy_BCGS(KSP ksp)173 PetscErrorCode KSPDestroy_BCGS(KSP ksp)
174 {
175 PetscFunctionBegin;
176 PetscCall(KSPReset_BCGS(ksp));
177 PetscCall(KSPDestroyDefault(ksp));
178 PetscFunctionReturn(PETSC_SUCCESS);
179 }
180
181 /*MC
182 KSPBCGS - Implements the BiCGStab (Stabilized version of Biconjugate Gradient) method {cite}`vorst92`
183
184 Level: beginner
185
186 Notes:
187 Supports left and right preconditioning but not symmetric
188
189 See `KSPBCGSL` for additional stabilization
190
191 See `KSPFBCGS`, `KSPFBCGSR`, and `KSPPIPEBCGS` for flexible and pipelined versions of the algorithm
192
193 This method can be a good alternative to `KSPGMRES` in that it does not require explicitly storing the Krylov space as `KSPGMRES` requires
194 and there is no acceptable restart value that can be set with `KSPGMRESSetRestart()` that balances cost per iteration and convergence rates with `KSPGMRES`.
195
196 .seealso: [](ch_ksp), `KSPFBCGS`, `KSPFBCGSR`, `KSPPIPEBCGS`, `KSPBCGSL`, `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPBICG`, `KSPBCGSL`, `KSPFBICG`, `KSPQMRCGS`, `KSPSetPCSide()`
197 M*/
KSPCreate_BCGS(KSP ksp)198 PETSC_EXTERN PetscErrorCode KSPCreate_BCGS(KSP ksp)
199 {
200 KSP_BCGS *bcgs;
201
202 PetscFunctionBegin;
203 PetscCall(PetscNew(&bcgs));
204
205 ksp->data = bcgs;
206 ksp->ops->setup = KSPSetUp_BCGS;
207 ksp->ops->solve = KSPSolve_BCGS;
208 ksp->ops->destroy = KSPDestroy_BCGS;
209 ksp->ops->reset = KSPReset_BCGS;
210 ksp->ops->buildsolution = KSPBuildSolution_BCGS;
211 ksp->ops->buildresidual = KSPBuildResidualDefault;
212 ksp->ops->setfromoptions = KSPSetFromOptions_BCGS;
213
214 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
215 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
216 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_LEFT, 1));
217 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
218 PetscFunctionReturn(PETSC_SUCCESS);
219 }
220