#define N 3  /* N為方程組系數矩陣的階數 */
void Lu(float a[N][N],float L[N][N],float U[N][N]) /* LU分解 */
{  int i,j,k;
   for(k=0;k<N;k++)
   {  for(i=k;i<N;i++)  /* 計算L矩陣的第k列元素 */
      {  L[i][k]=a[i][k];
         for(j=0;j<=k-1;j++)
	    L[i][k]-=(L[i][j]*U[j][k]);
      }
      for(j=k+1;j<N;j++)  /* 計算U矩陣的第k列元素 */
      {  U[k][j]=a[k][j];
         for(i=0;i<=k-1;i++)
           U[k][j]-=(L[k][i]*U[i][j]);
	 U[k][j]/=L[k][k]; 
      }
   }
}
void solve(float b[N],float L[N][N],float U[N][N],float x[N])
{  int j,k;
   float y[N];
   for(k=0;k<N;k++)  /* 計算Ly=b中的y */
   {  y[k]=b[k];
      for(j=0;j<=k-1;j++)
        y[k]=y[k]-L[k][j]*y[j];
      y[k]=y[k]/L[k][k];
   }
   for(k=N-1;k>=0;k--) /* 計算Ux=y中的x */
   {  x[k]=y[k];
      for(j=k+1;j<N;j++)
        x[k]=x[k]-U[k][j]*x[j];
   }
}
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");
}
main()
{  static float a[N][N]={{1,2,-1},{1,-1,5},{4,1,-2}};
   float b[N]={3,0,2};
   float x[N]={0,0,0};
   float U[N][N]={{1,0,0},{0,1,0},{0,0,1}};
   float L[N][N];
   int i;
   zg_matric(a,b);
   Lu(a,L,U);
   solve(b,L,U,x);
   for(i=0;i<N;i++)  /* 輸出方程組的解 */
     printf("  x%d=%11.7f\n",i+1,x[i]);
}