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