1 // This is core/vnl/algo/vnl_lsqr.h
2 #ifndef vnl_lsqr_h_
3 #define vnl_lsqr_h_
4 //:
5 // \file
6 // \brief Linear least squares
7 //
8 // vnl_lsqr implements an algorithm for large, sparse linear systems and
9 // sparse, linear least squares. It is a wrapper for the LSQR algorithm
10 // of Paige and Saunders (ACM TOMS 583). The sparse system is encapsulated
11 // by a vnl_linear_system.
12 //
13 // \author David Capel, capes@robots
14 // \date   July 2000
15 //
16 // \verbatim
17 //  Modifications
18 //   000705 capes@robots initial version.
19 //   4/4/01 LSB (Manchester) Documentation tidied
20 // \endverbatim
21 //-----------------------------------------------------------------------------
22 
23 #include <iosfwd>
24 #include <vnl/vnl_vector.h>
25 #include <vnl/vnl_linear_system.h>
26 #ifdef _MSC_VER
27 #  include <vcl_msvc_warnings.h>
28 #endif
29 
30 #include <vnl/algo/vnl_algo_export.h>
31 
32 //: Linear least squares
33 //  vnl_lsqr implements an algorithm for large, sparse linear systems and
34 //  sparse, linear least squares. It is a wrapper for the LSQR algorithm
35 //  of Paige and Saunders (ACM TOMS 583). The sparse system is encapsulated
36 //  by a vnl_linear_system.
37 
38 class VNL_ALGO_EXPORT vnl_lsqr
39 {
40  public:
vnl_lsqr(vnl_linear_system & ls)41   vnl_lsqr(vnl_linear_system& ls) :
42     ls_(&ls), max_iter_(4*ls.get_number_of_unknowns()) {}
43 
44   ~vnl_lsqr();
45 
set_max_iterations(long max_iter)46   void set_max_iterations(long max_iter) { max_iter_ = max_iter; }
47 
48   //: Perform the minimization starting from x=0 and putting the result into x.
49   // Return code may be translated with translate_return_code(), or the result of the
50   // minimization may be printed in more detail with diagnose_outcome()
51   int minimize(vnl_vector<double>& x);
52 
get_number_of_iterations()53   long get_number_of_iterations() const { return num_iter_; }
54 
55   //: Pontificate about the outcome of the last minimization.
56   void diagnose_outcome(std::ostream& os) const;
57 
58   static void translate_return_code(std::ostream& os, int return_code);
59 
60   //: Return the residual norm estimate:
get_resid_norm_estimate()61   double get_resid_norm_estimate() const { return resid_norm_estimate_; }
62 
63   //: Get the return code for the last minimization
return_code()64   inline long return_code() const { return return_code_; }
65 
66  protected:
67   vnl_linear_system* ls_;
68   long max_iter_;
69   long num_iter_;
70   double resid_norm_estimate_;
71   double result_norm_estimate_;
72   double A_condition_estimate_;
73   double result_norm_;
74   long return_code_;
75 
76   static int aprod_(const long* mode, const long* m, const long* n, double* x, double* y,
77                     long* leniw, long* lenrw, long* iw, double* rw,
78                     void* userdata);
79 };
80 
81 #endif // vnl_lsqr_h_
82