#include<stdio.h>

#include<math.h>

#define N 3

#define EPSILON 1.0E-5

#define TRUE 1

#define FALSE 0

/*** a[][]:[coefficient matrix,constant vector] ***/

double a[N][N+1]={{ 4.0, 2.0 , 1.0, 11.0 },

{ 3.0, 1.500001, 2.0, 12.000002 },

{ 5.0, 1.0 , 3.0, 16.0 }};

int sw=TRUE;

main(){

int i,j;

FILE *fp;

void sweep(int *);

fp=fopen("E:20060131.data","w");

printf("\n\n ====== Coefficient matrix and constant vector ====== \n\n");

fprintf(fp,"\n\n ====== Coefficient matrix and constant vector ====== \n\n");

for(i=0;i<N;i++){

for(j=0;j<N;j++){

printf(" %10.6f",a[i][j]);

fprintf(fp," %10.6f",a[i][j]);

}

printf(" %10.6f\n",a[i][N]);

fprintf(fp," %10.6f\n",a[i][N]);

}

sweep(&sw);

printf("\n\n ====== Computed matrix ====== \n\n");

fprintf(fp,"\n\n ====== Computed matrix ====== \n\n");

for(i=0;i<N;i++){

for(j=0;j<N;j++){

printf(" %10.6f",a[i][j]);

fprintf(fp," %10.6f",a[i][j]);

}

printf(" %10.6f\n",a[i][N]);

fprintf(fp," %10.6f\n",a[i][N]);

}

if(sw=TRUE){

printf("\n\n ====== Roots of simultaneous linear equations ====== \n\n");

fprintf(fp,"\n\n ====== Roots of simultaneous linear equations ====== \n\n");

for(i=0;i<N;i++){

printf(" x[%d] = %10.6f \n",i+1,a[i][N]);

fprintf(fp," x[%d] = %10.6f \n",i+1,a[i][N]);

}

}

else{

printf("\n\n Compute is not able,because pivot is very small or 0.\n\n");

fprintf(fp,"\n\n Compute is not able,because pivot is very small or 0.\n\n");

}

fclose(fp);

return(1);

}

void sweep(int *sw){

int i,j,k,ik;

double ak,aik;

void swap(double *,double *);

for(k=0;k<N;k++){

if(fabs(a[k][k])<=EPSILON){

ik=k+1;

while((ik<N) && (fabs(a[ik][k])<EPSILON))

ik++;

if(ik<N)

for(j=k;j<N+1;j++)

swap(&a[k][j],&a[ik][j]);

else{

*sw=FALSE;

goto end;

}

}

ak=a[k][k];

for(j=k;j<=N;j++)

a[k][j]=a[k][j]/ak;

for(i=0;i<N;i++)

if(i!=k){

aik=a[i][k];

for(j=k;j<=N;j++)

a[i][j]=a[i][j]-aik*a[k][j];

}

}

end:;

}

void swap(double *w1,double *w2){

double w;

w=*w1;*w1=*w2;*w2=w;

}