Fix documentation error and ordering bug in DLASD7

[MaartenBaert]

Description

For certain problems with many singular values that are close together, DLASD7 subtly mis-orders the resulting singular values, which causes cascading problems throughout the divide-and-conquer algorithm and eventually results in convergence failure.

There is also a documentation error: the documentation states that the deflated singular values are returned in D in increasing order, but in reality they are in decreasing order.

This fixes #1138.
Please refer to OpenMathLib/OpenBLAS#5333 for a detailed discussion of the issue and fix.

Checklist

• The documentation has been updated.
• If the PR solves a specific issue, it is set to be closed on merge.

[martin-frbg]

It is trivial to predict that this fix (if accepted) will also be needed in slasd7.f - probably slasd2.f/dlasd2.f as well. And there is broadly similar code in ?LAED2 and ?LAED8...

[langou]

(1) The change (increasing → decreasing) is great. Thanks @MaartenBaert

(2) As mentioned by @martin-frbg, it would better to report the same change to slasd7.f, slasd2.f and dlasd2.f. The same mistake is in these three files.

(3) Something to consider, I think "nonincreasing" would be better than "decreasing" in this context since we have D( I ) ≥ D(I+1), we do not necessarily have D( I ) > D(I+1). So I would suggest to change (increasing → nonincreasing).

[langou]

I looked at [s/d]laed2.f and and [s/d]laed8.f and these four routines use "increasing". They do not use the term "nondecreasing". This is to say that the proposed change (increasing → decreasing) is more consistent with is done in [s/d]laed2.f and and [s/d]laed8.f. So please discard my comment #3 above that starts with "(3) Something to consider, I think . . ."

lapack/SRC/dlaed2.f

Line 76 in efbd2fd

*> (those which were deflated) sorted into increasing order.

lapack/SRC/dlaed8.f

Line 86 in efbd2fd

*> (those which were deflated) sorted into increasing order.

[langou]

Most drivers routines use the terminology "ascending order" for HEEV (resp. "descending order" for SVD). Most drivers routines do not use the terminology "increasing order" for HEEV (resp. "decreasing order" for SVD). (This might be more due to comment-copy-paste than intentionality comment writing.)

Some examples for the HEEV drivers:

lapack/SRC/dsyev.f

Line 86 in efbd2fd

*> If INFO = 0, the eigenvalues in ascending order.

lapack/SRC/dsyevd.f

Line 92 in efbd2fd

*> If INFO = 0, the eigenvalues in ascending order.

lapack/SRC/dsyevr.f

Line 240 in efbd2fd

*> ascending order.

lapack/SRC/dsyevx.f

Line 167 in efbd2fd

*> eigenvalues in ascending order.

Some examples for the SVD drivers:

lapack/SRC/dgesvd.f

Line 47 in efbd2fd

*> are returned in descending order. The first min(m,n) columns of

lapack/SRC/dgesvdx.f

Line 51 in efbd2fd

*> are returned in descending order. The first min(m,n) columns of

GESVDQ avoids to use any term and

lapack/SRC/dgesvdq.f

Line 156 in efbd2fd

*> The singular values of A, ordered so that S(i) >= S(i+1).

lapack/SRC/dgesvd.f

Line 121 in efbd2fd

*> The singular values of A, sorted so that S(i) >= S(i+1).

[MaartenBaert]

Thanks for the feedback, I've pushed equivalent changes to SLASD7, SLASD2 and DLASD2.

Interestingly, SLAED8/DLAED8 appears to contain code specifically to reorder the indices, presumably in order to solve the same problem I encountered, though it does this in a far more convoluted way:

            IF( K2+I.LE.N ) THEN
               IF( D( JLAM ).LT.D( INDXP( K2+I ) ) ) THEN
                  INDXP( K2+I-1 ) = INDXP( K2+I )
                  INDXP( K2+I ) = JLAM
                  I = I + 1
                  GO TO 90
               ELSE
                  INDXP( K2+I-1 ) = JLAM
               END IF
            ELSE
               INDXP( K2+I-1 ) = JLAM
            END IF

This is basically doing an insertion sort step to find the correct insertion position into the output array to avoid mis-ordering D. I think this is functionally equivalent to what I would have done:

            DO 85 JP = JLAM, J - 1
               INDXP( K2 + J - 1 - JP ) = JP
   85       CONTINUE

The documentation of SLAED8/DLAED8 is still wrong though (increasing -> decreasing) so I will change that.

[MaartenBaert]

Regarding SLAED2/DLAED2/SLAED8/DLAED8, do you prefer to keep the existing (complicated) reordering code, or replace it with the new but simpler code?

[martin-frbg]

Perhaps simplify the ..ED codes in a separate PR, if desired ? (Their rewrite is only tangentially related to the bug this one fixes)