第28章 矩阵运算

这一章给出了求解线性方程组Ax=b(|A|！=0)和求矩阵的逆的伪代码，下面用c++实现。

线性方程组：

1： 先求出P,L,U，使得PA=LU成立

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];                }        }}

2：然后再用获得的L, U, P去求解AX=b;

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]<<" ";}