I’ve got some code where I need to perform interpolations at every step of a loop (it’s the integration in time of a PDE). While the code works, currently with `Dierckx`

, I was wondering if any of the interpolation packages can be used such that the memory allocation for the interpolant only needs to be done once, and I just update the interpolation data at each iterate (the number of data points is constnat).

I believe this is possible with https://github.com/JuliaMath/Interpolations.jl

Even if it isn’t possible with that package, it should still be faster.

1 Like

Im actually sure Dierckx can do evaluations into preallocated arrays.

Looks like I was wrong, but it would be very easy to fix; the evaluation routine

```
function _evaluate(t::Vector{Float64}, c::Vector{Float64}, k::Int,
x::Vector{Float64}, bc::Int)
bc in (0, 1, 2, 3) || error("bc = $bc not in (0, 1, 2, 3)")
m = length(x)
xin = convert(Vector{Float64}, x)
y = Vector{Float64}(undef, m)
ier = Ref{Int32}(0)
ccall((:splev_, libddierckx), Nothing,
(Ref{Float64}, Ref{Int32},
Ref{Float64}, Ref{Int32},
Ref{Float64}, Ref{Float64}, Ref{Int32},
Ref{Int32}, Ref{Int32}),
t, length(t), c, k, xin, y, m, bc, ier)
ier[] == 0 || error(_eval1d_messages[ier[]])
return y
end
```

could be rewritten, e.g., like this

```
function _evaluate(t::Vector{Float64}, c::Vector{Float64}, k::Int,
x::Vector{Float64}, bc::Int, y=Vector{Float64}(undef, length(x)))
bc in (0, 1, 2, 3) || error("bc = $bc not in (0, 1, 2, 3)")
m = length(x)
xin = convert(Vector{Float64}, x)
ier = Ref{Int32}(0)
ccall((:splev_, libddierckx), Nothing,
(Ref{Float64}, Ref{Int32},
Ref{Float64}, Ref{Int32},
Ref{Float64}, Ref{Float64}, Ref{Int32},
Ref{Int32}, Ref{Int32}),
t, length(t), c, k, xin, y, m, bc, ier)
ier[] == 0 || error(_eval1d_messages[ier[]])
return y
end
```

this - or something similar - should do what you want. Should be a straightforward PR.

Interpolations appears to not yet support better than linear interpolation when working with off grid data.