# Parallel computing on a subset of an array

**URL:** <https://discourse.julialang.org/t/parallel-computing-on-a-subset-of-an-array/11881>\
**Category:** Performance\
**Tags:** question\
**Created:** [June 22, 2018, 3:00pm UTC](https://discourse.julialang.org/t/parallel-computing-on-a-subset-of-an-array/11881 "2018-06-22T15:00:42Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Balinus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/balinus/32/243_2.png) [@Balinus](https://discourse.julialang.org/u/Balinus)\
**Post date:** [June 22, 2018, 3:00pm UTC](https://discourse.julialang.org/t/parallel-computing-on-a-subset-of-an-array/11881/1 "2018-06-22T15:00:42Z")

</div>

Hello!

I’m trying to do some parallel computation on the 2 first dimensions of a 3D array `A[j,k, t]`. To put things in context, the data are timeseries, where `j` and `k` are spatial grid points and `t` is the time dimension. There is also no possibility to permutedims in this context (Climate dataset can be huge).

My objective is to compute the same function for each grid points, using the whole timeserie for a given grid point `j, k`. Hence, going `eachindex` but only for the spatial dimensions.

The first attempt is to simply put `Threads.@threads` before one of the “spatial loop”. Here’s a MWE.

```julia
using BenchmarkTools

function simplethread()
    A = randn(3, 3, 10000)
    dataout = fill(NaN, 3, 3, 10000)
    
    for k = 1:size(A, 2)
        Threads.@threads for j = 1:size(A, 1)
            val = somefunction(A[j,k,:])
            dataout[j,k,:] = val
        end
    end
    return dataout
end

function somefunction(datain) # (in real case, I do `Polynomials.polyval` over the timeserie using a previous `polyfit`).
       dataout = datain .+ randn(1)
end

@benchmark simplethread()
BenchmarkTools.Trial: 
  memory estimate: 2.76 MiB
  allocs estimate: 442
  --------------
  minimum time: 1.414 ms (0.00% GC)
  median time: 1.576 ms (0.00% GC)
  mean time: 2.046 ms (16.58% GC)
  maximum time: 231.659 ms (98.22% GC)
  --------------
  samples: 2438
  evals/sample: 1

```

Now, I want to effectively do calculations in parallel for each grid point, not just one of the spatial dimension. I tried using CartesianRange, but it does not work. Here’s my attempt :

```julia
function cartesianthread()
    A = randn(3, 3, 10000)
    dataout = fill(NaN, 3, 3, 10000)
       
    R = CartesianRange(Base.front(indices(A)))
           
    Threads.@threads for r in R        
           val = somefunction(A[r,:])
           dataout[r,:] = val    
     end
     return dataout
end

cartesianthread()
ERROR: MethodError: no method matching unsafe_getindex(::CartesianRange{CartesianIndex{2}}, ::Int64)
Closest candidates are:
  unsafe_getindex(::StepRangeLen{T,#s45,#s44} where #s44<:Base.TwicePrecision where #s45<:Base.TwicePrecision, ::Integer) where T at twiceprecision.jl:193
  unsafe_getindex(::StepRangeLen{T,R,S} where S where R, ::Integer) where T at range.jl:505
  unsafe_getindex(::LinSpace, ::Integer) at range.jl:510
  ...
Stacktrace:
 [1] (::##342#threadsfor_fun#26{CartesianRange{CartesianIndex{2}},Array{Float64,3},Array{Float64,3}})(::Bool) at ./threadingconstructs.jl:63
 [2] macro expansion at ./threadingconstructs.jl:71 [inlined]
 [3] cartesianthread() at ./REPL[120]:7

```

Trying without the `Threads.@threads` macro inside the `cartesianthread()` function returns the correct values.

Thanks for any help or hints about what can be done!

---

<div class="post-metadata">

**Author:** ![Balinus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/balinus/32/243_2.png) [@Balinus](https://discourse.julialang.org/u/Balinus)\
**Post date:** [June 26, 2018, 1:29pm UTC](https://discourse.julialang.org/t/parallel-computing-on-a-subset-of-an-array/11881/2 "2018-06-26T13:29:31Z")

</div>

Anyone has an idea? 😄

---

<div class="post-metadata">

**Author:** ![favba](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/favba/32/2735_2.png) [@favba](https://discourse.julialang.org/u/favba)\
**Post date:** [June 27, 2018, 6:02pm UTC](https://discourse.julialang.org/t/parallel-computing-on-a-subset-of-an-array/11881/3 "2018-06-27T18:02:47Z")

</div>

Right now `@threads` only works with indexable iterables (more precisely, iterators with `Base.unsafe_index` defined).

You could try reshaping your array:

```julia
julia> function simplethread2()
                  A = randn(3, 3, 10000)
                  dataout = fill(NaN, 3, 3, 10000)
                  B = reshape(A,(size(A,1)*size(A,2),size(A,3)))
                  dataoutB = reshape(dataout,(size(A,1)*size(A,2),size(A,3)))
                  Threads.@threads for j = 1:size(B, 1)
                      val = somefunction(B[j,:])
                      dataoutB[j,:] = val
                  end
                  return dataout
              end
simplethread2 (generic function with 1 method)

julia> @benchmark simplethread2()
BenchmarkTools.Trial:
  memory estimate: 1.68 MiB
  allocs estimate: 43
  --------------
  minimum time: 1.766 ms (0.00% GC)
  median time: 2.372 ms (0.00% GC)
  mean time: 2.399 ms (0.00% GC)
  maximum time: 9.151 ms (0.00% GC)
  --------------
  samples: 2082
  evals/sample: 1

julia> @benchmark simplethread()
BenchmarkTools.Trial:
  memory estimate: 2.46 MiB
  allocs estimate: 470
  --------------
  minimum time: 2.051 ms (0.00% GC)
  median time: 2.333 ms (0.00% GC)
  mean time: 2.539 ms (0.00% GC)
  maximum time: 6.084 ms (0.00% GC)
  --------------
  samples: 1968
  evals/sample: 1

```

A and B share the same data.  
Ps: Note that I didn’t test if both functions provide the same results.

---

<div class="post-metadata">

**Author:** ![Balinus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/balinus/32/243_2.png) [@Balinus](https://discourse.julialang.org/u/Balinus)\
**Post date:** [June 27, 2018, 8:26pm UTC](https://discourse.julialang.org/t/parallel-computing-on-a-subset-of-an-array/11881/4 "2018-06-27T20:26:03Z")

</div>

Thanks! Looks like the rehaping does not cost too much in this case. I’ll have to try if it’s still the case with regular size array I’m used too (~ 200x200x20\_000). Will keep posted!

---

<div class="post-metadata">

**Author:** ![Balinus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/balinus/32/243_2.png) [@Balinus](https://discourse.julialang.org/u/Balinus)\
**Post date:** [June 28, 2018, 3:11pm UTC](https://discourse.julialang.org/t/parallel-computing-on-a-subset-of-an-array/11881/5 "2018-06-28T15:11:40Z")

</div>

The provided solution does work, thanks!

The speedup is not as important as I expected, due to the cost of reshaping I guess vs no reshape in the non-threaded function.

---

<div class="post-metadata">

**Author:** ![favba](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/favba/32/2735_2.png) [@favba](https://discourse.julialang.org/u/favba)\
**Post date:** [June 28, 2018, 3:47pm UTC](https://discourse.julialang.org/t/parallel-computing-on-a-subset-of-an-array/11881/6 "2018-06-28T15:47:05Z")

</div>

The way your data is laid-out is not the best for the type of work you are doing (using the whole time series of one point). Is reorganizing your data really out of question? You would see major speed improvement if your data was organized as `A[t,k,j]` instead of `A[j,k,t]`.  
`A[t,j,k]` would work too.

---

<div class="post-metadata">

**Author:** ![favba](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/favba/32/2735_2.png) [@favba](https://discourse.julialang.org/u/favba)\
**Post date:** [June 28, 2018, 3:51pm UTC](https://discourse.julialang.org/t/parallel-computing-on-a-subset-of-an-array/11881/7 "2018-06-28T15:51:00Z")

</div>

You could try transposing your array before this calculation, and then transposing it back to the original shape. You’d have to benchmark it to see whether the speed gain on the calculation is worth, considering the extra time spent transposing the data.

---

<div class="post-metadata">

**Author:** ![Balinus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/balinus/32/243_2.png) [@Balinus](https://discourse.julialang.org/u/Balinus)\
**Post date:** [June 28, 2018, 4:18pm UTC](https://discourse.julialang.org/t/parallel-computing-on-a-subset-of-an-array/11881/8 "2018-06-28T16:18:13Z")

</div>

The original datasets contained in netCDF files from climate models are `[k, j, t]` (standardized practice).

Initially, we were transposing it to `[t, k, j]` to get better performance on calculations. However, the bottleneck came from loading bigger dataset such as [Daymet](https://daymet.ornl.gov/) (e.g. 1km resolution → size of `(7814 × 8075 × 365)` for a given year. The transposition (using `permutedims`) was out of the question due to memory limits. Hence, our solution was to avoid transposing at the cost of slower calculations.

I’m open to suggestions! For example, is there a function that tests if a given operation would give a memory limit error? I could simply do a `if-else` perhaps.

---

<div class="post-metadata">

**Author:** ![greg\_plowman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/greg_plowman/32/8100_2.png) [@greg\_plowman](https://discourse.julialang.org/u/greg_plowman)\
**Post date:** [June 28, 2018, 11:57pm UTC](https://discourse.julialang.org/t/parallel-computing-on-a-subset-of-an-array/11881/9 "2018-06-28T23:57:18Z")

</div>

Instead of reshaping A, you could try linearly looping over a list of cartesian indices:

```julia
    R = CartesianRange((indices(A,1), indices(A,2)))
    cartesian_inds = collect(R)

    Threads.@threads for i in linearindices(cartesian_inds)
        j = cartesian_inds[i][1]
        k = cartesian_inds[i][2]
        val = somefunction(A[j,k,:])
        dataout[j,k,:] = val
    end

```

However, the biggest gains would probably come from:

1. cache-friendly layout with time series as the first dimension.

2. re-writing `somefunction` to accept `A`, `j`, `k` as separate arguments (to avoid `A[j,k,:]` slicing), and possibly to update `val` or `dataout` in-place, rather than returning an allocated array.

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [June 29, 2018, 7:31am UTC](https://discourse.julialang.org/t/parallel-computing-on-a-subset-of-an-array/11881/10 "2018-06-29T07:31:23Z")

</div>

Showing my utter ignorance here. The documentation for reshape says “The two arrays share the same underlying data, so that setting elements of R alters the values of A and vice versa.”  
So the reshaped array uses the same memory addresses. It is just indexed differenty.  
ie. you dont use twice the memory when reshape is used?

I guess I could Google, but at a low level how does an array “know” that it has dimensions A B C?  
I guess it is held in a data structure… Algorithms + Data Structures = Programs after all!

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [June 29, 2018, 7:49am UTC](https://discourse.julialang.org/t/parallel-computing-on-a-subset-of-an-array/11881/11 "2018-06-29T07:49:47Z")

</div>

> [@johnh](#):
>
> you dont use twice the memory when reshape is used?

Precisely. Eg

```julia
julia> A = zeros(4)
4-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0

julia> B = reshape(A, 2, :)
2×2 Array{Float64,2}:
 0.0 0.0
 0.0 0.0

julia> B[1, 1] = 1.0
1.0

julia> A
4-element Array{Float64,1}:
 1.0
 0.0
 0.0
 0.0

```

> [@johnh](#):
>
> at a low level how does an array “know” that it has dimensions A B C?

This is an implementation detail which you don’t need unless you want to contribute to Julia internals, but a good starting point is [the devdocs](https://docs.julialang.org/en/latest/devdocs/object/#Memory-layout-of-Julia-Objects-1).

---

<div class="post-metadata">

**Author:** ![johnh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnh/32/3615_2.png) [@johnh](https://discourse.julialang.org/u/johnh)\
**Post date:** [June 29, 2018, 8:03am UTC](https://discourse.julialang.org/t/parallel-computing-on-a-subset-of-an-array/11881/12 "2018-06-29T08:03:47Z")

</div>

Thankyou.  
Veering wildly off topic, at one point I had expertise in BOS memory management for FORTRAN  
[http://www.desy.de/~blobel/bosman.pdf](http://www.desy.de/~blobel/bosman.pdf)

Doing terrible terrible things to the insides of innocent COMMON blocks. The horror… the horror…

---

<div class="post-metadata">

**Author:** ![Balinus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/balinus/32/243_2.png) [@Balinus](https://discourse.julialang.org/u/Balinus)\
**Post date:** [June 29, 2018, 12:50pm UTC](https://discourse.julialang.org/t/parallel-computing-on-a-subset-of-an-array/11881/13 "2018-06-29T12:50:15Z")

</div>

Thanks for the suggestion! I’ll try it and see how it goes. (can’t test it before next week though).