使用 Lapack 和 Fortran 进行家庭 QR 分解

Question

我在使用 dgeqrf 进行 QR 分解时遇到问题代码。我需要创建 Q 和 R 矩阵，以便 A = Q * R 和然后计算 err=||A - Q*R|| / ||A||，其中|| ||是 Frobenius 范数（与二维数组的欧几里德范数相同）。我在此处 https://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga3766ea903391b5cf9008132f7440ec7b.html 阅读了有关 dgeqrf 的信息。错误肯定出在我对函数的使用中，因为我有一个类似的代码和我自己的 QR 分解，一切正常。请帮助me.I认为程序中有错误，因为err应该在10^(-15)左右（这里是10^(0)左右），那是因为A应该等于Q * R 有小误差。

program terra
    implicit none
    integer(4) :: n, diag = 1, h, wdth, u
    real(8), dimension(:, :), allocatable :: A, Q, R, M_aux
    real(8), dimension(:), allocatable :: arr, aux
    real(8) :: err
    n = 1024
    allocate (A(n,n), Q(n,n), R(n, n), M_aux(n, n), arr(n), aux(n))

    !filling matrices
    do h = 1, n
        do wdth = 1, n
            A(h, wdth) = 1 / real(h + wdth, 8)
            R(h,wdth) = A(h,wdth)
            Q(h, wdth) = 0
        end do
        Q(diag, diag) = 1
        diag = diag + 1
    end do

    !doing the factrorisation
    call dgeqrf(n, n, R, n , arr, aux, n, u)

    !getting Q matrix
    do h = 2, n
        diag = h - 1
        do wdth = 1, diag
            Q(h, wdth) = R(h, wdth)
            R(h, wdth) = 0
            !write(*, *) h, wdth, R(h,wdth), Q(h, wdth)
        end do
    end do

    !Transposing it
    do wdth = 1, n
        do h = 1, n
            M_aux(wdth, h) = Q(wdth, h)
        end do
    end do
    do wdth = 1, n
        do h = 1, n
            Q(wdth, h) = M_aux(h, wdth)
        end do
    end do

    !computing err
    do h = 1, n
        do wdth = 1, n
            M_aux(h, wdth) = 0
            do diag = 1, n
                M_aux(h, wdth) = M_aux(h, wdth) + Q(h, diag) * R(diag, wdth)
            end do
        end do
    end do
    M_aux = M_aux - A
    err = norm2(M_aux) / norm2(A)
    write(*,*) err
end program

Answer 1

program terra
    implicit none
    integer(4) :: n, diag = 1, h, wdth, u
    real(8), dimension(:, :), allocatable :: A, Q, R, M_aux
    real(8), dimension(:), allocatable :: arr, aux
    real(8) :: err
    n = 1024
    allocate (A(n,n), Q(n,n), R(n, n), M_aux(n, n), arr(n), aux(n), work(n))

    !filling matrices
    do h = 1, n
        do wdth = 1, n
            A(h, wdth) = 1 / real(h + wdth, 8)
            R(h,wdth) = A(h,wdth)
            Q(h, wdth) = 0
        end do
    end do

    !doing the factrorization
    call dgeqrf(n, n, R, n , arr, aux, n, u)
    Q = R

    !getting R matrix
    do h = 2, n
        diag = h - 1
        do wdth = 1, diag
            R(h, wdth) = 0
        end do
    end do

    !getting Q matrix
    call DORG2R(n, n, n, Q, n, arr, aux, u)

    !computing err
    do h = 1, n
        do wdth = 1, n
            M_aux(h, wdth) = 0
            do diag = 1, n
                M_aux(h, wdth) = M_aux(h, wdth) + Q(h, diag) * R(diag, wdth)
            end do
        end do
    end do
    M_aux = M_aux - A
    err = norm2(M_aux) / norm2(A)
    write(*,*) err
end program

现在可以了。错误大约是 10^(-16)。伊恩布什，我衷心感谢你！我可以在你的评论上点赞或类似的东西吗？我不明白怎么做。

Answer 2

这是我昨晚不想打字时想给出的答案。很高兴你解决了眼前的问题（需要使用 DORG2R），但我还想指出一些其他问题

Real( 8 ) 和 Integer( 4 ) 不可移植且不适用于所有编译器。下面的代码展示了一种可移植的方式
对于像 Fortran 这样的专栏主要语言，你的循环是错误的。按照下面代码中的顺序，内存访问模式更好，循环会更快
现在没有人应该使用循环来编写矩阵乘法。更简单更好的方法是使用 Fortran 内在函数 Matmul，但如果你真的想要性能，BLAS 可能仍然是你最好的选择。我在下面展示了两者。

不管怎样，这就是我所拥有的

ian@eris:~/work/stack$ cat qr.f90
Program terra

  Use, Intrinsic :: iso_fortran_env, Only :  wp => real64
  
  Implicit None

  Integer :: n, h, wdth, u
  Real(wp), Dimension(:, :), Allocatable :: A, Q, R, M_aux
  Real(wp), Dimension(:), Allocatable :: arr, aux, work
  Real(wp) :: err
  ! Smaller for quick testing
  n = 1024
!!$  n = 200
  Allocate (A(n,n), Q(n,n), R(n, n), M_aux(n, n), arr(n), aux(n), work(n))

  !filling matrices
  ! Simple array syntax to fill Q
  Q = 0.0_wp
  ! Order loops for best memory access order in a
  ! column major language like Fortran
  Do wdth = 1, n
     Do h = 1, n
        A(h, wdth) = 1.0_wp / Real(h + wdth, Kind( A ) )
        R(h,wdth) = A(h,wdth)
     End Do
  End Do

  !doing the factrorization
  Call dgeqrf(n, n, R, n , arr, aux, n, u)
  Q = R

  !getting R matrix
  ! Order loops for best memory access order in a
  ! column major language like Fortran
  Do wdth = 1, n
     Do h = wdth + 1, n
        R(h, wdth) = 0.0_wp
     End Do
  End Do

  !getting Q matrix
  Call DORG2R(n, n, n, Q, n, arr, aux, u)

  ! Find error using Matmul - probably quicker than simple loops
  M_aux = Matmul( Q, R )
  M_aux = M_aux - A
  err = Norm2(M_aux) / Norm2(A)
  Write(*,*) 'Error, matmul ', err

  ! Find error using BLAS - almost certainly quicker than simple loops
  M_aux = Matmul( Q, R )
  Call dgemm( 'N', 'N', n, n, n, 1.0_wp, Q    , Size( Q    , Dim = 1 ), &
                                         R    , Size( R    , Dim = 1 ), &
                                 0.0_wp, M_Aux, Size( M_aux, Dim = 1 ) )
  M_aux = M_aux - A
  err = Norm2(M_aux) / Norm2(A)
  Write(*,*) 'Error, Blas ', err

End Program terra
  
ian@eris:~/work/stack$ gfortran --version
GNU Fortran (Ubuntu 7.4.0-1ubuntu1~18.04.1) 7.4.0
Copyright (C) 2017 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
ian@eris:~/work/stack$ gfortran -Wall -Wextra -fcheck=all -g -std=f2008 -Wall qr.f90  -llapack -lblas
ian@eris:~/work/stack$ ./a.out
 Error, matmul    8.2474910291356388E-016
 Error, Blas    8.9420732364900517E-016

使用 Lapack 和 Fortran 进行家庭 QR 分解

Householder QR-factorization using Lapack and Fortran

fortran

lapack