rgauss.c

gauss顺序消去法及其全选主元的ＧＡＵＳＳ消去法的Ｃ语言实现

#include "stdio.h" 
#include "stdlib.h" 
#include "math.h" 
int rgauss(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[4][4]=
	{ {0.94455,-2.81716,-2.25806,-2.89224},
	{-6.20716,-3.85174,9.77783,3.49712},
	{-7.56157,-4.03668,9.38475,-8.64498},
	{-8.94467,-5.31114,-0.98666,-4.10932} };
	static double x[4], b[4]={1.51158,-6.46352,1.23264,0.61068};
i=rgauss(a,b,4,x); 
if(i!=0) 
  for(i=0;i<=3;i++) 
    printf("x(%d)=%lf\n",i,x[i]); 
}