1 // Copyright (c) 2018 ERGO-Code. See license.txt for license.
2
3 #include "conjugate_residuals.h"
4 #include <algorithm>
5 #include <cmath>
6 #include "timer.h"
7 #include "utils.h"
8
9 namespace ipx {
10
ConjugateResiduals(const Control & control)11 ConjugateResiduals::ConjugateResiduals(const Control& control) :
12 control_(control) {}
13
Solve(LinearOperator & C,const Vector & rhs,double tol,const double * resscale,Int maxiter,Vector & lhs)14 void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs,
15 double tol, const double* resscale, Int maxiter,
16 Vector& lhs) {
17 const Int m = rhs.size();
18 Vector residual(m); // rhs - C*lhs
19 Vector step(m); // update to lhs
20 Vector Cresidual(m); // C * residual
21 Vector Cstep(m); // C * step
22 double cdot = 0.0; // dot product from C.Apply
23 Timer timer;
24
25 errflag_ = 0;
26 iter_ = 0;
27 time_ = 0.0;
28 if (maxiter < 0)
29 maxiter = m+100;
30
31 // Initialize residual, step and Cstep.
32 if (Infnorm(lhs) == 0.0) {
33 residual = rhs; // saves a matrix-vector op
34 } else {
35 C.Apply(lhs, residual, nullptr);
36 residual = rhs-residual;
37 }
38 C.Apply(residual, Cresidual, &cdot);
39 step = residual;
40 Cstep = Cresidual;
41
42 while (true) {
43 // Termination check.
44 double resnorm = 0.0;
45 if (resscale)
46 for (Int i = 0; i < m; i++)
47 resnorm = std::max(resnorm, std::abs(resscale[i]*residual[i]));
48 else
49 resnorm = Infnorm(residual);
50 if (resnorm <= tol)
51 break;
52 if (iter_ == maxiter) {
53 control_.Debug(3)
54 << " CR method not converged in " << maxiter << " iterations."
55 << " residual = " << sci2(resnorm) << ','
56 << " tolerance = " << sci2(tol) << '\n';
57 errflag_ = IPX_ERROR_cr_iter_limit;
58 break;
59 }
60 if (cdot <= 0.0) {
61 errflag_ = IPX_ERROR_cr_matrix_not_posdef;
62 break;
63 }
64
65 // Update lhs, residual and Cresidual.
66 const double denom = Dot(Cstep,Cstep);
67 const double alpha = cdot/denom;
68 if (!std::isfinite(alpha)) {
69 errflag_ = IPX_ERROR_cr_inf_or_nan;
70 break;
71 }
72 lhs += alpha*step;
73 residual -= alpha*Cstep;
74 double cdotnew;
75 C.Apply(residual, Cresidual, &cdotnew);
76
77 // Update step and Cstep.
78 const double beta = cdotnew/cdot;
79 step = residual + beta*step;
80 Cstep = Cresidual + beta*Cstep;
81 cdot = cdotnew;
82
83 iter_++;
84 if ((errflag_ = control_.InterruptCheck()) != 0)
85 break;
86 }
87 time_ = timer.Elapsed();
88 }
89
Solve(LinearOperator & C,LinearOperator & P,const Vector & rhs,double tol,const double * resscale,Int maxiter,Vector & lhs)90 void ConjugateResiduals::Solve(LinearOperator& C, LinearOperator& P,
91 const Vector& rhs, double tol,
92 const double* resscale, Int maxiter, Vector& lhs){
93 const Int m = rhs.size();
94 Vector residual(m); // rhs - C*lhs
95 Vector sresidual(m); // preconditioned residual
96 Vector step(m); // update to lhs
97 Vector Csresidual(m); // C * sresidual
98 Vector Cstep(m); // C * step
99 double cdot = 0.0; // dot product from C.Apply
100 Timer timer;
101
102 // resnorm_precond_system = residual' * P * residual.
103 // This quantity is minimized by the preconditioned CR method over a Krylov
104 // subspace. Hence (in theory) must decrease strictly monotonically. If it
105 // does not, then the method stagnated due to round-off errors. This
106 // happened in a few cases with augmented diagonal preconditioning (i.e.
107 // diagonal preconditioning with dense columns treated as low-rank update)
108 // because operations with P were not sufficiently accurate.
109 double resnorm_precond_system = 0.0;
110
111 errflag_ = 0;
112 iter_ = 0;
113 time_ = 0.0;
114 if (maxiter < 0)
115 maxiter = m+100;
116
117 // Initialize residual, sresidual, step and Cstep.
118 if (Infnorm(lhs) == 0.0) {
119 residual = rhs; // saves a matrix-vector op
120 } else {
121 C.Apply(lhs, residual, nullptr);
122 residual = rhs-residual;
123 }
124 P.Apply(residual, sresidual, &resnorm_precond_system);
125 C.Apply(sresidual, Csresidual, &cdot);
126 step = sresidual;
127 Cstep = Csresidual;
128
129 while (true) {
130 // Termination check.
131 double resnorm = 0.0;
132 if (resscale)
133 for (Int i = 0; i < m; i++)
134 resnorm = std::max(resnorm, std::abs(resscale[i]*residual[i]));
135 else
136 resnorm = Infnorm(residual);
137 if (resnorm <= tol)
138 break;
139 if (iter_ == maxiter) {
140 control_.Debug(3)
141 << " PCR method not converged in " << maxiter << " iterations."
142 << " residual = " << sci2(resnorm) << ','
143 << " tolerance = " << sci2(tol) << '\n';
144 errflag_ = IPX_ERROR_cr_iter_limit;
145 break;
146 }
147 if (cdot <= 0.0) {
148 control_.Debug(3)
149 << " matrix in PCR method not posdef. cdot = " << sci2(cdot)
150 << ", infnorm(sresidual) = " << sci2(Infnorm(sresidual))
151 << ", infnorm(residual) = " << sci2(Infnorm(residual))
152 << '\n';
153 errflag_ = IPX_ERROR_cr_matrix_not_posdef;
154 break;
155 }
156
157 // Update lhs, residual, sresidual and Csresidual.
158 double cdotnew;
159 {
160 // Uses Csresidual as storage for preconditioned Cstep.
161 Vector& precond_Cstep = Csresidual;
162 double pdot;
163 P.Apply(Cstep, precond_Cstep, &pdot);
164 if (pdot <= 0.0) {
165 errflag_ = IPX_ERROR_cr_precond_not_posdef;
166 break;
167 }
168 const double alpha = cdot/pdot;
169 if (!std::isfinite(alpha)) {
170 errflag_ = IPX_ERROR_cr_inf_or_nan;
171 break;
172 }
173 lhs += alpha*step;
174 residual -= alpha*Cstep;
175 sresidual -= alpha*precond_Cstep;
176 C.Apply(sresidual, Csresidual, &cdotnew);
177 // Now Csresidual is restored and alias goes out of scope.
178 }
179
180 // Update step and Cstep.
181 const double beta = cdotnew/cdot;
182 step = sresidual + beta*step;
183 Cstep = Csresidual + beta*Cstep;
184 cdot = cdotnew;
185
186 iter_++;
187 if (iter_%5 == 0) {
188 // We recompute the preconditioned residual from its definition
189 // all 5 iterations. Otherwise it happened that the preconditioned
190 // residual approached zero but the true residual stagnated. This
191 // was caused by the update to the preconditioned residual being not
192 // sufficiently accurate.
193 // As a second safeguard, we check that resnorm_precond_system
194 // decreased during the last 5 iterations.
195 double rsdot;
196 P.Apply(residual, sresidual, &rsdot);
197 if (rsdot >= resnorm_precond_system) {
198 control_.Debug(3)
199 << " resnorm_precond_system old = "
200 << sci2(resnorm_precond_system) << '\n'
201 << " resnorm_precond_system new = "
202 << sci2(rsdot) << '\n';
203 errflag_ = IPX_ERROR_cr_no_progress;
204 break;
205 }
206 resnorm_precond_system = rsdot;
207 }
208
209 if ((errflag_ = control_.InterruptCheck()) != 0)
210 break;
211 }
212 time_ = timer.Elapsed();
213 }
214
errflag() const215 Int ConjugateResiduals::errflag() const { return errflag_; }
iter() const216 Int ConjugateResiduals::iter() const { return iter_; }
time() const217 double ConjugateResiduals::time() const { return time_; }
218
219 } // namespace ipx
220