JuliaLang
LinearAlgebra: LAPACK.gelsy!
I'm using julia 1.8.3 for Linux installed from https://julialang.org/downloads/.
julia> versioninfo()
Julia Version 1.8.3
Commit 0434deb161e (2022-11-14 20:14 UTC)
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 16 × Intel(R) Xeon(R) W-2145 CPU @ 3.70GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-13.0.1 (ORCJIT, skylake-avx512)
Threads: 1 on 16 virtual cores
Environment:
LD_LIBRARY_PATH = ...
I am using LAPACK.gelsy!(A,B) to solve AX=B. The docs say it overwrites B with the solution, but I do not observe this, e.g.,
julia> A = [1 1 1; 2 2 2; 1 0 1.]
3×3 Matrix{Float64}:
1.0 1.0 1.0
2.0 2.0 2.0
1.0 0.0 1.0
julia> b = [1,-3,-1.]
3-element Vector{Float64}:
1.0
-3.0
-1.0
julia> (x, rnk) = LAPACK.gelsy!(A, b)
([-0.5000000000000004, 1.0402585548239313e-15, -0.5000000000000007], 2)
julia> hcat(x,b)
3×2 Matrix{Float64}:
-0.5 1.0
1.04026e-15 -3.0
-0.5 -1.0
I observe that this also happens when B is a matrix, e.g.,
julia> A = [1 1 1; 2 2 2; 1 0 1.]
3×3 Matrix{Float64}:
1.0 1.0 1.0
2.0 2.0 2.0
1.0 0.0 1.0
julia> B = Matrix([1 -3 -1.]')
3×1 Matrix{Float64}:
1.0
-3.0
-1.0
julia> (X, rnk) = LAPACK.gelsy!(A, B)
([-0.5000000000000004; 1.0402585548239313e-15; -0.5000000000000007;;], 2)
julia> hcat(X,B)
3×2 Matrix{Float64}:
-0.5 1.0
1.04026e-15 -3.0
-0.5 -1.0
Perhaps it is related to: https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/lapack.jl#L1340.
@jacob-roth That is correct. XGELSY needs much more working space in B than its "correct" dimension. Not sure how else to handle this...
Perhaps the docstring should then be rewritten to reflect this. The line you were pointing at clearly allocates new memory, and then likely overwrites that.
Indeed.
However, while we're at this, would it be perhaps better to let the users manage the memory themselves (via new optional parameters to the method) instead of the method allocating newB and ipiv every call?
That way, users who repeatedly call gelsy! (and are possibly doing this for performance reasons) can avoid the allocations by reusing memory on their end...