高斯消元

//高斯消元法
#include <stdio.h>
#include <math.h>

#define n 3

void GSXiaoYuanFun(float c);

void main()
{
	GSXiaoYuanFun(1e-3);
}

void GSXiaoYuanFun(float c)
{
	float m,M=0,xx;
	float d,x[n];
	int k,i,j,l,q;
	float a[3][3]={{2,1,1},{1,3,2},{1,2,2}};
	float b[3]={4,6,5};
	printf("消元后的矩阵是：\n");
	for(k=0;k<n;k++)
	{
		d=a[k][k];
		for(i=k+1;i<n;i++)
		{
			m=a[i][k]/a[k][k];
			b[i]=b[i]-m*b[k];
			if(fabs(a[i][k])>fabs(d))
			{
				d=a[i][k];
				l=i;
			}
			if(a[k][k]<c)
			{
				for(q=0;q<n;q++)
				{
					xx=a[k][q];a[k][q]=a[l][q];a[l][q]=xx;
				}
			}/*列主元消元*/

			for(j=0;j<n;j++)
			{
				a[i][j]=a[i][j]-m*a[k][j];
			}
		}
	}/*消元成三角矩阵*/
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
			printf("%-10f",a[i][j]);
		printf("%f\n",b[i]);
	}

	/*解方程*/
	k=n-1;
	x[k]=b[k]/a[k][k];
	printf("\tx[%d]=%f\n",n,x[k]);
	for(i=k-1;i>=0;i--)
	{
		for(j=i+1;j<=k;j++)
		{
			M=M+(a[i][j]*x[j]);
		}
		x[i]=(b[i]-M)/a[i][i];
		M=0;
		printf("\tx[%d]=%f\n",i+1,x[i]);
	}
}