C语言高斯约旦法求矩阵逆时循环未迭代问题排查

我刚接触编程，尝试用高斯约旦法编写C语言程序求矩阵逆——通过行变换将给定矩阵转为单位矩阵，同时对初始单位矩阵执行相同变换。但程序未输出正确的逆矩阵，调试时发现k变量对应的循环仅执行k=0的情况，测试输入为3阶矩阵：1 2 3 4 5 6 7 8 8，无法定位问题。

原代码 #

#include <stdio.h>
void main()
{
    float a[5][5],b[5][5]={};
    int n,i,j,k,m,o,p,q;
    printf("Enter order n of square matrix : ");
    scanf("%d",&n);
    printf("Enter %d elements : ",n*n);
    for(i=0;i<n;i++)
        for(j=0;j<n;j++)
            scanf("%f",&a[i][j]);
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
            printf("%10f",a[i][j]);
        printf("\n");
    }
    if(n>=3)
    {
        for(i=0;i<n;i++)
            b[i][i]=1;
        for(i=0;i<n-1;i++)
        {
                if(a[i][i]==0)
                    for(m=i+1;m<n;m++)
                        if(a[m][i]!=0)
                        {
                            for(o=0;o<n;o++)
                            {
                                a[m][o]+=a[m][o]+a[i][o];
                                a[i][o]=a[m][o]-a[i][o];
                                a[m][o]=a[m][o]-a[i][o];
                                b[m][o]+=b[m][o]+b[i][o];
                                b[i][o]=b[m][o]-b[i][o];
                                b[m][o]=b[m][o]-b[i][o];
                            }
                            break;
                        }
            for(j=i+1;j<n;j++)
                {if(a[j][i]!=0)
                {
                    for(k=0;k<n;k++)
                    {
                        b[j][k]=b[j][k]*a[i][i]-a[j][i]*b[i][k];
                        a[j][k]=a[j][k]*a[i][i]-a[j][i]*a[i][k];
                    }
                }
                for(p=0;p<n;p++)
                {
                    for(q=0;q<n;q++)
                        printf("%10f",a[p][q]);
                    printf("\t\t");
                    for(q=0;q<n;q++)
                        printf("%10f",b[p][q]);
                    printf("\n");
                }}
        }
            for(i=n-1;i>0;i--)
                for(j=i-1;j>=0;j--)
                    if(a[j][i]!=0)
                    {
                        for(k=0;k<n;k++)
                        {
                            b[j][k]=b[j][k]*a[i][i]-a[j][i]*b[i][k];
                            a[j][k]=a[j][k]*a[i][i]-a[j][i]*a[i][k];
                        }
                    }
        for(i=0;i<n;i++)
            for(j=0;j<n;j++)
                b[i][j]=(b[i][j])/(a[i][i]);
        printf("\n\n\n");
        for(i=0;i<n;i++)
        {
        for(j=0;j<n;j++)
        printf("%10f",b[i][j]);
        printf("\n");
        }
    }
}

问题定位与修复步骤 #

1. 行交换逻辑完全错误 #

原代码中交换两行的代码逻辑混乱，导致行交换失败：

a[m][o]+=a[m][o]+a[i][o];
a[i][o]=a[m][o]-a[i][o];
a[m][o]=a[m][o]-a[i][o];

正确的行交换应该使用临时变量，同时对矩阵a和b执行相同操作：

float temp;
temp = a[m][o];
a[m][o] = a[i][o];
a[i][o] = temp;

temp = b[m][o];
b[m][o] = b[i][o];
b[i][o] = temp;

2. 高斯消元公式错误 #

原代码使用a[j][k] = a[j][k]*a[i][i] - a[j][i]*a[i][k]进行消元，这会错误放大主元行，导致后续计算全部偏离。高斯约旦法的正确消元逻辑是：
第j行 = 第j行 - (a[j][i]/a[i][i]) * 第i行
需要先计算缩放因子，再执行行变换：

float factor = a[j][i] / a[i][i];
for(k=0;k<n;k++)
{
    a[j][k] -= factor * a[i][k];
    b[j][k] -= factor * b[i][k];
}

3. 循环结构与缩进错误 #

原代码将打印矩阵的循环嵌套在j循环内，导致每次消元一行就打印一次，造成逻辑混乱。同时n>=3的限制完全没必要，高斯约旦法适用于所有可逆方阵，应去掉该判断。

4. 规范问题 #

void main()不符合C语言标准，应改为int main()并添加return 0;。

修复后的完整代码 #

#include <stdio.h>

int main()
{
    float a[5][5], b[5][5] = {0};
    int n, i, j, k, m;
    printf("Enter order n of square matrix : ");
    scanf("%d", &n);
    printf("Enter %d elements : ", n * n);
    for (i = 0; i < n; i++)
        for (j = 0; j < n; j++)
            scanf("%f", &a[i][j]);

    // 初始化单位矩阵b
    for (i = 0; i < n; i++)
        b[i][i] = 1.0f;

    // 高斯约旦消元：将原矩阵转为单位矩阵，同步变换b矩阵
    for (i = 0; i < n; i++)
    {
        // 处理主元为0的情况，交换行
        if (a[i][i] == 0)
        {
            for (m = i + 1; m < n; m++)
            {
                if (a[m][i] != 0)
                {
                    float temp;
                    for (k = 0; k < n; k++)
                    {
                        temp = a[i][k];
                        a[i][k] = a[m][k];
                        a[m][k] = temp;

                        temp = b[i][k];
                        b[i][k] = b[m][k];
                        b[m][k] = temp;
                    }
                    break;
                }
            }
        }

        // 将主元行缩放，使主元为1
        float pivot = a[i][i];
        for (k = 0; k < n; k++)
        {
            a[i][k] /= pivot;
            b[i][k] /= pivot;
        }

        // 消去其他行的第i列元素
        for (j = 0; j < n; j++)
        {
            if (j != i && a[j][i] != 0)
            {
                float factor = a[j][i];
                for (k = 0; k < n; k++)
                {
                    a[j][k] -= factor * a[i][k];
                    b[j][k] -= factor * b[i][k];
                }
            }
        }
    }

    // 输出逆矩阵
    printf("\nInverse matrix:\n");
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
            printf("%10.4f", b[i][j]);
        printf("\n");
    }

    return 0;
}

内容的提问来源于stack exchange，提问作者Lakshita