Go to the documentation of this file.00001 #ifndef _LINALG_TOOLS_HPP_
00002 #define _LINALG_TOOLS_HPP_
00003
00004 #include "types.hpp"
00005
00006
00007
00008 #if BOOST_VERSION < 103301 // missing includes in lu.hpp
00009 #include "boost/numeric/ublas/vector.hpp"
00010 #include "boost/numeric/ublas/vector_proxy.hpp"
00011 #include "boost/numeric/ublas/matrix.hpp"
00012 #include "boost/numeric/ublas/matrix_proxy.hpp"
00013 #include "boost/numeric/ublas/triangular.hpp"
00014 #include "boost/numeric/ublas/operation.hpp"
00015 #endif
00016
00017 #include <boost/numeric/ublas/lu.hpp>
00018
00019
00020 typedef boost::numeric::ublas::matrix<double> mat;
00021 typedef boost::numeric::ublas::identity_matrix<double> identity_mat;
00022
00023
00024
00028 template<class M1, class M2>
00029 void lu_inv(M1 const& m, M2& inv) {
00030
00031 assert( m.size1() == m.size2() );
00032
00033 assert( inv.size1() == m.size1() && inv.size1() == m.size2() );
00034
00035
00036 using namespace boost::numeric::ublas;
00037
00038 mat mLu(m);
00039
00040
00041 lu_factorize(mLu);
00042
00043
00044 inv.assign(identity_mat(m.size1()));
00045
00046
00047 lu_substitute<mat const, M2>(mLu, inv);
00048 };
00049
00050
00051
00055 template<class M>
00056 double lu_det(M const& m) {
00057
00058 using namespace boost::numeric::ublas;
00059
00060 assert( m.size1() == m.size2() );
00061
00062
00063
00064 mat mLu(m);
00065 permutation_matrix<std::size_t> pivots(m.size1());
00066
00067 lu_factorize(mLu, pivots);
00068
00069 double det = 1.0;
00070
00071 for (std::size_t i=0; i < pivots.size(); ++i) {
00072 if (pivots(i) != i)
00073 det *= -1.0;
00074 det *= mLu(i,i);
00075 }
00076 return det;
00077 };
00078
00079
00080
00081 #endif //_LINALG_TOOLS_HPP_