高斯消元模板(列主元)

这是在 人民邮电出版社 出版  <挑战程序设计竞赛>第二版  P287 抄录的 模板.  加了一些自己对该模板的思考.

</pre><pre name="code" class="cpp"><pre name="code" class="cpp">#include<stdio.h>
#include<vector>
using namespace std;
typedef vector<double> vec;
typedef vector<vec> mat;
const double EPS =1e-8;
int abb(int a)
{
	if(a<0)
		return -a;
	return a;
}
//方程无解时 或无穷多解时  返回长度为0的数组
vec gauss_jordan(const mat&A,const vec&b)//A是方阵   
{
	int n=A.size();//几行几列
	mat B(n,vec(n+1));//扩展矩阵  ,size=n 所以n行   n+1列,容器内每个单元大小都为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 pivot=i;  //关键,  也就是记录当前处理i行时 把i行以及其后的i列元素绝对值最大的哪一行,换到 i行
		for(int j=i;j<n;j++)
		{
			if(abb(B[j][i])>abb(B[pivot][i]))
				pivot=j;
		}
		swap(B[i],B[pivot]);
	    if(abb(B[i][i])<EPS) return vec();

		for(int j=i+1;j<=n;j++) B[i][j]/=B[i][i];//i行i列为1;  //i行i列一定会被处理为1  所以暂不处理 因为处理后面的数据还要用到它 
		for(int j=0;j<n;j++)// 在除了当前处理行外的每一行
		{
			if(i!=j)  //相当于 i行乘 j行i列 元素.j行减i行; i行再除回 原j行i列;  
			{
				for(int k=i+1;k<=n;k++) B[j][k]-=B[j][i]*B[i][k];//j行,i列及其之前以及化为零了 所以默认不处理.
                                                                                //而其后元素会影响到下面的计算所以要处理
			}
		}
	}
	/*---------------------------------①
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			if(i!=j)
				B[i][j]=0;
			else
				B[i][j]=1;
		}
	}
	*/
	vec x(n);
	for (int i=0;i<n;i++)
		x[i]=B[i][n];
	/*
	for(int i=0;i<n;i++)//输出处理后的B矩阵
	{
		for(int j=0;j<n+1;j++)
			printf("%.3lf\t",B[i][j]);
		puts("");
	}
	*/
	return x;
}
int main()
{
	vec ans;
	int n,i,j;
	double a;
	
	mat A;
	while(scanf("%d",&n)!=EOF)//AX=b
	{
		mat A(n,vec(n));
		vec b(n);
		for(i=0;i<n;i++)
    	{
			for(j=0;j<n;j++)
			{
				scanf("%lf",&a); //输入矩阵A
				A[i][j]=a;
			}
		}
	    for(i=0;i<n;i++)//输入矩阵b
		{
			scanf("%lf",&a);
			b[i]=a;
		}
		ans=gauss_jordan(A,b);
		for(int i=0;i<n;i++)//输出答案
			printf("%.3lf\t",ans[i]);
		puts("");
	}
	return 0;
}

/*
输入 
3
1 2 3
1 4 9
1 8 27
1 2 3  


如果①处不加 
B矩阵会被处理成 
1.000   2.000   -5.000  -0.500
1.000   6.000   4.000   1.000
1.000   2.000   -2.000  -0.167
最后一排即是答案  x1=-0.5  x2=1 x3=-0.167


如果加了
B矩阵会处理成
1.000   0.000   0.000   -0.500
0.000   1.000   0.000   1.000
0.000   0.000   1.000   -0.167
*/