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