# Performance of using CartesianIndices with MArray

**URL:** <https://discourse.julialang.org/t/performance-of-using-cartesianindices-with-marray/77167>\
**Category:** Performance\
**Tags:** staticarrays, cartesianindices\
**Created:** [February 27, 2022, 10:40pm UTC](https://discourse.julialang.org/t/performance-of-using-cartesianindices-with-marray/77167 "2022-02-27T22:40:07Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![nchisholm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nchisholm/32/37414_2.png) [@nchisholm](https://discourse.julialang.org/u/nchisholm)\
**Post date:** [February 27, 2022, 10:40pm UTC](https://discourse.julialang.org/t/performance-of-using-cartesianindices-with-marray/77167/1 "2022-02-27T22:40:07Z")

</div>

Multidimensional algorithms built on `CartesianIndices` seem to have significant performance overhead when applied to `StaticArrays.MArray`​s, but not regular `Array`​s. Writing several functions with nested loops for different dimensions (or leveraging the macros in `Base.Cartesian`) recovers performance, but I am wondering why the more elegant `CartesianIndices` solution is slower when used with `MArray`.

To demonstrate the problem, consider an `N`-dimensional algorithm that takes a function `f(i_1, ⋯, i_N)` of the `N` indices and uses it to fill the elements of an input array `arr`.

```julia
function fillfn_cartesianindices!(arr::AbstractArray, f::F) where F # force specialization
    @inbounds for I ∈ CartesianIndices(arr)
        arr[I] = f(I)
    end
    return arr
end

```

Here is a similar `@generated` function using the `Base.Cartesian` macros.

```julia
using Base.Cartesian

@generated function fillfn_basecartesian!(arr::AbstractArray{<:Any,N}, f::F) where {N,F}
    # Generate an expression of nested for loops:
    # for i_N ∈ axes(arr, N), ..., i_1 ∈ axes(arr, 1)
    # arr[i_1, ..., i_N] = f(i_1, ..., i_N)
    # end
    quote
        @inbounds @nloops $N i arr begin
            (@nref $N arr i) = (@ncall $N f i)
        end
        return arr
    end
end

```

Finally, here is function of indices to feed to either of the above functions (for the four dimensional case).

```julia
    indexfunc4d(i, j, k, l) = (i==k) * (j==l) - (i==l) * (j==k)
    @inline indexfunc4d(I::CartesianIndex{4}) = indexfunc4d(Tuple(I)...)

```

Benchmarking on a 3×3×3×3 `Array` shows that using `CartesianIndices` has no performance overhead.

```julia
using BenchmarkTools
using StaticArrays

let dims = (3,3,3,3)
    T = eltype(indexfunc4d(map(_ -> 1, dims)...))

    println("Using `CartesianIndices`")
    out1 = @btime(
        fillfn_cartesianindices!(A, $indexfunc4d),
        setup=(A = Array{$T}(undef, $dims))
    )
    println("Using `Base.Cartesian`")
    out2 = @btime(
        fillfn_basecartesian!(A, $indexfunc4d),
        setup=(A = Array{$T}(undef, $dims))
    )
    @assert out1 == out2
end;

```

```
Using `CartesianIndices`
  107.183 ns (0 allocations: 0 bytes)
Using `Base.Cartesian`
  114.619 ns (0 allocations: 0 bytes)

```

A similar benchmark where `arr isa MArray` shows that `fillfn_basecartesian!` 7x faster than `fillfn_cartesianindices!`.

```julia
    let SDims = NTuple{4,3}
        dims = fieldtypes(SDims)
        T = eltype(indexfunc4d(map(_ -> 1, dims)...)) # Element type
    
        println("Using `CartesianIndices`")
        out1 = @btime(
            fillfn_cartesianindices!(A, $indexfunc4d),
            setup=(A = MArray{$SDims,$T}(undef))
        )
        println("Using `Base.Cartesian`")
        out2 = @btime(
            fillfn_basecartesian!(A, $indexfunc4d),
            setup=(A = MArray{$SDims,$T}(undef))
        )
    end;

```

```
Using `CartesianIndices`
  94.878 ns (0 allocations: 0 bytes)
Using `Base.Cartesian`
  12.877 ns (0 allocations: 0 bytes)

```

Note that using `StaticArrays.sacollect` with `CartesianIndices` seems to be efficient.

```julia
@btime StaticArrays.sacollect(SArray{NTuple{4,3}},
                              $indexfunc4d(I) for I ∈ CartesianIndices(ntuple(_ -> SOneTo(3), Val(4))));
@btime StaticArrays.sacollect(MArray{NTuple{4,3}},
                              $indexfunc4d(I) for I ∈ CartesianIndices(ntuple(_ -> SOneTo(3), Val(4))));

```

```
9.858 ns (0 allocations: 0 bytes)
87.575 ns (1 allocation: 672 bytes)

```

However, AFAIK, `sacollect` does not perform a `for` loop over `CartesianIndices` at runtime but rather constructs an expression during compilation to construct the array “all-at-once”, as required by `SArray`.

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [March 1, 2022, 1:24pm UTC](https://discourse.julialang.org/t/performance-of-using-cartesianindices-with-marray/77167/2 "2022-03-01T13:24:48Z")

</div>

It looks to me like [`MArray`](https://github.com/JuliaArrays/StaticArrays.jl/blob/c1620f0595b02bab67e055294850b19dffafa549/src/MArray.jl) doesn’t have special support for `CartesianIndices` and I don’t know if this is even possible. Using `@code_llvm` I saw that

```julia
out2 = @btime(
            fillfn_basecartesian!(A, $indexfunc4d),
            setup=(A = MArray{$SDims,$T}(undef))
        )

```

compiles down to a series of `call void @llvm.memset.p0i8.i64`.

---

<div class="post-metadata">

**Author:** ![nchisholm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nchisholm/32/37414_2.png) [@nchisholm](https://discourse.julialang.org/u/nchisholm)\
**Post date:** [March 2, 2022, 3:29pm UTC](https://discourse.julialang.org/t/performance-of-using-cartesianindices-with-marray/77167/3 "2022-03-02T15:29:10Z")

</div>

Thanks for your reply and to looking into this. Perhaps this is one of those cases why `Base.Cartesian` is kept around. I opened an issue with `StaticArrays` in case others run into the same problem or know a way to fix the issue.

[https://github.com/JuliaArrays/StaticArrays.jl/issues/1010](https://github.com/JuliaArrays/StaticArrays.jl/issues/1010).