# Map Performance with CuArrays
**URL:** https://discourse.julialang.org/t/map-performance-with-cuarrays/33497
**Category:** GPU
**Tags:** question, fftw, cuda, broadcast
**Created:** [January 17, 2020, 5:54pm UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497 "2020-01-17T17:54:50Z")
**Posts on this page:** 16
**Page:** 1
### Author:  [@bmit](https://discourse.julialang.org/u/bmit)
#### Post date: [January 17, 2020, 5:54pm UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497/1 "2020-01-17T17:54:50Z")
I’m trying to benchmark different methods of calling the same in-place function on subsets of a large data block using CuArrays. I’m using fft as an example function because I can baseline against planning an fft over one dimension of a large array. Using this as a minimum working example, I created a testbench with the following methods:
1. Create an list of 1D CuArrays. map f(x) to each element of the list
2. Create a 2D CuArray. mapslices to f(x) for each column in array.
3. Baseline taking f(x) over the full 2D array
```julia
using FFTW
using BenchmarkTools
using CuArrays
CuArrays.allowscalar(false)
N = 100000
fftsz = 1024
# list of random vectors
A1_d = [CuArray(rand(ComplexF32, fftsz)) for i in 1:N]
# random 2D array
A2_d = CuArray(rand(ComplexF32, N, fftsz))
p1 = plan_fft!(A1_d[1])
p2 = plan_fft!(A2_d, 1)
function onefft!(data, plan)
plan * data
end
#1
@btime CuArrays.@sync map(x -> onefft!(x, p1), $A1_d)
#2
@btime CuArrays.@sync mapslices(x -> onefft!(x, p1), $A2_d, dims=2)
#3
@btime CuArrays.@sync onefft!($A2_d, p2)
```
The results are:
```julia
629.991 ms (200011 allocations: 3.81 MiB)
3.618 s (13198542 allocations: 440.96 MiB)
20.603 ms (10 allocations: 176 bytes)
```
I would think the approaches should be somewhat identical in terms of performance. However, the strangest things I notice are:
- Why do (1) and (2) allocate memory? All the data should already be on the GPU and all operations should be in-place
- Why are (1) and (2) so different from one another? Even if there is some penalty to map, it seems like (1) and (2) should be identical.
---
### Author:  [@maleadt](https://discourse.julialang.org/u/maleadt)
#### Post date: [January 20, 2020, 7:39am UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497/2 "2020-01-20T07:39:57Z")
> [@bmit](#):
>
> Why do (1) and (2) allocate memory? All the data should already be on the GPU and all operations should be in-place
These are CPU allocations; use `CuArrays.@time` to see if there are GPU allocations (likely not). There’s probably some inefficiency in the CuFFT wrappers leading to a couple of small CPU allocations per FFT here.
> [@bmit](#):
>
> Why are (1) and (2) so different from one another? Even if there is some penalty to map, it seems like (1) and (2) should be identical.
It looks like `mapslices` is badly implemented. I’m not surprised; we don’t have any GPU-optimized implementation in CuArrays/GPUArrays so this falls back to the Base implementation which surprisingly doesn’t use scalar iteration, but is still [know to be slow](https://github.com/JuliaLang/julia/issues/26868).
Any about the difference between (1) and (3), it’s always recommended to reduce the number of distinct calls into the GPU so that you can saturate the device with as much work as possible. There might also be some synchronizing behavior in how the wrappers are implemented (e.g., copying back some memory to inspect whether the operation succeeded), serializing the GPU at each individual FFT. So prefer batched interfaces whenever possible.
---
### Author:  [@bmit](https://discourse.julialang.org/u/bmit)
#### Post date: [January 23, 2020, 2:04pm UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497/3 "2020-01-23T14:04:15Z")
Thank you for the help. I’ll try and avoid `mapslices` for now.
> [@maleadt](#):
>
> Any about the difference between (1) and (3), it’s always recommended to reduce the number of distinct calls into the GPU so that you can saturate the device with as much work as possible. There might also be some synchronizing behavior in how the wrappers are implemented (e.g., copying back some memory to inspect whether the operation succeeded), serializing the GPU at each individual FFT. So prefer batched interfaces whenever possible.
This is the core of the question I was hoping to get some community advice on.
I used an FFT example here because I knew that I could explicitly create a batch example as a benchmark (3) with near-optimal performance. In a more general case, it seems to me like the batch approach is somewhat divergent from the Julia design philosophy of writing many simple functions. I’ll try and explain what I mean.
If I follow the [Cuda.jl Introduction](https://juliagapu.gitlab.io/CUDA.jl/tutorials/introduction/), I get extremely good performance with broadcasting simple functions onto the GPU (ex. writing a scalar addition function and broadcasting it to add N-dimensional arrays). From my limited CUDA knowledge, I think this means these simple atomic actions are being interpreted as “GPU kernels” and they distributed across the GPU in parallel. In other words this broadcast is a single “distinct call into the GPU.”
This doesn’t seem to scale up to slightly more complicated functions like the example I post here (1) (an arbitrary 1D function that can be broadcast across an N-dimensional array). From the Julia user perspective, it isn’t clear why the `map` call I use isn’t a single distinct call into the GPU the way the simple `broadcast` call is in the Cuda.jl example.
From a Julia perspective I’d rather write a general function `foo`
```julia
map(x -> foo(x, args...), bunch_of_datasets)
```
instead of equivalent macro-function function `bar`,
```julia
bar(bunch_of_datasets, args...)
```
Is this fundamentally a GPU and/or CUDA limitation? Is there anything we can do from the Julia user side to help inform when map/broadcast benefit performance and when it hurts? Are there design guides or principles from the Julia/GPU perspective that we can use? As you can tell, I’m still pretty new to all this.
---
### Author:  [@maleadt](https://discourse.julialang.org/u/maleadt)
#### Post date: [January 23, 2020, 2:11pm UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497/4 "2020-01-23T14:11:44Z")
> [@bmit](#):
>
> From a Julia perspective I’d rather write a general function `foo`
>
> ```julia
> map(x -> foo(x, args...), bunch_of_datasets)
>
> ```
But you were mapping across a CPU collection of GPU arrays, so isn’t it expected that this will result into multiple separate calls to the GPU? If you map across a GPU dataset, you _can_ use arbitrary functions and not be constrained to simple atomic actions as you mention. Furthermore, with broadcast fusion and using “dots” you can trivially create new and complex kernels that will be executed on the GPU in a single step. If you need more flexibility, you need to resort to writing your own kernels (which can be arbitrarily complex) or use something like GPUifyLoops for the vendor-agnostic equivalent.
---
### Author:  [@bmit](https://discourse.julialang.org/u/bmit)
#### Post date: [January 24, 2020, 3:07pm UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497/5 "2020-01-24T15:07:59Z")
> [@"maleadt](#):
>
> If you map across a GPU dataset, you _can_ use arbitrary functions and not be constrained to simple atomic actions as you mention.
This is exactly what I wanted to do originally, but I ran into technical difficulties. I’m probably missing something simple. If I create a multidimensional array (2), then it seems like I need to use `mapslices`, which is broken.
For the first approach I originally tried the following:
```julia
A1_d = CuArray([rand(ComplexF32, fftsz) for i in 1:N])
```
but got
`ERROR: type does not have a fixed size`
I also tried
`A1_d = CuArray([CuArray(rand(ComplexF32, fftsz)) for i in 1:N])`
but it crashed the REPL.
Do you have an approach in mind to fix (1) to be a GPU-only data structure I can map across?
---
### Author:  [@maleadt](https://discourse.julialang.org/u/maleadt)
#### Post date: [January 28, 2020, 7:43am UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497/6 "2020-01-28T07:43:02Z")
> [@bmit](#):
>
> If I create a multidimensional array (2), then it seems like I need to use `mapslices` , which is broken.
You are not creating a multidimensional array here, you are creating an array of arrays. That is not supported by CuArrays, and you need to create a 2D array instead.
---
### Author:  [@bmit](https://discourse.julialang.org/u/bmit)
#### Post date: [January 29, 2020, 1:40pm UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497/7 "2020-01-29T13:40:28Z")
I think my tests 2) and 3) create multidimensional arrays on the GPU. My test 2) didn’t work because of the issues you pointed out with mapslices. Is there another way to map or broadcast a function across rows, columns, etc of a multidimensional array that would make this work?
---
### Author:  [@maleadt](https://discourse.julialang.org/u/maleadt)
#### Post date: [January 29, 2020, 1:56pm UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497/8 "2020-01-29T13:56:37Z")
What kind of function do you want to map? `fft` calls into CUFFT, and functions like that need an API to explicitly support processing rows/columns of a dataset (either a batched interface like this, or a `dims` keyword as with `sum` or `reduce`). It would be easier to help if you have a more concrete use case, as the one you posted here (`fft`) is solved through use of its batched API. I can imagine a `maprows` (or similar) function, but that would be of limited use on the GPU because you would then be performing a coarse operation to process all elements of the dimension using a single thread. Typically, you need specifically-designed parallel kernels for workloads like that (e.g. a scan / prefix sum).
---
### Author:  [@bmit](https://discourse.julialang.org/u/bmit)
#### Post date: [January 29, 2020, 10:39pm UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497/9 "2020-01-29T22:39:51Z")
> If you map across a GPU dataset, you _can_ use arbitrary functions and not be constrained to simple atomic actions as you mention. Furthermore, with broadcast fusion and using “dots” you can trivially create new and complex kernels that will be executed on the GPU in a single step.
> I can imagine a `maprows` (or similar) function, but that would be of limited use on the GPU because you would then be performing a coarse operation to process all elements of the dimension using a single thread. Typically, you need specifically-designed parallel kernels for workloads like that (e.g. a scan / prefix sum).
I guess my confusion comes from the fact that these two comments seem a little contradictory. Your `maprows` idea is what I was trying to do with `mapslices`, but if that is executed serially then it isn’t helpful, as you say. If I need to design specific parallel kernels, rather than mapping/broadcasting arbitrary functions, then I must be missing what Julia adds vs CUDA.
Here’s the CPU code for something that is a little more representative:
```julia
using FFTW
using LinearAlgebra
N = 1000
M = 1024
L = 1024
# many datasets
datasets = [rand(M, L) for i in 1:N]
p = plan_fft(datasets[1], 1)
function process_dataset(A, p)
# the specific code here shouldn't be important
# some signal processing
B = p * A
# some linear algebra
B = B * Diagonal(B) * B'
# some custom function
B = some_custom_function.(A)
return B
end
function some_custom_function(A)
return sin(exp(A)) + 5
end
function combine_datasets(A)
B = reduce(+, A)
return B
end
# wrap static arguments for broadcasting
f(x) = process_dataset(x, p)
# independent processing on each dataset
data_tmp = f.(datasets)
# processing across datasets
data_out = combine_datasets(data_tmp)
```
Am I correct in concluding that I can’t broadcast on the GPU like I can on the CPU, at least not efficiently. Based on some of what has been suggested here, I would need to write a custom `process_many_datasets` that can batch process datasets. To me this seems more “MATLAB-like” and less “Julia-like.”
Do we understand where the limitation is (hardware, software, compiler)?
Am I thinking about this wrong?
---
### Author:  [@maleadt](https://discourse.julialang.org/u/maleadt)
#### Post date: [January 30, 2020, 7:25am UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497/10 "2020-01-30T07:25:02Z")
OK, that example clarifies a lot. We don’t efficiently support what you’re doing there: broadcasting a function that does its own broadcasting over the inner 1D dataset. If that inner broadcast is coarse enough, i.e. if it does enough work to saturate the GPU, a mapslices-like approach might work (assuming we optimize the implementation of that API). But generally it will be more efficient to create a function that you can broadcast over the entire 2D dataset (not an array of arrays) in one go, e.g., like the batched fft interface in your first post. Some of your operations are element-wise anyway, and the final `reduce` you can perform in a batched manner with the `dims` keyword.
> [@bmit](#):
>
> Am I correct in concluding that I can’t broadcast on the GPU like I can on the CPU, at least not efficiently.
Don’t expect magic from the GPU. We implement `broadcast` to compile to a single kernel; you’re broadcasting a function that does an FFT and further broadcasts; that’s never going to work as a single kernel on the GPU. If you stick to this approach, optimizing `mapslices` would be the way to go. Else, you’d need to think about restructuring your operations to be dataset-wide, which would have advantages on the CPU too.
> [@bmit](#):
>
> I must be missing what Julia adds vs CUDA.
Except for the fact you’re writing high-level Julia code instead of CUDA C, you mean? 🙂
---
### Author:  [@bmit](https://discourse.julialang.org/u/bmit)
#### Post date: [January 30, 2020, 3:43pm UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497/11 "2020-01-30T15:43:01Z")
This explanation was extremely helpful.
To summarize:
The issue was that a single call to a function `process_dataset` does not saturate the GPU, so I was trying to parallelize that call over many datasets using broadcast (or map). GPU parallel processing does not seem to be that generalized. The solution is to follow the common `dims=` style interface for a function that can explicitly processes many datasets in a single batch.
Also, `mapslices` [needs to be optimized](https://github.com/JuliaLang/julia/issues/26868) to be used effectively for even serial processing on the GPU - see example (2) above.
> [@maleadt](#):
>
> Except for the fact you’re writing high-level Julia code instead of CUDA C, you mean? 🙂
👍
---
### Author:  [@Ratingulate](https://discourse.julialang.org/u/Ratingulate)
#### Post date: [January 30, 2020, 5:11pm UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497/12 "2020-01-30T17:11:15Z")
Are you looking for something like this: [https://github.com/FluxML/Hydra.jl](https://github.com/FluxML/Hydra.jl) ?
---
### Author:  [@bjarthur](https://discourse.julialang.org/u/bjarthur)
#### Post date: [December 26, 2020, 3:57pm UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497/13 "2020-12-26T15:57:29Z")
re. options #1 and #2 with map and mapslices: shouldn’t those kernels run asynchronously on the GPU? if so, i would expect some inefficiency compared to a monolithic kernel which processes the entire data in a single batch, but not a \>30-fold slow down.
my specific problem is as follows. i have a 44x44 matrix of floats (call it P) that i need to multiply by a 44x1 vector (=R). i have 20,000 such products to compute. storing them as a vector of matrices, `map((x,y)->x*y, P, R)` takes 20,000 times longer than `P[1]*R[1]`. so seemingly no asynchronous parallelism. i tried pre-allocating the output and using `map!`, but no improvement. i tried `asyncmap` and got “GPU program failed to execute”. not sure how i could re-factor to multiply two 3D matrices in that way.
thanks for any advice.
---
### Author:  [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)
#### Post date: [December 26, 2020, 4:10pm UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497/14 "2020-12-26T16:10:20Z")
You probably want to call CUDA’s built-in [batched matrix multiplication](https://developer.nvidia.com/blog/cublas-strided-batched-matrix-multiply/) routine. You can get at the batched routines through NNLib.jl, specifically [`NNlib.batched_mul!`](https://fluxml.ai/Flux.jl/stable/models/nnlib/).
---
### Author:  [@bjarthur](https://discourse.julialang.org/u/bjarthur)
#### Post date: [December 28, 2020, 2:42pm UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497/15 "2020-12-28T14:42:42Z")
thanks for the pointer to NNlib!
still curious though that kernel calls do not seem to be asynchronous.
---
### Author:  [@maleadt](https://discourse.julialang.org/u/maleadt)
#### Post date: [January 4, 2021, 7:29am UTC](https://discourse.julialang.org/t/map-performance-with-cuarrays/33497/16 "2021-01-04T07:29:48Z")
> [@bjarthur](#):
>
> calls do not seem to be asynchronous.
They definitely are, why do you assume they aren’t? Kernel calls aren’t free, of course, so if your kernel execution time is shorter than the launch time you’ll just scale the time it takes to launch.