core/vnl/algo/vnl_lsqr.cxx
Go to the documentation of this file.
00001 // This is core/vnl/algo/vnl_lsqr.cxx
00002 #ifdef VCL_NEEDS_PRAGMA_INTERFACE
00003 #pragma implementation
00004 #endif
00005 //
00006 // vnl_lsqr
00007 // Author: David Capel
00008 // Created: July 2000
00009 //
00010 //-----------------------------------------------------------------------------
00011 
00012 #include "vnl_lsqr.h"
00013 #include <vcl_vector.h>
00014 #include <vcl_iostream.h>
00015 #include <vnl/vnl_vector_ref.h>
00016 
00017 #include <vnl/algo/vnl_netlib.h> // lsqr_()
00018 
00019 vnl_lsqr::~vnl_lsqr()
00020 {
00021 }
00022 
00023 // Requires number_of_residuals() of workspace in rw.
00024 int vnl_lsqr::aprod_(long* mode, long* m, long* n, double* x, double* y, long* /*leniw*/, long* /*lenrw*/, long* /*iw*/, double* rw, void* userdata)
00025 {
00026   vnl_lsqr* self = static_cast<vnl_lsqr*>(userdata);
00027 
00028   //  If MODE = 1, compute  y = y + A*x.
00029   //  If MODE = 2, compute  x = x + A(transpose)*y.
00030 
00031   vnl_vector_ref<double> x_ref(*n,x);
00032   vnl_vector_ref<double> y_ref(*m,y);
00033 
00034   if (*mode == 1) {
00035     vnl_vector_ref<double> tmp(*m,rw);
00036     self->ls_->multiply(x_ref, tmp);
00037     y_ref += tmp;
00038   }
00039   else {
00040     vnl_vector_ref<double> tmp(*n,rw);
00041     self->ls_->transpose_multiply(y_ref, tmp);
00042     x_ref += tmp;
00043   }
00044 
00045   return 0;
00046 }
00047 
00048 int vnl_lsqr::minimize(vnl_vector<double>& result)
00049 {
00050   long m = ls_->get_number_of_residuals();
00051   long n = ls_->get_number_of_unknowns();
00052   double damp = 0;
00053   long leniw = 1;
00054   long* iw = 0;
00055   long lenrw = m;
00056 #if defined __GNUC__ && !defined __STRICT_ANSI__
00057   double rw[m];
00058   double v[n];
00059   double w[n];
00060   double se[n];
00061 #else
00062   vcl_vector<double> rw(m);
00063   vcl_vector<double> v(n);
00064   vcl_vector<double> w(n);
00065   vcl_vector<double> se(n);
00066 #endif
00067   double atol = 0;
00068   double btol = 0;
00069   double conlim = 0;
00070   long nout = -1;
00071   double anorm, acond, rnorm, arnorm, xnorm;
00072 
00073   vnl_vector<double> rhs(m);
00074   ls_->get_rhs(rhs);
00075 
00076   v3p_netlib_lsqr_(
00077         &m, &n, aprod_, &damp, &leniw, &lenrw, iw, &rw[0],
00078         rhs.data_block(), &v[0], &w[0], result.data_block(), &se[0],
00079         &atol, &btol, &conlim, &max_iter_, &nout, &return_code_,
00080         &num_iter_, &anorm, &acond, &rnorm, &arnorm, &xnorm,
00081         this);
00082 
00083   resid_norm_estimate_ = rnorm;
00084   result_norm_estimate_ = xnorm;
00085   A_condition_estimate_ = acond;
00086 
00087 #if 0
00088   vcl_cerr << "A Fro norm estimate      = " << anorm << vcl_endl
00089            << "A condition estimate     = " << acond << vcl_endl
00090            << "Residual norm estimate   = " << rnorm << vcl_endl
00091            << "A'(Ax - b) norm estimate = " << arnorm << vcl_endl
00092            << "x norm estimate          = " << xnorm << vcl_endl;
00093 #endif
00094 
00095   // We should return the return code, as translate_return_code is public and
00096   // it is very misleading that the return code from this function can't be fed
00097   // into translate_return_code. (Brian Amberg)
00098   return return_code_;
00099 }
00100 
00101 void vnl_lsqr::diagnose_outcome(vcl_ostream& os) const
00102 {
00103   translate_return_code(os, return_code_);
00104   os << __FILE__ " : residual norm estimate = " << resid_norm_estimate_ << vcl_endl
00105      << __FILE__ " : result norm estimate   = " << result_norm_estimate_ << vcl_endl
00106      << __FILE__ " : condition no. estimate = " << A_condition_estimate_ << vcl_endl
00107      << __FILE__ " : iterations             = " << num_iter_ << vcl_endl;
00108 }
00109 
00110 void vnl_lsqr::translate_return_code(vcl_ostream& os, int rc)
00111 {
00112   const char* vnl_lsqr_reasons[] = {
00113    "x = 0  is the exact solution. No iterations were performed.",
00114    "The equations A*x = b are probably compatible.  "
00115        "Norm(A*x - b) is sufficiently small, given the "
00116        "values of ATOL and BTOL.",
00117    "The system A*x = b is probably not compatible.  "
00118        "A least-squares solution has been obtained that is "
00119        "sufficiently accurate, given the value of ATOL.",
00120    "An estimate of cond(Abar) has exceeded CONLIM.  "
00121        "The system A*x = b appears to be ill-conditioned.  "
00122        "Otherwise, there could be an error in subroutine APROD.",
00123    "The equations A*x = b are probably compatible.  "
00124        "Norm(A*x - b) is as small as seems reasonable on this machine.",
00125    "The system A*x = b is probably not compatible.  A least-squares "
00126        "solution has been obtained that is as accurate as seems "
00127        "reasonable on this machine.",
00128    "Cond(Abar) seems to be so large that there is no point in doing further "
00129        "iterations, given the precision of this machine. "
00130        "There could be an error in subroutine APROD.",
00131    "The iteration limit ITNLIM was reached."
00132   };
00133 
00134   if
00135     (rc < 0 || rc > 7) os << __FILE__ " : Illegal return code : " << rc << vcl_endl;
00136   else
00137     os << __FILE__ " : " << vnl_lsqr_reasons[rc] << vcl_endl;
00138 }