下面的代码在VS2005和VS2008里面可以编译运行,boost用的是1.35 // REMEMBER to update "lu.hpp" header includes from boost-CVS #include <boost/numeric/ublas/vector.hpp> #include <boost/numeric/ublas/vector_proxy.hpp> #include <boost/numeric/ublas/matrix.hpp> #include <boost/numeric/ublas/triangular.hpp> #include <boost/numeric/ublas/lu.hpp> #include <boost/numeric/ublas/io.hpp> #include <boost/numeric/ublas/vector.hpp> #include <iostream> namespace ublas = boost::numeric::ublas; 原来的地址:www.crystalclearsoftware.com/cgi-bin/boost_wiki/wiki.pl //! \brief 求矩阵的逆 //! \param input[in] 被求矩阵 //! \param inverse[out] 逆矩阵 /* Matrix inversion routine. Uses lu_factorize and lu_substitute in uBLAS to invert a matrix */ template<class T> bool InvertMatrix (const ublas::matrix<T>& input, ublas::matrix<T>& inverse) { using namespace boost::numeric::ublas; typedef permutation_matrix<std::size_t> pmatrix; // create a working copy of the input matrix<T> 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<T>(A.size1())); // backsubstitute to get the inverse lu_substitute(A, pm, inverse); return true; } int main() { namespace ublas = boost::numeric::ublas; //! \brief 三次Hermite矩阵 ublas::matrix<double> MHermite(4,4), MHInverse(4,4); MHermite(0,0) = 2; MHermite(0,1) = -2; MHermite(0,2) = 1; MHermite(0,3) = 1; MHermite(1,0) = -3; MHermite(1,1) = 3; MHermite(1,2) = -2; MHermite(1,3) = -1; MHermite(2,0) = 0; MHermite(2,1) = 0; MHermite(2,2) = 1; MHermite(2,3) = 0; MHermite(3,0) = 1; MHermite(3,1) = 0; MHermite(3,2) = 0; MHermite(3,3) = 0; std::cout <<"MHermit矩阵"<< MHermite <<std::endl; ublas::matrix<double> mTmp(4,4); InvertMatrix(MHermite, MHInverse); InvertMatrix(MHermite, mTmp); std::cout <<"MHermit逆矩阵"<< mTmp << std::endl; InvertMatrix(mTmp, mTmp); std::cout <<"MHermite逆矩阵的逆矩阵 " <<mTmp <<std::endl; std::cout <<"" <<std::endl; //! \brief 三次Bezier矩阵 ublas::matrix<double> MBezier(4,4),MBInverse(4,4); MBezier(0,0) = -1; MBezier(0,1) = 3; MBezier(0,2) = -3; MBezier(0,3) = 1; MBezier(1,0) = 3; MBezier(1,1) = -6; MBezier(1,2) = 3; MBezier(1,3) = 0; MBezier(2,0)= -3; MBezier(2,1) = 3; MBezier(2,2) = 0; MBezier(2,3) = 0; MBezier(3,0) = 1; MBezier(3,1) = 0; MBezier(3,2) = 0; MBezier(3,3) = 0; std::cout <<"MBezier矩阵"<< MBezier<<std::endl; InvertMatrix(MBezier, MBInverse); InvertMatrix(MBezier, mTmp); std::cout <<"MBezier逆矩阵 " <<mTmp <<std::endl; InvertMatrix(mTmp, mTmp); std::cout <<"MBezier逆矩阵的逆矩阵 " <<mTmp <<std::endl; std::cout <<"" <<std::endl; //! \brief 三次B样条矩阵 ublas::matrix<double> MBSpline(4,4),MBSInverse(4,4); MBSpline(0,0) = -1; MBSpline(0,1) = 3; MBSpline(0,2) = -3; MBSpline(0,3) = 1; MBSpline(1,0) = 3; MBSpline(1,1) = -6; MBSpline(1,2) = 3; MBSpline(1,3) = 0; MBSpline(2,0) = -3; MBSpline(2,1) = 0; MBSpline(2,2) = 3; MBSpline(2,3) = 0; MBSpline(3,0) = 1; MBSpline(3,1) = 4; MBSpline(3,2) = 1; MBSpline(3,3) = 0; std::cout <<"MBSpline矩阵"<< MBSpline<<std::endl; InvertMatrix(MBSpline, MBSInverse); InvertMatrix(MBSpline, mTmp); std::cout <<"MBSpline逆矩阵 " <<MBSInverse <<std::endl; ublas::matrix<double> mTmp1(4,4); InvertMatrix(mTmp, mTmp); std::cout <<"MBSpline逆矩阵的逆矩阵 " <<mTmp <<std::endl; std::cout <<"" <<std::endl; char a; std::cin >> a; return 0; }