svd for a 16x16 BigFloat matrix runs forever

[andreasvarga]

The following computation of a SVD decomposition runs forever on my Windows computer:

using GenericLinearAlgebra
using JLD
T = load("testt.jld","T");
F = svd(T)

These are the singular values for Float64 data:

julia> svdvals(Float64.(T))
16-element Vector{Float64}:
 3.226038366213185
 3.1574963242403324
 2.995684819873801
 1.592942650296512
 1.241696738163066
 1.0508426556228927
 3.716265036915672e-16
 3.0204507505017347e-16
 2.2093670317846095e-16
 1.2432064652599674e-16
 1.1473927285217181e-16
 7.799656720638606e-17
 5.375447706163234e-17
 1.4675517113065786e-32
 9.31581897714817e-33
 0.0

julia> versioninfo()
Julia Version 1.8.3
Commit 0434deb161 (2022-11-14 20:14 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 16 × Intel(R) Core(TM) i7-10700 CPU @ 2.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, skylake)
  Threads: 1 on 16 virtual cores

testt.zip

[andreasnoack]

Thanks for reporting this. I believe I've identified the root cause of this and the fix should be on master soon. I'll tag a patch release so it would be great if you could test it out.

[andreasvarga]

I encountered several similar failures with non-full rank matrices. The provided example was that one with the smallest dimensions. I will provide you a feedback to confirm that the new planed patch release correctly works with my examples.

[andreasnoack]

Thanks. The patch release is available in the registry. Hopefully it will work now. If not please let me know and I'll look into it.

[andreasvarga]

Here is a simple example which fails on my computer:

using GenericLinearAlgebra
a = BigFloat.([0 0 0 0; .1 .1 .1 .1; 2 2 2 2; 3 3 3 3; 0 0 0 0])
svd(a)

[andreasvarga]

I updated to the latest version. Here is an example which is still not working:

using GenericLinearAlgebra
using JLD
T = load("testt2.jld","T");
F = svd(T)

(@v1.8) pkg> status GenericLinearAlgebra
Status `C:\Users\Andreas\.julia\environments\v1.8\Project.toml`
  [14197337] GenericLinearAlgebra v0.3.7

testt2.zip

[andreasnoack]

I found some time tonight to reread https://netlib.org/lapack/lawnspdf/lawn03.pdf and I think I found the issue here. Each time you hit a zero diagonal and deflate, you should recompute an estimate of the smallest remaining singular value and I didn't do that. I've modified the algorithm with such a recalculation in #108 so hopefully things are now working. It would helpful if you could try to break it with more of your matrices.

[andreasvarga]

I performed successfully all tests, where previously error occured. Many thanks for your help to fix this issue.