gauss_jordan_elimination.c

使用高斯消去法Gauss_Jordan解方程式

/*CH06:Gauss-Jordan Elimination
  T.B. Chen  2002 08 03
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define dim 10
double a[dim+1][dim+1],b[dim+1];

main()
{
int 	n,i,j,k,pivot[dim+1],ipvt_store;
double 	temp,pvt;

printf("input the dimension of matrix\n");
scanf("%d",&n);
 if(n>dim){
  printf("dimension over preset value %d\n",dim);
  exit(0);
   }
printf("Input the a matrix\n\n");
 for(i=1; i<=n; i++){
   for(j=1;j<=n;j++){
     printf("input the element of a[%d,%d]\n",i,j);
     scanf("%lf",&a[i][j]);
    }
  }
printf("\n Input the b vector\n\n");
  for(i=1; i<=n ; i++){
    printf("input the value of b[%d]\n",i);
    scanf("%lf",&b[i]);
   }
printf("\n Print out Designed matrix:\n");
mprint(n,n,a,b);

/* start to the procedure of Gauss_Jordan Elimination */
  for(j=1; j<= n; j++){
/*Search for the maximum element */
  pvt=fabs(a[j][j]);
  pivot[j]=j;
  ipvt_store=j;
 for(i=j+1; i<=n; i++) {
    if(fabs(a[i][j])>pvt){
      pvt=fabs(a[i][j]);
      ipvt_store=i;
     }
   }
/* Row change */
  if(pivot[j]!=ipvt_store){
     pivot[j]=ipvt_store;
     pivot[ipvt_store]=j;
    for(k=1; k<=n;k++){
	temp=a[j][k];
	a[j][k]=a[pivot[j]][k];
	a[pivot[j]][k]=temp;
     }
	temp=b[j];
	b[j]=b[pivot[j]];
	b[pivot[j]]=temp;
  }
/*Make the diagonal element of key row unit */
temp=a[j][j];
   for(k=j;k<=n;k++)  a[j][k]=a[j][k]/temp;
	b[j]=b[j]/temp;
/*