linAlgTools.hpp
Go to the documentation of this file.
00001 #ifndef _LINALG_TOOLS_HPP_
00002 #define _LINALG_TOOLS_HPP_
00003 
00004 #include "types.hpp"
00005 
00006 // courtesy of Thomas Lemaire
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 // forget about JFR_PRECOND macros...
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   // ublasExtra::lu_inv(): input matrix must be squared
00033   assert( inv.size1() == m.size1() && inv.size1() == m.size2() );
00034   // ublasExtra::lu_inv(): invalid size for inverse matrix
00035 
00036   using namespace boost::numeric::ublas;
00037   // create a working copy of the input
00038   mat mLu(m);
00039 
00040   // perform LU-factorization
00041   lu_factorize(mLu);
00042 
00043   // create identity matrix of "inverse"
00044   inv.assign(identity_mat(m.size1()));
00045 
00046   // backsubstitute to get the inverse
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   // ublasExtra::lu_det: matrix must be square
00062 
00063   // create a working copy of the input
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_


gaussian_process
Author(s): Maintained by Juergen Sturm
autogenerated on Mon Oct 6 2014 00:09:34