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