Shared Memory Parallelizaton (Question)

[psv2]

Does Interpolations.jl support shared memory parallelization in Julia, as is present in the SharedArray construct for standard arrays? I ask because I want to interpolate an array of size 300753007520 over all dimensions but the last (of size 20); this array already has 10^10 elements, which would certainly require shared memory even leaving the issue of interpolation aside.

[timholy]

I believe an interpolation object can be constructed from any kind of AbstractArray, so in principle you should be fine. However, keep in mind that interpolation of order higher than Linear requires inversion of a tridiagonal system at "construction" time, so you may be in for a big surprise in terms of memory allocation and runtime. If you do want higher-order interpolation, you'll likely want to consider interpolate! and the InPlace() or InPlaceQ() boundary conditions. You may also want to created a distributed version of prefilter.

For such a demanding application, you'll basically have to dive into the code to see how things work. The doc folder might have some useful resources, especially https://github.com/JuliaMath/Interpolations.jl/blob/master/doc/devdocs.md.

[psv2]

Thanks for the reply. I'll give it a shot; I only need to interpolate this large array once (and then evaluate it many times, passing it as a function argument), so I'm not too concerned about runtime increasing from that single interpolation. Also, just to clarify: what exactly does the InPlace() boundary condition do? (There doesn't seem to be much documentation regarding its function, as opposed to the other popular boundary conditions.)

[psv2]

An additional clarification that I forgot to ask in the last comment: if I use interpolate! on a SharedArray, will the resulting object still be stored in shared memory, suitable for parallelization in Julia?

[timholy]

Also, just to clarify: what exactly does the InPlace() boundary condition do? (There doesn't seem to be much documentation regarding its function, as opposed to the other popular boundary conditions.)

It was documented once (in LaTeX), but I don't think that documentation is hanging around anymore. Your best bet is to look at the source code for details. Suffice it to say that it matches the values at the grid points while also not adding any padding.

if I use interpolate! on a SharedArray, will the resulting object still be stored in shared memory

Yes, because it's doing everything in place. Note that if the SharedArray is backed by a disk file, that involves writing data to disk.

[psv2]

Thanks for the information. I'll give this a go, and only reopen this issue if I find a feature related to this to be problematic or lacking in some way.

[psv2]

I'm reopening this issue because it seems like there is an issue with interpolate() versus interpolate!(), specifically with regard to how types are handled. I attach the following code as an illustration.
julia> C = SharedArray(Float64, 6, 6); julia> for ii=1:6 for jj=1:6 C[ii,jj] = sqrt((ii-1)^2 + (jj-1)^2) end end
produces
julia> C 6×6 SharedArray{Float64,2}: 0.0 1.0 2.0 3.0 4.0 5.0 1.0 1.41421 2.23607 3.16228 4.12311 5.09902 2.0 2.23607 2.82843 3.60555 4.47214 5.38516 3.0 3.16228 3.60555 4.24264 5.0 5.83095 4.0 4.12311 4.47214 5.0 5.65685 6.40312 5.0 5.09902 5.38516 5.83095 6.40312 7.07107
confirmed by
julia> typeof(C) SharedArray{Float64,2}.
However, if I then do
julia> D = interpolate(C, (BSpline(Cubic(Natural())), NoInterp()), OnGrid()) 6×6 Interpolations.BSplineInterpolation{Float64,2,Array{Float64,2},Tuple{Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line}},Interpolations.NoInterp},Interpolations.OnGrid,(1,0)}: -2.77556e-17 1.0 2.0 3.0 4.0 5.0 1.0 1.41421 2.23607 3.16228 4.12311 5.09902 2.0 2.23607 2.82843 3.60555 4.47214 5.38516 3.0 3.16228 3.60555 4.24264 5.0 5.83095 4.0 4.12311 4.47214 5.0 5.65685 6.40312 5.0 5.09902 5.38516 5.83095 6.40312 7.07107
while the output is correct, the type is wrong:
julia> typeof(D) Interpolations.BSplineInterpolation{Float64,2,Array{Float64,2},Tuple{Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line}},Interpolations.NoInterp},Interpolations.OnGrid,(1,0)}.
Likewise,
julia> E = interpolate(deepcopy(C), (BSpline(Cubic(Natural())), NoInterp()), OnGrid()) 6×6 Interpolations.BSplineInterpolation{Float64,2,Array{Float64,2},Tuple{Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line}},Interpolations.NoInterp},Interpolations.OnGrid,(1,0)}: -2.77556e-17 1.0 2.0 3.0 4.0 5.0 1.0 1.41421 2.23607 3.16228 4.12311 5.09902 2.0 2.23607 2.82843 3.60555 4.47214 5.38516 3.0 3.16228 3.60555 4.24264 5.0 5.83095 4.0 4.12311 4.47214 5.0 5.65685 6.40312 5.0 5.09902 5.38516 5.83095 6.40312 7.07107
(using deepcopy() instead of copy() or the object itself)
produces the wrong type
julia> typeof(E) Interpolations.BSplineInterpolation{Float64,2,Array{Float64,2},Tuple{Interpolations.BSpline{Interpolations.Cubic{Interpolations.Line}},Interpolations.NoInterp},Interpolations.OnGrid,(1,0)}
with no hint of SharedArray-ness (as with D).
If I now do
julia> F = interpolate!(deepcopy(C), (BSpline(Cubic(InPlace())), NoInterp()), OnGrid()) 6×6 Interpolations.BSplineInterpolation{Float64,2,SharedArray{Float64,2},Tuple{Interpolations.BSpline{Interpolations.Cubic{Interpolations.InPlace}},Interpolations.NoInterp},Interpolations.OnGrid,0}: -1.11022e-16 1.11849 2.08797 3.06136 4.0455 5.03577 1.0 1.48215 2.29545 3.20911 4.16076 5.13021 2.0 2.25346 2.85922 3.63785 4.50194 5.41177 3.0 3.16805 3.62046 4.26269 5.0215 5.85202 4.0 4.12562 4.47988 5.01227 5.67176 6.41909 5.0 5.09828 5.38574 5.83544 6.41143 7.08193
the type is correct
julia> typeof(F) Interpolations.BSplineInterpolation{Float64,2,SharedArray{Float64,2},Tuple{Interpolations.BSpline{Interpolations.Cubic{Interpolations.InPlace}},Interpolations.NoInterp},Interpolations.OnGrid,0}
but the values are wrong. What's going on?

[psv2]

For the record, this was tested on Julia 0.5.1.

[andrewsteck]

Did anything ever come of this? Or else can someone point me in the direction to start digging?

I am seeing similar issues in Julia 0.6. I have a different use-case (need many, smaller interpolations), but the SharedArray construction would be very helpful in speeding things up pre- and post-interpolation.

[psv2]

In Julia 0.4.6, I can say that it is possible to use interpolate() on a multidimensional SharedArray object; in fact, I've been able to create arrays of such interpolation objects and manipulate them without issues. I haven't tried on newer versions.

[timholy]

Presumably this is just the behavior of similar; see JuliaLang/julia#12964 for why this happens. It looks like you can allocate the array manually and call prefilter! directly, though.

[andrewsteck]

Great, thanks. I'll push forward then and report back if I have any new issues...