Skip to content

Speed up Thiele construction by using array evaluatoin, and add barycentric array evaluation #29

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
tobydriscoll merged 2 commits into from
Jun 16, 2026
Merged

Speed up Thiele construction by using array evaluatoin, and add barycentric array evaluation #29

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
317d8fe
use array evaluation for Thiele construction
tobydriscoll Jun 16, 2026
4930bb1
faster evaluation of Barycentric for array inputs
tobydriscoll Jun 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 25 additions & 1 deletion src/barycentric.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -154,12 +154,36 @@ function evaluate(r::Barycentric, z::Number)
end
end

# faster version for evaluation on an array
function evaluate(r::Barycentric, z::AbstractArray{<:Number})
u = similar(z, eltype(r.values))
return evaluate!(u, r, z)
end

function evaluate!(u::AbstractArray, r::Barycentric, z::AbstractArray{<:Number})
num = @. r.w_times_f[1] / (z - r.nodes[1])
den = @. r.weights[1] / (z - r.nodes[1])
safe = trues(size(z))
idx = isinf.(num)
safe[idx] .= false
u[idx] .= r.values[1]
@inbounds for k in 2:length(r.nodes)
c = @. 1 / (z - r.nodes[k])
idx = isinf.(c)
safe[idx] .= false
u[idx] .= r.values[k]
axpy!(r.w_times_f[k], c, num)
axpy!(r.weights[k], c, den)
end
@views u[safe] .= num[safe] ./ den[safe]
return u
end

function _evaluate_numden(r::Barycentric, z::Number)
C = @. 1 / (z - r.nodes)
return sum(C .* r.w_times_f), sum(C .* r.weights)
end


# Evaluate when given Cauchy matrix
function evaluate!(u::AbstractArray, r::Barycentric, C::AbstractMatrix)
n = length(r.nodes)
Expand Down
33 changes: 17 additions & 16 deletions src/thiele.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,7 +107,7 @@ function _evaluate_onediv(r::Thiele, z::Number)
else
numer / denom
end
end
end

_evaluator = _evaluate_onediv # default choice

Expand Down Expand Up @@ -338,6 +338,7 @@ function approximate(::Type{Thiele},
idx_test = CartesianIndices((1:1, 2:num_ref+1))
idx_new_test = idx_test
fτ = Matrix{eltype(fσ)}(undef, size(τ)) # f at test points
rτ = similar(complex(fτ)) # r at test points
fτ[idx_test] .= f.(τ[idx_test])
fmax = maximum(abs, view(fτ, idx_test)) # scale of f

Expand All @@ -349,7 +350,8 @@ function approximate(::Type{Thiele},
n = 1 # iteration counter
while true
# test_actual = view(fτ, idx_test) # f at test points
err = @. abs(fτ[idx_test] - r(τ[idx_test]))
evaluate!(view(rτ, idx_test), r, view(τ, idx_test)) # r at test points
err = abs.(view(fτ, idx_test) - view(rτ, idx_test))
err_max, idx_max = findmax(err)
history[n].error = err_max

Expand Down Expand Up @@ -417,21 +419,18 @@ function approximate(::Type{Thiele},
fmax = maximum(abs, y) # scale of f
_, idx_min = findmin(abs, y)
r = Thiele([z[idx_min]], [y[idx_min]])
idx_test[idx_min] = false
y_test = copy(collect(y))
z_test = copy(collect(z))
deleteat!(y_test, idx_min)
deleteat!(z_test, idx_min)
r_test = similar(y_test)

history = [IterationRecord(r, NaN, missing)]
n = 1 # iteration counter
while length(z) > 0
idx_max, err_max = 0, -Inf
for i in eachindex(y)
if idx_test[i]
err = abs(y[i] - r(z[i]))
if err > err_max
err_max = err
idx_max = i
end
end
end
evaluate!(r_test, r, z_test)
r_test .-= y_test
err_max, idx_max = findmax(abs(e) for e in r_test)
history[n].error = err_max

status = quitting_check(history, stagnation, tol, fmax, max_iter, allowed)
Expand All @@ -448,14 +447,16 @@ function approximate(::Type{Thiele},

# Add new node:
try
add_node!(r, z[idx_max], y[idx_max])
add_node!(r, z_test[idx_max], y_test[idx_max])
push!(history, IterationRecord(r, NaN, missing))
idx_test[idx_max] = false
deleteat!(y_test, idx_max)
deleteat!(z_test, idx_max)
deleteat!(r_test, idx_max)
catch(e)
# look for the best acceptable case
status = quitting_check(history, stagnation, tol, fmax, 1, allowed)
r = history[status].interpolant
@warn("NaN weight encountered; stopping with estimated error $(round(history[status].error, sigdigits=4))")
@warn("Adding node failed; stopping with estimated error $(round(history[status].error, sigdigits=4))")
@debug("Error $e")
break
end
Expand Down
Loading