[Boost-users] matrix inversion in boost.

Subject: [Boost-users] matrix inversion in boost.
From: Mohit Singh (mohit1007_at_[hidden])
Date: 2013-06-20 20:25:20

Hi,
 I am trying to do a simple matrix inversion operation using boost. But I
am getting an error.

Basically what I am trying to find is inversted_matrix  =
inverse(trans(matrix) * matrix)

But I am getting an error
Check failed in file boost_1_53_0/boost/numeric/ublas/lu.hpp at line 299:
detail::expression_type_check (prod (triangular_adaptor<const_matrix_type,
upper> (m), e), cm2)
terminate called after throwing an instance of
'boost::numeric::ublas::internal_logic'
  what():  internal logic
Aborted (core dumped)

My attempt:

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/io.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>
namespace ublas = boost::numeric::ublas;
template<class T>
bool InvertMatrix (const ublas::matrix<T>& input, ublas::matrix<T>&
inverse) {
  using namespace boost::numeric::ublas;
  typedef permutation_matrix<std::size_t> pmatrix;
  // create a working copy of the input
  matrix<T> A(input);
  // create a permutation matrix for the LU-factorization
  pmatrix pm(A.size1());
  // perform LU-factorization
  int res = lu_factorize(A,pm);
        if( res != 0 ) return false;
  // create identity matrix of "inverse"
  inverse.assign(ublas::identity_matrix<T>(A.size1()));
  // backsubstitute to get the inverse
  lu_substitute(A, pm, inverse);
  return true;
 }

int main(){
using namespace boost::numeric::ublas;
matrix<double> m(4,5);
 vector<double> v(4);
vector<double> thetas;

m(0,0) = 1; m(0,1) = 2104; m(0,2) = 5; m(0,3) = 1;m(0,4) = 45;
 m(1,0) = 1; m(1,1) = 1416; m(1,2) = 3; m(1,3) = 2;m(1,4) = 40;
m(2,0) = 1; m(2,1) = 1534; m(2,2) = 3; m(2,3) = 2;m(2,4) = 30;
 m(3,0) = 1; m(3,1) = 852;  m(3,2) = 2; m(3,3) = 1;m(3,4) = 36;
std::cout<<m<<std::endl;
 matrix<double> product = prod(trans(m), m);
std::cout<<product<<std::endl;
 matrix<double> inversion(5,5);
bool inverted;
inverted = InvertMatrix(product, inversion);
 std::cout<<inversion<<std::endl;

}

-- 
Mohit
"When you want success as badly as you want the air, then you will get it.
There is no other secret of success."
-Socrates