Memory Allocation and Interpolation

General Usage
1

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).

2

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
3

Im actually sure Dierckx can do evaluations into preallocated arrays.

4

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.

5

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