基于Gauss-Jordan消元法的Julia矩阵求逆算法报错求助

问题根源 #

你遇到的InexactError是因为输入矩阵A是整数类型（Matrix{Int64}），构造的增广矩阵Inv默认继承了整数类型，但消元过程中除法会产生浮点数（比如1.25），整数数组无法存储浮点数，导致类型不匹配。此外代码还有两处逻辑问题：选主元时用的是未更新的原矩阵A，以及缺少置换矩阵P的同步更新。

修正后的完整代码 #

using LinearAlgebra

function GJinv(A)
    # 转换为浮点类型，避免整数运算的类型冲突
    A_float = float(A)
    n = size(A_float, 1)
    # 构造浮点类型的增广矩阵
    Inv = [A_float I(n)]
    # 初始化置换矩阵
    P = Matrix{Float64}(I(n))
    
    for k in 1:n
        # 选取当前列从k行开始的最大绝对值主元
        col_segment = Inv[k:n, k]
        max_abs_val, row_idx = findmax(abs.(col_segment))
        # 转换为原矩阵的行号
        pivot_row = k + row_idx - 1
        
        # 交换增广矩阵的行
        Inv[[k, pivot_row], :] = Inv[[pivot_row, k], :]
        # 同步交换置换矩阵的行
        P[[k, pivot_row], :] = P[[pivot_row, k], :]
        
        # 主元行归一化
        pivot_val = Inv[k, k]
        Inv[k, :] /= pivot_val
        
        # 消去其他行的当前列元素
        for i in 1:n
            if i != k
                factor = Inv[i, k]
                Inv[i, :] -= factor * Inv[k, :]
            end
        end
    end
    
    # 提取逆矩阵部分
    inv_matrix = Inv[:, n+1:end]
    return inv_matrix, P
end

# 测试用例
A = [2 2 3; 4 5 6; 1 2 4]
inv_A, P = GJinv(A)
println("计算得到的逆矩阵：")
display(inv_A)
println("\n置换矩阵P：")
display(P)

核心修改说明 #

• 类型处理：一开始就把输入矩阵转成浮点类型，确保整个消元过程的运算和存储都支持浮点数，从根源解决InexactError。
• 主元选择逻辑修正：直接从增广矩阵的当前列子段找主元，而不是原矩阵——因为行交换后原矩阵和增广矩阵已经不同步了，用原矩阵选主元会导致错误。同时用findmax一步获取最大值和索引，比原代码的maximum+findfirst更可靠，避免浮点精度匹配问题。
• 置换矩阵同步更新：补上了Matlab代码中P矩阵的逻辑，每次行交换时同步更新P，保证输出和原Matlab行为一致。
• 索引简化：把原代码中复杂的行号计算n-(n-k+1)+r改成更直观的k + row_idx -1，减少出错概率。

内容的提问来源于stack exchange，提问作者I Like Algebra