使用Lapack和Fortran的Householder QR因子分解

Question

我的fortran代码。我需要创建Q和R矩阵，这样A=Q*R和则计算err=||A-Q*R||/||A||，其中||||是Frobenius范数(与二维数组的欧几里得范数相同(。我在这里读过关于dgeqrf的文章https://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga3766ea903391b5cf9008132f7440ec7b.html.错误肯定是在我使用该函数的某个地方，因为我有一个类似的代码，带有我自己的QR factrorization，一切都在那里工作。请帮帮我。我认为程序中有一个错误，因为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

这是我昨晚打字不多时想给出的答案。很高兴你已经解决了眼前的问题(需要使用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

Answer 2

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

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