Add Chbv

[mforets]

No description provided.

[mforets]

in commit bfe6cee, the inner accumulations are done inplace. notice that the auxiliary vector y should be complex in both cases.

do we want to have an inplace version similar to expmv!, like padm!{T}(w::Vector{T}, A, v::Vector{T})?

[acroy]

I guess this could be p=length(α) (or even min(length(α),length(θ))) to make sure that we stay in the bounds. One could also put θ and α into one matrix ...

[mforets]

ok, changed.

[acroy]

If one wants to keep y real, one could just sum the real parts.

[mforets]

yes, thanks for the observation; this allows me to get rid of w.

[acroy]

I think we don't need to construct an explicit identity matrix, we can rely on UniformScaling (I) below?

[mforets]

true, this saves some allocation, changed.

[mforets]

FYI for this function i collected some benchmark data, it's here.

Now i've started to do the chbv! version, locally, but i get some errors.

[acroy]

I was wondering if we should add those by hand to the definitions of θ and α above. It is not saving a lot of memory, but will reduce garbage collection if called frequently.

[mforets]

you mean to define a table for theta_conj and alpha_conj as well, right? yes, we can add those defs on top.

[acroy]

Yes, or you can just add it to the vectors you already have (but then p=length(α)/2).

[mforets]

i'm struggling with the inplace version, let me put the broken version here 🤕

actually if i use the three input inplace versions, chbv!(Hv, H, v), work properly, but the two-component one is wrong.

function chbv{T}(A, vec::Vector{T})
    result = convert(Vector{promote_type(eltype(A), T)}, copy(vec))
    chbv!(A, result)
    return result
end

# wrong, use: chbv!(H, v);
chbv!{T}(A, vec::Vector{T}) = chbv!(vec, A, vec)

# OK, use: Hv = similar(v); chbv!(Hv, H, v); norm(Hv - expm(H) * v) ~= 1e-12
function chbv!{T<:Real}(w::Vector{T}, A, vec::Vector{T})
    p = min(length(θ), length(α))
    scale!(copy!(w, vec), α0)

    @inbounds for i = 1:p
        w .+= real((A - θ[i]*I) \ (α[i] * vec))
    end
    return w
end

# OK, use: Hv = similar(v); chbv!(Hv, H, v); norm(Hv - expm(H) * v) ~= 1e-12
function chbv!{T<:Complex}(w::Vector{T}, A, vec::Vector{T})
    p = min(length(θ), length(α))
    scale!(copy!(w, vec), α0)
    t = [θ; conj(θ)]
    a = 0.5 * [α; conj(α)]

    @inbounds for i = 1:2*p
        w .+= (A - t[i]*I) \ (a[i] * vec)
    end
    return w
end

EDIT: some comments and removed vec=w in the end of each function.

[acroy]

The two argument version of chbv! cannot work since you overwrite the input vec which you need further down. I think you have to do chbv!{T}(A, vec::Vector{T}) = chbv!(similar(vec), A, vec) to make it work (or even using similar(vec, promote_type(eltype(A), T))).

[mforets]

but the vec doesn't get updated: chbv!(H, v) returns the correct value, although it is not attached to v after the call.

[acroy]

True, then maybe chbv!{T}(A, vec::Vector{T}) = chbv!(vec, A, copy(vec)) will do?

[mforets]

Well done! That does it. I've updated the branch.

[acroy]

I think this will result in an unnecessary copy. It's better to use the three argument version here.

[acroy]

Do you have an idea why the tests for 0.4 fail? I will try to test locally ...

[mforets]

No idea. I'm in v0.6 and for the failing test in v0.4 i get that the residuum is smaller than 1e-11 for a set of 10.000 tests.

I saw that JuliaBox proposes older kernels too, i can check there..

[mforets]

i can confirm that there is an issue with v0.4.7.

with the 3 argument variant, it does produce the correct result, but it doesn't modify its first argument when the function is returned. like a different handling of variables scope?

[mforets]

to be precise, this happens:

function f(w)
    w .+= ones(length(w))
    return w
end
w = 2.0 * ones(4);
f(w), w
# ([3.0,3.0,3.0,3.0],[2.0,2.0,2.0,2.0]) in v0.4.7
# ([3.0, 3.0, 3.0, 3.0], [3.0, 3.0, 3.0, 3.0]) in v0.5.1 and later

i dunno if Compat is supposed to "fix this".

[acroy]

Ah, good catch! This is probably because .+= was not in place pre 0.5. We should be able to "fix" this by returning chbv!(result, A, vec) instead of result. At least the tests will pass.

[mforets]

Done. Shall i squash commits into a single one?

[acroy]

Yes, that would be great. We will also have to add the information to the README file, or do you want to generate that from the docstrings?

[mforets]

ok, done.

regarding documentation, yes, in case you haven't tried it yet with Documenter.jl the setup is super easy and plays very well with github. kudos to Documenter devs :)
in this case we can enclose the function's name in a @docs macro, in a file that can be e.g. docs/src/lib/chbv.md.

[acroy]

Cool. I had not seen Documenter.jl, but it looks really nice and useful!