# CuArray + GLMakie

**URL:** <https://discourse.julialang.org/t/cuarray-glmakie/52461>\
**Category:** Visualization\
**Tags:** plotting\
**Created:** [December 27, 2020, 3:37pm UTC](https://discourse.julialang.org/t/cuarray-glmakie/52461 "2020-12-27T15:37:13Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [December 27, 2020, 3:37pm UTC](https://discourse.julialang.org/t/cuarray-glmakie/52461/1 "2020-12-27T15:37:13Z")

</div>

Hi,  
I wonder if it would make sense (from the performance point of view) to be able to use `GLMakie.jl` with `CuArray`s from `CUDA.jl` ?

---

<div class="post-metadata">

**Author:** ![sdanisch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdanisch/32/1406_2.png) [@sdanisch](https://discourse.julialang.org/u/sdanisch)\
**Post date:** [December 27, 2020, 6:58pm UTC](https://discourse.julialang.org/t/cuarray-glmakie/52461/2 "2020-12-27T18:58:42Z")

</div>

Yes would be lovely and should be very efficient - but isn’t currently easy.  
First step is to figure out if the API parts for OpenGL interop are wrapped in CUDA.jl ( I think not, we should open an issue), and then we need to port the samples in [CUDA Samples :: CUDA Toolkit Documentation](https://docs.nvidia.com/cuda/cuda-samples/index.html#simple-opengl) to Julia + CUDA… Makie should mostly work with GLBuffers or GLTextures directly, so after they got created and hooked up with CUDA.jl, there shouldn’t be a big problem to plot those 😉

---

<div class="post-metadata">

**Author:** ![sdanisch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdanisch/32/1406_2.png) [@sdanisch](https://discourse.julialang.org/u/sdanisch)\
**Post date:** [December 28, 2020, 12:35pm UTC](https://discourse.julialang.org/t/cuarray-glmakie/52461/3 "2020-12-28T12:35:09Z")

</div>

Actually they’re wrapped! Let me see if I can create a simple example. What do you actually want to plot? Sadly, surfaces, meshes and scatter will need quite the different setup to share the resources with CUDA.

---

<div class="post-metadata">

**Author:** ![LaurentPlagne](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/laurentplagne/32/10103_2.png) [@LaurentPlagne](https://discourse.julialang.org/u/LaurentPlagne)\
**Post date:** [December 28, 2020, 4:41pm UTC](https://discourse.julialang.org/t/cuarray-glmakie/52461/4 "2020-12-28T16:41:37Z")

</div>

Hi Simon,  
Here is an example which is unfortunately a bit too long for the 7-lines thread…  
Anyway, It takes about 5s on my machine to run and show a rather convincing combination of `CUDA.jl` and `Makie.jl`.  
In the plot (can be `heatmap` or `surface`), I copy the `CuArray zs` to a standard `Array czs`…

![toto](https://global.discourse-cdn.com/julialang/original/3X/7/6/76ae1d2c297b7968245a3435f83275cf3feb8bd2.gif)

```julia
using CUDA,GLMakie,AbstractPlotting
const T=Float32
const V=CuArray{T,2}
function go()
    CUDA.allowscalar(false)
    ls,ns,ms,λ,dt,sr,k,nt,two,eight,fr =1.e-2,800,1.e-4,0.1,2.e-7,0.1,5e7,5000,T(2),T(8),10
    dt2sm=dt^2/ms
    r,r0,r1,r2=1:ns,1:ns-2,2:ns-1,3:ns
    cxs = [T(j-1)*ls for i in r, j in r]
    dx(i,j) = cxs[i,j]-cxs[ns÷2,ns÷2]
    cxc = [cxs[i,j]+sr*(dx(i,j))*exp(-0.5(dx(i,j)^2+dx(j,i)^2)/λ^2)/λ for i in r, j in r]
    xs,xc,ys,yc=V(cxs),V(cxc),V(cxs'),V(cxc')
    # ys,yc=collect(xs'),collect(xc')
    xp,yp,xt,yt,fx,fy,zs=copy(xc),copy(yc),copy(xc),copy(yc),zero(xc),zero(xc),zero(xc)
    @. zs = sqrt((xc-xs)^2+(yc-ys)^2) 
    czs=Array(zs) ; zn=Node(czs) ; zl=sr*0.05
    scene = Scene(resolution = (400,400)) ; scale!(scene, 1, 1, 10000)
    heatmap!(scene,1:ns,1:ns,lift(z->z,zn),colorrange = (-zl,zl))
    # surface!(scene,1:ns,1:ns,lift(z->z,zn),colorrange = (-zl,zl))
    zlims!((-zl,zl))
    display(scene)
    GLMakie.record(scene, "output.gif", 1:nt÷fr, framerate=30) do j
    # for j in 1:nt÷fr
        for i in 1:fr
            @views @. fx[r1,r1] = -k*(eight*xc[r1,r1]-xc[r0,r0]-xc[r1,r0]-xc[r2,r0]-xc[r0,r1]-xc[r2,r1]-xc[r0,r2]-xc[r1,r2]-xc[r2,r2])
            @views @. fy[r1,r1] = -k*(eight*yc[r1,r1]-yc[r0,r0]-yc[r1,r0]-yc[r2,r0]-yc[r0,r1]-yc[r2,r1]-yc[r0,r2]-yc[r1,r2]-yc[r2,r2])
            @. xt = two*xc-xp+fx*dt2sm
            @. yt = two*yc-yp+fy*dt2sm
            xc,xp,xt = xt,xc,xp
            yc,yp,yt = yt,yc,yp
        end
        @. zs = sqrt((xc-xs)^2+(yc-ys)^2)
        czs=Array(zs)
        zn[] = czs
        yield()
    end
end
@time go()

```

![toto](https://global.discourse-cdn.com/julialang/original/3X/1/c/1c0a8283dd2e18237210ae53fd6b400de8cffaa6.gif)

---

<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 1, 2022, 1:17pm UTC](https://discourse.julialang.org/t/cuarray-glmakie/52461/5 "2022-04-01T13:17:09Z")

</div>

Here’s an example on how to keep the data on the GPU:

```julia
using GLMakie, CUDA
using AbstractFFTs

using GLMakie: GLAbstraction

function plot(T=Float32, N=1024*1024) # 1024*1024 Float32's = 8MiB
    # generate dummy data: a simple gaussian
    # NOTE: for demonstration purposes, we're keeping the time series on the CPU
    t = range(-5, 5, length=N)
    f = CuArray{T}(exp.(-t.^2))

    ## initialization
    #
    # this should be only done once, keeping the buffer and resources across frames

    # get a buffer object and register it with CUDA
    buffer = GLAbstraction.GLBuffer(eltype(f), length(f))
    resource = let
        ref = Ref{CUDA.CUgraphicsResource}()
        CUDA.cuGraphicsGLRegisterBuffer(ref, buffer.id,
                                        CUDA.CU_GRAPHICS_MAP_RESOURCE_FLAGS_WRITE_DISCARD)
        ref[]
    end

    ## main processing
    #
    # this needs to be done every time we get new data and want to plot it

    plot = NVTX.@range "main" begin
        NVTX.@range "CUDA" CUDA.@sync begin # doesn't actually need @sync, just for timings
            # map OpenGL buffer object for writing from CUDA
            CUDA.cuGraphicsMapResources(1, [resource], stream())

            # get a CuArray object that we can work with
            array = let
                ptr_ref = Ref{CUDA.CUdeviceptr}()
                numbytes_ref = Ref{Csize_t}()
                CUDA.cuGraphicsResourceGetMappedPointer_v2(ptr_ref, numbytes_ref, resource)

                ptr = reinterpret(CuPtr{T}, ptr_ref[])
                len = Int(numbytes_ref[] ÷ sizeof(T))

                unsafe_wrap(CuArray, ptr, len)
            end

            # example processing: compute the FFT
            # NOTE: real applications will want to perform a pre-planned in-place FFT
            F = fft(f)
            # shift the zero frequency component to the center
            array[1:N÷2] .= abs.(view(F, N÷2+1:N))
            array[N÷2+1:N] .= abs.(view(F, 1:N÷2))

            CUDA.cuGraphicsUnmapResources(1, [resource], stream())
        end

        NVTX.@range "Makie" lines(t, buffer)
    end

    ## clean-up

    CUDA.cuGraphicsUnregisterResource(resource)

    return plot
end

function main()
    # XXX: we need create a screen, which initializes a GL Context,
    # so that we can create a GLBuffer before having rendered anything.
    GLMakie.Screen()

    save("plot.png", plot())
    return
end

```

I was experimenting with this in order to create a waterfall plot for data that resides on the GPU (hence the `fft`), but performance isn’t great: Processing the 8MB dummy data data using CUDA takes 500us, rendering that using GLMakie takes 500ms… Anything obvious I’m doing wrong here?

---

<div class="post-metadata">

**Author:** ![sdanisch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdanisch/32/1406_2.png) [@sdanisch](https://discourse.julialang.org/u/sdanisch)\
**Post date:** [April 1, 2022, 1:27pm UTC](https://discourse.julialang.org/t/cuarray-glmakie/52461/6 "2022-04-01T13:27:52Z")

</div>

Multiple things that look problematic:

1. I’m not sure how you time the rendering, but `save("plot.png", plot())` does quite a bit more than just rendering one frame, so I guess 500ms could be (almost) plausible. Rendering exactly one frame takes a bit more manual setup
2. `GLMakie.Screen()` + `save("plot.png", plot())` isn’t guaranteed to use the same context
3. `lines(t, buffer)` will run into the conversion machine and will just collect the buffer I’m afraid. I will need to double check, but `GLBuffer(Point2f[...])` may get to OpenGL without getting converted

I can take a look later if I can update the example to correctly render one frame and use the same context and not convert the GLBuffer to an array

---

<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 1, 2022, 1:43pm UTC](https://discourse.julialang.org/t/cuarray-glmakie/52461/7 "2022-04-01T13:43:01Z")

</div>

I was timing the NVTX ranges, so it should exclude the saving. But yeah, I’m new to Makie, so probably doing some silly mistakes 🙂 If this ends up working properly I’ll add the necessary high-level wrappers to CUDA.jl so that we can avoid some of the boilerplate here.

---

<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 1, 2022, 2:49pm UTC](https://discourse.julialang.org/t/cuarray-glmakie/52461/8 "2022-04-01T14:49:14Z")

</div>

> [@sdanisch](#):
>
> `lines(t, buffer)` will run into the conversion machine and will just collect the buffer I’m afraid. I will need to double check, but `GLBuffer(Point2f[...])` may get to OpenGL without getting converted

When tracing with NSight, I didn’t see any copies to CPU memory (i.e., no DtoH copies), so I don’t think it converted back.

---

<div class="post-metadata">

**Author:** ![sdanisch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdanisch/32/1406_2.png) [@sdanisch](https://discourse.julialang.org/u/sdanisch)\
**Post date:** [April 1, 2022, 4:14pm UTC](https://discourse.julialang.org/t/cuarray-glmakie/52461/9 "2022-04-01T16:14:45Z")

</div>

> When tracing with NSight, I didn’t see any copies to CPU memory

Does that trace OpenGL Buffer copies? The OpenGL backend can’t render separated x, y buffers, so if the plot produces an output, it should copy it to CPU…

---

<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 5, 2022, 9:01am UTC](https://discourse.julialang.org/t/cuarray-glmakie/52461/10 "2022-04-05T09:01:52Z")

</div>

Turns out there was some “scalar indexing” going on, where individual elements of the GPU buffers were being copied to the CPU for processing. By using `GLBuffer{Point2f}` directly, and using low-level APIs that avoid some of the automatic processing (e.g. `center!` inspects data to determine limits, `lines!` does a `map` to detect invalid vertices) I have a working example that keeps all data on the GPU:

```julia
function plot(; T=Float32, N=1024, resolution=(800,600))
    t = CUDA.rand(T, N)
    X = CUDA.rand(T, N)

    ## initialization
    #
    # this should be only done once, keeping the buffer and resources across frames

    # XXX: we need create a screen, which initializes a GL Context,
    # so that we can create a GLBuffer before having rendered anything.
    screen = GLMakie.global_gl_screen(resolution, true)

    # get a buffer object and register it with CUDA
    buffer = GLAbstraction.GLBuffer(Point2f, N)
    resource = let
        ref = Ref{CUDA.CUgraphicsResource}()
        CUDA.cuGraphicsGLRegisterBuffer(ref, buffer.id,
                                        CUDA.CU_GRAPHICS_MAP_RESOURCE_FLAGS_WRITE_DISCARD)
        ref[]
    end

    # NOTE: Makie's out-of-place API (lines, scatter) performs may iterating operations,
    # like determining the range of the data, so we use a manual scene instead.
    scene = Scene(; resolution)
    cam2d!(scene)

    # XXX: manually position the cameral (`center!` would iterate data)
    cam = Makie.camera(scene)
    cam.projection[] = Makie.orthographicprojection(
        #= x =# 0f0, 1f0,
        #= y =# 0f0, 1f0,
        #= z =# 0f0, 1f0)

    ## main processing
    #
    # this needs to be done every time we get new data and want to plot it

    NVTX.@range "main" begin
        # process data, generate points
        NVTX.@range "CUDA" begin
            # map OpenGL buffer object for writing from CUDA
            CUDA.cuGraphicsMapResources(1, [resource], stream())

            # get a CuArray object that we can work with
            array = let
                ptr_ref = Ref{CUDA.CUdeviceptr}()
                numbytes_ref = Ref{Csize_t}()
                CUDA.cuGraphicsResourceGetMappedPointer_v2(ptr_ref, numbytes_ref, resource)

                ptr = reinterpret(CuPtr{Point2f}, ptr_ref[])
                len = Int(numbytes_ref[] ÷ sizeof(Point2f))

                unsafe_wrap(CuArray, ptr, len)
            end

            # generate points
            broadcast!(array, t, X) do x, y
                Point2f(x, y)
            end

            # wait for the GPU to finish
            synchronize()

            CUDA.cuGraphicsUnmapResources(1, [resource], stream())
        end

        # generate and render plot
        NVTX.@range "Makie" begin
            scatter!(scene, buffer)

            # force everything to render (for benchmarking purposes)
            GLMakie.render_frame(screen, resize_buffers=false)
            GLMakie.glFinish()
        end

    end

    save("plot.png", scene)

    ## clean-up

    CUDA.cuGraphicsUnregisterResource(resource)

    return
end

```

This performs well, doing all the rendering in a couple of 100s of us. It requires [https://github.com/JuliaPlots/Makie.jl/pull/1803](https://github.com/JuliaPlots/Makie.jl/pull/1803), and I’d also recommend to disable `gpu_getindex` so that scalar iteration of GLBuffer (a performance trap) is disallowed:

```julia
@eval GLMakie.GLAbstraction begin
    # XXX: to make scalar iteration error
    function gpu_getindex(b::GLBuffer{T}, range::UnitRange) where T
        error("GLBuffer getindex") # XXX: for development
        multiplicator = sizeof(T)
        offset = first(range)-1
        value = Vector{T}(undef, length(range))
        bind(b)
        glGetBufferSubData(b.buffertype, multiplicator*offset, sizeof(value), value)
        bind(b, 0)
        return value
    end
end

```

It’s too bad that GLBuffer doesn’t support many array operations to make, e.g., `lines!` work properly. I wonder if it wouldn’t be better to pass `CuArray`s into Makie and have it do the necessary GL Interop calls automatically (which would make is possible to keep the `CuArray` around longer, and perform array operations on it, instead of eagerly converting it to a `GLBuffer`).

---

<div class="post-metadata">

**Author:** ![duanestorti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/duanestorti/32/37374_2.png) [@duanestorti](https://discourse.julialang.org/u/duanestorti)\
**Post date:** [July 26, 2022, 12:54am UTC](https://discourse.julialang.org/t/cuarray-glmakie/52461/11 "2022-07-26T00:54:39Z")

</div>

Any chance there has been further progress along this path? I am trying to implement CUDA/OpenGL graphics interop (i.e. to compute an array using CUDA and then render the array of data as a texture drawn on a rectangle in an OpenGL - or Makie - window).

Any chance you can provide more specific pointers on how to proceed?

For specific comparison code, a C version is available at:  
[Live Display via Graphics Interop | CUDA for Engineers: 2D Grids and Interactive Graphics | InformIT](https://www.informit.com/articles/article.aspx?p=2455391&seqNum=2)

---

<div class="post-metadata">

**Author:** ![waywardpidgeon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/waywardpidgeon/32/212368_2.png) [@waywardpidgeon](https://discourse.julialang.org/u/waywardpidgeon)\
**Post date:** [October 1, 2024, 9:01pm UTC](https://discourse.julialang.org/t/cuarray-glmakie/52461/12 "2024-10-01T21:01:08Z")

</div>

When attempting to run the nice example of LaurentPlagne on 2nd Oct 24 I first added AbstractPlotting and then rm it having see it is deprecated and not needed. I obtained the undefined error for Node with and without the call to CUDA allowscalar.

I am using Julia 1.10.5 in a win 11 Power Shell with NVIDIA+CUDA.jl  
Thanks - Kevin

---

<div class="post-metadata">

**Author:** ![sdanisch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdanisch/32/1406_2.png) [@sdanisch](https://discourse.julialang.org/u/sdanisch)\
**Post date:** [October 1, 2024, 9:07pm UTC](https://discourse.julialang.org/t/cuarray-glmakie/52461/13 "2024-10-01T21:07:33Z")

</div>

Can you try it with a new pkg env?  
And if that doesn’t work, show what’s in the environment?  
[https://pkgdocs.julialang.org/v1/environments/](https://pkgdocs.julialang.org/v1/environments/)

---

<div class="post-metadata">

**Author:** ![waywardpidgeon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/waywardpidgeon/32/212368_2.png) [@waywardpidgeon](https://discourse.julialang.org/u/waywardpidgeon)\
**Post date:** [October 2, 2024, 8:53pm UTC](https://discourse.julialang.org/t/cuarray-glmakie/52461/14 "2024-10-02T20:53:16Z")

</div>

Ok, I did that starting from scratch in the environments subfolder. Then ran the LaurentPlagne example and obtained the undefined “Node” error again.  
Here is the error and status:

julia\> @time go()  
ERROR: UndefVarError: `Node` not defined  
Stacktrace:  
[1] go()  
@ Main .\REPL[42]:13  
[2] macro expansion  
@ .\timing.jl:279 [inlined]  
[3] top-level scope  
@ .\REPL[48]:1

(@v1.10) pkg\> status  
Status `C:\Users\kab\.julia\environments\v1.10\Project.toml`  
[052768ef] CUDA v5.5.2  
[0376089a] ClimaOcean v0.1.3  
[e9467ef8] GLMakie v0.10.12  
⌅ [9e8cae18] Oceananigans v0.90.14  
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

I did not of course add AbstractPlotting. Actually examples from the Oceananigans set with CUDA and GLMakie, using GLMakie.record in place of record and surroundng the plotting instructions with “CUDA.allowscalars() do … end” did work (that was before my going astray with AbstractPlotting ).

---

<div class="post-metadata">

**Author:** ![sdanisch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sdanisch/32/1406_2.png) [@sdanisch](https://discourse.julialang.org/u/sdanisch)\
**Post date:** [October 3, 2024, 8:55am UTC](https://discourse.julialang.org/t/cuarray-glmakie/52461/15 "2024-10-03T08:55:57Z")

</div>

Ahh, sorry, when I read your post it looked like a new topic, so I didn’t realize you’re referring to an example just above.  
Node is now `Observable`.  
And there are a few other changes, but actually surprisingly little:

```julia
using CUDA,GLMakie
const T=Float32
const V=CuArray{T,2}
CUDA.allowscalar(false)
ls,ns,ms,λ,dt,sr,k,nt,two,eight,fr =1.e-2,800,1.e-4,0.1,2.e-7,0.1,5e7,5000,T(2),T(8),10
dt2sm=dt^2/ms
r,r0,r1,r2=1:ns,1:ns-2,2:ns-1,3:ns
cxs = [T(j-1)*ls for i in r, j in r]
dx(i,j) = cxs[i,j]-cxs[ns÷2,ns÷2]
cxc = [cxs[i,j]+sr*(dx(i,j))*exp(-0.5(dx(i,j)^2+dx(j,i)^2)/λ^2)/λ for i in r, j in r]
xs,xc,ys,yc=V(cxs),V(cxc),V(cxs'),V(cxc')
# ys,yc=collect(xs'),collect(xc')
xp,yp,xt,yt,fx,fy,zs=copy(xc),copy(yc),copy(xc),copy(yc),zero(xc),zero(xc),zero(xc)
@. zs = sqrt((xc-xs)^2+(yc-ys)^2) 
czs=Array(zs) ; zn=Observable(czs) ; zl=sr*0.05
scene = Scene(resolution = (400,400)) ; scale!(scene, 1, 1, 10000)
cam3d!(scene)
surface!(scene, 1:ns, 1:ns, zn, colorrange=(-zl, zl))

center!(scene)
display(scene)
for j in 1:nt÷fr
# for j in 1:nt÷fr
    for i in 1:fr
        @views @. fx[r1,r1] = -k*(eight*xc[r1,r1]-xc[r0,r0]-xc[r1,r0]-xc[r2,r0]-xc[r0,r1]-xc[r2,r1]-xc[r0,r2]-xc[r1,r2]-xc[r2,r2])
        @views @. fy[r1,r1] = -k*(eight*yc[r1,r1]-yc[r0,r0]-yc[r1,r0]-yc[r2,r0]-yc[r0,r1]-yc[r2,r1]-yc[r0,r2]-yc[r1,r2]-yc[r2,r2])
        @. xt = two*xc-xp+fx*dt2sm
        @. yt = two*yc-yp+fy*dt2sm
        xc,xp,xt = xt,xc,xp
        yc,yp,yt = yt,yc,yp
    end
    @. zs = sqrt((xc-xs)^2+(yc-ys)^2)
    czs=Array(zs)
    zn[] = czs
    yield()
end

```

And here a more modern version:

```julia
czs=Array(zs) ; zn=Observable(czs) ; zl=sr*0.05
f, ax, pl = surface(1:ns, 1:ns, zn, colorrange=(-zl, zl), axis=(; type=Axis3))

```

[https://github.com/user-attachments/assets/28fa8b53-0c07-4e7b-b8ae-53c7ca024028(image larger than 4 MB)](https://github.com/user-attachments/assets/28fa8b53-0c07-4e7b-b8ae-53c7ca024028)