# Illegal memory indexing for broadcast of a recursive function

**URL:** <https://discourse.julialang.org/t/illegal-memory-indexing-for-broadcast-of-a-recursive-function/99534>\
**Category:** GPU\
**Tags:** gpu, indexing, cuda, array, memory\
**Created:** [May 29, 2023, 12:29am UTC](https://discourse.julialang.org/t/illegal-memory-indexing-for-broadcast-of-a-recursive-function/99534 "2023-05-29T00:29:19Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [May 29, 2023, 12:29am UTC](https://discourse.julialang.org/t/illegal-memory-indexing-for-broadcast-of-a-recursive-function/99534/1 "2023-05-29T00:29:19Z")

</div>

Hello,

I’m trying to implement the generalized Laguerre function acting element-wise to the elements of a CuArray.

The basic code uses a recursive function to evaluate the laguerre polynomial

```julia
using CUDA
using LinearAlgebra: BlasFloat
CUDA.allowscalar(false)

function gen_laguerre(n::Integer, alpha::Number, x::T)::T where T <: BlasFloat
    n == 0 && return 1
    n == 1 && return 1 + alpha - x
    return ((2*n-1+alpha-x)*gen_laguerre(n-1,alpha,x)-(n-1+alpha)*gen_laguerre(n-2,alpha,x))/n
end

A = CUDA.rand(1000)
gen_laguerre.(10, 3, A)

```

works, but if I increase the degree of the polynomial

```julia
gen_laguerre.(30, 3, A)

CUDA error: an illegal memory access was encountered (code 700, ERROR_ILLEGAL_ADDRESS)

Stacktrace:
  [1] throw_api_error(res::CUDA.cudaError_enum)
    @ CUDA ~/.julia/packages/CUDA/HzJ3E/lib/cudadrv/libcuda.jl:27
  [2] check
    @ ~/.julia/packages/CUDA/HzJ3E/lib/cudadrv/libcuda.jl:34 [inlined]
  [3] cuMemAllocAsync
    @ ~/.julia/packages/CUDA/HzJ3E/lib/utils/call.jl:26 [inlined]
  [4] alloc(::Type{CUDA.Mem.DeviceBuffer}, bytesize::Int64; async::Bool, stream::CuStream, pool::Nothing)
    @ CUDA.Mem ~/.julia/packages/CUDA/HzJ3E/lib/cudadrv/memory.jl:83
  [5] alloc
    @ ~/.julia/packages/CUDA/HzJ3E/lib/cudadrv/memory.jl:71 [inlined]
  [6] actual_alloc(bytes::Int64; async::Bool, stream::CuStream)
    @ CUDA ~/.julia/packages/CUDA/HzJ3E/src/pool.jl:65
  [7] actual_alloc
    @ ~/.julia/packages/CUDA/HzJ3E/src/pool.jl:59 [inlined]
  [8] #993
    @ ~/.julia/packages/CUDA/HzJ3E/src/pool.jl:417 [inlined]
  [9] retry_reclaim(f::CUDA.var"#993#995"{CuStream, Int64}, isfailed::typeof(isnothing))
    @ CUDA ~/.julia/packages/CUDA/HzJ3E/src/pool.jl:337
 [10] macro expansion
    @ ~/.julia/packages/CUDA/HzJ3E/src/pool.jl:416 [inlined]
 [11] macro expansion
    @ ./timing.jl:393 [inlined]
 [12] #_alloc#992
    @ ~/.julia/packages/CUDA/HzJ3E/src/pool.jl:413 [inlined]
 [13] _alloc
    @ ~/.julia/packages/CUDA/HzJ3E/src/pool.jl:408 [inlined]
 [14] #alloc#991
    @ ~/.julia/packages/CUDA/HzJ3E/src/pool.jl:398 [inlined]
 [15] alloc
    @ ~/.julia/packages/CUDA/HzJ3E/src/pool.jl:392 [inlined]
 [16] CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}(#unused#::UndefInitializer, dims::Tuple{Int64})
    @ CUDA ~/.julia/packages/CUDA/HzJ3E/src/array.jl:93
 [17] CuArray
    @ ~/.julia/packages/CUDA/HzJ3E/src/array.jl:176 [inlined]
 [18] CuArray
    @ ~/.julia/packages/CUDA/HzJ3E/src/array.jl:187 [inlined]
 [19] similar
    @ ./abstractarray.jl:882 [inlined]
 [20] similar
    @ ./abstractarray.jl:881 [inlined]
 [21] similar
    @ ~/.julia/packages/CUDA/HzJ3E/src/broadcast.jl:11 [inlined]
 [22] copy
    @ ~/.julia/packages/GPUArrays/TnEpb/src/host/broadcast.jl:37 [inlined]
 [23] materialize(bc::Base.Broadcast.Broadcasted{CUDA.CuArrayStyle{1}, Nothing, typeof(gen_laguerre), Tuple{Int64, Int64, CuArray{Float32, 1, CUDA.Mem.DeviceBuffer}}})
    @ Base.Broadcast ./broadcast.jl:873
 [24] top-level scope
    @ In[25]:1

```

And then it is impossible to calculate again a lower order Laguerre polynomial.

I can’t figure out the source of the problem.

**EDIT**  
I noticed that on the CPU it tooks a lot of time when n\>10. Maybe this is the reason of the fail. Howerver I can’t understand why it fails, while I understand the long calculation time on the CPU.

---

<div class="post-metadata">

**Author:** ![jling](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jling/32/212909_2.png) [@jling](https://discourse.julialang.org/u/jling)\
**Post date:** [May 29, 2023, 4:20am UTC](https://discourse.julialang.org/t/illegal-memory-indexing-for-broadcast-of-a-recursive-function/99534/2 "2023-05-29T04:20:25Z")

</div>

Yeah, avoid recursive is probably the easiest solution regardless of the exact underlying problem

---

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [May 29, 2023, 3:42pm UTC](https://discourse.julialang.org/t/illegal-memory-indexing-for-broadcast-of-a-recursive-function/99534/3 "2023-05-29T15:42:28Z")

</div>

But why CUDA doesn’t support recursive functions?

---

<div class="post-metadata">

**Author:** ![mcabbott](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mcabbott/32/6603_2.png) [@mcabbott](https://discourse.julialang.org/u/mcabbott)\
**Post date:** [May 29, 2023, 4:13pm UTC](https://discourse.julialang.org/t/illegal-memory-indexing-for-broadcast-of-a-recursive-function/99534/4 "2023-05-29T16:13:06Z")

</div>

> [@albertomercurio](#):
>
> ```julia
> n == 0 && return 1
> n == 1 && return 1 + alpha - x
> 
> ```

This looks type-unstable. Maybe the assertion `::T` fixes it but maybe not enough?

I would fix that, and try inlining everything. Then try using `Val(n)` not `n` so that the recursion can be inferred.

> ****
>
> ```julia
> julia> function gen_laguerre_2(n::Integer, alpha::Number, x::Number)
> n == 0 && return 1
> n == 1 && return 1 + alpha - x
> return ((2*n-1+alpha-x)*gen_laguerre(n-1,alpha,x)-(n-1+alpha)*gen_laguerre(n-2,alpha,x))/n
> end
> gen_laguerre_2 (generic function with 1 method)
> 
> julia> @code_warntype gen_laguerre_2(3, 0.1, 0.1)
> MethodInstance for gen_laguerre_2(::Int64, ::Float64, ::Float64)
> from gen_laguerre_2(n::Integer, alpha::Number, x::Number) @ Main REPL[96]:1
> Arguments
> #self#::Core.Const(gen_laguerre_2)
> n::Int64
> alpha::Float64
> x::Float64
> Body::Union{Float64, Int64}
> 
> ```

---

<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:** [May 30, 2023, 10:12am UTC](https://discourse.julialang.org/t/illegal-memory-indexing-for-broadcast-of-a-recursive-function/99534/5 "2023-05-30T10:12:33Z")

</div>

It’s likely overflowing the stack, indeed. If you’re bent on using this implementation, you can try querying and increasing the stack size by calling `CUDA.limit!` with `CU_LIMIT_STACK_SIZE`.

---

<div class="post-metadata">

**Author:** ![albertomercurio](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertomercurio/32/27051_2.png) [@albertomercurio](https://discourse.julialang.org/u/albertomercurio)\
**Post date:** [May 30, 2023, 11:14pm UTC](https://discourse.julialang.org/t/illegal-memory-indexing-for-broadcast-of-a-recursive-function/99534/6 "2023-05-30T23:14:25Z")

</div>

I can’t see an example on the documentation. Will the following code work?

```julia
CUDA.limit!(CUDA.CU_LIMIT_STACK_SIZE, 1000)

```

where 1000 is the maximum number of recursion I expect.

---

<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:** [May 31, 2023, 7:20am UTC](https://discourse.julialang.org/t/illegal-memory-indexing-for-broadcast-of-a-recursive-function/99534/7 "2023-05-31T07:20:16Z")

</div>

> [@albertomercurio](#):
>
> I can’t see an example on the documentation.

We don’t re-document all of CUDA, so refer to the NVIDIA documentation for low-level details: [CUDA Driver API :: CUDA Toolkit Documentation](https://docs.nvidia.com/cuda/cuda-driver-api/group __CUDA__ CTX.html#group __CUDA__ CTX_1g9f2d47d1745752aa16da7ed0d111b6a8).

Stack size is expressed in bytes, not in number of recursions.