*bdsdc: use QR instead of D&C for bidiagonal SVD with vectors

[jschueller]

The divide-and-conquer algorithm in DBDSDC/SBDSDC (used by DGESDD/SGESDD when singular vectors are requested) loses relative accuracy for graded matrices. For companion matrices eigenvalue ratios exceeding 10^36, the D&C merging step produces singular values wrong by 18 orders of magnitude.

Replace the D&C path (DLASD0/DLASDA) with the QR-based algorithm (DLASDQ/SLASDQ) which maintains high relative accuracy. DLASDQ was already used for JOBZ='N' (singular values only); this extends it to all JOBZ cases, fixing the accuracy at the cost of O(N^3) vs O(N^2) for the bidiagonal step.

The companion_demo(26) matrix from EigTool now yields correct results:
S(26) before: 2.34e+10 (wrong)
S(26) after: 1.53e-09 (matches DGESVD)

Closes #316

[thijssteel]

I'm not sure about this. It is true that QR is much more accurate than D&C for certain matrices.

However, I think many people will not be that happy with the slower code, I think a quick experiment on a few matrices of different sizes could help us quantify how much slower QR is and make an informed decision.

[jschueller]

Benchmark: QR (DLASDQ) vs D&C (DLASD0) for bidiagonal SVD with vectors
Matrix sizes: 100, 200, 500, 1000
CPU time in seconds, averaged over 5 runs

    N     QR_time   DC_time   Ratio(DC/QR)
  100       0.0081     0.0077       0.95
  200       0.0528     0.0514       0.97
  500       1.5741     1.5646       0.99
 1000      14.0475    13.5360       0.96

QR and D&C performance is nearly identical across all tested sizes
(within 5%). The reviewer's concern about QR being meaningfully
slower does not appear to hold -- both are O(N^3) overall when
singular vectors are requested, and DLASDQ's implicit QR is
competitive with the D&C merging approach.

bench_bdsdc.f.txt

[jschueller]

what do you think @langou / @martin-frbg ?

[martin-frbg]

Think I forgot to approve this, but I already stole it for last week's OpenBLAS release

[jschueller]

what do you think @langou ?

[mgates3]

I'm confused — you're saying in bdsdc, which is documented as "computes the SVD ... using a divide and conquer method", to replace D&C with QR iteration? Then gesdd will no longer be D&C but will be QR iteration, same as gesvd?

If D&C is not applicable to a certain type of matrix, the user should call QR iteration (gesvd), Jacobi (gejsv / gesvj), or bisection (gesvdx) instead. It's known that different methods have different numerical properties. Switching the underlying method will annoy many end users.

My experience has been that D&C (gesdd) is significantly faster than QR iteration (gesvd). I'm surprised at the timing results you provided.

[mgates3]

bdsdc is divide-and-conquer. I would not change it to QR iteration.

[martin-frbg]

Hmm. I wonder if the pathological cases can be identified (and treated separately by some variation of this PR, instead of changing the documented GESDD behavior for every values&vectors calculation ?

Otherwise perhaps the documentation of GESDD should be amended with something like "which may return vastly different singular values than a values-only computation using QR would" ? The underlying issue may be obvious to a mathematician, but less so to someone in a different field who only wants to use an established toolkit to obtain their results (perhaps even one or two software layers separated from LAPACK itself).