www.pudn.com > cghost.rar > GAGAUS.C

/****************************高斯消去法******************/ 
/a存放系数矩阵;b存放右端常向量;x返回方程组的解向量.*****/ 
#include "stdio.h" 
#include "stdlib.h" 
#include "math.h" 
int cagaus(a,b,n,x) 
int n; 
double a[],b[],x[]; 
{int *js,l,k,i,j,is,p,q; 
double d,t; 
js=malloc(n*sizeof(int)); 
l=1; 
for(k=0;k<=n-2;k++) 
{d=0.0; 
for(i=k;i<=n-1;i++) 
for(j=k;j<=n-1;j++) 
{t=fabs(a[i*n+j]); 
if(t>d){d=t;js[k]=j;is=i;} 
} 
if(d+1.0==1.0)l=0; 
else 
  {if(js[k]!=k) 
    for(i=0;i<=n-1;i++) 
    {p=i*n+k;q=i*n+js[k]; 
    t=a[p];a[p]=a[q];a[q]=t; 
    } 
   if(is!=k) 
     {for(j=k;j<=n-1;j++) 
       {p=k*n+j;q=is*n+j; 
       t=a[p];a[p]=a[q];a[q]=t; 
       } 
       t=b[k];b[k]=b[is];b[is]=t; 
       } 
   } 
if(l==0) 
  {free(js);printf("fail\n"); 
   return(0); 
   } 
d=a[k*n+k]; 
for(j=k+1;j<=n-1;j++) 
  {p=k*n+j;a[p]=a[p]/d;} 
b[k]=b[k]/d; 
for(i=k+1;i<=n-1;i++) 
  {for(j=k+1;j<=n-1;j++) 
    {p=i*n+j; 
     a[p]=a[p]-a[i*n+k]*a[k*n+j]; 
     } 
   b[i]=b[i]-a[i*n+k]*b[k]; 
  } 
} 
d=a[(n-1)*n+n-1]; 
if(fabs(d)+1.0==1.0) 
  {free(js);printf("fail\n"); 
    return(0); 
  } 
x[n-1]=b[n-1]/d; 
for(i=n-2;i>=0;i--) 
  {t=0.0; 
   for(j=i+1;j<=n-1;j++) 
     t=t+a[i*n+j]*x[j]; 
     x[i]=b[i]-t; 
  } 
js[n-1]=n-1; 
for(k=n-1;k>=0;k--) 
  if(js[k]!=k) 
    {t=x[k];x[k]=x[js[k]];x[js[k]]=t;} 
free(js); 
return(1); 
} 
/**************************************************/ 
main() 
{int i; 
static double a[3][3]= 
       {{1.,-1.,1.},{5.,-4.,3.},{2.,1.,1.}}; 
static double x[3],b[3]= 
       {-4.,-12.,11}; 
i=cagaus(a,b,3,x); 
if(i!=0) 
  for(i=0;i<=2;i++) 
    printf("x(%d)=%e\n",i+1,x[i]); 
}