/* gcpelimm.c:use Gauss elimination method to solve linear system.這是主函數*/
#define DIM 5
#include "gcpelim.c"
extern int gelim(int process,double A[DIM][DIM],double b[DIM]);
void main(void)
{
int i,j;
static double A[DIM][DIM]={{2.0,-1.0,4.0,-3.0,1.0},
                           {-1.0,1.0,2.0,1.0,3.0},
                           {4.0,2.0,3.0,3.0,-1.0},
                           {-3.0,1.0,3.0,2.0,4.0},
                           {1.0,3.0,1.0,4.0,4.0}};
static double b[DIM]={11.0,14.0,4.0,16.0,18.0};
system("cls");
if(gcpelim(1,A,b)==1)
  {
   printf("The linear system has't solution!\n");
   printf("Strike key to exit !\n"); getch();  exit(1);
   }
   printf("Column principle elimination for the sulotion:\n");
for(i=0;i<DIM;i++)
  printf("  %10.6f\n",b[i]);
  getch();
}


--------------------------------------------
/*gcpelim.c:gauss列主元消去法的子函數*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#ifndef EPSILON
#define EPSILON 0.000000001
#endif
int gcpelim(int process,double A[DIM][DIM],double xx[DIM])
{
int k,i,j,i0;
double pelement;
if(process==1)  printf("The process of elimination\n");
/*  elimination step  */
for(k=0;k<DIM;k++)
{ 
   /*for principal element */
pelement=fabs(A[k][k]);i0=k;
for(i=k;i<DIM;i++)
   if(fabs(A[i][k])>pelement) {pelement=fabs(A[i][k]);i0=i;}
if(i0!=k)
{
  for(j=0;j<DIM;j++)
   {
    pelement=A[k][j];A[k][j]=A[i0][j];A[i0][j]=pelement;}
    pelement=xx[k];xx[k]=xx[i0];xx[i0]=pelement;
    }
if(fabs(A[k][k])<EPSILON) return(1);
for(i=k+1;i<DIM;i++)
  {
     A[i][k]=A[i][k]/A[k][k];
     for(j=k+1;j<DIM;j++)  A[i][j]=A[i][j]-A[i][k]*A[k][j];
     xx[i]=xx[i]-A[i][k]*xx[k];
      }
     if(process==1)
       {
       for(i=0;i<DIM;i++)
       {
       for(j=0;j<DIM;j++)  printf("%10.6f",A[i][j]);
       printf("  |%10.6f\n",xx[i]);
        }
       printf("\n");
        }
       }
  /* backward step  */

  for(i=DIM-1;i>=0;i--)
    {
     for(j=i+1;j<DIM;j++)  xx[i]=xx[i]-A[i][j]*xx[j];
     xx[i]=xx[i]/A[i][i];
     }
   return(0);
   }