I want to solve a linear system, A*x = b for a fixed A but different b's.
Further, I'd like to overwrite the solution (i.e. x) in b. I found a LAPACK
command that does this:
julia> import Base.LinAlg.LAPACK: gesv!
julia> A = randn(100,100);
julia> B = randn(100);
julia> x1 = A \ B;
julia> gesv!(A,B);
julia> isapprox(B,x1)
true
It makes a lot of sense for me to cache the LU factorization of A, since
I'll be solving this multiple times, but gesv! doesn't seem to accept an LU
factorization as input:
julia> A_lu = lufact(A);
julia> gesv!(A_lu,B)
ERROR: MethodError: no method matching gesv!(::Base.LinAlg.LU{Float64,Array{
Float64,2}}, ::Array{Float64,1})
My guess is that there is a different function that will operate like
gesv!, but accepts an LU factorization as input? Any pointers would be
greatly appreciated.
Apologies if this is documented somewhere obvious. I dug around but
couldn't find much, and I'm somewhat new to this LAPACK business.
Thanks in advance.
-- Alex
P.S. The documentation
<http://docs.julialang.org/en/release-0.4/stdlib/linalg/#Base.LinAlg.LAPACK.gesv!>
says
that gesv! should overwrite A with its LU factorization. I don't see this.
typeof(A) == Array{Float64,2} after I call gesv!(A,B).
