Numerical issues with svd

[olof3]

I've been running into some problems with svd when multiple singular values have the same magnitude. Consider for example the following.

import GenericLinearAlgebra

# Construct a matrix with singular values [1; 1; 1; 0]
U0, _, V0 = svd(big.(reshape(0:15,4,4)))
A = U0[:, 1:3] * V0[:, 1:3]'

U, S, V = svd(A)

norm(A - U*Diagonal(S)*V')

This gives an error of 2.07 (!). With GenericSVD the same error is 5.6e-77. Casting A to Float64 gives an error of 1.9e-16.

Note that the singular values themselves seem to be computed alright.

[andreasnoack]

Sorry for the slow response here. I've identified the issue but I'm still working on the fix. I think the issue is that a shift is causing catastrophic cancellation which introduces an inaccurate Givens rotation. I'll probably have to adjust the convergence criterion a bit.