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>
35 using namespace boost::numeric::ublas;
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,
73 using namespace boost::numeric::ublas;
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_