线性方程组的直接解

#include <bits/stdc++.h>
using namespace std;

double cut_double(double n){       //降低double精度
    double ans=floor((n*pow(10,6)+0.5)) / pow(10,6);
    return ans;
}

vector<double>Gauss(vector<vector<double>>A){
    int n=A.size();
    for(int i=0;i<n;++i){
        if(A[i][i]==0){   
            int j=i;
            while(A[j][i]==0){
                if(j==n) {       //到最后一行都没有非0
                    cout<<"错误"<<endl;
                    break;
                }
                ++j;
            }
            swap(A[i],A[j]);   //A[j]!=0,交换行
        }
        for(int j=0;j<n;++j){      //原本行A[i],消去A[0]-A[n-1]行的第i列元素
            if(j==i) continue;
            double r=A[j][i]/A[i][i];
            r=cut_double(r);
            for(int k=i;k<=n;++k){   //要消去第i列，i列之前已经为0
                A[j][k]-=r*A[i][k];
                A[j][k]=cut_double(A[j][k]);
            }
        }
    }
    vector<double> ans(n);
    for(int i=0;i<n;++i){
        ans[i]=cut_double(A[i][n]/A[i][i]);
    }
    return ans;
}

vector<double> Gauss_ColPivot(vector<vector<double>> A)
{
    int n=A.size();
    for(int i=0;i<n;++i)
    {
        int col_max=A[i][i], index=i;
        for(int j=i+1;j<n;++j)         //从i+1行开始选列主元
        {
            if(fabs(A[j][i]>col_max))      //相比高斯消元增加了选择列元和交换两行元素
            {
                col_max=fabs(A[j][i]);
                index=j;
            }
        }
        swap(A[i],A[index]);
        for(int j=0;j<n;++j)           
        {
            if(j==i) continue;
            double r=cut_double(A[j][i]/A[i][i]);
            for(int k=i;k<=n;++k)
            {
                A[j][k]-=r*A[i][k];
                A[j][k]=cut_double(A[j][k]);
            }
        }
    }
    vector<double> ans(n);
    for(int i=0;i<n;++i){
        ans[i]=cut_double(A[i][n]/A[i][i]);
    }
    return ans;
}

vector<double> Doolittle(vector<vector<double>> A)
{
    int n=A.size();
    vector<vector<double>> L(n,vector<double>(n)),U(n,vector<double>(n));
    for(int i=0;i<n;++i){
        L[i][i]=1;              //下三角对角元素为1        
    }
    for(int i=0;i<n;++i)        //计算U第i行 和 L第i列 各个元素
    {
        for(int j=i;j<n;++j)    //U[i]的第j列前面的都是0，从第j列开始算,U[i][j]中i<=j
        {
            double s=0;              //s=L[i][]*U[][j]
            for(int k=0;k<i;++k)    //算s的时候U的第j列的i行之后为0
                s+=L[i][k]*U[k][j];
            U[i][j]=A[i][j]-s; 
        }
        for(int j=i+1;j<n;++j)   //L[][i]的j行前面是1或0，不用计算
        {
            double s=0;
            for(int k=0;k<j;++k)      //L[i][j]中i>=j,由j控制相乘长度
                s+=L[i][k]*U[k][j];
            L[i][j]=(A[i][j]-s)/U[j][j];
        }
    }
    vector<double>Y(n),X(n);
    for(int i=0;i<n;++i){
        double b=A[i][n];
        for(int j=0;j<i;++j)
            b-=L[i][j]*Y[j];
        Y[i]=b;
    }
    for(int i=n-1;i>=0;--i){
        double y=Y[i];
        for(int j=n-1;j>i;--j)
            y-=U[i][j]*X[j];
        X[i]=y/U[i][i];
    }
}


vector<double> Crout(vector<vector<double>> A)
{
    int n=A.size();
    vector<vector<double>> L(n,vector<double>(n)),U(n,vector<double>(n));
    for(int i=0;i<n;++i){
        U[i][i]=1;              //上三角对角元素为1        
    }
    for(int i=0;i<n;++i)        //计算 L第i列 和 U第i行 各个元素
    {
        for(int j=i;j<n;++j)    //L的 i列 的 第i行前面的都是0，j从i开始，L[i][j]中i>=j
        {
            double s=0;              //s=L[j][]*U[][i]
            for(int k=0;k<i;++k)    //算s的时候算到L[j][i-1]相乘即可
                s+=L[j][k]*U[k][i];
            L[j][i]=A[j][i]-s; 
        }
        for(int j=i+1;j<n;++j)   //U i行的i+1列 的前面是1或0，不用计算
        {
            double s=0;
            for(int k=0;k<i;++k)      //算U[i][j]乘到U[][j-1]
                s+=L[i][k]*U[k][j];
            U[i][j]=(A[i][j]-s)/L[i][i];
        }
    }
    vector<double> Y(n),X(n);
    for(int i=0;i<n;++i)
    {
        double b=A[i][n];
        for(int j=0;j<i;++j) 
            b-=L[i][j]*Y[j];
        Y[i]=b/L[i][i];
    }
    for(int i=n-1;i>=0;--i){
        double y=Y[i];
        for(int j=n-1;j>i;--j)
            y-=U[i][j]*X[j];
        X[i]=y;
    }
}

signed main()
{
    vector<vector<double>> A={{1,2,1,1},{3,1,1,2},{4,0,5,3}};
    vector<double> ans =Gauss(A);
    for(int i=0;i<(int)ans.size();++i) 
        cout<<ans[i]<<" ";
    cout<<endl;
    system("pause");
    return 0;
}