BLAS Matrix multiplication NaN

Question

I am performing an SVD of a matrix using the LAPACK library and then multiplying the matrices to double check that they are correct. See the code below

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

The output I get is

 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

So I am not sure where is this NaN coming from. The odd thing is that if before printing D in such way I print it as follows:

  print*,'D'
  print*, d

Then I don't get the NaN anymore, so the output for D is

 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

Any idea why is this happening?

PS: Information for the dgesvd (LAPACK) and dgemm (BLAS) subroutines.

Which error checking compiler flags do you use? How do you compile it at all? — Vladimir F
This is a subroutine within a big program, but basically looks like this: gfortran -cpp -Ofast -funroll-loops --param max-unroll-times=4 -o lotus lotus.o -L/home/bfg/Workspace/solver/oop/. -L/home/bfg/Workspace/solver/geom_lib/. -L/usr/lib/lapack/. -L/usr/lib/libblas/. -lfluid -lgeom -llapack -lblas — b-fg
As an update: I have compiled the subroutine in another test program that only contains this subroutine and now the NaN has disappeared. I used the same subroutine, same compiler and flags. It's very weird! — b-fg
Try to compile with -fbounds-check (checks for out of bounds array calls) or -Wuninitialized (checks for variables that are not initialized with a value before they are used). If that doesn't throw errors/warnings, check very carefully if the subroutines have the right arguments and if you can use them as you think you can (for example, twice n in dgemm). — StefanS
Checking the previous flags and -fcheck=all I don't get any error. The thing is: why if I print D before printing it component by component the result changes? — b-fg

Jonatan Öström Jonatan Öström · Accepted Answer · 2016-07-10T16:51:23

So from our comments dialogue it seems like you have a problem that stems from not initializing an array. It is always good practice to do this, and in situations where you do operations like var = var +1 it is required. If you are unlucky your program will work fine anyway. But then strange things happen once in a while when some garbage happened to reside in memory where the array gets allocated.

A double should be initialized like this

array = 0.0d0 ! for double precision

or

array = 0  ! ok for single,double and integer

Initialize single precision but not double precision like this:

array = 0.0 ! single (not double) precision.

or this

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

I recommend the page fortran90.org.

BLAS Matrix multiplication NaN

1 Answers