subtle bug in SGESDD lapack ilp64 for large matrices

[Basden__Alastair]

Hi,
just incase anyone else runs into this problem, there is a subtle bug in SGESDD in MKL LAPACK when linking the ilp64 version.

When performing a workspace size query, the size of the work needed is returned as the first element of the 32-bit float array work. However, internally to sgesdd, this is computed as a long integer (64 bit) before being cast to a 32 bit float. This cast results in a loss of precision when the A matrix size is large, and so the work size estimate is wrong, meaning that if you use this value to allocate the work array, subsequent calls to sgesdd fail when using it, as the array is too small. The way round this (which I think should be done internally in sgesdd) is to convert this floating point value to 64bit integer, and increment it until a cast back to 32 bit float yields a different number. This integer value is then the value that should be used for the work size.

Hope this helps someone...

[Michael_C_Intel4]

Hi,

you're right, it's a potential issue. It can be fixed on MKL side. At the moment I would suggest several workarounds:

1) apply the algorithm you described when you get the value exceeding 2^32;

2) don't query the workspace, but use the formula suggested in the documentation explicitly:

For real flavors

:

If

jobz = 'N', lwork = 3*min(m, n) + max (max(m,n), 6*min(m, n));

If

jobz = 'O', lwork = 3*(min(m, n))2 + max (max(m, n), 5*(min(m, n))2 + 4*min(m, n));

If

jobz = 'S' or 'A', lwork = 3*(min(m, n))2 + max (max(m, n), 4*(min(m, n))2 + 4*min(m, n))

3) query DGESDD for the workspace - it should be the same as for sgesdd and contain all digits.

Michael.