使用Fortran LAPACK dgesv求解多右端项方程组结果异常求助

问题背景 #

尝试用LAPACK的dgesv子程序求解线性方程组，系统矩阵B定义如下：

B(1,:) = [1.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000]
B(2,:) = [0.11111111111111113, 0.66666666666666674, 0.22222222222222221, 0.0000000000000000]
B(3,:) = [0.0000000000000000, 0.22222222222222227, 0.66666666666666674, 0.11111111111111106]
B(4,:) = [0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 1.0000000000000000]

多右端项矩阵P定义如下：

P(1,:) = [0.0000000000000000, 1.0000000000000000, 2.0000000000000000, 3.0000000000000000]
P(2,:) = [0.0000000000000000, 0.0000000000000000, 0.0000000000000000, 0.0000000000000000]
P(3,:) = [0.0000000000000000, 1.0000000000000000, 2.0000000000000000, 3.0000000000000000]

异常现象 #

1. 调用多右端项求解：call dgesv(4, 3, B, 4, pivot, P, 4, info)，得到的结果矩阵所有行均不正确。
2. 逐个调用单右端项求解时，仅第一次调用能得到正确结果，后续调用结果均错误（零右端项除外）。

错误原因分析 #

1. 系数矩阵被dgesv修改 #

dgesv子程序会直接覆盖输入的系数矩阵B，将其替换为LU分解的结果。第一次调用后，B已经不再是初始化时的原始矩阵，后续调用使用的是分解后的LU矩阵，自然无法得到正确结果。

2. 右端项矩阵维度/存储顺序错误 #

Fortran采用列主序存储，LAPACK所有子程序均遵循此规则：

• dgesv的第二个参数NRHS是右端项的数量（此处为3）
• 第六个参数（右端项矩阵）要求维度为(N, NRHS)，即每一列对应一个右端项
  你的代码中定义P(3,4)（3行4列），与要求的(4,3)（4行3列）完全相反，导致数据读取错位，多右端项调用结果完全错误。

修正方案 #

方案一：保存原始系数矩阵副本，每次调用前恢复 #

program main

    integer :: pivot(4) = 0
    integer :: info = 0
    doubleprecision :: B(4,4), B_orig(4,4)
    doubleprecision :: P(4,3)  ! 修正维度：4行3列，每列对应一个右端项

    ! 初始化原始系数矩阵并保存副本
    B_orig(1,:) = [1.0d0, 0.0d0, 0.0d0, 0.0d0]
    B_orig(2,:) = [1.0d0/9.0d0, 2.0d0/3.0d0, 2.0d0/9.0d0, 0.0d0]
    B_orig(3,:) = [0.0d0, 2.0d0/9.0d0, 2.0d0/3.0d0, 1.0d0/9.0d0]
    B_orig(4,:) = [0.0d0, 0.0d0, 0.0d0, 1.0d0]

    ! 初始化右端项矩阵（每列对应一个右端项）
    P(:,1) = [0.0d0, 1.0d0, 2.0d0, 3.0d0]   ! 原P(1,:)
    P(:,2) = [0.0d0, 0.0d0, 0.0d0, 0.0d0]   ! 原P(2,:)
    P(:,3) = [0.0d0, 1.0d0, 2.0d0, 3.0d0]   ! 原P(3,:)

    ! 逐个求解单右端项，每次调用前恢复原始矩阵
    B = B_orig
    call dgesv(4, 1, B, 4, pivot, P(:,1), 4, info)
    
    B = B_orig
    call dgesv(4, 1, B, 4, pivot, P(:,2), 4, info)
    
    B = B_orig
    call dgesv(4, 1, B, 4, pivot, P(:,3), 4, info)

end program main

方案二：分解与求解分离（更高效） #

先调用dgetrf做LU分解，再用dgetrs求解所有右端项，仅需分解一次：

program main

    integer :: pivot(4) = 0
    integer :: info = 0
    doubleprecision :: B(4,4)
    doubleprecision :: P(4,3)

    ! 初始化系数矩阵和右端项矩阵（同方案一）
    B(1,:) = [1.0d0, 0.0d0, 0.0d0, 0.0d0]
    B(2,:) = [1.0d0/9.0d0, 2.0d0/3.0d0, 2.0d0/9.0d0, 0.0d0]
    B(3,:) = [0.0d0, 2.0d0/9.0d0, 2.0d0/3.0d0, 1.0d0/9.0d0]
    B(4,:) = [0.0d0, 0.0d0, 0.0d0, 1.0d0]

    P(:,1) = [0.0d0, 1.0d0, 2.0d0, 3.0d0]
    P(:,2) = [0.0d0, 0.0d0, 0.0d0, 0.0d0]
    P(:,3) = [0.0d0, 1.0d0, 2.0d0, 3.0d0]

    ! LU分解（修改B为分解结果）
    call dgetrf(4, 4, B, 4, pivot, info)
    
    if (info == 0) then
        ! 求解所有右端项
        call dgetrs('N', 4, 3, B, 4, pivot, P, 4, info)
    end if

end program main

方案三：修正维度后直接调用多右端项dgesv #

program main

    integer :: pivot(4) = 0
    integer :: info = 0
    doubleprecision :: B(4,4), B_temp(4,4)
    doubleprecision :: P(4,3)

    ! 初始化原始矩阵和右端项
    B(1,:) = [1.0d0, 0.0d0, 0.0d0, 0.0d0]
    B(2,:) = [1.0d0/9.0d0, 2.0d0/3.0d0, 2.0d0/9.0d0, 0.0d0]
    B(3,:) = [0.0d0, 2.0d0/9.0d0, 2.0d0/3.0d0, 1.0d0/9.0d0]
    B(4,:) = [0.0d0, 0.0d0, 0.0d0, 1.0d0]

    P(:,1) = [0.0d0, 1.0d0, 2.0d0, 3.0d0]
    P(:,2) = [0.0d0, 0.0d0, 0.0d0, 0.0d0]
    P(:,3) = [0.0d0, 1.0d0, 2.0d0, 3.0d0]

    ! 使用临时矩阵避免修改原始B
    B_temp = B
    call dgesv(4, 3, B_temp, 4, pivot, P, 4, info)

end program main

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