# Why functions in SpecialFunctions package work on CUDA arrays?

**URL:** <https://discourse.julialang.org/t/why-functions-in-specialfunctions-package-work-on-cuda-arrays/128194>\
**Category:** Performance\
**Tags:** question, cuda, specialfunctions\
**Created:** [April 18, 2025, 10:54pm UTC](https://discourse.julialang.org/t/why-functions-in-specialfunctions-package-work-on-cuda-arrays/128194 "2025-04-18T22:54:55Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![singularity](https://avatars.discourse-cdn.com/v4/letter/s/898d66/32.png) [@singularity](https://discourse.julialang.org/u/singularity)\
**Post date:** [April 18, 2025, 10:54pm UTC](https://discourse.julialang.org/t/why-functions-in-specialfunctions-package-work-on-cuda-arrays/128194/1 "2025-04-18T22:54:56Z")

</div>

HI. I am confused how the behavior of functions inside the SpecialFunctions package on CUDA arrays. I am confused why the code below works and doesn’t throw an error at the last line.

```julia

using CUDA
using SpecialFunctions
x = rand(Float32, 64, 64, 25)
x_gpu = cu(x)
z = loggamma.(x_gpu)

```

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [April 19, 2025, 2:21am UTC](https://discourse.julialang.org/t/why-functions-in-specialfunctions-package-work-on-cuda-arrays/128194/2 "2025-04-19T02:21:30Z")

</div>

_ **EDIT:** I realized after writing this post that CUDA comes with its own `loggamma` implementation, so what I wrote below is inaccurate: the implementation details of `SpecialFunctions.loggamma` don’t matter. See the next post for more details. I’m leaving the remainder of this post unchanged, so just imagine that we’re discussing some generic Julia function, and then read the next post to learn how this simplifies in the case of `loggamma`._

* * *

The first thing to realize is that many functions in SpecialFunctions.jl, including `loggamma`, are implemented in Julia and do not call out to a compiled C library. You can see the guts of `loggamma(::Float32)` here: [SpecialFunctions.jl/src/logabsgamma/e\_lgammaf\_r.jl at v2.5.1 · JuliaMath/SpecialFunctions.jl · GitHub](https://github.com/JuliaMath/SpecialFunctions.jl/blob/v2.5.1/src/logabsgamma/e_lgammaf_r.jl).

Another important detail is that `loggamma` does not allocate any intermediate arrays. It’s just a sequence of arithmetic and binary operations applied to plain numbers.

So your code works for exactly the same reason that the following works:

```julia
using CUDA
x = rand(Float32, 64, 64, 25)
x_gpu = cu(x)
f(x) = log(2cosh(x))
z = f.(x_gpu)

```

You’re just using a slightly more complicated `f`.

* * *

The obvious follow-up question is, well, so why _does_ that work? I can give you a rough idea, but will have to defer to the experts for details.

- Think about what broadcasting does on a regular CPU array. Essentially, `z_cpu = loggamma.(x)` would be translated into a loop like the following:

```julia
function loggamma_broadcast!(output, input)
    for i in eachindex(output, input)
        output[i] = loggamma(input[i])
    end
end
z_cpu = similar(x)
loggamma_broadcast!(z_cpu, x)
z_cpu

```

- GPU arrays have their own implementation of the broadcasting machinery, which writes GPU kernels to execute the broadcasted operations. Your example `z = loggamma.(x_gpu)` is essentially translated into something like the following (don’t worry about the details here; the point is that this also looks like a regular Julia function that calls `loggamma` on elements of `input`):

```julia
function loggamma_broadcast_cuda!(output, input)
    index = threadIdx().x
    stride = blockDim().x
    for i in index:stride:length(output)
        @inbounds output[i] = loggamma(input[i])
    end
    return nothing
end
z = similar(x_gpu)
@cuda threads=256 loggamma_broadcast_cuda!(z, x_gpu)
z

```

To read more about how to write and launch CUDA kernels in Julia, see [Introduction · CUDA.jl](https://cuda.juliagpu.org/stable/tutorials/introduction/)
- The magic of CUDA.jl (with the help of GPUCompiler.jl) is that regular Julia code such as this, which calls other regular Julia functions such as `loggamma`, can be compiled to CUDA kernels that run on the GPU. Someone else will have to explain the details of how that works. (I’ve seen terms like _method overlay tables_ thrown around, and I know that LLVM has a backend for emitting Nvidia PTX assembly: [https://llvm.org/docs/NVPTXUsage.html.](https://llvm.org/docs/NVPTXUsage.html.))

* * *

**P.S.:** In reality, GPU array broadcasting uses KernelAbstractions.jl to produce device-agnostic kernels that work on many different GPUs, not just Nvidia ones. KernelAbstractions.jl then produces CUDA.jl kernel code like the above. I skipped this extra layer of magic to show more directly how `loggamma` is used in the code that CUDA.jl compiles to the GPU. However, the KernelAbstractions.jl version is much simpler; it looks something like this:

```julia
@kernel function loggamma_broadcast_KA!(output, input)
    I = @index(Global)
    @inbounds output[I] = loggamma(input[I])
end

```

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [April 19, 2025, 3:38am UTC](https://discourse.julialang.org/t/why-functions-in-specialfunctions-package-work-on-cuda-arrays/128194/3 "2025-04-19T03:38:03Z")

</div>

> [@danielwe](#):
>
> The first thing to realize is that all the functions in [SpecialFunctions.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/SpecialFunctions) are implemented in Julia.

Oops! In the case of `loggamma`, what I said here is actually beside the point. When compiling for the GPU, CUDA.jl obviously needs to replace Julia functions with GPU native functions at some point—if not before, then at least at the level of hardware instructions like `+` and `*`. But Nvidia provides its own library of math functions for their GPUs, including things like `sin` and `exp` and, you guessed it, `loggamma`. So when CUDA.jl is compiling `loggamma_kernel_cuda!` and encounters `SpecialFunctions.loggamma`, it simply replaces this call with the GPU native `loggamma`. It never looks at the implementation of `SpecialFunctions.loggamma`.

This is an example of the method overlay table in action. You can see how the override for `SpecialFunctions.loggamma` is registered here: [CUDA.jl/ext/SpecialFunctionsExt.jl at v5.7.3 · JuliaGPU/CUDA.jl · GitHub](https://github.com/JuliaGPU/CUDA.jl/blob/v5.7.3/ext/SpecialFunctionsExt.jl#L31-L32)

In the example I provided for comparison, with `f(x) = log(2cosh(x))`, there are overrides for `log` and `cosh`, but obviously not for the function `f` that I defined myself, so this is actually a better example of CUDA.jl being able to compile GPU kernels from Julia functions that call other Julia functions.

---

<div class="post-metadata">

**Author:** ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)\
**Post date:** [April 19, 2025, 4:18am UTC](https://discourse.julialang.org/t/why-functions-in-specialfunctions-package-work-on-cuda-arrays/128194/4 "2025-04-19T04:18:04Z")

</div>

Halfway through the first comment, I was about to say wasn’t this [an issue](https://github.com/JuliaGPU/CUDA.jl/issues/1528) when `SpecialFunctions` deprecated `lgamma` for `loggamma` and CUDA.jl didn’t know to use the CUDA C intrinsics for and failed to compile `loggamma`. But looking at it more closely, it failed because `lgamma`/`loggamma` was `ccall`-ing routines for `Float32`/`Float64`/`BigFloat` and converting other `Real`s to and from those. Currently `loggamma` for `Float32`/`Float64` uses [pure Julia ports of OpenLibm’s algorithms](https://github.com/JuliaMath/SpecialFunctions.jl/pull/413), I wonder if CUDA.jl could compile those if it weren’t using CUDA C intrinsics.

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [April 19, 2025, 5:57am UTC](https://discourse.julialang.org/t/why-functions-in-specialfunctions-package-work-on-cuda-arrays/128194/5 "2025-04-19T05:57:39Z")

</div>

> [@Benny](#):
>
> I wonder if [CUDA.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/CUDA) could compile those if it weren’t using CUDA C intrinsics.

Indeed, it can:

```julia-repl
julia> using CUDA

julia> using IrrationalConstants:
           twoπ,
           halfπ,
           sqrtπ,
           sqrt2π,
           invπ,
           inv2π,
           invsqrt2,
           invsqrt2π,
           logtwo,
           logπ,
           log2π

julia> function logabsgamma(x::Float32)
           # [...]
           # implementation copied from https://github.com/JuliaMath/SpecialFunctions.jl/blob/v2.5.1/src/logabsgamma/e_lgammaf_r.jl
           # [...]
       end
logabsgamma (generic function with 1 method)

julia> function loggamma(x::Float32)
           # just skipping the f(x::Real) = f(float(x)) wrappers
           # otherwise identical to SpecialFunctions.loggamma
           (y, s) = logabsgamma(x)
           s < 0 && throw(DomainError(x, "`gamma(x)` must be non-negative"))
           return y
       end
loggamma (generic function with 1 method)

julia> x = rand(Float32, 64, 64, 25);

julia> x_gpu = cu(x);

julia> z = loggamma.(x_gpu) # success!
64×64×25 CuArray{Float32, 3, CUDA.DeviceMemory}:
[...]

```

More surprisingly, it’s only 10 % slower than the intrinsic version (tested on an NVIDIA A100):

```julia-repl
julia> using BenchmarkTools

julia> @btime CUDA.@sync loggamma.($x_gpu);
  40.833 μs (98 allocations: 3.41 KiB)

julia> import SpecialFunctions

julia> @btime CUDA.@sync SpecialFunctions.loggamma.($x_gpu); # overridden in CUDA.jl
  36.625 μs (98 allocations: 3.41 KiB)

```

The `SpecialFunctions` implementation doesn’t look very GPU-friendly with all its branches, but maybe it isn’t possible to do all that much better.

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [April 19, 2025, 6:23am UTC](https://discourse.julialang.org/t/why-functions-in-specialfunctions-package-work-on-cuda-arrays/128194/6 "2025-04-19T06:23:04Z")

</div>

those branches are fairly necessary. you can’t cover the domain with a single implementation.

---

<div class="post-metadata">

**Author:** ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)\
**Post date:** [April 19, 2025, 1:50pm UTC](https://discourse.julialang.org/t/why-functions-in-specialfunctions-package-work-on-cuda-arrays/128194/7 "2025-04-19T13:50:39Z")

</div>

CUDA C’s [lgamma and lgammaf](https://docs.nvidia.com/cuda/cuda-c-programming-guide/#standard-functions) have 4 and 6 maximum ULPs error (actually more if inside some intervals), while the `SpecialFunction`’s `loggamma` were tested to usually max out at 1.5 ULPs (2.5 for a miniscule fraction of positives, though the underlying `logsabsgamma` can get up to 3.6e6 ULPs for a small fraction of negatives). I’m however not confident lgamma and lgammaf are the same as CUDA’s libdevice functions [\_\_nv\_lgamma and \_\_nv\_lgammaf](https://docs.nvidia.com/cuda/libdevice-users-guide/ __nv_lgamma.html#__ nv_lgamma), which CUDA.jl uses in SpecialFunctionsExt.

---

<div class="post-metadata">

**Author:** ![maleadt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maleadt/32/10097_2.png) [@maleadt](https://discourse.julialang.org/u/maleadt)\
**Post date:** [April 20, 2025, 7:04am UTC](https://discourse.julialang.org/t/why-functions-in-specialfunctions-package-work-on-cuda-arrays/128194/8 "2025-04-20T07:04:45Z")

</div>

> [@Benny](#):
>
> I’m however not confident lgamma and lgammaf are the same as CUDA’s libdevice functions [\_\_nv\_lgamma and \_\_nv\_lgammaf](https://docs.nvidia.com/cuda/libdevice-users-guide/ __nv_lgamma.html#__ nv_lgamma), which [CUDA.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/CUDA) uses in SpecialFunctionsExt.

Yes, those are from libdevice.

If the additional precision is important, I’d be happy to take a PR that switches from calling `libdevice` to a native Julia implementation. We generally strive to be compatible with Julia code on the CPU, so that can even come at a slight performance cost (we could always keep the `libdevice` calls as a `_fast` version for use with `@fastmath`, if SpecialFunctions.jl supports that).