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