ublas-invert.h

Go to the documentation of this file.

#ifndef _RHEOLEF_UBLAS_INVERT_H
#define _RHEOLEF_UBLAS_INVERT_H
// http://www.crystalclearsoftware.com/cgi-bin/boost_wiki/wiki.pl?action=browse&diff=1&id=LU_Matrix_Inversion
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/io.hpp>
namespace ublas = boost::numeric::ublas;

namespace rheolef { 
template<class T>
bool invert (const ublas::matrix<T>& input, ublas::matrix<T>& inverse) {
  using namespace boost::numeric::ublas;
  typedef permutation_matrix<std::size_t> pmatrix;
  matrix<T> A(input);
  pmatrix pm(A.size1());
 
  int res = lu_factorize(A,pm);
  if( res != 0 ) return false;

  inverse.assign(ublas::identity_matrix<T>(A.size1()));

  lu_substitute(A, pm, inverse);

  return true;
}
template<class T>
boost::numeric::ublas::matrix<T>
invert(const boost::numeric::ublas::matrix<T> &m, bool &is_singular)
{ 
   ublas::matrix<T> inv_m (m.size1(), m.size2());
   is_singular = ! invert (m, inv_m);
   return inv_m;
}
// http://archives.free.net.ph/message/20080909.064313.59c122c4.fr.html 
template<class matrix_T>
double determinant(ublas::matrix_expression<matrix_T> const& mat_r) {
  double det = 1.0;
  matrix_T mLu(mat_r());
  ublas::permutation_matrix<std::size_t> pivots(mat_r().size1() );
  bool is_singular = lu_factorize(mLu, pivots);
  if (!is_singular) {
    for (std::size_t i = 0; i < pivots.size(); ++i) {
      if (pivots(i) != i) {
        det *= -1.0;
      }
      det *= mLu(i,i);
    }
  } else {
    det = 0.0;
  }
  return det;
} 
}// namespace rheolef
#endif // _RHEOLEF_UBLAS_INVERT_H