第28章 矩阵运算

线性方程组：

void LUPDecomposition(vector< vector<double> >& array,vector<int>& P)
{
int dimension=array.size();

P.resize(dimension);
for(int i=0;i!=P.size();++i)
P[i]=i;

for(int k=0;k!=array.size();++k)
{
double maxAbsoluteValue=0;
int rowLabel=k;
for(int i=k;i!=array.size();++i)
if(maxAbsoluteValue<abs(array[i][k])){
maxAbsoluteValue=abs(array[i][k]);
rowLabel=i;
}

if(maxAbsoluteValue==0)
throw runtime_error("singular matrix");

swap(P[k],P[rowLabel]);

for(int i=0;i!=array.size();++i)
swap(array[k][i],array[rowLabel][i]);

for(int i=k+1;i!=array.size();++i)
{
array[i][k]=array[i][k]/array[k][k];
for(int j=k+1;j!=array.size();++j)
array[i][j]=array[i][j]-array[i][k]*array[k][j];
}
}
}

void LUPSolution(const vector< vector<double> >& array,const vector<int>& P,const vector<double>& b,vector<double>& solution)
{
int dimension=array.size();

vector<double> tmp(dimension);

tmp[0]=b[P[0]];

for(int i=1;i!=dimension;++i)
{
double tmpSum=0;
for(int j=0;j<=i-1;++j)
tmpSum+=array[i][j]*tmp[j];

tmp[i]=b[P[i]]-tmpSum;
}

solution.resize(dimension);
solution[dimension-1]=tmp[dimension-1]/array[dimension-1][dimension-1];

for(int i=dimension-2;i>=0;--i)
{
double tmpSum=0;
for(int j=i+1;j<=dimension-1;++j)
tmpSum+=array[i][j]*solution[j];

solution[i]=(tmp[i]-tmpSum)/array[i][i];
}

}

void linearSystem(vector< vector<double> >& array,const vector<double>& b)
{
vector<int> P;
LUPDecomposition(array,P);

vector<double> solution;

LUPSolution(array,P,b,solution);

for(int i=0;i!=solution.size();++i)
cout<<solution[i]<<" ";
cout<<endl;
}

求解矩阵的逆的代码：

void inverseMatrix(vector< vector<double> >& array)
{
vector<int> P;
LUPDecomposition(array,P);

int dimension=array.size();
vector< vector<double> > inverseMat;
inverseMat.resize(dimension);
for(int i=0;i!=inverseMat.size();++i)
inverseMat[i].resize(dimension);

vector<double> b;
b.resize(dimension,0);
b[0]=1;

for(int i=0;i!=dimension;++i)
{
vector<double> solution;
LUPSolution(array,P,b,solution);

for(int j=0;j!=dimension;++j)
inverseMat[j][i]=solution[j];

if(i<dimension-1)
swap(b[i],b[i+1]);
}

for(int i=0;i!=inverseMat.size();++i)
for(int j=0;j!=inverseMat[i].size();++j)
cout<<inverseMat[i][j]<<" ";
}