# Access arrays in memory order, along columns

**URL:** <https://discourse.julialang.org/t/access-arrays-in-memory-order-along-columns/77073>\
**Category:** New to Julia\
**Tags:** array, memory-allocation\
**Created:** [February 25, 2022, 2:26pm UTC](https://discourse.julialang.org/t/access-arrays-in-memory-order-along-columns/77073 "2022-02-25T14:26:47Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![fredrikpaues](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikpaues/32/34080_2.png) [@fredrikpaues](https://discourse.julialang.org/u/fredrikpaues)\
**Post date:** [February 25, 2022, 2:26pm UTC](https://discourse.julialang.org/t/access-arrays-in-memory-order-along-columns/77073/1 "2022-02-25T14:26:47Z")

</div>

I’m trying to wrap my head around how I best should iterate over arrays. I’m coming to Julia after primarily using MATLAB and as both MATLAB and Julia are column major arrays should be accessed in the exact same order. My understanding from MATLAB is that I below have ordered the loops for optimal speed.

```nohighlight
a = 1:25*50*75;
a = reshape(a, 25, 50, 75);
for k = 1:size(a, 3)
    for j = 1:size(a, 2)
        for i = 1:size(a, 1)
            a(i, j, k) = i + j + k;
        end
    end
end

```

Should the loops be ordered the same way in Julia?

```julia
a = collect(1:25*50*75)
a = reshape(a, (25, 50, 75))
for k = 1:size(a, 3)
    for j = 1:size(a, 2)
        for i = 1:size(a, 1)
            a[i, j, k] = i + j + k
        end
    end
end

```

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [February 25, 2022, 2:28pm UTC](https://discourse.julialang.org/t/access-arrays-in-memory-order-along-columns/77073/2 "2022-02-25T14:28:33Z")

</div>

Yes.

ps: you can use `for i in axes(a,1)` for example, such that the loop is independent on how the array is indexed (for example, not starting at 1 - not common, but possible).

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [February 25, 2022, 3:37pm UTC](https://discourse.julialang.org/t/access-arrays-in-memory-order-along-columns/77073/3 "2022-02-25T15:37:16Z")

</div>

You can also use `eachindex` or `CartesianIndices` or `pairs` for this sort of thing, depending on the calculations you want to perform. They will choose the right order for you. You shouldn’t have to think about this too often.

In most cases, do

```julia
for i in eachindex(a)
    a[i] = ...
end

```

If you need the actual indices, like in your example,

```julia
for ind in CartesianIndices(a)
    a[ind] = sum(Tuple(ind))
end

```

These work for arrays of arbitrary dimensions.

There are several indexing functions that make working with arrays much nicer than in Matlab.

---

<div class="post-metadata">

**Author:** ![Seif\_Shebl](https://avatars.discourse-cdn.com/v4/letter/s/eada6e/32.png) [@Seif\_Shebl](https://discourse.julialang.org/u/Seif_Shebl)\
**Post date:** [February 26, 2022, 12:20pm UTC](https://discourse.julialang.org/t/access-arrays-in-memory-order-along-columns/77073/4 "2022-02-26T12:20:31Z")

</div>

Besides what others have suggested, in Julia you can write this in a single line without worrying about column or row order and get the best perormance. This is possible using packages like `Tullio`:

```julia
using Tullio
a = collect(1:25*50*75)
a = reshape(a, 25, 50, 75)
@tullio a[i, j, k] = i + j + k

```

---

<div class="post-metadata">

**Author:** ![photor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/photor/32/14343_2.png) [@photor](https://discourse.julialang.org/u/photor)\
**Post date:** [February 26, 2022, 1:23pm UTC](https://discourse.julialang.org/t/access-arrays-in-memory-order-along-columns/77073/5 "2022-02-26T13:23:50Z")

</div>

This Tullio version takes only 26.500 μs on my machine, while the multi-threaded version

```julia
Threads.@threads for k in axes(a, 3)
    for j in axes(a, 2)
        for i in axes(a, 1)
            a[i, j, k] = i + j + k
        end
    end
end

```

takes 2.954 ms. How can Tullio achieve that speed?

---

<div class="post-metadata">

**Author:** ![Seif\_Shebl](https://avatars.discourse-cdn.com/v4/letter/s/eada6e/32.png) [@Seif\_Shebl](https://discourse.julialang.org/u/Seif_Shebl)\
**Post date:** [February 26, 2022, 1:41pm UTC](https://discourse.julialang.org/t/access-arrays-in-memory-order-along-columns/77073/6 "2022-02-26T13:41:21Z")

</div>

Tullio doesn’t do black magic, it should be equivalent to properly ordered loops, be sure you added `-t auto` as compiler option. Multihtreading should not do much here because this problem is essentially memory-bound rather than cpu-bound. I observe only 2.5x speedup on a 4-core/8-threads machine.

```julia
using Tullio

function do_loops!(a) 
    Threads.@threads for k in axes(a, 3)
        for j in axes(a, 2)
            for i in axes(a, 1)
                a[i, j, k] = i + j + k
            end
        end
    end
end

a = collect(1:25*50*75)
a = reshape(a, (25, 50, 75))
@btime do_loops!($a) 
  # 24.900 μs (0 allocations: 0 bytes) - Serial
  # 10.600 μs (48 allocations: 5.08 KiB) - Parallel
@btime @tullio $a[i, j, k] = i + j + k
  # 25.800 μs (0 allocations: 0 bytes)

```

---

<div class="post-metadata">

**Author:** ![photor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/photor/32/14343_2.png) [@photor](https://discourse.julialang.org/u/photor)\
**Post date:** [February 26, 2022, 1:46pm UTC](https://discourse.julialang.org/t/access-arrays-in-memory-order-along-columns/77073/7 "2022-02-26T13:46:09Z")

</div>

I see. I forgot to put the for loops in a function.

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [February 26, 2022, 2:31pm UTC](https://discourse.julialang.org/t/access-arrays-in-memory-order-along-columns/77073/8 "2022-02-26T14:31:48Z")

</div>

`@tturbo` from `LoopVectorization` sometimes does magic in situations like these, because it has a allocation-free and fast threading mechanism (which can be found in `Polyester`).

`@inbounds` may also help, although the use `axes` guarantees that already, and thus it may not make any difference here.

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [February 26, 2022, 5:56pm UTC](https://discourse.julialang.org/t/access-arrays-in-memory-order-along-columns/77073/9 "2022-02-26T17:56:48Z")

</div>

> [@fredrikpaues](#):
>
> ```julia
> a = collect(1:25*50*75)
> a = reshape(a, (25, 50, 75))
> 
> ```

I would also suggest that since you aren’t using the original values stored in `a`, you change the above Matlabish code to a more Julian form:

```julia
a = zeros(Int, 25, 50, 75)

```

or

```julia
a = Array{Int, 3}(undef, 25, 50, 75)

```

The first form is more convenient to write and initializes `a` to all zeros. The second form may be a bit faster and leaves the contents of `a` unspecified.

---

<div class="post-metadata">

**Author:** ![Klaas\_Pauly](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaas_pauly/32/14819_2.png) [@Klaas\_Pauly](https://discourse.julialang.org/u/Klaas_Pauly)\
**Post date:** [March 1, 2022, 11:42am UTC](https://discourse.julialang.org/t/access-arrays-in-memory-order-along-columns/77073/10 "2022-03-01T11:42:14Z")

</div>

> Multithreading should not do much here

The strange thing is that when I run your code on my computer, I don’t get any improvement from the parallel version with `Threads.@threads` - in fact, sometimes worse - but I do get a significant improvement `using Distributed` `@distributed`. On my 8 core, 16 threads Windows laptop, Tullio is slowest, with time for `@distributed` \< `@tturbo` \< `@tullio`, each with about a factor 10 difference.

---

<div class="post-metadata">

**Author:** ![Nathan\_Boyer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nathan_boyer/32/14825_2.png) [@Nathan\_Boyer](https://discourse.julialang.org/u/Nathan_Boyer)\
**Post date:** [March 1, 2022, 4:10pm UTC](https://discourse.julialang.org/t/access-arrays-in-memory-order-along-columns/77073/11 "2022-03-01T16:10:02Z")

</div>

If you need `i+1`, `j-1`, etc. on the right-hand side, is this the suggested method?

```julia
for ind in CartesianIndices(a)
    (i, j, k) = Tuple(ind)
    a[ind] = (a[i+1, j, k] + a[i, j-1, k])/2
end

```

I don’t see `pairs` discussed in the array indexing [documentation](https://docs.julialang.org/en/v1/manual/arrays/#Iteration) or [blog post](https://julialang.org/blog/2016/02/iteration/). I interpret the `pairs` docstring to mean it should be used if you need indices and values at the same time and that `pairs` is always preferred over `enumerate`. Is that correct?

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [March 1, 2022, 4:38pm UTC](https://discourse.julialang.org/t/access-arrays-in-memory-order-along-columns/77073/12 "2022-03-01T16:38:39Z")

</div>

This is going to fail, since `i-1` and `i+1` must necessarily go out of bounds at the edges. You’ll have to figure out some way to peel of the edge values.

> [@Nathan\_Boyer](#):
>
> `pairs` is always preferred over `enumerate` . Is that correct?

`pairs` will give you the value and the index/indices. Enumerate will just count from one, so it depends on what you need. If you are looking for indices, use `pairs`.

---

<div class="post-metadata">

**Author:** ![fredrikpaues](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fredrikpaues/32/34080_2.png) [@fredrikpaues](https://discourse.julialang.org/u/fredrikpaues)\
**Post date:** [March 20, 2022, 5:58am UTC](https://discourse.julialang.org/t/access-arrays-in-memory-order-along-columns/77073/13 "2022-03-20T05:58:56Z")

</div>

`eachindex` and `CartesianIndices` are great, but often I have some computations that are unique only to one or two dimensions. And, so as not to repeat the computations on each iteration of the remaining dimensions, I need to explicitly loop over the array myself.

---

<div class="post-metadata">

**Author:** ![photor](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/photor/32/14343_2.png) [@photor](https://discourse.julialang.org/u/photor)\
**Post date:** [March 20, 2022, 8:45am UTC](https://discourse.julialang.org/t/access-arrays-in-memory-order-along-columns/77073/14 "2022-03-20T08:45:22Z")

</div>

I am new to Distributed, so could you give an example code about the problem in the OP? Thanks.