# Laguerre polynomials on GPU

**URL:** <https://discourse.julialang.org/t/laguerre-polynomials-on-gpu/99571>\
**Category:** GPU\
**Tags:** gpu, cuda, polynomials\
**Created:** [May 29, 2023, 3:42pm UTC](https://discourse.julialang.org/t/laguerre-polynomials-on-gpu/99571 "2023-05-29T15:42:54Z")\
**Posts on this page:** 8\
**Page:** 2

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [June 8, 2023, 12:32pm UTC](https://discourse.julialang.org/t/laguerre-polynomials-on-gpu/99571/21 "2023-06-08T12:32:52Z")

</div>

This looks rally efficient indeed. Thanks for the link i learned something!

---

<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:** [June 8, 2023, 6:42pm UTC](https://discourse.julialang.org/t/laguerre-polynomials-on-gpu/99571/22 "2023-06-08T18:42:48Z")

</div>

> [@dlfivefifty](#):
>
> But you can make your own implementation in 3 lines that can easily be used on a GPU:
> 
> ```julia
> function laguerre(n, x)
> p0,p1 = 0.0,1.0
> for k = 0:n-1
> p1,p0 = ((2k+1)/(k+1) - x/(k+1))*p1 - k/(k+1)*p0,p1
> end
> p1
> end
> 
> ```

Thanks! it definitely worked! It is slightly slower than the other implementation for higher degrees of the polynomial, but it is still quick. And it works with the GPU.

```julia
function laguerre(n::Int, x::T) where {T <: BlasFloat}
    p0, p1 = 0.0, 1.0
    for k = 0:n-1
        p1, p0 = ((2k+1)/(k+1) - x/(k+1))*p1 - k/(k+1)*p0, p1
    end
    p1
end

function _genlaguerre(n::Int, α::Int, x::T) where {T<:BlasFloat}
    t = binomial(n+α,n)
    L = t
    for k = 1:n
        t *= (-x) * (n-k+1) / (k * (k+α))
        L += t
    end
    return L
end

A = rand(500, 500);
B = cu(A);

@benchmark laguerre.(80, $A) seconds=10

BenchmarkTools.Trial: 214 samples with 1 evaluation.
 Range (min … max): 46.647 ms … 49.016 ms ┊ GC (min … max): 0.00% … 3.79%
 Time (median): 46.703 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 46.839 ms ± 414.309 μs ┊ GC (mean ± σ): 0.17% ± 0.77%

  ▂█▄▁                                                          
  ████▆▁▁▁▁▄▁▁▁▁▁█▇▆▆▅▆▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▅
  46.6 ms Histogram: log(frequency) by time 48.6 ms <

 Memory estimate: 1.91 MiB, allocs estimate: 2.

@benchmark _genlaguerre.(80, 0, $A) seconds=10

BenchmarkTools.Trial: 342 samples with 1 evaluation.
 Range (min … max): 28.836 ms … 31.543 ms ┊ GC (min … max): 0.00% … 5.92%
 Time (median): 29.149 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 29.253 ms ± 399.656 μs ┊ GC (mean ± σ): 0.26% ± 1.19%

         ▆█▅▁ ▁                                       
  ▄▁▁▁▁▄▆████▅▁▁▁▁▁▁▁▁▁▆█▆▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▆█ ▆
  28.8 ms Histogram: log(frequency) by time 31 ms <

 Memory estimate: 1.91 MiB, allocs estimate: 2.

@benchmark laguerre.(80, $B) seconds=10

BenchmarkTools.Trial: 2167 samples with 5 evaluations.
 Range (min … max): 6.682 μs … 1.078 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 1.019 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 922.869 μs ± 296.625 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

                                                              █  
  ▃▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▁▁█ ▂
  6.68 μs Histogram: frequency by time 1.02 ms <

 Memory estimate: 1.28 KiB, allocs estimate: 26.

```

It would be better to extend it also for the generalized case. But it is what I needed for the moment.

I have only one question. The input array `B` is a Float32, while the output is Float64. I guess that it comes from the initialization of the variables `p0, p1 = 0.0, 1.0`. I tried something like `p0, p1 = zero(x), one(x)`, but it doesn’t work.

```julia
GPU broadcast resulted in non-concrete element type AbstractFloat.

```

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [June 8, 2023, 7:16pm UTC](https://discourse.julialang.org/t/laguerre-polynomials-on-gpu/99571/23 "2023-06-08T19:16:20Z")

</div>

You also need to make sure the rational expressions have the right types since dividing Int will return a Float64

---

<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:** [June 8, 2023, 7:16pm UTC](https://discourse.julialang.org/t/laguerre-polynomials-on-gpu/99571/24 "2023-06-08T19:16:50Z")

</div>

I previously said, together with @lrnv, that when alpha is a positive integer, the generalized Laguerre can be seen as the case with alpha=0 and different n. However I don’t remember if it is effectively true.

---

<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:** [June 8, 2023, 7:59pm UTC](https://discourse.julialang.org/t/laguerre-polynomials-on-gpu/99571/25 "2023-06-08T19:59:15Z")

</div>

> [@dlfivefifty](#):
>
> You also need to make sure the rational expressions have the right types since dividing Int will return a Float64

Ok the following code works

```julia
function laguerre(n::Int, x::T) where {T <: BlasFloat}
    p0, p1 = zero(T), one(T)
    for k = zero(T):n-1
        p1, p0 = ((2k+1)/(k+1) - x/(k+1))*p1 - k/(k+1)*p0, p1
    end
    p1
end

```

In my implementation I use integer alphas, which I remeber that they can be mapped to the case of zero alpha and different n of the polynomial. But I don’t remeber the relation, if any.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [June 8, 2023, 8:03pm UTC](https://discourse.julialang.org/t/laguerre-polynomials-on-gpu/99571/26 "2023-06-08T20:03:13Z")

</div>

Adding general α is a very simple modification, see the DLMF link I sent for how to change the rationals

---

<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:** [June 8, 2023, 11:46pm UTC](https://discourse.julialang.org/t/laguerre-polynomials-on-gpu/99571/27 "2023-06-08T23:46:06Z")

</div>

Thanks! The generalized Laguerre Polynomial now works also on the GPU!

```julia
function laguerre(n::Int, α::Number, x::T) where {T <: BlasFloat}
    α = convert(T, α)
    p0, p1 = one(T), -x+(α+1)
    n == 0 && return p0
    for k = one(T):n-1
        p1, p0 = ((2k+α+1)/(k+1) - x/(k+1))*p1 - (k+α)/(k+1)*p0, p1
    end
    p1
end

A = CUDA.rand(1000, 1000)

@benchmark laguerre.(10, 2, $A)

BenchmarkTools.Trial: 4153 samples with 5 evaluations.
 Range (min … max): 6.432 μs … 54.395 ms ┊ GC (min … max): 0.00% … 1.86%
 Time (median): 250.764 μs ┊ GC (median): 0.00%
 Time (mean ± σ): 240.418 μs ± 1.456 ms ┊ GC (mean ± σ): 0.30% ± 0.05%

                                                             █  
  ▆▁▁▁▁▂▁▁▁▁▁▁▁▁▂▁▁▁▁▁▂▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂█ ▂
  6.43 μs Histogram: frequency by time 251 μs <

 Memory estimate: 1.28 KiB, allocs estimate: 26.

```

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [June 9, 2023, 11:35am UTC](https://discourse.julialang.org/t/laguerre-polynomials-on-gpu/99571/28 "2023-06-09T11:35:19Z")

</div>

> [@albertomercurio](#):
>
> ```julia
> for k = one(T):n-1
> p1, p0 = ((2k+α+1)/(k+1) - x/(k+1))*p1 - (k+α)/(k+1)*p0, p1
> end
> 
> ```

There is no need for `one(T)` here, since all the numerators are of type `T`. And iterating over a float range can be significantly slower than over an `Int` range. On the CPU the difference is dramatic, but on the GPU it seems small.

Anyway, change `one(T)` to `1`.

[Previous page](https://discourse.julialang.org/t/laguerre-polynomials-on-gpu/99571.md?page=1)