Workspace queries return too small values for xGESDD

[rmlarsen]

We have found that the workspace query for xGESDD returns values that are too small for certain matrix sizes, causing the subsequent call to xGESDD to fail with INFO=-12. We have verified that this happens for JOBZ='A' and JOBZ='S' but not for JOBZ='N' or JOBZ='O'. Example program:

program sgesddbug
implicit none
integer, parameter :: m = 6163, n = 2967, min_lwork = n*(6+4*n)+m
integer :: lwork, info
real :: A(m,n), U(m,n), VT(n,n), sigma(n), iwork(8*n)
real, allocatable :: work(:)
lwork = -1
allocate(work(1))
call sgesdd('s', m, n, A, m, sigma, U, m, VT, n, work, lwork, iwork, info)
write (*,*) 'lwork = ',work(1), ', minimum lwork = ', min_lwork

lwork = nint(work(1))
deallocate(work)
allocate(work(lwork))

call sgesdd('s', m, n, A, m, sigma, U, m, VT, n, work, lwork, iwork, info)
write (*,*) 'info = ', info
end

prints:

lwork = 35233124.0 , minimum lwork = 35236321
<some BLAS error string>
info = -12


Best regards,
Rasmus Munk Larsen,
Google Knowledge & Research, Mountain View.

[admin]

Thank you for the example code.
We listed the Bug on our Errata page.
We hope to fix the issue very soon.
Julie

[admin]

Dear Rasmus,
We fixed the issue you reported in revision 1503.
Actually, there were two issues in the code you reported.
The second issue is more tricky.
The number corresponding to the workzise 35233125 is not represented in real.
The two REAL around the integer 35,233,125 are 35,233,124 et 35,233,128.
the closest REAL of 35,233,125 est 35,233,124 which in binary is represented by:
( 1.00001100110011101011001 )_2 * 2^25

I would encourage to use double precision. That would solve the issue.