void ColGauss(float a[N][N],float b[N])
{  float t,max_fab;
   int i,j,k,l;
   for(i=0;i<N-1;i++)
   {  l=i;
      max_fab=fabs(a[i][i]);
      for(j=i+1;j<N;j++)  /* 找主元素 */
	     if(fabs(a[j][i]>max_fab)) 
        {  l=j;  max_fab=fabs(a[j][i]);   }
      if(i<l)  /* 交換i、l兩行 */
      {  t=b[i];
         b[i]=b[l];
         b[l]=t;
         for(j=i;j<N;j++)
         {  t=a[i][j]; a[i][j]=a[l][j]; a[l][j]=t; }
      }
      for(j=i+1;j<N;j++) /* 消元過(guò)程開始 */
      {  t=-a[j][i]/a[i][i];
         b[j]=b[j]+t*b[i];
         for(k=i;k<N;k++)
	        a[j][k]=a[j][k]+t*a[i][k];  }   }
}
void zg_matric(float a[N][N],float b[N]) /* 輸出增廣矩陣 */
{  int i,j;
   for(i=0;i<N;i++)
   {  for(j=0;j<N;j++)
	      printf("%10f",a[i][j]);
      printf("%10f",b[i]);
      printf("\n");
   }
   printf("\n");
} 
#include <math.h>
#define N 4  /* N為方程組系數(shù)矩陣的階數(shù) */
main()
{  float a[N][N]={{1.003,0.333,1.504,-0.333},{-2.011,1.455,0.506,2.956},
                  {4.329,-1.952,0.006,2.087},{5.113,-4.004,3.332,-1.112}};
   float b[N]={3.005,5.407,0.136,3.772};
   float x[N]={0,0,0,0};
   int i,j;
   zg_matric(a,b);
   ColGauss(a,b); 
   x[N-1]=b[N-1]/a[N-1][N-1]; /* 回代過(guò)程開始 */
   for(i=N-2;i>=0;i--)
   {  x[i]=b[i];
      for(j=i+1;j<N;j++)
	     x[i]=x[i]-a[i][j]*x[j];
      x[i]=x[i]/a[i][i];
   }
   for(i=0;i<N;i++) /* 輸出方程組的解 */
     printf("  x%d=%11.7f\n",i+1,x[i]);
}