Skip to content

Commit

Permalink
changes functions finite_differences
Browse files Browse the repository at this point in the history
improve allocations in computing
  • Loading branch information
rveltz committed Dec 21, 2024
1 parent aa95a45 commit 18b286a
Show file tree
Hide file tree
Showing 4 changed files with 15 additions and 30 deletions.
33 changes: 9 additions & 24 deletions src/Utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -84,25 +84,6 @@ function finite_differences(F, x::AbstractVector; δ = 1e-9)
N = length(x)
Nf = length(F(x))
J = zeros(eltype(x), Nf, N)
finite_differences!(F, J, x; δ)
return J
end

function finite_differences(F, x::AbstractVector, p; δ = 1e-9)
N = length(x)
Nf = length(F(x, p))
J = zeros(eltype(x), Nf, N)
finite_differences!(F, J, x, p; δ = δ)
return J
end

"""
$(SIGNATURES)
Compute a Jacobian by Finite Differences, update J. Use the centered formula (f(x+δ)-f(x-δ))/2δ.
"""
function finite_differences!(F, J, x::AbstractVector; δ = 1e-9)
f = F(x)
x1 = copy(x)
@inbounds for i in eachindex(x)
x1[i] += δ
Expand All @@ -115,15 +96,19 @@ function finite_differences!(F, J, x::AbstractVector; δ = 1e-9)
return J
end

function finite_differences!(F, J, x::AbstractVector, p; δ = 1e-9)
f = F(x,p)
"""
$(SIGNATURES)
Same as finite_differences but with inplace `F`
"""
@views function finite_differences!(F, J, x::AbstractVector; δ = 1e-9, tmp = copy(x))
x1 = copy(x)
@inbounds for i in eachindex(x)
x1[i] += δ
J[:, i] .= F(x1,p)
F(J[:, i], x1)
x1[i] -= 2δ
J[:, i] .-= F(x1,p)
J[:, i] ./= 2δ
F(tmp, x1)
J[:, i] .= @. (J[:, i] - tmp) / (2δ)
x1[i] += δ
end
return J
Expand Down
8 changes: 4 additions & 4 deletions src/periodicorbit/Floquet.jl
Original file line number Diff line number Diff line change
Expand Up @@ -420,7 +420,7 @@ end
vals = values[indvalid]
# these are the Floquet multipliers
μ = @. Complex(1 / (1 + vals))
vp0 = minimum(abslog, μ)
vp0 = minimum(abs log, μ)
if vp0 > 1e-8
@warn "The precision on the Floquet multipliers is $vp0. Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetCollGEV`"
end
Expand Down Expand Up @@ -455,7 +455,7 @@ end
_eig_floquet_col(J, n, m, Ntst, nev)
end

function _eig_floquet_col(J::AbstractMatrix{𝒯}, n, m, Ntst, nev) where 𝒯
@views function _eig_floquet_col(J::AbstractMatrix{𝒯}, n, m, Ntst, nev) where 𝒯
nbcoll = n * m
N = n
In = LinearAlgebra.I(N)
Expand All @@ -478,7 +478,7 @@ function _eig_floquet_col(J::AbstractMatrix{𝒯}, n, m, Ntst, nev) where 𝒯

Jcop = copy(J)
Jcop[end-N:end-1,end-N:end-1] .= In
Jcop[end-N:end-1,1:N] .= -In
Jcop[end-N:end-1,1:N] .= (-1) .* In
Jcop[end, end] = J[end,end]
Lₜ = LowerTriangular(copy(blockⱼ))
for 𝐢 in 1:Ntst
Expand Down Expand Up @@ -532,7 +532,7 @@ function _eig_floquet_col(J::AbstractMatrix{𝒯}, n, m, Ntst, nev) where 𝒯
# give indications on the precision on the Floquet coefficients
vp0 = minimum(abs, σ)
if vp0 > 1e-9
@warn "The precision on the Floquet multipliers is $vp0.\n It may be not enough to allow for precise bifurcation detection.\n Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetColl` for computing Floquet coefficients."
@debug "The precision on the Floquet multipliers is $vp0.\n It may be not enough to allow for precise bifurcation detection.\n Either decrease `tol_stability` in the option ContinuationPar or use a different method than `FloquetColl` for computing Floquet coefficients."
end
return σ, Complex.(vecs[I, :]), true, 1
end
2 changes: 1 addition & 1 deletion test/problems.jl
Original file line number Diff line number Diff line change
Expand Up @@ -31,4 +31,4 @@ BK.re_make(prob, J = (x,p)->zeros(2,2), Jᵗ = (x,p)->zeros(2,2), d2F=(x,p,dx1,d
######################################################################
# test finite differences
BK.finite_differences(identity, zeros(2))
BK.finite_differences((x,p)->x, zeros(2), nothing)
BK.finite_differences!((o,x)->o.=x, zeros(2, 2), zeros(2))
2 changes: 1 addition & 1 deletion test/simple_continuation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -130,7 +130,7 @@ br0 = @time continuation(prob,
###############
br0 = @time continuation(prob, PALC(), opts; callback_newton = BK.cbMaxNormAndΔp(10,10)) #(6.20 k allocations: 409.469 KiB)
try
continuation(prob, PALC(), opts; callback_newton = BK.cbMaxNormAndΔp(10,10), bla = 1)
continuation(prob, PALC(), opts; callback_newton = BK.cbMaxNormAndΔp(10,10))
catch
end

Expand Down

0 comments on commit 18b286a

Please sign in to comment.