1.cpp

求解线性方程组的高斯主元消去法程序

//求解线性方程组的高斯主元消去法子程序
#include<stdlib.h>
#include<math.h>
#include<stdio.h>
#include<conio.h>
#define brow 4 //方程组未知数的个数
int n;
double a[],b[]; //ax=b,a为方程组的系数用一维数组表示，b为方程右端数据
int agaus(double *a,double *b,int n)//全选主元高斯消去法解方程
{int *js,l,k,i,j,is,p,q;
double d,t;
js=(int *)malloc(n*sizeof(int));//堆内存分配，指针指向该内存的地址返回，该内存的内容未知。
l=1;

for(k=0;k<=n-2;k++) //方程组系数数组a中选取绝对值最大的元素作为主元素，
{d=0.0;            // 交换行使得具有最大绝对值的元素被交换到第k行，k列。
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=1;
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]; //进行归一化计算，a(kj)=a(kj)/a(kk),j=k+1,k+2,...,n-1
for(j=k+1;j<=n-1;j++)//       b(k)=b(k)/a(kk)
{p=k*n+j;a[p]=a[p]/d;}
b[k]=b[k]/d;

for(i=k+1;i<=n-1;i++)//进行消去计算：a(ij)=a(ij)-a(ik)a(kj)   i,j=k+1,k+2,...,n-1
{for(j=k+1;j<=n-1;j++)//             b(i)=b(i)-a(ik)b(k)
{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);
}

b[n-1]=b[n-1]/d;//回代过程:x(n-1)=b(n-1)/a(n-1,n-1)
for(i=n-2;i>=0;i--)
{t=0.0;
for(j=i+1;j<=n-1;j++)
t=t+a[i*n+j]*b[j];//x(i)=b(i)-sum(j=i+1,...,n-1)a(ij)x(j)
b[i]=b[i]-t;}

js[n-1]=n-1;
for(k=n-1;k>=0;k--)
if(js[k]!=k)
{t=b[k];b[k]=b[js[k]];b[js[k]]=t;}
free(js);
return(1);
}

//主函数选择了一个四元方程组为例来验证算法。
//0.2368x0+0.2471x1+0.2568x2+1.2671x3=1.8471
//0.1968x0+0.2071x1+1.2168x2+0.2271x3=1.7471
//0.1581x0+1.1675x1+0.1768x2+0.1871x3=1.6471
//1.1161x0+0.1254x1+0.1397x2+0.1490x3=1.5471
//x(0)=1.040577
//x(1)=0.987051
//x(2)=0.935040
//x(3)=0.881282
int main() 
{int i;
static double a[16]={0.2368,0.2471,0.2568,1.2671,0.1968,0.2071,1.2168,0.2271,
0.1581,1.1675,0.1768,0.1871,1.1161,0.1254,0.1397,0.149};
static double b[4]={1.8471,1.7471,1.6471,1.5471};
if(agaus(a,b,brow)!=0)
for(i=0;i<brow;i++)  //输出求解结果
printf("x(%d)=%f\n",i,b[i]);
getch();
return 0;
}