ldiv!(A::SVD, B) and ldiv!(Y, A::SVD, B) do not store the result in place

[mzilhao]

The documentation for ldiv! states

  ldiv!(Y, A, B) -> Y                                   
                            
  Compute A \ B in-place and store the result in Y, returning the result.

  ldiv!(A, B)

  Compute A \ B in-place and overwriting B to store the result.

But this does not seem to be honored when A is an SVD object.

julia> A = [-1.0 -1.0 1.0 ; 1.0 3.0 3.0 ; -1 -1 5.0 ];b = [ -1.0 ; 23 ; 15];

The solution x to A x = b is given by

julia> A \ b
3-element Vector{Float64}:
 2.0
 3.0
 4.0

When doing this via an SVD decomposition:

julia> Asvd = svd(A);

julia> x = copy(b)
3-element Vector{Float64}:
 -1.0
 23.0
 15.0

julia> ldiv!(x, Asvd, b)
3-element Vector{Float64}:
 1.9999999999999996
 3.0000000000000004
 4.000000000000001

julia> x
3-element Vector{Float64}:
 -1.0
 23.0
 15.0

Note that ldiv! returned the correct solution, but this was not written into x, which it should according to the documentation. The same happens with the 2-argument version:

julia> ldiv!(Asvd, b)
3-element Vector{Float64}:
 1.9999999999999996
 3.0000000000000004
 4.000000000000001

julia> b
3-element Vector{Float64}:
 -1.0
 23.0
 15.0

When using LU factorization everything is as expected:

julia> Alu = lu(A);

julia> ldiv!(x, Alu, b)
3-element Vector{Float64}:
 2.0
 3.0
 4.0

julia> x
3-element Vector{Float64}:
 2.0
 3.0
 4.0

julia> ldiv!(Alu, b)
3-element Vector{Float64}:
 2.0
 3.0
 4.0

julia> b
3-element Vector{Float64}:
 2.0
 3.0
 4.0

julia> versioninfo()
Julia Version 1.6.1
Commit 6aaedecc44 (2021-04-23 05:59 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-8565U CPU @ 1.80GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)

[andreasnoack]

I hit that one yesterday as well so it's being fixed as part of JuliaLang/julia#40899

[aviks]

Close, now that JuliaLang/julia#40899 is merged?