1 /*
2 This file contains an implementation of Tony Chan's transpose-free QMR.
3
4 Note: The vector dot products in the code have not been checked for the
5 complex numbers version, so most probably some are incorrect.
6 */
7
8 #include <../src/ksp/ksp/impls/tcqmr/tcqmrimpl.h>
9
KSPSolve_TCQMR(KSP ksp)10 static PetscErrorCode KSPSolve_TCQMR(KSP ksp)
11 {
12 PetscReal rnorm0, rnorm, dp1, Gamma;
13 PetscScalar theta, ep, cl1, sl1, cl, sl, sprod, tau_n1, f;
14 PetscScalar deltmp, rho, beta, eptmp, ta, s, c, tau_n, delta;
15 PetscScalar dp11, dp2, rhom1, alpha, tmp;
16
17 PetscFunctionBegin;
18 ksp->its = 0;
19
20 PetscCall(KSPInitialResidual(ksp, x, u, v, r, b));
21 PetscCall(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 PetscCall((*ksp->converged)(ksp, 0, ksp->rnorm, &ksp->reason, ksp->cnvP));
26 if (ksp->reason) PetscFunctionReturn(PETSC_SUCCESS);
27
28 PetscCall(VecSet(um1, 0.0));
29 PetscCall(VecCopy(r, u));
30 rnorm = rnorm0;
31 tmp = 1.0 / rnorm;
32 PetscCall(VecScale(u, tmp));
33 PetscCall(VecSet(vm1, 0.0));
34 PetscCall(VecCopy(u, v));
35 PetscCall(VecCopy(u, v0));
36 PetscCall(VecSet(pvec1, 0.0));
37 PetscCall(VecSet(pvec2, 0.0));
38 PetscCall(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;
50
51 /*
52 CALCULATE SQUARED LANCZOS vectors
53 */
54 if (ksp->normtype != KSP_NORM_NONE) ksp->rnorm = rnorm;
55 else ksp->rnorm = 0;
56 PetscCall((*ksp->converged)(ksp, ksp->its, ksp->rnorm, &ksp->reason, ksp->cnvP));
57 while (!ksp->reason) {
58 PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
59 ksp->its++;
60
61 PetscCall(KSP_PCApplyBAorAB(ksp, u, y, vtmp)); /* y = A*u */
62 PetscCall(VecDot(y, v0, &dp11));
63 KSPCheckDot(ksp, dp11);
64 PetscCall(VecDot(u, v0, &dp2));
65 alpha = dp11 / dp2; /* alpha = v0'*y/v0'*u */
66 deltmp = alpha;
67 PetscCall(VecCopy(y, z));
68 PetscCall(VecAXPY(z, -alpha, u)); /* z = y - alpha u */
69 PetscCall(VecDot(u, v0, &rho));
70 beta = rho / (f * rhom1);
71 rhom1 = rho;
72 PetscCall(VecCopy(z, utmp)); /* up1 = (A-alpha*I)*
73 (z-2*beta*p) + f*beta*
74 beta*um1 */
75 PetscCall(VecAXPY(utmp, -2.0 * beta, p));
76 PetscCall(KSP_PCApplyBAorAB(ksp, utmp, up1, vtmp));
77 PetscCall(VecAXPY(up1, -alpha, utmp));
78 PetscCall(VecAXPY(up1, f * beta * beta, um1));
79 PetscCall(VecNorm(up1, NORM_2, &dp1));
80 KSPCheckNorm(ksp, dp1);
81 f = 1.0 / dp1;
82 PetscCall(VecScale(up1, f));
83 PetscCall(VecAYPX(p, -beta, z)); /* p = f*(z-beta*p) */
84 PetscCall(VecScale(p, f));
85 PetscCall(VecCopy(u, um1));
86 PetscCall(VecCopy(up1, u));
87 beta = beta / Gamma;
88 eptmp = beta;
89 PetscCall(KSP_PCApplyBAorAB(ksp, v, vp1, vtmp));
90 PetscCall(VecAXPY(vp1, -alpha, v));
91 PetscCall(VecAXPY(vp1, -beta, vm1));
92 PetscCall(VecNorm(vp1, NORM_2, &Gamma));
93 KSPCheckNorm(ksp, Gamma);
94 PetscCall(VecScale(vp1, 1.0 / Gamma));
95 PetscCall(VecCopy(v, vm1));
96 PetscCall(VecCopy(vp1, v));
97
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 }
119
120 delta = -c * deltmp + s * Gamma;
121 tau_n = -c * tau_n1;
122 tau_n1 = -s * tau_n1;
123 PetscCall(VecCopy(vm1, pvec));
124 PetscCall(VecAXPY(pvec, -theta, pvec2));
125 PetscCall(VecAXPY(pvec, -ep, pvec1));
126 PetscCall(VecScale(pvec, 1.0 / delta));
127 PetscCall(VecAXPY(x, tau_n, pvec));
128 cl1 = cl;
129 sl1 = sl;
130 cl = c;
131 sl = s;
132
133 PetscCall(VecCopy(pvec1, pvec2));
134 PetscCall(VecCopy(pvec, pvec1));
135
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 PetscCall((*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 PetscCall(KSPMonitor(ksp, ksp->its, ksp->rnorm));
148 PetscCall(KSPUnwindPreconditioner(ksp, x, vtmp));
149 PetscFunctionReturn(PETSC_SUCCESS);
150 }
151
KSPSetUp_TCQMR(KSP ksp)152 static PetscErrorCode KSPSetUp_TCQMR(KSP ksp)
153 {
154 PetscFunctionBegin;
155 PetscCall(KSPSetWorkVecs(ksp, TCQMR_VECS));
156 PetscFunctionReturn(PETSC_SUCCESS);
157 }
158
159 /*MC
160 KSPTCQMR - A variant of a tranpose-free QMR (quasi minimal residual) {cite}`chan1998transpose` algorithm
161
162 Level: beginner
163
164 Notes:
165 Supports either left or right preconditioning, but not symmetric
166
167 The "residual norm" computed in this algorithm is actually just an upper bound on the actual residual norm.
168 For left preconditioning it is a bound on the preconditioned residual norm and for right preconditioning
169 it is a bound on the true residual norm.
170
171 This algorithm is related to the standard transpose-free QMR implemented by `KSPTFQMR`.
172
173 .seealso: [](ch_ksp), `KSPCreate()`, `KSPSetType()`, `KSPType`, `KSP`, `KSPTFQMR`
174 M*/
175
KSPCreate_TCQMR(KSP ksp)176 PETSC_EXTERN PetscErrorCode KSPCreate_TCQMR(KSP ksp)
177 {
178 PetscFunctionBegin;
179 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_PRECONDITIONED, PC_LEFT, 3));
180 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_UNPRECONDITIONED, PC_RIGHT, 2));
181 PetscCall(KSPSetSupportedNorm(ksp, KSP_NORM_NONE, PC_RIGHT, 1));
182
183 ksp->data = NULL;
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 PetscFunctionReturn(PETSC_SUCCESS);
192 }
193