DGELSD fails to converge for some matrices

[MaartenBaert]

DGELSD fails to converge in some rare cases. I have attached a set of matrices (as binary files) that reproduce the problem. The LHS matrix is 1602x1602, the RHS matrix is 1602x1. I am only able to reliably reproduce the problem with OMP_NUM_THREADS=1.

matrices.zip

Here's a minimal Fortran program that reproduces the problem:

program solve_least_squares
    implicit none
    integer, parameter :: m = 1602, n = 1602, nrhs = 1
    double precision, dimension(m, n) :: lhs
    double precision, dimension(m, nrhs) :: rhs
    double precision, dimension(min(m, n)) :: s
    double precision :: rcond
    integer :: rank, info, lwork, liwork
    double precision, allocatable :: work(:)
    integer, allocatable :: iwork(:)

    ! Declare the LAPACK routine
    interface
        subroutine DGELSD(m, n, nrhs, a, lda, b, ldb, s, rcond, rank, work, lwork, iwork, info)
            integer, intent(in) :: m, n, nrhs, lda, ldb, lwork
            double precision, intent(inout) :: a(lda, *), b(ldb, *), s(*), work(*)
            double precision, intent(in) :: rcond
            integer, intent(out) :: rank, iwork(*), info
        end subroutine DGELSD
    end interface

    ! Use machine precision as the default
    rcond = -1.0d0

    ! Open and read lhs.bin
    print *, "Reading lhs.bin"
    open(unit=10, file="lhs.bin", form="unformatted", access="stream", status="old")
    read(10) lhs
    close(10)

    ! Open and read rhs.bin
    print *, "Reading rhs.bin"
    open(unit=11, file="rhs.bin", form="unformatted", access="stream", status="old")
    read(11) rhs
    close(11)

    ! Workspace query
    print *, "Querying workspace"
    lwork = -1
    allocate(work(1))
    allocate(iwork(1))

    ! Query optimal workspace size
    print *, "Calling DGELSD to get workspace size"
    call DGELSD(m, n, nrhs, lhs, m, rhs, m, s, rcond, rank, work, lwork, iwork, info)
    if (info /= 0) then
        print *, "Error during workspace query: DGELSD failed with info =", info
        stop
    end if

    ! Allocate optimal workspace
    lwork = int(work(1))
    liwork = iwork(1)
    deallocate(work, iwork)
    allocate(work(lwork), iwork(liwork))

    ! Solve the least-squares problem
    print *, "Calling DGELSD to solve the least-squares problem"
    call DGELSD(m, n, nrhs, lhs, m, rhs, m, s, rcond, rank, work, lwork, iwork, info)
    if (info /= 0) then
        print *, "Error: DGELSD failed with info =", info
        stop
    end if

    ! Write the result to res.bin
    print *, "Writing res.bin"
    open(unit=12, file="res.bin", form="unformatted", access="stream", status="replace")
    write(12) rhs(1:n, 1)  ! Write only the solution vector
    close(12)

    ! Deallocate workspace
    deallocate(work, iwork)

    print *, "Solution written to res.bin"
end program solve_least_squares

When the problem occurs, it produces the following error:

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

... which indicates that the SVD failed to converge.

I can reproduce the problem with OpenBLAS versions v0.3.18 through v0.3.30. I bisected it down to this commit:

47ba85f314808476c8254779389607f9af60231f is the first bad commit
commit 47ba85f314808476c8254779389607f9af60231f
Author: Martin Kroeker <martin@ruby.chemie.uni-freiburg.de>
Date:   Thu Jul 22 17:24:15 2021 +0200

    Fix regex to match kernels suffixed with cpuname too

 cmake/utils.cmake | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

... which had the effect of enabling FMA instructions for many files.

This problem seems to be highly sensitive to floating point rounding errors. Even tiny changes to the input matrix make the problem go away, as does setting OMP_NUM_THREADS to anything other than 1. I don't know whether the introduction of FMA instructions actually caused this problem, or merely changed the rounding errors slightly such that my example matrix no longer works to trigger the problem.

I have compiled OpenBLAS with make TARGET=HASWELL on Red Hat Enterprise Linux release 8.10 (Ootpa). The pre-compiled OpenBLAS library included with Scipy shares the same issue. The corresponding Scipy issue can be found here: scipy/scipy#23220

[MaartenBaert]

On closer inspection, it looks like the first commit that caused the problem was actually this one:

49bbf330ca592f439a07f24f137e61af1cc9c616 is the first bad commit
commit 49bbf330ca592f439a07f24f137e61af1cc9c616
Author: Martin Kroeker <martin@ruby.chemie.uni-freiburg.de>
Date:   Sun Jul 18 22:19:19 2021 +0200

    Empirical workaround for numpy SVD NaN problem from issue 3318

 kernel/Makefile.L2 | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

So perhaps this problem is related to issue #3318?

It looks like the only file affected by this change is kernel/x86_64/dgemv_t_4.c. When this file is compiled with -mfma, the SVD fails to converge. When compiled without -mfma, it works. Unfortunately this is the exact opposite of issue #3318.

[martin-frbg]

Just to clarify - are you building with CMake or plain make ? The file from your first bisect is only needed in one, and the file from your second only in the other...

[MaartenBaert]

I am using plain make. I have also tried manually recompiling just the one file kernel/x86_64/dgemv_t_4.c with and without -mfma and confirmed that this is enough to make the problem appear or disappear.

I may have had issues with incomplete rebuilds the first time I ran the bisect. The second time I completely cleaned the repository each time, so it should be accurate.

[MaartenBaert]

Apparently the FMA flag is a red herring - I managed to reproduce the problem with or without FMA enabled by writing a bidiagonalized matrix generated by an affected version to a file, and then using that matrix as the input of either DGELSD or DLALSD in other (previously unaffected) versions of OpenBLAS, which reproduces the problem more reliably. It also makes reproducing the problem independent of OMP_NUM_THREADS.

Here's the bidiagonal matrix that triggers the problem (D = main diagonal, E = super-diagonal):
bidiag_de.zip

[martin-frbg]

Agree that it has nothing to do with the FMA flag, the original "effect" is not reproducible at -O1 -mfma instead of the default -O2 (or conventional -O3 for CMake "Release" builds). I'm currently trying to find the exact GCC option switched on at -O2 that triggers it, but it might be equally enlightening to see if the problem occurs with the original Reference-LAPACK (and its BLAS) as well - this could mean that either the matrix is ill conditioned, or the DGELSD algorithm ill equipped to handle it.

[angsch]

Blind shot: In reference-LAPACK, the algorithm swap of NRM2 changed the convergence behavior of some matrices (in both directions,some previously converged went to diverged and some previously diverged went to converged). If the matrix is sensitive to small roundoff errors, you can try swapping out NRM2.

Also, I second Martin, it would be good to know if the problem reproduces with reference-LAPACK + reference-BLAS. If so, is the 1602x1602 problem the smallest size (the title mentions "some matrices")?

[MaartenBaert]

It turns out that with the bidiagonal form of the matrix I can reproduce the problem with the latest reference LAPACK+BLAS as well.

I have not yet found a smaller test case that triggers the problem, only larger ones. Also, I found that only the LHS is relevant, the RHS can be zero and it will still trigger the problem.

I haven't been able to reproduce the problem with MKL so far.

[martin-frbg]

It appears that the decisive difference in compiler options between the default -O2 and a build at -O1 that passes the test is -fexpensive-optimizations , but I have not checked yet if this matters specifically for compilation of GEMV or one of the other functions in the call graph (like NRM2 as @angsch mentioned) - I think we've seen this effect before, somewhere in Reference-LAPACK code.

[martin-frbg]

hmm.. the -fexpensive-optimizations-related problem in Reference-LAPACK was with NRM2 on Apple M1 (issue 1073 there), despite BLAS being Fortran code there as well, and also using a revised algorithm for NRM2...

[MaartenBaert]

@martin-frbg Did you test this with the original example, or with the bidiagonal matrix I shared later? For reference, here's the updated example that uses the bidiagonal matrix, which triggers the problem reliably for me:
example.zip
This example fails for me with reference LAPACK with or without optimization.

[martin-frbg]

That was still with the original testcase. The bidiagonal case seems to fail no matter what variations in compiler options or NRM2 code I throw at it. (If I replace the Haswell DNRM2 kernel with a simple C version, the testcase produces an additional IEEE_DIVIDE_BY_ZERO exception that is raised by "ETA = B / A" in line 350 of LAPACK's DLAED6.f , but that may be an unrelated oddity. (Though it could perhaps explain why DGELSD seems to be vulnerable to DNRM2's precision - the next line after this division checks if A is less or equal zero...)

[martin-frbg]

(fixing that division in DLAED6 doesn't help, the IEEE_DENORMAL message doesn't provide a lead either - just DLAMCH computing the smallest "safely" usable number - and the IEEE_UNDERFLOW is scaling a 1.e-298 matrix element by 1.e-14 in DLASR line 273 )

[MaartenBaert]

I have narrowed down the problem to a convergence failure in DLASD4, which is a root finding algorithm for calculating updated singular values. It is failing for a subproblem with N=233 on the 4th singular value. The reason for the failure seems to be that the caller (DLASD6) has violated the precondition that the original singular values D must be in sorted order. The first 10 values of D are:

0.0
0.9999999970963765
0.9999999970963985
0.99999999709642
0.9999999970964196
0.999999997096431
0.9999999970964446
0.9999999970964584
0.9999999970964952
0.9999999970965138

Note that the 5th entry is smaller than the 4th, which is not allowed. These values come directly from DLASD7, which is responsible for merging and sorting singular values. So I suspect that's the real source of the problem. The documentation states that the output of DLASD7 will be in sorted order, but in this case it isn't.

[martin-frbg]

That then would seem to be a failure of DLAMRG (which is supposed to generate the index array to turn D into a sorted list, if I understand the comments correctly ? At least this function does not make use of any BLAS calls...