2019-05-19 线性方程组-高斯消元法

高斯消元法：（列主元高斯消元法）
1.正常的情况下，我们将矩阵化成上三角，然后化成行最简形式，这里算法略微不同；
2.这里处理的时候，每循环一次，消除一个未知数，没有意义的计算直接省略；

const double EPS=1e-8;
typedef vector<double>vec;
typedef vector<vec>mat;

//计算Ax=B
//当方程组无解或者无穷解时候,返回一个长度为0的数组 
vec gauss_jordan(const mat&A,const vec&b)
{
    int n=A.size();
    mat B(n,vec(n+1));//增广矩阵，只能对增广矩阵做初等行变换，得到的新的方程组同解
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++) B[i][j]=A[i][j];
    for(int i=0;i<n;i++) B[i][n]=b[i];
    
    for(int i=0;i<n;i++){
        int privot=i;
        for(int j=i;j<n;j++){
            if(abs(B[j][i])>abs(B[pivot][i])) pivot=j;
        }
        swap(B[i],B[pivot]);
        //无解或者无穷解 
        if(abs(B[i][i])<EPS) return vec();
        
        //这里只想得到方程组的解，因此很多后面用不到的结果直接不去计算！！ 
        for(int j=i+1;j<=n;j++) B[i][j]/=B[i][j];//省略了B[i][i]/=B[i][j];
        for(int j=0;j<n;j++){
            if(i!=j){
                for(int k=i+1;k<=n;k++) B[j][k]-=B[i][k]*B[j][i];
            } 
        } 
    }
    
    vec x(n);
    for(int i=0;i<n;i++) x[i]=B[i][n]; 
    return x;
}