# For dimensions >= 2
@inline function reduce_forces!(Fij, Dijd, Pij, ::Val{D}) where {D}
    @tullio Fij[d, i] = Dijd[i, j, d] * Pij[i, j]
    return nothing
end

using Tullio, KernelAbstractions, Metal
mul(A, B, C) = @tullio C[k] = A[k] * B[k];  # run macro after loading packages

n = 2^9
A = rand(n, n)
B = rand(n, n)
C = similar(A)

# Works on CPU
mul(A, B, C)

# Throws "scalar indexing" error on Metal backend, not on CUDA though
mul(mtl(A), mtl(B), mtl(C))

julia> mul(mtl(A), mtl(B), mtl(C))
ERROR: Scalar indexing is disallowed.
Invocation of getindex resulted in scalar indexing of a GPU array.
This is typically caused by calling an iterating implementation of a method.
Such implementations *do not* execute on the GPU, but very slowly on the CPU,
and therefore should be avoided.

If you want to allow scalar iteration, use `allowscalar` or `@allowscalar`
to enable scalar iteration globally or for the operations in question.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] errorscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:151
  [3] _assertscalar(op::String, behavior::GPUArraysCore.ScalarIndexing)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:124
  [4] assertscalar(op::String)
    @ GPUArraysCore ~/.julia/packages/GPUArraysCore/aNaXo/src/GPUArraysCore.jl:112
  [5] getindex
    @ ~/.julia/packages/GPUArrays/ZRk7Q/src/host/indexing.jl:50 [inlined]
  [6] 𝒜𝒸𝓉!
    @ ~/.julia/packages/Tullio/2zyFP/src/macro.jl:1036 [inlined]
  [7] tile_halves(fun!::var"#𝒜𝒸𝓉!#1", ::Type{…}, As::Tuple{…}, Is::Tuple{…}, Js::Tuple{}, keep::Nothing, final::Bool)
    @ Tullio ~/.julia/packages/Tullio/2zyFP/src/threads.jl:139
  [8] tile_halves(fun!::var"#𝒜𝒸𝓉!#1", ::Type{…}, As::Tuple{…}, Is::Tuple{…}, Js::Tuple{}, keep::Nothing, final::Bool) (repeats 12 times)
    @ Tullio ~/.julia/packages/Tullio/2zyFP/src/threads.jl:142
  [9] tile_halves
    @ ~/.julia/packages/Tullio/2zyFP/src/threads.jl:136 [inlined]
 [10] threader
    @ ~/.julia/packages/Tullio/2zyFP/src/threads.jl:65 [inlined]
 [11] macro expansion
    @ ~/.julia/packages/Tullio/2zyFP/src/macro.jl:1004 [inlined]
 [12] mul(A::MtlMatrix{…}, B::MtlMatrix{…}, C::MtlMatrix{…})
    @ Main ./REPL[2]:1
 [13] top-level scope
    @ REPL[8]:1
 [14] top-level scope
    @ ~/.julia/packages/Metal/UahNC/src/initialization.jl:80
Some type information was truncated. Use `show(err)` to see complete types.