Support SharedArrays

[ranocha]

SharedArrays are currently not supported

julia> using OrdinaryDiffEq
prob
julia> prob = ODEProblem((t,u,du)-> begin du .= u; nothing end, SharedArray{Float64}(3,2), (0., 1.))
DiffEqBase.ODEProblem with uType SharedArray{Float64,2} and tType Float64. In-place: true
timespan: (0.0, 1.0)
u0: [0.0 0.0; 0.0 0.0; 0.0 0.0]

julia> solve(prob, Euler(), dt=1.e-2)
ERROR: MethodError: no method matching OrdinaryDiffEq.EulerCache(::SharedArray{Float64,2}, ::SharedArray{Float64,2}, ::Array{Float64,2}, ::Array{Float64,2}, ::Array{Float64,2})
Closest candidates are:
  OrdinaryDiffEq.EulerCache(::uType, ::uType, ::uType, ::rateType, ::rateType) where {uType, rateType} at .../.julia/v0.6/OrdinaryDiffEq/src/caches/low_order_rk_caches.jl:2
Stacktrace:
 [1] alg_cache(::OrdinaryDiffEq.Euler, ::SharedArray{Float64,2}, ::Array{Float64,2}, ::Type{T} where T, ::Type{T} where T, ::SharedArray{Float64,2}, ::SharedArray{Float64,2}, ::Function, ::Float64, ::Float64, ::Float64, ::Type{Val{true}}) at .../.julia/v0.6/OrdinaryDiffEq/src/caches/low_order_rk_caches.jl:29
 [2] #init#1831(::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::Void, ::Bool, ::Void, ::Bool, ::Bool, ::Bool, ::Bool, ::Float64, ::Bool, ::Rational{Int64}, ::Void, ::Void, ::Int64, ::Rational{Int64}, ::Int64, ::Int64, ::Rational{Int64}, ::Bool, ::Int64, ::Rational{Int64}, ::Rational{Int64}, ::Int64, ::Float64, ::Float64, ::DiffEqBase.#ODE_DEFAULT_NORM, ::DiffEqBase.#ODE_DEFAULT_ISOUTOFDOMAIN, ::DiffEqBase.#ODE_DEFAULT_UNSTABLE_CHECK, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Int64, ::String, ::DiffEqBase.#ODE_DEFAULT_PROG_MESSAGE, ::Void, ::Void, ::Bool, ::Bool, ::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.ODEProblem{SharedArray{Float64,2},Float64,true,##1#2,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Euler, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at .../.julia/v0.6/OrdinaryDiffEq/src/solve.jl:241
 [3] (::DiffEqBase.#kw##init)(::Array{Any,1}, ::DiffEqBase.#init, ::DiffEqBase.ODEProblem{SharedArray{Float64,2},Float64,true,##1#2,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Euler, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0
 [4] #solve#1830(::Array{Any,1}, ::Function, ::DiffEqBase.ODEProblem{SharedArray{Float64,2},Float64,true,##1#2,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Euler, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at .../.julia/v0.6/OrdinaryDiffEq/src/solve.jl:6
 [5] (::DiffEqBase.#kw##solve)(::Array{Any,1}, ::DiffEqBase.#solve, ::DiffEqBase.ODEProblem{SharedArray{Float64,2},Float64,true,##1#2,Void,UniformScaling{Int64},DiffEqBase.StandardODEProblem}, ::OrdinaryDiffEq.Euler, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at ./<missing>:0 (repeats 2 times)

I'm using julia v0.6.1

julia> versioninfo()
Julia Version 0.6.1-pre.1
Commit 4b967d4a9d* (2017-08-21 21:11 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: Intel(R) Core(TM) i5-4590S CPU @ 3.00GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.9.1 (ORCJIT, haswell)

This is due to the problem that copy(::SharedArray) and similar(::SharedArray) do not return SharedArrays

julia> u0 = SharedArray{Float64}(3,2)
3×2 SharedArray{Float64,2}:
 0.0  0.0
 0.0  0.0
 0.0  0.0

julia> copy(u0)
3×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0
 0.0  0.0

julia> deepcopy(u0)
3×2 SharedArray{Float64,2}:
 0.0  0.0
 0.0  0.0
 0.0  0.0

See also JuliaLang/julia#12964, JuliaLang/julia#22218, https://discourse.julialang.org/t/why-does-similar-sharedarray-create-an-array/6673/3.

Shall we use deepcopy instead of copy for SharedArrays?

[ChrisRackauckas]

I think we just try to get SharedArrays fixed since the current behavior seems broken in many ways. If we need it right now, since it's a Base type we can add some special handling. What's the use case?

[ranocha]

I'm experimenting with parallelisation strategies for (hyperbolic) PDE discretisations. I've tried using threads and standard Arrays and wanted to test SharedArrays.

[ranocha]

Shall I implement such a check in the initialisation part?

[ChrisRackauckas]

Shall I implement such a check in the initialisation part?

If that fixes it, go for it.

[ChrisRackauckas]

I've tried using threads and standard Arrays and wanted to test SharedArrays.

@dextorious found before (in a different case) that SharedArrays have quite a bit of overhead. I think the more interesting case is trying DistributedArrays and a strategy for having it automatically handle the operator's actions at the boundaries between processes.

[ranocha]

Sounds also interesting. I will write tests concerning the usage of both types of parallel arrays.

I found a very strange behaviour of threads. I wanted to construct types with field parallel::Parallel where Parallel==Val{:serial} or Parallel==Val{:threads} or something similar in order to be able to turn threading on and off and compare the performance. I found the very strange behaviour that I have to run the serial code once bafore I run the parallel code. Otherwise, the parallel code has a very bad performance. However, I've not been able to build a minimal example. Therefore, I wanted to know whether something similar happens to SharedArrays.

Does this approach with parallel::Parallel as part of the semidiscretisation (i.e. the f) seem to be useable or is it just stupid?

[ChrisRackauckas]

Does this approach with parallel::Parallel as part of the semidiscretisation (i.e. the f) seem to be useable or is it just stupid?

You mean auto-parallelism here? I am not sure what you're asking. A user can always parallelize in their f.

[ranocha]

I have something like

struct Semidiscretisation{OtherStuff, Parallel}
  other_stuff::OtherStuff
  parallel::Parallel
end

function (semidiscretisation::Semidiscretisation)(t, u, du)
  foo!(du, u, t, semidiscretisation.other_stuff, semidiscretisation.parallel)
end

where a serial version is used for semidiscretisation.parallel == Val{:serial}() and threading is used if semidiscretisation.parallel == Val{:threads}(). Does this sound acceptable or just stupid and there is a better alternative to switch between parallel and serial verisons?

[ChrisRackauckas]

Nah, that seems fine. Normally I think you can assume that operators will faster when threaded though, so DiffEqOperators.jl implicitly multithreads its A_mul_B!. However, I am not sure we have Clik-style multithreading yet (core keeps saying stuff about it though?) so putting multithreaded codes on top of each other could mess with it, so the ability to turn it off like that could be useful.

[ranocha]

Thanks for sharing your opinion.

so the ability to turn it off like that could be useful.

That's what I thought. Additionally, threading prevents the normal error messages and makes debugging much more complicated if some errors occur. And it seems to be complicated to get all the performance from threading.

[ChrisRackauckas]

And it seems to be complicated to get all the performance from threading.

This frequently gets in the way: JuliaLang/julia#15276

[ranocha]

This frequently gets in the way: JuliaLang/julia#15276

Yeah, that's it.

It seems to be possible to use SharedArrays with the following approach:

same(u::AbstractArray) = similar(u)
same(u::SharedArray) = typeof(u)(size(u)...)

Then, we would have to use same(u) instead of similar(u) in many places, especially the functions alg_cache(...). A typical example is alg_cache(alg::BS3, ...), where tmp = similar(u) has to be changed to tmp = same(u). This seems to work with broadcasting and indexing, probably without parallelisation.
@ChrisRackauckas Do you think we should implement it this way? Or should we just stay with standard Arrays etc. and ignore this part at the moment? Maybe there will be some better alternatives in the future, especially concerning broadcasting and automatic parallelisation.

[ChrisRackauckas]

Do you think we should implement it this way? Or should we just stay with standard Arrays etc. and ignore this part at the moment?

I think we should keep to similar since then the code is much easier for newcomers to understand. If we go too far down in our own direction it can get harder and harder to follow. The issue here is that SharedArrays should fix its similar output, and changing the similar dispatch is a quick hack to support it if you need it now.

Maybe there will be some better alternatives in the future, especially concerning broadcasting and automatic parallelisation.

I think that's the real answer. An array type that multithreads broadcast. I thought GPUArrays' JLArray did that, but I guess not. We can go bug Simon to make it happen 👍. In general I don't think SharedArrays are that great for this kind of usage.

[ranocha]

Okay, then I'll close this issue now. Thanks for your thoughts!