列主元高斯消去法

#include<stdio.h>
#include<math.h>
#define N 4 /*定义方程组系数矩阵的阶数N*/
static float a[N][N]={{1.5,0.5,2.5,-3},{6,0.6,0,-5},{1,1.2,14,1},{-3,-5,1,11.1}};/*输入被求解线性方程组*/
static float b[N]={14.4095,16.166,24.628,-5.8256};
float x[N]={0,0,0,0}; /*定义初始值*/
void ColGauss(float a[N][N],float b[N])/*列主元高斯消去法计算核心*/
{
       float t,max_fab;
       int i,j,k,l;
       for(i=0;i<N-1;i++)
       {
              l=i;
              max_fab=fabs(a[i][i]);
              for(j=i+1;j<N;j++) /* 找主元素*/
                     if(fabs(a[j][i]>max_fab))
                     {
                            l=j;
                            max_fab=fabs(a[j][i]);
                     }                      if(i<l) /* 交换i、l两行*/
                     {
                            t=b[i];
                            b[i]=b[l];
                            b[l]=t;
                            for(j=i;j<N;j++)
                            {
                                   t=a[i][j];
                                   a[i][j]=a[l][j];
                                   a[l][j]=t;
                            }
                     }
                     for(j=i+1;j<N;j++) /* 消元过程开始*/
                     {
                            t=-a[j][i]/a[i][i];
                            b[j]=b[j]+t*b[i];
                            for(k=i;k<N;k++)
                                   a[j][k]=a[j][k]+t*a[i][k];
                     }
       }
} void augm_matrix(float a[N][N],float b[N]) /* 输出增广矩阵*/
{
       int i,j;
       printf("Linear equations:\n");
       for(i=0;i<N;i++)
       {
              for(j=0;j<N;j++)
                     printf("%10f",a[i][j]);
              printf("%10f",b[i]);
              printf("\n");
       }
       printf("\n");
} main()
{
       int i,j;
       augm_matrix(a,b);
       ColGauss(a,b);
       x[N-1]=b[N-1]/a[N-1][N-1]; /* 回代过程*/
       for(i=N-2;i>=0;i--)
       {
              x[i]=b[i];
              for(j=i+1;j<N;j++)
                     x[i]=x[i]-a[i][j]*x[j];
              x[i]=x[i]/a[i][i];
       }
       printf("Results:\n");
       for(i=0;i<N;i++) /* 输出方程组的解*/
              printf(" x%d=%f\n",i+1,x[i]);
       system("pause");
}

参考资料:
[1] 王汉青. 暖通空调流体流动数值计算方法与应用[M]. 北京: 科学出版社, 2013.10.