Currently, `svd`

defaults to using `LinearAlgebra.LAPACK.gesdd!`

under the hood. I vaguely mentioned some time ago that `gesdd!`

gives incorrect results in some cases which can be circumvented by explicitly calling `LinearAlgebra.LAPACK.gesvd!`

instead. However, for large enough matrices, `gesvd!`

is much slower than `gesdd!`

, which uses a divide-and-conquer approach. Whether the inaccuracies (which can be huge) are due to a bug in the implementation or a shortcoming of this different method I don’t know.

In any case, it naturally raises the question which method we should default to. It seems to me that we should prefer accuracy over speed and default to `gesvd!`

in `svd`

. Of course, we could introduce a keyword flag `divide`

or similar that allows one to switch to the faster method. Just getting wrong results without a warning or anything, as it is the case right now, seems bad.

Note that Octave, which was also defaulting to `gesdd`

has recently made precisely this switch as discussed in this issue.

Sort of taking the example from this thread, here is a MWE that demonstrates the issue:

```
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_gesdd = LinearAlgebra.LAPACK.gesdd!('A', copy(A))[2];
svals_gesvd = LinearAlgebra.LAPACK.gesvd!('A', 'A', copy(A))[2];
return maximum(svals_gesdd - svals_gesvd)
end
function gesdd_gesvd_comparison(Nmax=40)
for N in 1:Nmax
println("N = $(N), max diff svdvals: ", svdvaldiff(N))
end
end
```

Which gives the following output:

```
julia> gesdd_gesvd_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
```