Feature Request: `mapreduce` over cartesian indices

`mapreduce` over `CartesianIndices` now fallbacks to the one in Base (as it should be).
I find this is also a useful interface for CUDA programming, is there some existing approach to realize similar functionality?
If not, where can I start to support this kind of interface?

Here is an minimal example of calculating means distance
```julia console
julia> using CuArrays
julia> using CUDAnative
julia> function expect(f, xs::CuVector, ys::CuVector)
M = length(xs)
N = length(ys)
ci = CartesianIndices((M, N))
mapreduce(inds->f(xs[inds[1]], ys[inds[2]]), +, ci)/M/N
end

julia> @benchmark expect((x,y)->abs2(x-y), randn(100)|>cu, randn(100)|>cu)
BenchmarkTools.Trial:
memory estimate: 3.97 MiB
allocs estimate: 100033
--------------
minimum time: 182.878 ms (0.00% GC)
median time: 185.717 ms (0.00% GC)
mean time: 190.299 ms (1.61% GC)
maximum time: 246.243 ms (24.11% GC)
--------------
samples: 27
evals/sample: 1
```