acroy / Expokit.jl

Julia implementation of EXPOKIT routines
Other
26 stars 11 forks source link

Rename chbv -> chbmv; improved in-place treatment #23

Open acroy opened 6 years ago

acroy commented 6 years ago

As partly discussed in #20. Might still need some smoothing.

mforets commented 6 years ago

i was thinking if it's more efficient to directly replace the diagonal element of B from the previous one, avoiding copy!(B, A) each time. this is the idea concretely, but i haven't benchmarked it (yet):

function chbmv!{T<:Complex}(w::Vector{T}, A, vec::Vector{T}, fac=factorize)
    p = min(length(θ), length(α))
    t = [θ; θconj]
    a = 0.5 * [α; αconj]

    tmp_v = similar(vec)
    B = A - θ[1]*I
    scale!(copy!(tmp_v, vec), a[1])
    A_ldiv_B!(fac(B), tmp_v)
    copy!(scale!(w, α0), tmp_v)

    @inbounds for i = 2:2*p
        scale!(copy!(tmp_v, vec), a[i])
        for n=1:size(B,1)
            B[n,n] = B[n,n] + t[i-1] - t[i]
        end
        A_ldiv_B!(fac(B), tmp_v)
        w .+=  tmp_v
    end

    return w
end # chbmv!
acroy commented 6 years ago

You are basically right. However if fac is changing B (like lufact!) we need a copy.

mforets commented 6 years ago

ok, i see. one would need to do A_ldiv_B!(fac(copy(B)), tmp_v), and in the end it's pretty much the same speed as doing a copy of A to B and then B[n,n] = A[n,n] - t[i].