Out of bounds write in dgeev(x) when matrix has NaNs

[andreasnoackjensen]

The Julia collaborator Yichao Yu has debugged the following bug. When a matrix has Infs or NaNs, it is possible that dgeev(x) writes out of bounds. The reason is that e.g. the lines 449 and 479 assume that the array wi has been correctly set by dhseqr such that wi[i]>0 implies wi[i+1]<0 and therefore that i < n. However, with Infs and NaN, dgebal returns with info != 0 and in consequence dhseqr returns early without modifying wi. Maybe dgeev should return early when info != 0 or maybe dhseqr should set wi(:) = 0 before returning early. With NaN input, the program below shows the bug on Mac with reference LAPACK provided by OpenBLAS and on a Ubuntu server with a recent MKL

Code: Select all

program test

    double precision :: a(2,2), wr(2), wi(2), vl(0,0), vr(2,3), work(68)
    integer*8 :: info, lwork

    read(*,*) a(1,1)
    a(2,1) = a(1,1)
    a(1,2) = a(1,1)
    a(2,2) = a(1,1)
    wi(1) = 0.0
    wi(2) = 0.0
    lwork = 68

    write(*,*) 'vr:', vr

    call dgeev('N', 'V', 2_8, a, 2_8, wr, wi, vl, 1_8, vr, 2_8, work, lwork, info)

    write(*,*) 'vr:', vr
    wi(2) = 0.1

    call dgeev('N', 'V', 2_8, a, 2_8, wr, wi, vl, 1_8, vr, 2_8, work, lwork, info)

    write(*,*) 'vr:', vr

end program


In addition: Intel has a warning in their documentation, https://software.intel.com/en-us/node/520866, stating that

LAPACK routines assume that input matrices do not contain IEEE 754 special values such as INF or NaN values. Using these special values may cause LAPACK to return unexpected results or become unstable.

Is that really true that I should expect segfaults from LAPACK when having NaNs and Infs or is Intel wrong here?