为何基于LAPACK的简单Fortran矩阵求逆代码未返回预期值？

我尝试使用LAPACK包计算矩阵的逆，调用了sgetrf和sgetri例程，但测试代码时逆矩阵和LU分解结果均不符合预期。测试用的是4×4矩阵，代码基于Fortran 90编写，在Code::blocks中使用最新GNU编译器和LAPACK库，运行环境为Windows 10。

代码如下：

program test
    implicit none

    real, dimension(4,4)   :: R,R_inv

    integer                              :: info, ipiv
    real, dimension(4)                   :: lol

    R = reshape([real :: 1,0,0,0,          &
                         0,1.0/34.0,0,0,   &
                         0,0,1,2,          &
                         0,0,2,1]          &
                         ,shape(R), order = [2,1]                              )

R_inv = R
call sgetrf(4,4,R_inv,4,ipiv,info)

print *,info
print *, "_"
do pp = 1,4
print *, R_inv(pp,:)
end do
print *, "_"

call sgetri(2,R_inv,4,ipiv,lol,4,info)

print *,info
print *, "_"
print *, lol
print *, "_"
do pp = 1,4
print *, R_inv(pp,:)
end do
R_inv = matmul(R_inv,R)
print *, "_"
do pp = 1,4
print *, R_inv(pp,:)
end do
print *, "_"

end program 

观察到的现象：

• 第一个info输出为2，根据sgetrf文档，这表示因子U(2,2)完全奇异，但Matlab和在线LU分解计算器均给出有效结果；
• sgetri例程似乎运行正常；
• 若矩阵为对角矩阵，代码可正常运行；
• LU分解后的矩阵前半部分对角线正常，后半部分异常；
• 在线计算器和Matlab均显示矩阵R的逆存在；
• 改为双精度并调用dgetrf、dgetri后结果无变化。

请问这段代码存在什么问题？如何用Fortran结合LAPACK正确计算矩阵逆？

核心错误点 #

• ipiv变量维度错误：LAPACK的sgetrf例程中，ipiv参数需要是长度等于矩阵阶数（此处为4）的整数数组，用于存储LU分解的行交换枢轴索引。你定义的ipiv是单个integer变量，会导致内存越界，破坏矩阵数据，这是info=2的根本原因。
• sgetri的阶数参数错误：调用sgetri时第一个参数是矩阵的阶数，你传入了2而非4，导致仅计算前2×2子矩阵的逆，完全不符合需求。

修正后的代码 #

program test
    implicit none

    real, dimension(4,4) :: R, R_inv
    integer, dimension(4) :: ipiv  ! 修正：ipiv改为长度4的数组
    integer :: info
    real, dimension(4) :: lol

    ! 初始化矩阵R
    R = reshape([real :: 1,0,0,0, &
                     0,1.0/34.0,0,0, &
                     0,0,1,2, &
                     0,0,2,1], shape(R), order = [2,1])

    R_inv = R
    ! 调用LU分解，修正ipiv参数
    call sgetrf(4, 4, R_inv, 4, ipiv, info)

    print *, "sgetrf info: ", info
    print *, "_"
    do pp = 1,4
        print *, R_inv(pp,:)
    end do
    print *, "_"

    ! 调用矩阵求逆，修正第一个参数为4
    call sgetri(4, R_inv, 4, ipiv, lol, 4, info)

    print *, "sgetri info: ", info
    print *, "_"
    print *, "work array: ", lol
    print *, "_"
    do pp = 1,4
        print *, R_inv(pp,:)
    end do

    ! 验证逆矩阵：R_inv * R 应该接近单位矩阵
    R_inv = matmul(R_inv, R)
    print *, "_"
    print *, "验证结果（单位矩阵）："
    do pp = 1,4
        print *, R_inv(pp,:)
    end do
    print *, "_"

end program 

运行结果说明 #

修正后：

1. sgetrf的info会返回0，表示LU分解成功；
2. sgetri的info也返回0，矩阵求逆成功；
3. 最后的验证结果会输出接近4×4单位矩阵的内容，误差在浮点精度允许范围内，证明逆矩阵计算正确。

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