BLAS 矩阵乘法 NaN

Question

我正在使用 LAPACK 库执行矩阵的 SVD，然后将矩阵相乘以仔细检查它们是否正确。见下方代码

subroutine svd_and_dgemm() ! -- Matrix decomp: A = USV^t
  implicit none

  integer,parameter :: m = 2
  integer,parameter :: n = 3
  integer i,info,lda,ldu,ldv,lwork,l,lds,ldc,ldvt,ldd
  real*8 :: a(m,n),a_copy(m,n),sdiag(min(m,n)),s(m,n),u(m,m),vt(n,n),alpha,beta,c(m,n),d(m,n)
  character jobu, jobv, transu, transs
  real*8, allocatable, dimension ( : ) :: work

  lwork = max(1,3*min(m,n) + max(m,n), 5*min(m,n))
  allocate (work(lwork))

  a = reshape((/3,1,1,-1,3,1/),shape(a),order=(/2, 1/)) !column-wise
  print*,'A'
  print*, a(1,1), a(1,2), a(1,3)
  print*, a(2,1), a(2,2), a(2,3)

  jobu = 'A'
  jobv = 'A'
  lda = m
  ldu = m
  ldv = n
  a_copy = a

  call dgesvd (jobu, jobv, m, n, a_copy, lda, sdiag, u, ldu, vt, ldv, work, lwork, info)
  if ( info /= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'R8MAT_SVD_LAPACK - Failure!'
    write ( *, '(a)' ) '  The SVD could not be calculated.'
    write ( *, '(a)' ) '  LAPACK routine DGESVD returned a nonzero'
    write ( *, '(a,i8)' ) '  value of the error flag, INFO = ', info
    return
  end if
!
!  Make the MxN matrix S from the diagonal values in SDIAG.
  s(1:m,1:n) = 0.0D+00
  do i = 1, min ( m, n )
    s(i,i) = sdiag(i)
  end do

  print*,'U'
  print*, u(1,1), u(1,2)
  print*, u(2,1), u(2,2)
  print*,'S'
  print*, s(1,1), s(1,2), s(1,3)
  print*, s(2,1), s(2,2), s(2,3)
  print*,'Vt'
  print*, vt(1,1), vt(1,2), vt(1,3)
  print*, vt(2,1), vt(2,2), vt(2,3)
  print*, vt(3,1), vt(3,2), vt(3,3)
  deallocate (work)

  ! -- Verify SVD: A = USV^t
  ! -- Compute C = US
  transu = 'N'
  transs = 'N'
  ldu = m;     lds = m;   ldc = m
  alpha = 1.;  beta = 1.
  call dgemm(transu,transs,m,n,m,alpha,u,ldu,s,lds,beta,c,ldc)
  ! -- Compute A = D = CV^t
  l = m ! nrows C
  ldvt = n;   ldd = m
  call dgemm(transu,transs,m,n,n,alpha,c,ldc,vt,ldvt,beta,d,ldd)
  print*,'C'
  print*, c(1,1), c(1,2), c(1,3)
  print*, c(2,1), c(2,2), c(2,3)
  print*,'D'
  print*, d(1,1), d(1,2), d(1,3)
  print*, d(2,1), d(2,2), d(2,3)    
end subroutine svd_and_dgemm

我得到的输出是

 A
   3.0000000000000000        1.0000000000000000        1.0000000000000000     
  -1.0000000000000000        3.0000000000000000        1.0000000000000000     
 U
 -0.70710678118654835      -0.70710678118654657     
 -0.70710678118654668       0.70710678118654846     
 S
   3.4641016151377553        0.0000000000000000        0.0000000000000000     
   0.0000000000000000        3.1622776601683795        0.0000000000000000     
 Vt
 -0.40824829046386402      -0.81649658092772526      -0.40824829046386291     
 -0.89442719099991508       0.44721359549995882        5.2735593669694936E-016
 -0.18257418583505536      -0.36514837167011066       0.91287092917527679     
 C
  -2.4494897427831814       -2.2360679774997867        0.0000000000000000     
  -2.4494897427831757        2.2360679774997929        0.0000000000000000     
 D
   2.9999999999999991        1.0000000000000002       0.99999999999999989     
                       NaN   2.9999999999999991        1.0000000000000000

所以我不确定这个 NaN 是从哪里来的。奇怪的是，如果在以这种方式打印 D 之前，我按如下方式打印它：

  print*,'D'
  print*, d

然后我不再得到 NaN，所以 D 的输出是

 D
   2.9999999999999991      -0.99999999999999933        1.0000000000000002        2.9999999999999991       0.99999999999999989        1.0000000000000000     
 D
   2.9999999999999991        1.0000000000000002       0.99999999999999989     
 -0.99999999999999933        2.9999999999999991        1.0000000000000000

知道为什么会这样吗？

PS：dgesvd (LAPACK) and dgemm (BLAS) 子例程的信息。

Answer 1

因此，从我们的评论对话来看，您的问题似乎源于未初始化数组。这样做总是好的做法，并且在您执行 var = var +1 之类的操作的情况下，这是必需的。如果你不幸你的程序无论如何都可以正常工作。但是当一些垃圾恰好驻留在分配数组的内存中时，偶尔会发生奇怪的事情。

double 应该这样初始化

array = 0.0d0 ! for double precision

或

array = 0  ! ok for single,double and integer

初始化单精度但不是双精度像这样：

array = 0.0 ! single (not double) precision.

或这个

array = 0.0e0 ! single (not double) precision.

我推荐页面 fortran90.org。

BLAS 矩阵乘法 NaN

BLAS Matrix multiplication NaN

fortran

blas

svd

lapack

matrix-multiplication