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