Closed mforets closed 6 years ago
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.
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.
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))
).
but the vec
doesn't get updated: chbv!(H, v)
returns the correct value, although it is not attached to v
after the call.
True, then maybe chbv!{T}(A, vec::Vector{T}) = chbv!(vec, A, copy(vec))
will do?
Well done! That does it. I've updated the branch.
Do you have an idea why the tests for 0.4 fail? I will try to test locally ...
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..
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?
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".
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.
Done. Shall i squash commits into a single one?
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?
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
.
Cool. I had not seen Documenter.jl, but it looks really nice and useful!
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!
, likepadm!{T}(w::Vector{T}, A, v::Vector{T})
?