Fortran实现QR分解结果异常求助：求无LAPACK的正确代码

以下是基于Householder变换实现QR分解的纯Fortran代码，使用real64高精度类型，无LAPACK依赖：

module qr_mod
    use, intrinsic :: iso_fortran_env, only: real64
    implicit none

contains

    ! 计算Householder变换向量v和系数beta，将向量x变换为平行于单位向量e1
    subroutine householder(x, v, beta)
        real(real64), intent(inout) :: x(:)
        real(real64), intent(out) :: v(:), beta
        real(real64) :: sigma, mu

        sigma = dot_product(x(2:), x(2:))
        if (sigma == 0.0_real64) then
            beta = 0.0_real64
            v = 0.0_real64
            return
        end if

        mu = sqrt(x(1)**2 + sigma)
        ! 选择符号避免数值抵消，提升稳定性
        if (x(1) <= 0.0_real64) then
            x(1) = x(1) - mu
        else
            x(1) = -sigma / (x(1) + mu)
        end if

        v = x / x(1)
        beta = 2.0_real64 * x(1)**2 / (sigma + x(1)**2)
        x(1) = mu  ! 变换后x的第一个元素为mu，其余为0
    end subroutine householder

    ! QR分解：输入矩阵A，输出正交矩阵Q和上三角矩阵R，满足A = Q*R
    subroutine qr_householder(A, Q, R)
        real(real64), intent(inout) :: A(:,:)
        real(real64), intent(out) :: Q(:,:), R(:,:)
        integer :: m, n, k
        real(real64), allocatable :: v(:)
        real(real64) :: beta

        m = size(A, 1)
        n = size(A, 2)
        if (size(Q,1)/=m .or. size(Q,2)/=m) error stop "Q的维度必须与A的行数匹配"
        if (size(R,1)/=m .or. size(R,2)/=n) error stop "R的维度必须与A完全匹配"

        ! 初始化Q为单位矩阵，R为A的副本
        Q = 0.0_real64
        do k = 1, m
            Q(k,k) = 1.0_real64
        end do
        R = A

        allocate(v(m))

        do k = 1, n
            ! 对R的第k列，从第k行开始执行Householder变换
            call householder(R(k:m, k), v(k:m), beta)

            ! 更新R：R = H_k * R，H_k = I - beta*v*v^T
            R(k:m, k+1:n) = R(k:m, k+1:n) - beta * matmul(spread(v(k:m), 2, n - k), &
                                                          spread(v(k:m), 1, n - k)) * R(k:m, k+1:n)

            ! 更新Q：Q = Q * H_k^T（H_k是正交矩阵，H_k^T = H_k）
            Q(:,k:m) = Q(:,k:m) - beta * matmul(Q(:,k:m), spread(v(k:m), 2, m - k + 1)) * transpose(v(k:m))
        end do

        deallocate(v)
    end subroutine qr_householder

end module qr_mod

program test_qr_decomposition
    use qr_mod
    implicit none
    real(real64) :: A(3,3), Q(3,3), R(3,3), recon_A(3,3)
    integer :: i

    ! 测试用3x3矩阵
    A = reshape([1.0_real64, 2.0_real64, 3.0_real64, &
                 4.0_real64, 5.0_real64, 6.0_real64, &
                 7.0_real64, 8.0_real64, 10.0_real64], [3,3])

    print *, "=== 原始矩阵 A ==="
    do i = 1, 3
        print "(3F12.6)", A(i,:)
    end do

    ! 执行QR分解
    call qr_householder(A, Q, R)

    print *, new_line('a')//"=== 正交矩阵 Q ==="
    do i = 1, 3
        print "(3F12.6)", Q(i,:)
    end do

    print *, new_line('a')//"=== 上三角矩阵 R ==="
    do i = 1, 3
        print "(3F12.6)", R(i,:)
    end do

    ! 验证分解正确性：A 应等于 Q*R
    recon_A = matmul(Q, R)
    print *, new_line('a')//"=== 重构矩阵 Q*R ==="
    do i = 1, 3
        print "(3F12.6)", recon_A(i,:)
    end do

    ! 验证Q的正交性：Q^T*Q 应接近单位矩阵
    print *, new_line('a')//"=== Q^T * Q（验证正交性） ==="
    do i = 1, 3
        print "(3F12.6)", matmul(transpose(Q), Q)(i,:)
    end do

end program test_qr_decomposition

关键说明 #

• 使用real64类型确保数值计算精度，避免单精度下的累积误差
• householder子程序通过选择合适的符号处理数值稳定性问题，避免减法抵消
• 分解过程中逐列构造Householder矩阵，同时更新上三角矩阵R和正交矩阵Q
• 测试程序包含分解正确性验证和正交性验证，可直接编译运行

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