全选主元Gauss 消元法

这个算法是健壮的（数学术语）。我们通过选主元的方法来弥补顺序Gauss消元法的弊端，增加了pivrow和pivcol来记录包含主元素的行和列。由线性代数行列变换的知识知道，增广益矩阵行变换时，方程的解不变，但是列变换时，未知数的顺序发生了交换，因此需要记录列交换的信息在pivcol中，以便求解后还原。

全选主元高斯消元法的代码如下：

---

 

#include<iostream.h> #include<fstream.h> #include<math.h> #include<process.h> class gauss { private: int i,j,k,n,*pivrow,**pivcol; double eps,pivot,aik,al,sum,*x,**a; public: void gauss_input(); void gauss_elimination(); void gauss_output(); ~gauss() { delete [] x; delete [] pivrow; for(i = 0;i < n; i++) { delete [] a[i]; } delete [] a; for(i=0;i<n;i++) { delete [] pivcol[i]; } } }; void main() { gauss solution; solution.gauss_input(); solution.gauss_elimination(); solution.gauss_output(); } void gauss::gauss_input() { ifstream fin("gauss.dat"); fin>>n; x=new double[n]; pivrow = new int[n]; pivcol = new int *[n]; for(i=0;i<n;i++) { pivcol[i] = new int[2]; } a=new double*[n]; for(i=0;i<n;i++) { a[i] = new double[n+1]; } for(i=0;i<n;i++) for(j=0;j<n;j++) { fin>>a[i][j]; } for(i=0;i<n;i++) { fin>>a[i][n]; } fin>>eps; //输入段结束 fin.close(); } void gauss::gauss_elimination() //gauss 消去函数 { for(k=0;k<(n-1);k++) { pivot = fabs(a[k][k]); pivrow[k] = k; pivcol[k][0] = k; pivcol[k][1] = k; //下面是扫描所有行的最大值，并保存位置信息 for(i=k;i<n;i++) for(j=k;j<n;j++) { if(pivot<fabs(a[i][j])) { pivot = fabs(a[i][j]); pivrow[k] = i; pivcol[k][1] = j; } } if(pivot<eps) { cout<<"/n矩阵可能奇异.失败..."<<endl; exit(0); } //交换行 if(pivrow[k]!=k) { for(j=k;j<(n+1);j++) { al = a[pivrow[k]][j]; a[pivrow[k]][j] = a[k][j]; a[k][j] = al; } } if(pivcol[k][1]!=k)//交换列 { for(i=0;i<n;i++) { al = a[i][pivcol[k][1]]; a[i][pivcol[k][1]] = a[i][k]; a[i][k] = al; } } if(k!=(n-1)) { for(i=(k+1);i<n;i++) { aik = a[i][k]; for(j=k;j<(n+1);j++) { a[i][j] -= aik*a[k][j]/a[k][k]; } } } } x[n-1] = a[n-1][n]/a[n-1][n-1];//回代 for(i=(n-2);i>=0;i--) { sum = 0.0; for(j=(i+1);j<n;j++) { sum += a[i][j] * x[j]; } x[i] = (a[i][n] - sum)/a[i][i]; } for(k=(n-2);k>=0;k--) { al = x[pivcol[k][1]]; x[pivcol[k][1]]=x[pivcol[k][0]]; x[pivcol[k][0]] = al; } } void gauss::gauss_output() //结果输出函数 { ofstream fout("gauss.out"); fout<<"/n结果是："<<endl; for(i=0;i<n;i++) { fout<<"/nx["<<i<<"]="<<x[i]<<endl; } fout.close(); for(i=0;i<n;i++) { cout<<"/nx["<<i<<"]="<<x[i]<<endl; } }