1 /*
2 This file implements QMRCGS (QMRCGStab).
3 */
4 #include <../src/ksp/ksp/impls/bcgs/bcgsimpl.h> /*I "petscksp.h" I*/
5
KSPSetUp_QMRCGS(KSP ksp)6 static PetscErrorCode KSPSetUp_QMRCGS(KSP ksp)
7 {
8 PetscFunctionBegin;
9 PetscCall(KSPSetWorkVecs(ksp, 14));
10 PetscFunctionReturn(PETSC_SUCCESS);
11 }
12
13 /* Only need a few hacks from KSPSolve_BCGS */
14
KSPSolve_QMRCGS(KSP ksp)15 static PetscErrorCode KSPSolve_QMRCGS(KSP ksp)
16 {
17 PetscInt i;
18 PetscScalar eta, rho1, rho2, alpha, eta2, omega, beta, cf, cf1, uu;
19 Vec X, B, R, P, PH, V, D2, X2, S, SH, T, D, S2, RP, AX, Z;
20 PetscReal dp = 0.0, final, tau, tau2, theta, theta2, c, F, NV, vv;
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 P = ksp->work[1];
30 PH = ksp->work[2];
31 V = ksp->work[3];
32 D2 = ksp->work[4];
33 X2 = ksp->work[5];
34 S = ksp->work[6];
35 SH = ksp->work[7];
36 T = ksp->work[8];
37 D = ksp->work[9];
38 S2 = ksp->work[10];
39 RP = ksp->work[11];
40 AX = ksp->work[12];
41 Z = ksp->work[13];
42
43 /* Only supports right preconditioning */
44 PetscCheck(ksp->pc_side == PC_RIGHT, PetscObjectComm((PetscObject)ksp), PETSC_ERR_SUP, "KSP qmrcgs does not support %s", PCSides[ksp->pc_side]);
45 if (!ksp->guess_zero) {
46 if (!bcgs->guess) PetscCall(VecDuplicate(X, &bcgs->guess));
47 PetscCall(VecCopy(X, bcgs->guess));
48 } else {
49 PetscCall(VecSet(X, 0.0));
50 }
51
52 /* Compute initial residual */
53 PetscCall(KSPGetPC(ksp, &pc));
54 PetscCall(PCGetOperators(pc, &mat, NULL));
55 if (!ksp->guess_zero) {
56 PetscCall(KSP_MatMult(ksp, mat, X, S2));
57 PetscCall(VecCopy(B, R));
58 PetscCall(VecAXPY(R, -1.0, S2));
59 } else {
60 PetscCall(VecCopy(B, R));
61 }
62
63 /* Test for nothing to do */
64 if (ksp->normtype != KSP_NORM_NONE) PetscCall(VecNorm(R, NORM_2, &dp));
65 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
66 ksp->its = 0;
67 ksp->rnorm = dp;
68 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
69 PetscCall(KSPLogResidualHistory(ksp, dp));
70 PetscCall(KSPMonitor(ksp, 0, dp));
71 PetscCall((*ksp->converged)(ksp, 0, dp, &ksp->reason, ksp->cnvP));
72 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
73
74 /* Make the initial Rp == R */
75 PetscCall(VecCopy(R, RP));
76
77 eta = 1.0;
78 theta = 1.0;
79 if (dp == 0.0) {
80 PetscCall(VecNorm(R, NORM_2, &tau));
81 } else {
82 tau = dp;
83 }
84
85 PetscCall(VecDot(RP, RP, &rho1));
86 PetscCall(VecCopy(R, P));
87
88 i = 0;
89 do {
90 PetscCall(KSP_PCApply(ksp, P, PH)); /* ph <- K p */
91 PetscCall(KSP_MatMult(ksp, mat, PH, V)); /* v <- A ph */
92
93 PetscCall(VecDot(V, RP, &rho2)); /* rho2 <- (v,rp) */
94 if (rho2 == 0.0) {
95 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to division by zero");
96 ksp->reason = KSP_DIVERGED_NANORINF;
97 break;
98 }
99
100 if (rho1 == 0) {
101 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has stagnated");
102 ksp->reason = KSP_DIVERGED_BREAKDOWN; /* Stagnation */
103 break;
104 }
105
106 alpha = rho1 / rho2;
107 PetscCall(VecWAXPY(S, -alpha, V, R)); /* s <- r - alpha v */
108
109 /* First quasi-minimization step */
110 PetscCall(VecNorm(S, NORM_2, &F)); /* f <- norm(s) */
111 theta2 = F / tau;
112
113 c = 1.0 / PetscSqrtReal(1.0 + theta2 * theta2);
114
115 tau2 = tau * theta2 * c;
116 eta2 = c * c * alpha;
117 cf = eta * theta * theta / alpha;
118 PetscCall(VecWAXPY(D2, cf, D, PH)); /* d2 <- ph + cf d */
119 PetscCall(VecWAXPY(X2, eta2, D2, X)); /* x2 <- x + eta2 d2 */
120
121 /* Apply the right preconditioner again */
122 PetscCall(KSP_PCApply(ksp, S, SH)); /* sh <- K s */
123 PetscCall(KSP_MatMult(ksp, mat, SH, T)); /* t <- A sh */
124
125 PetscCall(VecDotNorm2(S, T, &uu, &vv));
126 if (vv == 0.0) {
127 PetscCall(VecDot(S, S, &uu));
128 if (uu != 0.0) {
129 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to division by zero");
130 ksp->reason = KSP_DIVERGED_NANORINF;
131 break;
132 }
133 PetscCall(VecAXPY(X, alpha, SH)); /* x <- x + alpha sh */
134 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
135 ksp->its++;
136 ksp->rnorm = 0.0;
137 ksp->reason = KSP_CONVERGED_RTOL;
138 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
139 PetscCall(KSPLogResidualHistory(ksp, dp));
140 PetscCall(KSPMonitor(ksp, i + 1, 0.0));
141 break;
142 }
143 PetscCall(VecNorm(V, NORM_2, &NV)); /* nv <- norm(v) */
144
145 if (NV == 0) {
146 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to singular matrix");
147 ksp->reason = KSP_DIVERGED_BREAKDOWN;
148 break;
149 }
150
151 if (uu == 0) {
152 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has stagnated");
153 ksp->reason = KSP_DIVERGED_BREAKDOWN; /* Stagnation */
154 break;
155 }
156 omega = uu / vv; /* omega <- uu/vv; */
157
158 /* Computing the residual */
159 PetscCall(VecWAXPY(R, -omega, T, S)); /* r <- s - omega t */
160
161 /* Second quasi-minimization step */
162 if (ksp->normtype != KSP_NORM_NONE && ksp->chknorm < i + 2) PetscCall(VecNorm(R, NORM_2, &dp));
163
164 if (tau2 == 0) {
165 PetscCheck(!ksp->errorifnotconverged, PetscObjectComm((PetscObject)ksp), PETSC_ERR_NOT_CONVERGED, "KSPSolve has not converged due to division by zero");
166 ksp->reason = KSP_DIVERGED_NANORINF;
167 break;
168 }
169 theta = dp / tau2;
170 c = 1.0 / PetscSqrtReal(1.0 + theta * theta);
171 if (dp == 0.0) {
172 PetscCall(VecNorm(R, NORM_2, &tau));
173 } else {
174 tau = dp;
175 }
176 tau = tau * c;
177 eta = c * c * omega;
178
179 cf1 = eta2 * theta2 * theta2 / omega;
180 PetscCall(VecWAXPY(D, cf1, D2, SH)); /* d <- sh + cf1 d2 */
181 PetscCall(VecWAXPY(X, eta, D, X2)); /* x <- x2 + eta d */
182
183 PetscCall(VecDot(R, RP, &rho2));
184 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)ksp));
185 ksp->its++;
186 ksp->rnorm = dp;
187 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)ksp));
188 PetscCall(KSPLogResidualHistory(ksp, dp));
189 PetscCall(KSPMonitor(ksp, i + 1, dp));
190
191 beta = (alpha * rho2) / (omega * rho1);
192 PetscCall(VecAXPBYPCZ(P, 1.0, -omega * beta, beta, R, V)); /* p <- r - omega * beta* v + beta * p */
193 rho1 = rho2;
194 PetscCall(KSP_MatMult(ksp, mat, X, AX)); /* Ax <- A x */
195 PetscCall(VecWAXPY(Z, -1.0, AX, B)); /* r <- b - Ax */
196 PetscCall(VecNorm(Z, NORM_2, &final));
197 PetscCall((*ksp->converged)(ksp, i + 1, dp, &ksp->reason, ksp->cnvP));
198 if (ksp->reason) break;
199 i++;
200 } while (i < ksp->max_it);
201
202 /* mark lack of convergence */
203 if (ksp->its >= ksp->max_it && !ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
204 PetscFunctionReturn(PETSC_SUCCESS);
205 }
206
207 /*MC
208 KSPQMRCGS - Implements the QMRCGStab method {cite}`chan1994qmrcgs` and Ghai, Lu, and Jiao (NLAA 2019).
209
210 Level: beginner
211
212 Note:
213 Only right preconditioning is supported.
214
215 Contributed by:
216 Xiangmin Jiao (xiangmin.jiao@stonybrook.edu)
217
218 .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPBICG`, `KSPFBICGS`, `KSPBCGSL`, `KSPSetPCSide()`
219 M*/
KSPCreate_QMRCGS(KSP ksp)220 PETSC_EXTERN PetscErrorCode KSPCreate_QMRCGS(KSP ksp)
221 {
222 KSP_BCGS *bcgs;
223 static const char citations[] = "@article{chan1994qmrcgs,\n"
224 " title={A quasi-minimal residual variant of the {Bi-CGSTAB} algorithm for nonsymmetric systems},\n"
225 " author={Chan, Tony F and Gallopoulos, Efstratios and Simoncini, Valeria and Szeto, Tedd and Tong, Charles H},\n"
226 " journal={SIAM Journal on Scientific Computing},\n"
227 " volume={15},\n"
228 " number={2},\n"
229 " pages={338--347},\n"
230 " year={1994},\n"
231 " publisher={SIAM}\n"
232 "}\n"
233 "@article{ghai2019comparison,\n"
234 " title={A comparison of preconditioned {K}rylov subspace methods for large-scale nonsymmetric linear systems},\n"
235 " author={Ghai, Aditi and Lu, Cao and Jiao, Xiangmin},\n"
236 " journal={Numerical Linear Algebra with Applications},\n"
237 " volume={26},\n"
238 " number={1},\n"
239 " pages={e2215},\n"
240 " year={2019},\n"
241 " publisher={Wiley Online Library}\n"
242 "}\n";
243 PetscBool cite = PETSC_FALSE;
244
245 PetscFunctionBegin;
246 PetscCall(PetscCitationsRegister(citations, &cite));
247 PetscCall(PetscNew(&bcgs));
248
249 ksp->data = bcgs;
250 ksp->ops->setup = KSPSetUp_QMRCGS;
251 ksp->ops->solve = KSPSolve_QMRCGS;
252 ksp->ops->destroy = KSPDestroy_BCGS;
253 ksp->ops->reset = KSPReset_BCGS;
254 ksp->ops->buildresidual = KSPBuildResidualDefault;
255 ksp->ops->setfromoptions = KSPSetFromOptions_BCGS;
256 ksp->pc_side = PC_RIGHT; /* set default PC side */
257
258 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
259 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
260 PetscFunctionReturn(PETSC_SUCCESS);
261 }
262