ublas求矩阵的逆

下面的代码在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;

}