全选主元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;}}