`svd` gives wrong results (upstream LAPACK bug)

[carstenbauer]

See this discourse thread.

There is an upstream bug in LAPACK's gesdd when run in full mode ('A'). Since our svd, only call LinearAlgebra.LAPACK.gesdd! under the hood, it currently gives wrong results:

using LinearAlgebra

function compan(c::AbstractVector{T}) where T
    N = length(c)
    A = diagm(-1 => ones(T, N-1))
    @inbounds for i in 1:N
        A[1,i] = - c[i]/c[1]
    end
    A
end

function svdvaldiff(N)
  A = compan(reverse( vcat(1, 1 ./ cumprod(1:N, dims=1)) ));
  svals_DnC = svd(A).S
  svals_QR  = svdvals(A)

  return maximum(svals_DnC - svals_QR)
end

function DnC_QR_comparison(Nmax=40)
    for N in 1:Nmax
        println("N = $(N), max diff svdvals: ", svdvaldiff(N)) 
    end
end

DnC_QR_comparison()

which produces

julia> DnC_QR_comparison()
N = 1, max diff svdvals: 0.0
N = 2, max diff svdvals: 0.0
N = 3, max diff svdvals: 0.0
N = 4, max diff svdvals: 0.0
N = 5, max diff svdvals: 0.0
N = 6, max diff svdvals: 0.0
N = 7, max diff svdvals: 0.0
N = 8, max diff svdvals: 0.0
N = 9, max diff svdvals: 0.0
N = 10, max diff svdvals: 0.0
N = 11, max diff svdvals: 0.0
N = 12, max diff svdvals: 0.0
N = 13, max diff svdvals: 0.0
N = 14, max diff svdvals: 0.0
N = 15, max diff svdvals: 0.0
N = 16, max diff svdvals: 0.0
N = 17, max diff svdvals: 0.0
N = 18, max diff svdvals: 0.0
N = 19, max diff svdvals: 0.0
N = 20, max diff svdvals: 0.0
N = 21, max diff svdvals: 0.0
N = 22, max diff svdvals: 0.0
N = 23, max diff svdvals: 0.0
N = 24, max diff svdvals: 0.0
N = 25, max diff svdvals: 749.7516745517077
N = 26, max diff svdvals: 166.767736950586
N = 27, max diff svdvals: 584.3097624800283
N = 28, max diff svdvals: 635.9955873150014
N = 29, max diff svdvals: 752.0372910883041
N = 30, max diff svdvals: 934.3056077045787
N = 31, max diff svdvals: 532.6321090336445
N = 32, max diff svdvals: 644.2861698543643
N = 33, max diff svdvals: 362.0563181932318
N = 34, max diff svdvals: 524.8519471908377
N = 35, max diff svdvals: 681.9577012165146
N = 36, max diff svdvals: 959.8179240905338
N = 37, max diff svdvals: 116.35110997021529
N = 38, max diff svdvals: 510.2091542800206
N = 39, max diff svdvals: 245.12683184534578
N = 40, max diff svdvals: 6.869915826281117

Note that if one replace 'A' by 'N' in the gesdd! call everything is fine.

Since it doesn't seem to be clear whether this is a bug in the implementation of gesdd or a shortcoming of the divide-and-conquer approach, I suggest that we default to gesvd for full=true until this is resolved. (Note that this comes with a performance penalty for large matrices.)

[simonbyrne]

Note that if one replace 'A' by 'N' in the gesdd! call everything is fine.

The LAPACK docs state that when just computing singular values both DGESVD and DGESDD use the same algorithm.

The best solution here would probably be to add a keyword argument, say method=, to svd and svd! which lets users choose between the simple, divide-and-conquer (and Jacobi methods if they're also wanted).

[carstenbauer]

Makes sense to me. I'll make a PR some time soon.

[andreasnoack]

I've made a few updates of the original comment since the "issue" isn't related to the use of the full keyword. The difference is between svd and svdvals because LAPACK doesn't use the DnC algorhtm when only the values are requested (as Simon also pointed out in
#607).