#include<stdio.h>
#include<math.h>
#define M 4

void gauss(double a[][M]); /*高斯變換*/
void slove(double [][M],double [M-1]);/*回代*/
int colmax(double [][M],int ); /*行變換*/
void print(double [][M]);
void print(double [M-1]);

void main()
{double a[3][M]={{0.00000001,2,3,1},{-1,3.712,4.623,2},{-2,1.072,5.643,3}};
 double x[M-1];
 printf("原來的矩陣為：\n");
 print(a);
 gauss(a);
 printf("用Guass列主元素消去法后的矩陣為：\n");
 print(a);
 slove(a,x); 
}

void gauss(double a[][M])/*高斯變換*/
{double x;
 int i,j,k;
 for (i=0;i<M-1;i++)
    {k=colmax(a,i);
     if (k!=i)
	for (j=i;j<M;j++)
	    {x=a[i][j];a[i][j]=a[k][j];a[k][j]=x;}
     for (j=i+1;j<M-1;j++)
	 {x=a[j][i]/a[i][i];
	  for (k=i;k<M;k++)
	      a[j][k]=a[j][k]-x*a[i][k];	   
	 }
    }
}

void slove(double a[][M],double x[M-1])/*回代*/
{double sum;
 int i,j;
 x[M-2]=a[M-2][M-1]/a[M-2][M-2];
 for (i=M-3;i>=0;i--)
     {sum=0;
      for (j=i+1;j<M-1;j++) sum=sum+a[i][j]*x[j];
      x[i]=(a[i][M-1]-sum)/a[i][i];
     }
 printf("距陣的解為:\n");
 for (i=0;i<M-1;i++)
 {printf("x[%d]=%.9f",i+1,x[i]);
 printf("\n");}
}

int colmax(double a[][M],int row)/*行變換*/
{double x;
 int i,k;
 x=fabs(a[row][row]); k=row;
 for (i=row+1;i<M-1;i++)
     if (fabs(a[i][row])>x)
       {x=fabs(a[i][row]);k=i;}
 return k;
}

void print(double a[][M])
{int i,j;
 for (i=0;i<M-1;i++)
     {for (j=0;j<M;j++)
      printf(" %12.9f",a[i][j]);
      printf("\n");
     }
 printf("\n");
}