# Reordering for memory contiguity/memory access

**URL:** <https://discourse.julialang.org/t/reordering-for-memory-contiguity-memory-access/43040>\
**Category:** Numerics\
**Tags:** question\
**Created:** [July 14, 2020, 12:42pm UTC](https://discourse.julialang.org/t/reordering-for-memory-contiguity-memory-access/43040 "2020-07-14T12:42:21Z")\
**Posts on this page:** 4\
**Page:** 1

<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:** [July 14, 2020, 12:42pm UTC](https://discourse.julialang.org/t/reordering-for-memory-contiguity-memory-access/43040/1 "2020-07-14T12:42:21Z")

</div>

I am simulating a “mass” that obeys a law of motion. The algorithm I use to discretize to a finite grid has the property that each cell is distributed to `K` others after being multiplied by something else. Formally, each x z\_k is added to some y\_{p\_k} according to “patterns” `p` (ie collections of indices).

There is an MWE below using `SVector`, but that’s just for the purposes of an example, generally my type supports the relevant methods for `*`, `zero`, `+`. `K` is known statically.

This is not unlike a generalization of a sparse matrix multiplication, but with a very special index structure. I am looking for suggestions to improve it; specifically, is there a simple way to calculate a permutation of indexes in `xs`/`ys` so that the memory access for `ys` is improved?

```julia
using StatsBase, StaticArrays, BenchmarkTools

N = 100
K = 4
M = 3

patterns = [SVector{K}(sample(1:N, 4; replace = false)...) for _ in 1:N]
xs = rand(SVector{M,Float64}, N)
ys = similar(xs)
zs = rand(SVector{K,Float64})

function combine_forward!(ys, xs::Vector{T},
                          patterns::Vector{SVector{K,Int}},
                          zs::SVector{K}) where {T,K}
    N = length(xs)
    fill!(ys, zero(T))
    @inbounds for (x, pattern) in zip(xs, patterns)
        for (j, z) in zip(pattern, zs)
            ys[j] += x * z
        end
    end
    ys
end

@btime combine_forward!(ys, xs, patterns, zs);

```

---

<div class="post-metadata">

**Author:** ![jlchan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlchan/32/10958_2.png) [@jlchan](https://discourse.julialang.org/u/jlchan)\
**Post date:** [July 14, 2020, 2:50pm UTC](https://discourse.julialang.org/t/reordering-for-memory-contiguity-memory-access/43040/2 "2020-07-14T14:50:57Z")

</div>

What about [CuthillMcKee ordering](https://github.com/rleegates/CuthillMcKee.jl)?

If you ignore the scaling and represent the accumulation as a sparse matvec, Cuthill McKee gives a bandwidth-reducing permutation of the sparse matrix. It might improve data locality (though this isn’t what it’s designed for).

---

<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:** [July 14, 2020, 3:44pm UTC](https://discourse.julialang.org/t/reordering-for-memory-contiguity-memory-access/43040/3 "2020-07-14T15:44:19Z")

</div>

Thanks, that could improve things (also, I did not realize that the proper term is _data locality_). I found [this topic](https://discourse.julialang.org/t/implementing-cuthill-mckee/9268/) which mentions some packages implementing Cluthill-McKee, but I am unsure which one of those would make it easiest to extract the permutation itself.

Also, I am wondering if instead of specifying targets in `patterns`, I should instead specify “sources” — ie effectively the sparsity pattern of the transpose. That would allow me to parallelize the outer loop without worrying about concurrent writes.

---

<div class="post-metadata">

**Author:** ![jlchan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlchan/32/10958_2.png) [@jlchan](https://discourse.julialang.org/u/jlchan)\
**Post date:** [July 14, 2020, 3:54pm UTC](https://discourse.julialang.org/t/reordering-for-memory-contiguity-memory-access/43040/4 "2020-07-14T15:54:39Z")

</div>

The “source” pattern would certainly another option. Without knowing more about the sparsity of the access pattern, it’s hard to say which would be more effective for your problem.

I just experimented briefly with CuthillMcKee.jl, but the permutation is pretty easy to extract (and is only a small precomputation if the pattern is fixed in time). .