C语言：选列主元约当消元法求线性方程组唯一解

采用选列主元约当消元法求线性方程组唯一解：输入4*5的线性方程组的满秩增广矩阵，输出方程组的解（使用双精度double类型精确到小数点后9位有效数）。

输入格式

m(整数，实验数据总数)

a11 a12 a13 a14 a15(实数，第1个增广矩阵的第1行5个数据)

a21 a22 a23 a24 a25(实数，第1个增广矩阵的第2行5个数据)

a31 a32 a33 a34 a35(实数，第1个增广矩阵的第3行5个数据)

a41 a42 a43 a44 a45(实数，第1个增广矩阵的第4行5个数据)

...

输出格式

 x1 x2 x3 x4(输出第1个方程组的解，精确到小数点后9位有效数)

 x1 x2 x3 x4(输出第2个方程组的解，精确到小数点后9位有效数)

...

 x1 x2 x3 x4(输出第m个方程组的解，精确到小数点后9位有效数)

输入样例

2

5.96912 -2.83547 -9.63744 8.13715 -0.46236

-9.31394 -7.37968 -3.62896 -0.62227 8.82260

-8.62545 2.22327 -9.32371 2.77932 2.70363

-1.63610 2.69204 -7.04398 -9.11191 5.51378

-9.97497 1.27171 -6.13392 6.17481 1.70019

-0.40254 -2.99417 7.91925 6.45680 4.93210

-6.51784 7.17887 4.21003 0.27070 -3.92010

-9.70031 -8.17194 -2.71096 -7.05374 -6.68203

对于本次实验主要是选择主元，然后采用约当消元法求线性方程组唯一解。首先说一下就是如何选择，我是采用选择排序的方式来选取出主元，就是在该列下一个个比较，如果数大那么就把该数所在的行与第k行交换，以得到主元，选择主元需要传递参数k，就是为第k个方程选取列主元，以克服当对角线上的值akk的绝对值很小对，后续计算产生影响。所以需要选取主元，在选主元是因为列主元相对于全主元简单易操作，所以采取选取列主元。如下：

void getm(int k)//选取最大主列
{
    double temp;
    for(int i=k+1; i<4; i++)
    {
        if(fabs(b[k][k])<fabs(b[i][k]))
        {
            for(int a=k; a<5; a++)
            {
                temp = b[k][a];
                b[k][a] = b[i][a];
                b[i][a] = temp;
            }
        }
    }
  return ;
}

选取列主元后需要进行约当消元，求出增广矩阵的解，计算线性方程组的解主要分为两步，先把其增广矩阵化为上三角矩阵，然后再求解上三角矩阵。首先化为上三角矩阵，既对K=1,2,…,M-1；将第K个方程两边同时除以xK项的系数aKK；对I=K+1,…,M ；方程(I)=方程(I)-[方程(K)]* aI K。最外层for循环是为了把前k（k<m）个方程全部计算出上三角方程组，内部嵌套循环第一个则是把第k个方程的第k个系数化为1，就是把对角线下面变成0，然后其项目两个循环则是把第k列的第k个后面全部变成0，既化上三角。

代码设计如：

void Uppertriangle()//化上三角
{
    double x;
    for(int i=0; i<3; i++)
    {
        getm(i);
        x = b[i][i];
        for(int a=i; a<5; a++)
        {

            b[i][a] = b[i][a]/x;//化1
        }
        for(int j=i+1; j<4; j++)
        {
            x = b[j][i];
            for(int a=i; a<5; a++)
            {
                b[j][a] = b[j][a]-b[i][a]*x;//化0
            }

        }
}
}

然后解上三角，主要是采用高斯消元法来计算。既第k个方程-（减去第k+1个方程*第k个方程的第k个系数,既b[k][i] = b[k][i]-b[k+1][i]*b[k][k+1],这里需要注意的是在实际上解完第k个未知数是，既对滴k个方程运算结束之后，其增广矩阵的第k+1列以及除了第k个系数为1其他已经全部为0了，既此时常数b[k][n]为对应的第k个未知数的解。代码设计如：

void Lowertriangle()//化下三角就是把对角线上面变成0
{
    double x;
    for(int i=3; i>0; i--)
    {
        x = b[i][i];
        for(int a=i; a<5; a++)
        {
            b[i][a] = b[i][a]/x;//把对角线上值系数变成1
        }
        for(int j=i-1; j>=0; j--)
        {
            x = b[j][i];
            for(int a=j; a<5; a++)
            {
                b[j][a] = b[j][a]-b[i][a]*x;//在上三角基础上把第i列的第i项之外其他变成0
            }
        }
    }
}

完整代码：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double b[4][5];

void getm(int k)//选取最大主列
{
    double temp;
    for(int i=k+1; i<4; i++)
    {
        if(fabs(b[k][k])<fabs(b[i][k]))
        {
            for(int a=k; a<5; a++)
            {
                temp = b[k][a];
                b[k][a] = b[i][a];
                b[i][a] = temp;
            }
        }
    }
  return ;
}


void Uppertriangle()//化上三角就是把对角线下面变成0
{
    double x;
    for(int i=0; i<3; i++)
    {
        getm(i);
        x = b[i][i];
        for(int a=i; a<5; a++)
        {

            b[i][a] = b[i][a]/x;//化1
        }
        for(int j=i+1; j<4; j++)
        {
            x = b[j][i];
            for(int a=i; a<5; a++)
            {
                b[j][a] = b[j][a]-b[i][a]*x;//化0
            }

        }
    }
}

void Lowertriangle()//化下三角就是把对角线上面变成0
{
    double x;
    for(int i=3; i>0; i--)
    {
        x = b[i][i];
        for(int a=i; a<5; a++)
        {
            b[i][a] = b[i][a]/x;//把对角线上值系数变成1
        }
        for(int j=i-1; j>=0; j--)
        {
            x = b[j][i];
            for(int a=j; a<5; a++)
            {
                b[j][a] = b[j][a]-b[i][a]*x;//在上三角基础上把第i列的第i项之外其他变成0
            }
        }
    }
}


int main()
{
    int n;
    scanf("%d",&n);
    while(n>0)
    {
        for(int i=0; i<4; i++)
        {
            for(int j=0; j<5; j++)
            {
                scanf("%lf",&b[i][j]);
            }
        }
        Uppertriangle();
        Lowertriangle();
printf("%0.9lf %0.9lf %0.9lf %0.9lf \n",
       b[0][4],b[1][4],b[2][4],b[3][4]);
        n--;
    }
    return 0;
}

主要是计算线性方程组的两个步骤：化上三角，解上三角，因为考虑到对角线上的值的绝对值会很小，所以根据题目提议选择选列主元高斯消元法进行接线性方程组。先对其选列主元，我这里采用选择排序算法对其进行选取列主元，这里的话对比课本例程时间复杂度同一个级数但是运算量大，不过也相对结构简单。然后是化上三角，这里我在开始之初，没意识到每一次对一个方程的第一个系数（b[k][0]）运算之后它会发生变化，所以要么最后对其化0、化1操作，或者做一个临时变量来记录初始的第一系数值。避免后面扰乱计算。我主要就是在这里出问题，还有一个小问题就是代码提交系统oj上面提交不过是因为，每次输出结果后需要换行。