[^0]: https://github.com/JuliaGPU/KernelAbstractions.jl/blob/maste...
https://juliagpu.gitlab.io/KernelAbstractions.jl/examples/na...
and you get ~60% of the best known CUDA implementation on that test:
https://juliagpu.gitlab.io/KernelAbstractions.jl/examples/pe...
Not perfect but not awful.
What I showed was the the essentially no-knowledge pure Julia code gets to ~60% of the optimized CUDA kernel that libraries tend to wrap. Take that as you will for whether that means that means that a user requires to write a "fast kernel": it's a data point and if you define that as fast enough then Julia does it, or if you require 100% as fast as the fastest CUDA kernel ever found :shrug: it doesn't get there quite yet because there's a lot of people working on getting it there.
If you could express the problem in generalized index notation, tullio.jl is an even higher level abstraction that uses KA
I don't know about the numpy function you mentioned, but here's a KA matmul. Definitely doesn't look like it has any CUDA specific things. Just a regular loop with some restrictions
@kernel function matmul_kernel!(a, b, c) i, j = @index(Global, NTuple)
# creating a temporary sum variable for matrix multiplication
tmp_sum = zero(eltype(c))
for k = 1:size(a)[2]
tmp_sum += a[i,k] * b[k, j]
end
c[i,j] = tmp_sumWish I had time to write the code and compare to pytorch for you though
using Tullio, OffsetArrays
# A convolution with cyclic indices mat = zeros(10,10,1); mat[2,2] = 101; mat[10,10] = 1; @tullio kern[i,j] := 1/(1+i^2+j^2) (i in -3:3, j in -3:3)
@tullio out[x,y,c] := begin xi = mod(x+i, axes(mat,1)) # xi = ... means that it won't be summed, yj = mod(y+j, axes(mat,2)) @inbounds trunc(Int, mat[xi, yj, c] * kern[i,j]) # and disables automatic @inbounds, end (x in 1:10, y in 1:10) # and prevents range of x from being inferred.
It's mostly just math!
Here's an example: not long ago I needed a fast GPU implementation of numpy.unpackbits. I can do it in Pytorch, but it's slow and it wastes a ton of memory. So, my question is: can I write it in Julia as efficiently as I can do it in CUDA? If yes, will it be easier to do so?