高斯消元法:(列主元高斯消元法)
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;
}