LCOV - code coverage report
Current view: directory - src/vxl/core/vnl/algo - vnl_sparse_lu.cxx (source / functions) Found Hit Coverage
Test: vnl.info Lines: 141 1 0.7 %
Date: 2011-09-02 Functions: 15 2 13.3 %

       1                 : // This is core/vnl/algo/vnl_sparse_lu.cxx
       2                 : #ifdef VCL_NEEDS_PRAGMA_INTERFACE
       3                 : #pragma implementation
       4                 : #endif
       5                 : //:
       6                 : // \file
       7                 : #include "vnl_sparse_lu.h"
       8                 : #include <vcl_cassert.h>
       9                 : #include <vcl_iostream.h>
      10                 : 
      11                 : #include <sparse/spMatrix.h>
      12                 : 
      13                 : // destructor - undo the spCreate() from the constructor(s)
      14                 : // (memory leak fix of 7 Feb. 2008 by Toon Huysmans)
      15               0 : vnl_sparse_lu::~vnl_sparse_lu()
      16                 : {
      17               0 :   spDestroy( pmatrix_ );
      18               0 : }
      19                 : 
      20                 : //: constructor - controls if condition information is computed
      21               0 : vnl_sparse_lu::vnl_sparse_lu(vnl_sparse_matrix<double> const & M, operation mode):
      22               0 :   A_(M), factored_(false),condition_computed_(false), mode_(mode),norm_(0), rcond_(0), largest_(0), pivot_thresh_(0),absolute_thresh_(0),diag_pivoting_(1),pmatrix_(0)
      23                 : {
      24               0 :   int n = (int)M.columns();
      25               0 :   assert(n == (int)(M.rows()));
      26               0 :   int error = 0;
      27               0 :   pmatrix_ = spCreate(n, 0, &error);
      28               0 :   if (error!=spOKAY)
      29                 :   {
      30               0 :     vcl_cout << "In vnl_sparse_lu::vnl_sparse_lu - error in creating matrix\n";
      31               0 :     return;
      32                 :   }
      33                 :   // fill the internal sparse matrix from A_
      34               0 :   spElement* pelement = 0;
      35               0 :   for (A_.reset(); A_.next();)
      36                 :   {
      37               0 :     int r = A_.getrow();
      38               0 :     int c = A_.getcolumn();
      39               0 :     double v = A_.value();
      40               0 :     pelement = spGetElement(pmatrix_, r+1, c+1);
      41               0 :     if (pelement == 0)
      42                 :     {
      43               0 :       vcl_cout<< "In vnl_sparse_lu::vnl_sparse_lu - error in getting element\n";
      44               0 :       return;
      45                 :     }
      46               0 :     *pelement = v;
      47                 :   }
      48               0 :   if (mode==estimate_condition || mode_==estimate_condition_verbose)
      49                 :   {
      50               0 :     largest_ = spLargestElement(pmatrix_);
      51               0 :     if (mode_==estimate_condition_verbose)
      52               0 :       vcl_cout << " Largest element in matrix = " << largest_ << '\n';
      53               0 :     norm_ = spNorm(pmatrix_);
      54                 :   }
      55               0 : }
      56                 : 
      57                 : //: estimate the condition number.
      58               0 : bool vnl_sparse_lu::est_condition()
      59                 : {
      60               0 :   int error = spOKAY;
      61               0 :   rcond_ = spCondition(pmatrix_, norm_, &error);
      62               0 :   if (error!=spOKAY)
      63                 :   {
      64               0 :     vcl_cout << "In vnl_sparse_lu::est_condition(..) - error in condition number calculation\n";
      65               0 :     return false;
      66                 :   }
      67               0 :   condition_computed_ = true;
      68               0 :   return true;
      69                 : }
      70                 : 
      71                 : //: Solve least squares problem M x = b.
      72               0 : void vnl_sparse_lu::solve(vnl_vector<double> const& b, vnl_vector<double>* x)
      73                 : {
      74               0 :   if (!pmatrix_)
      75                 :   {
      76               0 :     vcl_cout << "In vnl_sparse_lu::solve(..) - matrix not defined\n";
      77               0 :     return;
      78                 :   }
      79               0 :   unsigned n = b.size();
      80               0 :   assert(n == A_.columns());
      81               0 :   spREAL* rhs = new spREAL[n+1];
      82               0 :   for (unsigned i = 0; i<n; ++i)
      83               0 :     rhs[i+1]=b[i];
      84               0 :   if (mode_==verbose || mode_==estimate_condition_verbose)
      85                 :   {
      86               0 :     vcl_cout << "Matrix before ordering\n";
      87               0 :     spPrint(pmatrix_,1,1,1);
      88                 :   }
      89                 : 
      90               0 :   int error = 0;
      91                 :   //check if the matrix needs ordering
      92                 :   //if not, run the decomposition
      93               0 :   if (!factored_)
      94                 :   {
      95                 :     error = spOrderAndFactor(pmatrix_, rhs, pivot_thresh_,
      96               0 :                              absolute_thresh_, diag_pivoting_);
      97               0 :     if (error != spOKAY)
      98                 :     {
      99               0 :       vcl_cout << "In vnl_sparse_lu::solve(..) - error in factoring\n";
     100               0 :       return;
     101                 :     }
     102               0 :     if (mode_==estimate_condition || mode_==estimate_condition_verbose)
     103               0 :       if (!est_condition())
     104               0 :         return;
     105               0 :     factored_ = true;
     106                 :   }
     107                 : 
     108               0 :   if (mode_==verbose || mode_==estimate_condition_verbose)
     109                 :   {
     110               0 :     vcl_cout << "Matrix after ordering\n";
     111               0 :     spPrint(pmatrix_,1,1,1);
     112                 :   }
     113                 : 
     114               0 :   spSolve(pmatrix_, rhs, rhs);
     115                 : 
     116               0 :   for (unsigned i = 0; i<n; ++i)
     117               0 :     (*x)[i]=rhs[i+1];
     118                 : 
     119               0 :   delete [] rhs;
     120                 : }
     121                 : 
     122                 : //: Solve least squares problem M x = b.
     123               0 : vnl_vector<double> vnl_sparse_lu::solve(vnl_vector<double> const& b)
     124                 : {
     125               0 :   vnl_vector<double> ret(b.size());
     126               0 :   this->solve(b, &ret);
     127               0 :   return ret;
     128                 : }
     129                 : 
     130                 : //: Solve problem M^t x = b
     131               0 : void vnl_sparse_lu::solve_transpose(vnl_vector<double> const& b, vnl_vector<double>* x)
     132                 : {
     133               0 :   if (!pmatrix_)
     134                 :   {
     135               0 :     vcl_cout << "In vnl_sparse_lu::solve(..) - matrix not defined\n";
     136               0 :     return;
     137                 :   }
     138               0 :   unsigned n = b.size();
     139               0 :   assert(n == A_.columns());
     140               0 :   spREAL* rhs = new spREAL[n+1];
     141               0 :   for (unsigned i = 0; i<n; ++i)
     142               0 :     rhs[i+1]=b[i];
     143               0 :   int error = 0;
     144               0 :   if (mode_== verbose || mode_== estimate_condition_verbose)
     145                 :   {
     146               0 :     vcl_cout << "Matrix before ordering\n";
     147               0 :     spPrint(pmatrix_,1,1,1);
     148                 :   }
     149                 : 
     150                 :   //check if the matrix needs ordering
     151                 :   //if not, run the decomposition
     152               0 :   if (!factored_)
     153                 :   {
     154                 :     error = spOrderAndFactor(pmatrix_, rhs, pivot_thresh_,
     155               0 :                              absolute_thresh_, diag_pivoting_);
     156               0 :     if (error != spOKAY)
     157                 :     {
     158               0 :       vcl_cout << "In vnl_sparse_lu::solve(..) - error in factoring\n";
     159               0 :       return;
     160                 :     }
     161               0 :     if (mode_==estimate_condition || mode_==estimate_condition_verbose)
     162               0 :       if (!est_condition())
     163               0 :         return;
     164               0 :     factored_ = true;
     165                 :   }
     166                 : 
     167               0 :   if (mode_==verbose || mode_== estimate_condition_verbose)
     168                 :   {
     169               0 :     vcl_cout << "Matrix after ordering\n";
     170               0 :     spPrint(pmatrix_,1,1,1);
     171                 :   }
     172                 : 
     173               0 :   spSolveTransposed(pmatrix_, rhs, rhs);
     174                 : 
     175               0 :   for (unsigned i = 0; i<n; ++i)
     176               0 :     (*x)[i]=rhs[i+1];
     177                 : 
     178               0 :   delete [] rhs;
     179                 : }
     180                 : 
     181                 : //: Solve problem M^t x = b
     182               0 : vnl_vector<double> vnl_sparse_lu::solve_transpose(vnl_vector<double> const& b)
     183                 : {
     184               0 :   vnl_vector<double> ret(b.size());
     185               0 :   this->solve_transpose(b, &ret);
     186               0 :   return ret;
     187                 : }
     188                 : 
     189                 : // This routine might be eventually useful if other operations are to be done
     190                 : // after the factoring. Note that the routine spRowColOrder was
     191                 : // added to the sparse library (in spoutput.c).
     192                 : #if 0
     193                 : //: copy the matrix into a vnl_sparse_matrix
     194                 : vnl_sparse_matrix<double> vnl_sparse_lu::lu_matrix()
     195                 : {
     196                 :   unsigned n = A_.rows();
     197                 :   vnl_sparse_matrix<double> temp(n, n);
     198                 :   int error = 0;
     199                 :   spElement* pelement = 0;
     200                 :   if (!factored_)
     201                 :   {
     202                 :     error = spFactor(pmatrix_);
     203                 :     if (error != spOKAY)
     204                 :     {
     205                 :       vcl_cout << "In vnl_sparse_lu::lu_matrix() - factoring failed\n";
     206                 :       return temp;
     207                 :     }
     208                 :     factored_=true;
     209                 :   }
     210                 :   // get the row and column maps
     211                 :   int* row_map = new int[n+1];
     212                 :   int* col_map = new int[n+1];
     213                 :   spRowColOrder(pmatrix_, row_map, col_map);
     214                 : 
     215                 :   // create inverse maps
     216                 :   int* inv_row_map = new int[n+1];
     217                 :   int* inv_col_map = new int[n+1];
     218                 :   for (unsigned i = 1; i<=n; ++i)
     219                 :   {
     220                 :     inv_row_map[row_map[i]]=i;
     221                 :     inv_col_map[col_map[i]]=i;
     222                 :   }
     223                 : 
     224                 :   if (mode_==verbose || mode_==estimate_condition_verbose)
     225                 :     for (unsigned i = 1; i<=n; ++i)
     226                 :       vcl_cout << "inv_row_map[" << i << "]= " << inv_row_map[i]
     227                 :                << "    inv_col_map[" << i << "]= " << inv_col_map[i] << '\n';
     228                 :   for (unsigned r = 0; r<n; r++)
     229                 :     for (unsigned c = 0; c<n; c++)
     230                 :     {
     231                 :       pelement = spGetElement(pmatrix_, inv_row_map[r+1], inv_col_map[c+1]);
     232                 :       if (pelement)
     233                 :       {
     234                 :         double v = *pelement;
     235                 :         temp(r,c)= v;
     236                 :       }
     237                 :     }
     238                 :   delete [] row_map;
     239                 :   delete [] col_map;
     240                 :   return temp;
     241                 : }
     242                 : #endif
     243                 : 
     244                 : //: Compute determinant.
     245               0 : double vnl_sparse_lu::determinant()
     246                 : {
     247                 :   int exponent;
     248                 :   double determ;
     249               0 :   if (!factored_)
     250                 :   {
     251               0 :     spFactor(pmatrix_);
     252               0 :     if (mode_==estimate_condition || mode_==estimate_condition_verbose)
     253               0 :       if (!est_condition())
     254               0 :         return 0;
     255               0 :     factored_ = true;
     256                 :   }
     257               0 :   spDeterminant(pmatrix_, &exponent, &determ);
     258               0 :   if (exponent<0)
     259                 :   {
     260               0 :     while (exponent<0)
     261                 :     {
     262               0 :       determ *= 0.1;
     263               0 :       exponent++;
     264                 :     }
     265               0 :     return determ;
     266                 :   }
     267               0 :   else if (exponent>0)
     268               0 :     while (exponent>0)
     269                 :     {
     270               0 :       determ *= 10.0;
     271               0 :       exponent--;
     272                 :     }
     273                 : 
     274               0 :   return determ;
     275                 : }
     276                 : 
     277                 : //: the reciprocal of the condition number
     278               0 : double vnl_sparse_lu::rcond()
     279                 : {
     280               0 :   if (!factored_)
     281                 :   {
     282               0 :     spFactor(pmatrix_);
     283               0 :     if (mode_==estimate_condition || mode_==estimate_condition_verbose)
     284               0 :       if (!est_condition())
     285               0 :         return 0;
     286               0 :     factored_ = true;
     287                 :   }
     288               0 :   return rcond_;
     289                 : }
     290                 : 
     291                 : //:Estimated upper bound of error in solution
     292               0 : double vnl_sparse_lu::max_error_bound()
     293                 : {
     294               0 :   if (mode_!=estimate_condition && mode_ != estimate_condition_verbose)
     295               0 :     return 0;
     296                 : 
     297               0 :   if (!factored_)
     298                 :   {
     299               0 :     spFactor(pmatrix_);
     300               0 :     if (!est_condition())
     301               0 :       return 0;
     302               0 :     factored_ = true;
     303                 :   }
     304               0 :   double roundoff = spRoundoff(pmatrix_, largest_);
     305               0 :   if (rcond_>0)
     306               0 :     return roundoff/rcond_;
     307               0 :   return 0;
     308            6312 : }

Generated by: LCOV version 1.7