boost库

ThefollowingcodeinvertsthematrixinputusingLU-decompositionwith

nce:NumericalRecipiesinC,2nd

ed.,byPress,Teukolsky,Vetterling&Flannery.

#ifndefINVERT_MATRIX_HPP

#defineINVERT_MATRIX_HPP

//REMEMBERtoupdate""headerincludesfromboost-CVS

#include

#include

#include

#include

#include

#include

namespaceublas=boost::numeric::ublas;

/*Matrixinversionroutine.

Useslu_factorizeandlu_substituteinuBLAStoinvertamatrix*/

template

boolInvertMatrix?(constublas::matrix&input,ublas::matrix&

inverse){

usingnamespaceboost::numeric::ublas;

typedefpermutation_matrixpmatrix;

//createaworkingcopyoftheinput

matrixA(input);

//createapermutationmatrixfortheLU-factorization

pmatrixpm(1());

//performLU-factorization

intres=lu_factorize(A,pm);

if(res!=0)returnfalse;

//createidentitymatrixof"inverse"

(ublas::identity_matrix(1()));

//backsubstitutetogettheinverse

lu_substitute(A,pm,inverse);

returntrue;

}

#endif//INVERT_MATRIX_HPP

s,FredrikOrderud.

MatrixinversionwithSingularityDetection

omplexandprobably

inefficientversionisatthesecondsectionofEffectiveUBLAS/Matrix

Inversion.

#include

#include

#include

#include

#include

int_tmain(intargc,_TCHAR*argv[])

{

//可不要与std::vector混在一起啊.

//可考虑usingnamespaceboost::numeric;

//然后用ublas::vector

usingnamespaceboost::numeric::ublas;

constintn=4;

permutation_matrixP(n);

matrixm(n,n);

vectorx(n);//用来放解向量

vectorv(n);//常量向量

v(0)=5;

v(1)=3;

v(2)=2;

v(3)=2;

m(0,0)=1;

m(0,1)=2;

m(0,2)=1;

m(0,3)=1;

m(1,0)=2;

m(1,1)=1;

m(1,2)=0;

m(1,3)=0;

m(2,0)=1;

m(2,1)=0;

m(2,2)=1;

m(2,3)=0;

m(3,0)=1;

m(3,1)=0;

m(3,2)=0;

m(3,3)=1;

try{

lu_factorize(m,P);

x=v;

lu_substitute(m,P,x);

std::cout<<"解向量:"<

}

catch(...){

std::cout<<"奇异矩阵(无解或无穷解)."<

}

getchar();

return0;

}

m=[1,2,1,1;2,1,0,0;1,0,1,0;1,0,0,1];

lu_factorize(m)的结果如下：

m=[1,2,1,1;2,-3,-2,-2;1,2/3,4/3,1/3;1,2/3,1/4,5/4];

其中L,U阵为

L=[1,0,0,0;2,1,0,0;1,2/3,1,0;1,2/3,1/4,1];即L=tril(m,-1)+eye(1);m的下三角阵与单位矩阵之和

U=[1,2,1,1;0,-3,-2,-2;0,0,4/3,1/3;0,0,0,4/5];即U=tril(m);m的上三角阵

可得m=L*U;

lu_factorize(m,P)的结果

P=[1,1,2,3];注意P是一permutation_matrix

m=[2,1,0,0;1/2,3/2,1,1;1/2,-1/3,4/3,1/3;1/2,-1/3,1/4,5/4];

L=[1,0,0,0;1/2,1,0,0;1/2,-1/3,1,0;1/2,-1/3,1/4,1];

U=[2,1,0,0;0,3/2,1,1;0,0,4/3,1/3;0,0,0,5/4];

Pm=LU，P=[1,1,2,3]

Pm=A，A的0行为m的1行，A的1行为新的1行(m的第0行)，A的2行为m的2行，

A的3的行为m的3行。