Actual source code: tcqmr.c
2: /*
3: This file contains an implementation of Tony Chan's transpose-free QMR.
5: Note: The vector dot products in the code have not been checked for the
6: complex numbers version, so most probably some are incorrect.
7: */
9: #include <../src/ksp/ksp/impls/tcqmr/tcqmrimpl.h>
11: static PetscErrorCode KSPSolve_TCQMR(KSP ksp)
12: {
13: PetscReal rnorm0, rnorm, dp1, Gamma;
14: PetscScalar theta, ep, cl1, sl1, cl, sl, sprod, tau_n1, f;
15: PetscScalar deltmp, rho, beta, eptmp, ta, s, c, tau_n, delta;
16: PetscScalar dp11, dp2, rhom1, alpha, tmp;
18: ksp->its = 0;
20: KSPInitialResidual(ksp, x, u, v, r, b);
21: VecNorm(r, NORM_2, &rnorm0); /* rnorm0 = ||r|| */
22: KSPCheckNorm(ksp, rnorm0);
23: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm0;
24: else ksp->rnorm = 0;
25: (*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP);
26: if (ksp->reason) return 0;
28: VecSet(um1, 0.0);
29: VecCopy(r, u);
30: rnorm = rnorm0;
31: tmp = 1.0 / rnorm;
32: VecScale(u, tmp);
33: VecSet(vm1, 0.0);
34: VecCopy(u, v);
35: VecCopy(u, v0);
36: VecSet(pvec1, 0.0);
37: VecSet(pvec2, 0.0);
38: VecSet(p, 0.0);
39: theta = 0.0;
40: ep = 0.0;
41: cl1 = 0.0;
42: sl1 = 0.0;
43: cl = 0.0;
44: sl = 0.0;
45: sprod = 1.0;
46: tau_n1 = rnorm0;
47: f = 1.0;
48: Gamma = 1.0;
49: rhom1 = 1.0;
51: /*
52: CALCULATE SQUARED LANCZOS vectors
53: */
54: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm;
55: else ksp->rnorm = 0;
56: (*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP);
57: while (!ksp->reason) {
58: KSPMonitor(ksp, ksp->its, ksp->rnorm);
59: ksp->its++;
61: KSP_PCApplyBAorAB(ksp, u, y, vtmp); /* y = A*u */
62: VecDot(y, v0, &dp11);
63: KSPCheckDot(ksp, dp11);
64: VecDot(u, v0, &dp2);
65: alpha = dp11 / dp2; /* alpha = v0'*y/v0'*u */
66: deltmp = alpha;
67: VecCopy(y, z);
68: VecAXPY(z, -alpha, u); /* z = y - alpha u */
69: VecDot(u, v0, &rho);
70: beta = rho / (f * rhom1);
71: rhom1 = rho;
72: VecCopy(z, utmp); /* up1 = (A-alpha*I)*
73: (z-2*beta*p) + f*beta*
74: beta*um1 */
75: VecAXPY(utmp, -2.0 * beta, p);
76: KSP_PCApplyBAorAB(ksp, utmp, up1, vtmp);
77: VecAXPY(up1, -alpha, utmp);
78: VecAXPY(up1, f * beta * beta, um1);
79: VecNorm(up1, NORM_2, &dp1);
80: KSPCheckNorm(ksp, dp1);
81: f = 1.0 / dp1;
82: VecScale(up1, f);
83: VecAYPX(p, -beta, z); /* p = f*(z-beta*p) */
84: VecScale(p, f);
85: VecCopy(u, um1);
86: VecCopy(up1, u);
87: beta = beta / Gamma;
88: eptmp = beta;
89: KSP_PCApplyBAorAB(ksp, v, vp1, vtmp);
90: VecAXPY(vp1, -alpha, v);
91: VecAXPY(vp1, -beta, vm1);
92: VecNorm(vp1, NORM_2, &Gamma);
93: KSPCheckNorm(ksp, Gamma);
94: VecScale(vp1, 1.0 / Gamma);
95: VecCopy(v, vm1);
96: VecCopy(vp1, v);
98: /*
99: SOLVE Ax = b
100: */
101: /* Apply last two Given's (Gl-1 and Gl) rotations to (beta,alpha,Gamma) */
102: if (ksp->its > 2) {
103: theta = sl1 * beta;
104: eptmp = -cl1 * beta;
105: }
106: if (ksp->its > 1) {
107: ep = -cl * eptmp + sl * alpha;
108: deltmp = -sl * eptmp - cl * alpha;
109: }
110: if (PetscAbsReal(Gamma) > PetscAbsScalar(deltmp)) {
111: ta = -deltmp / Gamma;
112: s = 1.0 / PetscSqrtScalar(1.0 + ta * ta);
113: c = s * ta;
114: } else {
115: ta = -Gamma / deltmp;
116: c = 1.0 / PetscSqrtScalar(1.0 + ta * ta);
117: s = c * ta;
118: }
120: delta = -c * deltmp + s * Gamma;
121: tau_n = -c * tau_n1;
122: tau_n1 = -s * tau_n1;
123: VecCopy(vm1, pvec);
124: VecAXPY(pvec, -theta, pvec2);
125: VecAXPY(pvec, -ep, pvec1);
126: VecScale(pvec, 1.0 / delta);
127: VecAXPY(x, tau_n, pvec);
128: cl1 = cl;
129: sl1 = sl;
130: cl = c;
131: sl = s;
133: VecCopy(pvec1, pvec2);
134: VecCopy(pvec, pvec1);
136: /* Compute the upper bound on the residual norm r (See QMR paper p. 13) */
137: sprod = sprod * PetscAbsScalar(s);
138: rnorm = rnorm0 * PetscSqrtReal((PetscReal)ksp->its + 2.0) * PetscRealPart(sprod);
139: if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm;
140: else ksp->rnorm = 0;
141: (*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP);
142: if (ksp->its >= ksp->max_it) {
143: if (!ksp->reason) ksp->reason = KSP_DIVERGED_ITS;
144: break;
145: }
146: }
147: KSPMonitor(ksp, ksp->its, ksp->rnorm);
148: KSPUnwindPreconditioner(ksp, x, vtmp);
149: return 0;
150: }
152: static PetscErrorCode KSPSetUp_TCQMR(KSP ksp)
153: {
154: KSPSetWorkVecs(ksp, TCQMR_VECS);
155: return 0;
156: }
158: /*MC
159: KSPTCQMR - A variant of QMR (quasi minimal residual) [1]
161: Level: beginner
163: Notes:
164: Supports either left or right preconditioning, but not symmetric
166: The "residual norm" computed in this algorithm is actually just an upper bound on the actual residual norm.
167: That is for left preconditioning it is a bound on the preconditioned residual and for right preconditioning
168: it is a bound on the true residual.
170: References:
171: . [1] - Tony F. Chan, Lisette de Pillis, and Henk van der Vorst, Transpose free formulations of Lanczos type methods for nonsymmetric linear systems,
172: Numerical Algorithms, Volume 17, 1998.
174: .seealso: [](chapter_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPTFQMR`
175: M*/
177: PETSC_EXTERN PetscErrorCode KSPCreate_TCQMR(KSP ksp)
178: {
179: KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3);
180: KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2);
181: KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1);
183: ksp->data = (void *)0;
184: ksp->ops->buildsolution = KSPBuildSolutionDefault;
185: ksp->ops->buildresidual = KSPBuildResidualDefault;
186: ksp->ops->setup = KSPSetUp_TCQMR;
187: ksp->ops->solve = KSPSolve_TCQMR;
188: ksp->ops->destroy = KSPDestroyDefault;
189: ksp->ops->setfromoptions = NULL;
190: ksp->ops->view = NULL;
191: return 0;
192: }