1 // This is core/vnl/algo/vnl_lsqr.cxx
2 //
3 // vnl_lsqr
4 // Author: David Capel
5 // Created: July 2000
6 //
7 //-----------------------------------------------------------------------------
8
9 #include <vector>
10 #include <iostream>
11 #include <algorithm>
12 #include "vnl_lsqr.h"
13 #include <vnl/vnl_vector_ref.h>
14
15 #include "lsqrBase.h"
16
17 class lsqrVNL : public lsqrBase
18 {
19 public:
20
lsqrVNL()21 lsqrVNL()
22 {
23 this->ls_ = nullptr;
24 }
25
26 ~lsqrVNL() override = default;
27
28 /**
29 * computes y = y + A*x without altering x,
30 * where A is a matrix of dimensions A[m][n].
31 * The size of the vector x is n.
32 * The size of the vector y is m.
33 */
Aprod1(unsigned int m,unsigned int n,const double * x,double * y) const34 void Aprod1(unsigned int m, unsigned int n, const double * x, double * y ) const override
35 {
36 vnl_vector_ref<double> x_ref(n, const_cast<double*>(x) );
37 vnl_vector_ref<double> y_ref(m,y);
38
39 vnl_vector_ref<double> tmp(m,rw);
40 this->ls_->multiply(x_ref, tmp);
41 y_ref += tmp;
42 }
43
44
45 /**
46 * computes x = x + A'*y without altering y,
47 * where A is a matrix of dimensions A[m][n].
48 * The size of the vector x is n.
49 * The size of the vector y is m.
50 */
Aprod2(unsigned int m,unsigned int n,double * x,const double * y) const51 void Aprod2(unsigned int m, unsigned int n, double * x, const double * y ) const override
52 {
53 vnl_vector_ref<double> x_ref(n,x);
54 vnl_vector_ref<double> y_ref(m, const_cast<double*>( y ) );
55
56 vnl_vector_ref<double> tmp(n,this->rw);
57 this->ls_->transpose_multiply(y_ref, tmp);
58 x_ref += tmp;
59 }
60
61 /** Set the linear system to be solved A*x = b. */
SetLinearSystem(vnl_linear_system * inls)62 void SetLinearSystem( vnl_linear_system * inls )
63 {
64 this->ls_ = inls;
65 }
66
SetWorkingSpace(double * inrw)67 void SetWorkingSpace( double * inrw )
68 {
69 this->rw = inrw;
70 }
71
72 private:
73
74 vnl_linear_system * ls_;
75
76 double * rw;
77 };
78
79
80 vnl_lsqr::~vnl_lsqr() = default;
81
82 // Requires number_of_residuals() of workspace in rw.
aprod_(const long * mode,const long * m,const long * n,double * x,double * y,long *,long *,long *,double * rw,void * userdata)83 int vnl_lsqr::aprod_(const long* mode, const long* m, const long* n, double* x, double* y, long* /*leniw*/, long* /*lenrw*/, long* /*iw*/, double* rw, void* userdata)
84 {
85 //
86 // THIS CODE IS DEPRECATED
87 // THE FUNCTIONALITY HAS BEEN MOVED TO THE lsqrVNL class above.
88 // THE FUNCTIONS IS CONSERVED HERE ONLY FOR BACKWARD COMPATIBILITY.
89 //
90 auto* self = static_cast<vnl_lsqr*>(userdata);
91
92 // If MODE = 1, compute y = y + A*x.
93 // If MODE = 2, compute x = x + A(transpose)*y.
94
95 vnl_vector_ref<double> x_ref(*n,x);
96 vnl_vector_ref<double> y_ref(*m,y);
97
98 if (*mode == 1) {
99 vnl_vector_ref<double> tmp(*m,rw);
100 self->ls_->multiply(x_ref, tmp);
101 y_ref += tmp;
102 }
103 else {
104 vnl_vector_ref<double> tmp(*n,rw);
105 self->ls_->transpose_multiply(y_ref, tmp);
106 x_ref += tmp;
107 }
108
109 return 0;
110 }
111
minimize(vnl_vector<double> & result)112 int vnl_lsqr::minimize(vnl_vector<double>& result)
113 {
114 long m = ls_->get_number_of_residuals();
115 long n = ls_->get_number_of_unknowns();
116 double damp = 0;
117
118 //NOTE: rw is a scratch space used for both intermediate residual and unknown computations
119 std::vector<double> rw(std::max(m,n));
120 std::vector<double> v(n);
121 std::vector<double> se(n);
122
123 double atol = 0;
124 double btol = 0;
125 #ifdef THIS_CODE_IS_DISABLED_BECAUSE_THE_LSQR_CODE_FROM_NETLIB_WAS_COPYRIGHTED_BY_ACM
126 double conlim = 0;
127 long nout = -1;
128 double acond, rnorm, xnorm;
129 #endif
130 double anorm, arnorm;
131
132 vnl_vector<double> rhs(m);
133 ls_->get_rhs(rhs);
134
135 lsqrVNL solver;
136
137 solver.SetDamp( damp );
138 solver.SetLinearSystem( this->ls_ );
139 solver.SetWorkingSpace( rw.data() );
140 solver.SetMaximumNumberOfIterations( max_iter_ );
141 solver.SetStandardErrorEstimates( se.data() );
142 solver.SetToleranceA( atol );
143 solver.SetToleranceB( btol );
144
145 solver.Solve( m, n, rhs.data_block(), result.data_block() );
146
147 #ifdef THIS_CODE_IS_DISABLED_BECAUSE_THE_LSQR_CODE_FROM_NETLIB_WAS_COPYRIGHTED_BY_ACM
148 v3p_netlib_lsqr_(
149 &m, &n, aprod_, &damp, &leniw, &lenrw, iw, &rw[0],
150 rhs.data_block(), &v[0], &w[0], result.data_block(), &se[0],
151 &atol, &btol, &conlim, &max_iter_, &nout, &return_code_,
152 &num_iter_, &anorm, &acond, &rnorm, &arnorm, &xnorm,
153 this);
154 #endif
155
156 resid_norm_estimate_ = solver.GetFinalEstimateOfNormRbar();
157 result_norm_estimate_ = solver.GetFinalEstimateOfNormOfX();
158 A_condition_estimate_ = solver.GetConditionNumberEstimateOfAbar();
159 return_code_ = solver.GetStoppingReason();
160 num_iter_ = solver.GetNumberOfIterationsPerformed();
161 anorm = solver.GetFrobeniusNormEstimateOfAbar();
162 arnorm = solver.GetFinalEstimateOfNormOfResiduals();
163
164 #ifdef THIS_CODE_IS_DISABLED_BECAUSE_THE_LSQR_CODE_FROM_NETLIB_WAS_COPYRIGHTED_BY_ACM
165 std::cerr << "A Fro norm estimate = " << anorm << std::endl
166 << "A condition estimate = " << acond << std::endl
167 << "Residual norm estimate = " << rnorm << std::endl
168 << "A'(Ax - b) norm estimate = " << arnorm << std::endl
169 << "x norm estimate = " << xnorm << std::endl;
170 #endif
171
172 // We should return the return code, as translate_return_code is public and
173 // it is very misleading that the return code from this function can't be fed
174 // into translate_return_code. (Brian Amberg)
175 return return_code_;
176 }
177
diagnose_outcome(std::ostream & os) const178 void vnl_lsqr::diagnose_outcome(std::ostream& os) const
179 {
180 translate_return_code(os, return_code_);
181 os << __FILE__ " : residual norm estimate = " << resid_norm_estimate_ << std::endl
182 << __FILE__ " : result norm estimate = " << result_norm_estimate_ << std::endl
183 << __FILE__ " : condition no. estimate = " << A_condition_estimate_ << std::endl
184 << __FILE__ " : iterations = " << num_iter_ << std::endl;
185 }
186
translate_return_code(std::ostream & os,int rc)187 void vnl_lsqr::translate_return_code(std::ostream& os, int rc)
188 {
189 const char* vnl_lsqr_reasons[] = {
190 "x = 0 is the exact solution. No iterations were performed.",
191 "The equations A*x = b are probably compatible. "
192 "Norm(A*x - b) is sufficiently small, given the "
193 "values of ATOL and BTOL.",
194 "The system A*x = b is probably not compatible. "
195 "A least-squares solution has been obtained that is "
196 "sufficiently accurate, given the value of ATOL.",
197 "An estimate of cond(Abar) has exceeded CONLIM. "
198 "The system A*x = b appears to be ill-conditioned. "
199 "Otherwise, there could be an error in subroutine APROD.",
200 "The equations A*x = b are probably compatible. "
201 "Norm(A*x - b) is as small as seems reasonable on this machine.",
202 "The system A*x = b is probably not compatible. A least-squares "
203 "solution has been obtained that is as accurate as seems "
204 "reasonable on this machine.",
205 "Cond(Abar) seems to be so large that there is no point in doing further "
206 "iterations, given the precision of this machine. "
207 "There could be an error in subroutine APROD.",
208 "The iteration limit ITNLIM was reached."
209 };
210
211 if
212 (rc < 0 || rc > 7) os << __FILE__ " : Illegal return code : " << rc << std::endl;
213 else
214 os << __FILE__ " : " << vnl_lsqr_reasons[rc] << std::endl;
215 }
216