Fortran的MKL dgemm给出零结果

Question

在以下fortran程序中，我使用Intel的MKL库使用dgemm执行矩阵乘法。最初，我使用了matmul子例程并获得了正确的结果。当我在下面的循环中将matmul转换为dgemm时，我得到了所有零向量而不是正确的输出。感谢您的帮助。

program spectral_norm    
implicit none
!
integer, parameter :: n = 5500, dp = kind(0.0d0)
real(dp), allocatable :: A(:, :), u(:), v(:), Au(:), Av(:)
integer :: i, j
allocate(u(n), v(n), A(n, n), Au(n), Av(n))
do j = 1, n
    do i = 1, n
        A(i, j) = Ac(i, j)
    end do
end do
u = 1
do i = 1, 10
    call dgemm('N','N', n, 1, n, 1.0, A,  n, u, n, 0.0, Au, n) 
    call dgemm('N','N', n, 1, n, 1.0, Au, n, A, n, 0.0, v,  n) 
    call dgemm('N','N', n, 1, n, 1.0, A,  n, v, n, 0.0, Av, n) 
    call dgemm('N','N', n, 1, n, 1.0, Av, n, A, n, 0.0, u,  n) 
    !v = matmul(matmul(A, u), A)
    !u = matmul(matmul(A, v), A)
end do
write(*, "(f0.9)") sqrt(dot_product(u, v) / dot_product(v, v))
contains
pure real(dp) function Ac(i, j) result(r)
integer, intent(in) :: i, j
r = 1._dp / ((i+j-2) * (i+j-1)/2 + i)
end function
end program spectral_norm

这给出了NaN，而matmul的正确输出为1.274224153。

Answer 1

好吧，谢谢大家的建议。我想我找出了错误的根源。在两种情况下，乘法的顺序反转，应该是A * Au和A * Av。这是因为A具有n x n订单，并且Au ANS Av都有订单n x 1。因此，由于尺寸不匹配，我们不能乘以Au * A或Av * A。我在下面发布了更正的版本。

program spectral_norm    
implicit none
!
integer, parameter :: n = 5500, dp = kind(0.d0)
real(dp), allocatable :: A(:,:), u(:), v(:), Au(:), Av(:)
integer :: i, j
allocate(u(n), v(n), A(n, n), Au(n), Av(n))
do j = 1, n
    do i = 1, n
        A(i, j) = Ac(i, j)
    end do
end do
u = 1
do i = 1, 10
    call dgemm('N', 'N', n, 1, n, 1._dp, A, n, u,  n, 0._dp, Au, n)
    call dgemm('T', 'N', n, 1, n, 1._dp, A, n, Au, n, 0._dp, v,  n)
    call dgemm('N', 'N', n, 1, n, 1._dp, A, n, v,  n, 0._dp, Av, n)
    call dgemm('T', 'N', n, 1, n, 1._dp, A, n, Av, n, 0._dp, u,  n)
end do
write(*, "(f0.9)") sqrt(dot_product(u, v) / dot_product(v, v))
contains
pure real(dp) function Ac(i, j) result(r)
    integer, intent(in) :: i, j
    r = 1._dp / ((i+j-2) * (i+j-1)/2 + i)
end function
end program spectral_norm

这给出了正确的结果：

1.274224153
 Elapsed time   0.5150000     seconds