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