Description

I have found that DGELSD sometimes fails for certain matrices, even though the matrix is well-conditioned and does not contain any unreasonably large or small values. I have attached a 1602x1062 bidiagonal matrix that triggers the problem as an example (this is the smallest test case I have been able to find so far).

Minimal example

This includes a minimal Fortran example that triggers the problem, and the required data files with matrix data.
example.zip

When the problem occurs, this example will produce the following error:

Error: DGELSD failed with info =           1
Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL

Tested with the latest reference blas+lapack as well as various older versions, on Linux (RHEL8) with gcc/gfortran 8.5.0. Originally found in OpenBLAS: OpenMathLib/OpenBLAS#5333

Checklist

• I've included a minimal example to reproduce the issue
• I'd be willing to make a PR to solve this issue