// 帶有列主元的高斯消元法
// 功能：　求解線性方程組 Ax = b
// 參數：　A － 指向n*n系數矩陣的指針
// 　　　　b － 常數向量的指針
// 　　　　n － 方程組的維數
// 返回值：0 － 如果成功。線性方程組的解保存在 b 中
// 　　　　1 － 求解失敗
// ----------------------------------------------------------------------------
# include "math.h"
# include "stdio.h"
int LEquations_Gauss(double A[], double b[], int n)
{
	int i, j, k;
	double s, t;

	// 高斯列主元消元,使 Ax=b -> Ux=b'
	for(k=0; k<n; k++)
	{
		// 選取列主元素
		j = k;  t = fabs(A[k*n+k]);
		for(i=k+1; i<n; i++)
		{
			if((s = fabs(A[i*n+k])) > t)
			{
				t = s;
				j = i;
			}
		}
		if(t < 1.0e-30) return 1;	// 方程為病態,或無窮多解

		if(j != k)			// 交換方程的第 j 行和 k 行
		{
			for(i=k; i<n; i++)
			{
				t = A[j*n+i];
				A[j*n+i] = A[k*n+i];
				A[k*n+i] = t;
			}
			t = b[j];    b[j] = b[k];    b[k] = t;
		}

		t = 1.0 / A[k*n+k];		// 先將對角線元素變為 1.0
		for(i=k+1; i<n; i++) A[k*n+i] *= t;
		b[k] *= t;

		for(i=k+1; i<n; i++)		// 依次消元第k+1到n-1行的第k列元素
		{
			t = A[i*n+k];
			for(j=k+1; j<n; j++) A[i*n+j] -= A[k*n+j] * t;
			b[i] -= b[k] * t;
		}
	}

	// 回代,求解 Ux=b'
	for(i=n-2; i>=0; i--)
		for(j=i+1; j<n; j++) b[i] -= A[i*n+j] * b[j];

	return 0;
}




main()
{double A[]={1,2,1,-2,2,5,3,-2,-2,-2,3,5,1,3,2,3};
 double b[]={4,7,-1,0}; 
 int n=4,i;
  LEquations_Gauss(A,b,n);
	 for(i=0;i<4;i++)
		 printf("%f\n",b[i]);
}