title: matrix inversion in boost
body: <p>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 </p>
  
  <pre><code>Check failed in file boost_1_53_0/boost/numeric/ublas/lu.hpp at line 299: 
  detail::expression_type_check (prod (triangular_adaptor&lt;const_matrix_type, 
  upper&gt; (m), e), cm2) 
  terminate called after throwing an instance of 
  'boost::numeric::ublas::internal_logic' 
    what(): internal logic 
  Aborted (core dumped) 
  </code></pre>
  
  <p>My attempt: </p>
  
  <pre><code>#include &lt;boost/numeric/ublas/matrix.hpp&gt; 
  #include &lt;boost/numeric/ublas/vector.hpp&gt; 
  #include &lt;boost/numeric/ublas/io.hpp&gt; 
  #include &lt;boost/numeric/ublas/vector_proxy.hpp&gt; 
  #include &lt;boost/numeric/ublas/matrix.hpp&gt; 
  #include &lt;boost/numeric/ublas/triangular.hpp&gt; 
  #include &lt;boost/numeric/ublas/lu.hpp&gt; 
  
  namespace ublas = boost::numeric::ublas; 
  template&lt;class T&gt; 
  bool InvertMatrix (const ublas::matrix&lt;T&gt;&amp; input, ublas::matrix&lt;T&gt;&amp; inverse) { 
      using namespace boost::numeric::ublas; 
      typedef permutation_matrix&lt;std::size_t&gt; pmatrix; 
      // create a working copy of the input 
      matrix&lt;T&gt; 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&lt;T&gt;(A.size1())); 
      // backsubstitute to get the inverse 
      lu_substitute(A, pm, inverse); 
      return true; 
  }
  
  int main(){ 
      using namespace boost::numeric::ublas; 
      matrix&lt;double&gt; m(4,5); 
      vector&lt;double&gt; v(4); 
      vector&lt;double&gt; 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&lt;&lt;m&lt;&lt;std::endl; 
      matrix&lt;double&gt; product = prod(trans(m), m); 
      std::cout&lt;&lt;product&lt;&lt;std::endl; 
      matrix&lt;double&gt; inversion(5,5); 
      bool inverted; 
      inverted = InvertMatrix(product, inversion); 
      std::cout &lt;&lt; inversion &lt;&lt; std::endl; 
  } 
  </code></pre>