12 #ifndef _MATRIXINVERSE_HPP_ 13 #define _MATRIXINVERSE_HPP_ 15 #include <boost/numeric/ublas/matrix.hpp> 16 #include <boost/numeric/ublas/matrix_proxy.hpp> 17 #include <boost/numeric/ublas/io.hpp> 18 #include <boost/numeric/ublas/matrix_expression.hpp> 24 #include <boost/numeric/ublas/vector.hpp> 25 #include <boost/numeric/ublas/vector_proxy.hpp> 26 #include <boost/numeric/ublas/triangular.hpp> 27 #include <boost/numeric/ublas/lu.hpp> 37 typedef permutation_matrix<std::size_t> pmatrix;
42 pmatrix pm(A.size1());
45 int res = lu_factorize(A,pm);
46 if( res != 0 )
return false;
49 inverse.assign(identity_matrix<typename M::value_type>(A.size1()));
52 lu_substitute(A, pm, inverse);
69 boost::numeric::ublas::matrix<T>
70 gjinverse(
const boost::numeric::ublas::matrix<T> &m,
75 const size_t size = m.size1();
80 if (size != m.size2() || size == 0)
105 matrix<T>
A(size, 2*size);
106 matrix_range<matrix<T> > Aleft(A,
110 matrix_range<matrix<T> > Aright(A,
112 range(size, 2*size));
113 Aright = identity_matrix<T>(size);
116 for (
size_t k = 0;
k < size;
k++)
119 for (
size_t kk = 0; kk < size; kk++)
125 for (
size_t i = kk+1; i < size; i++)
137 std::cerr <<
"Error:" << __FUNCTION__ <<
":" 138 <<
"Input matrix is singular, because cannot find" 139 <<
" a row to swap while eliminating zero-diagonal.";
145 matrix_row<matrix<T> > rowk(A, kk);
146 matrix_row<matrix<T> > rowl(A, l);
160 for (
size_t j =
k+1; j < 2*size; j++)
165 for (
size_t i = 0; i < size; i++)
171 for (
size_t j =
k+1; j < 2*size; j++)
172 A(i,j) -=
A(
k,j) *
A(i,
k);
188 #endif // _MATRIXINVERSE_HPP_ boost::numeric::ublas::matrix< T > gjinverse(const boost::numeric::ublas::matrix< T > &m, bool &singular)
Invert a matrix via gauss-jordan algorithm (PARTIAL PIVOT)
bool InvertMatrix(const M &input, M &inverse)