Author Archives: Tim Besard

cuTile.jl 1.0: Tile windows, atomics, and sparse views

By: Tim Besard

Re-posted from: https://juliagpu.org/post/2026-08-19-cutile_1.0/index.html

cuTile.jl has reached 1.0! The release adds tile windows via eachtile, masked and view-based atomics, an @atomic macro, sparse views, compiler remarks, and support for Tile IR 13.4, alongside a dedicated documentation site.

The package started as an experiment in expressing NVIDIA's tile-based programming model in Julia. Six months and three releases later, we're confident tagging an initial stable release. This comes with a dedicated documentation site.

Compared to v0.3, the 1.0 release adds the following features.

Tile windows

Previously, walking an array tile by tile meant passing the array, index, and shape to every ct.load and ct.store. The new ct.eachtile instead returns an indexable collection of fixed-shape windows. A blocked matrix multiplication shows the difference:

using CUDA, cuTile
import cuTile as ctfunction matmul!(C, A, B)
    a_tiles = ct.eachtile(A, (64, 32))
    b_tiles = ct.eachtile(B, (32, 64))
    c_tiles = ct.eachtile(C, (64, 64))    m, n = ct.bid(1), ct.bid(2)
    acc = zeros(Float32, (64, 64))
    for k in Int32(1):Int32(size(a_tiles, 2))
        acc = muladd(a_tiles[m, k], b_tiles[k, n], acc)
    end
    c_tiles[m, n] = acc
    return
endA = CUDA.rand(Float16, 256, 128)
B = CUDA.rand(Float16, 128, 256)
C = CUDA.zeros(Float32, 256, 256)
@cuda backend=cuTile blocks=(4, 4) matmul!(C, A, B)

size(a_tiles, 2) returns the number of windows along that dimension, avoiding a separate trip-count calculation. step controls the distance between window origins: a smaller value produces overlap, while a larger one leaves gaps.

adjacent = ct.eachtile(a, (8, 8))               # step defaults to the shape
overlap  = ct.eachtile(a, (8, 8); step=(4, 8))  # neighboring windows overlap

Partial edge windows are handled by the padding mode, as with a normal load. Unequal shape and step require Tile IR bytecode v13.3 or newer.

Atomics

cuTile now has three atomic-operation families, with different return values and ordering options.

The read-modify-write functions (ct.atomic_add and friends) return the old value and take a configurable memory order. In 1.0 they also accept a mask, useful for the tail block of a grid that does not divide the data evenly:

function histogram!(counts, data, n::Int32)
    pid = ct.bid(1)
    offs = (pid - Int32(1)) * Int32(128) .+ ct.arange(128)
    vals = ct.load(data; index=pid, shape=(128,))
    active = offs .<= n                    # mask off the tail
    ct.atomic_add(counts, vals, Int32(1); mask=active)
    return
end

The new ct.atomic_store_* family lowers to Tile IR's view-based atomic reductions. These reduce a tile into an array or an eachtile window and return nothing, using relaxed device-wide ordering:

function accumulate_tiles!(out, src)
    tiles = ct.eachtile(out, (128,))
    pid = ct.bid(1)
    ct.atomic_store_add(tiles, 1, ct.load(src; index=pid, shape=(128,)))
    return
end

ct.@atomic provides Base-style statement and value forms:

ct.@atomic counters[i] += update
ct.@atomic counters[i] = max(counters[i], value)
old_new = ct.@atomic counters[i] + value      # returns old => new

Statement forms default to relaxed ordering, while value forms default to acquire-release. View-based reductions and ct.@atomic require Tile IR 13.3.

Sparse views

view and @view on a TileArray now accept positive step ranges. On arrays with two or more dimensions, one dimension may instead use a 1D integer tile, creating a sparse view that ct.load and ct.store lower to a Tile IR gather/scatter view:

function pick_rows!(dst, src)
    rows = ct.arange(4; start=1, step=2)      # rows 1, 3, 5, 7
    selected = @view src[rows, 1:8]
    tile = ct.load(selected, (4, 8))
    ct.store(dst, (1, 1), tile)
    return
end

The load shape is explicit and static, while the range starts may be runtime values. Sparse loads apply the requested padding and stores clip partially out-of-bounds elements; repeated indices are fine for loads, but conflicting stores are undefined. Step ranges and sparse views require Tile IR 13.3.

Compiler remarks

tileiras can report whether it selected tensor cores, vector loads, and other optimizations. code_tiled and @device_code_tiled now print those diagnostics with remarks=true.

Compile the matmul above with Float32 inputs:

A = CUDA.rand(Float32, 256, 128)
B = CUDA.rand(Float32, 128, 256)
C = CUDA.zeros(Float32, 256, 256)
ct.@device_code_tiled remarks=true @cuda backend=cuTile blocks=(4, 4) matmul!(C, A, B)
// tileiras optimization remarks
// Name:            RemarkMemoryLoadInstructionSelected
//   - RemarkId:    3
//   - Remark:      Load instruction selected
// Name:            RemarkTensorCoreMMA
//   - RemarkId:    1
//   - Remark:      MMA operation failed to optimize to use Tensor Cores, it is using FMA instructions instead

For this kernel, the Float32 multiply uses FMA instructions. With Float16 inputs and a Float32 accumulator, the compiler instead reports:

// Name:            RemarkTensorCoreMMA
//   - Remark:      MMA operation successfully optimized to use Tensor Cores

Remarks require tileiras 13.4 or newer, which is still in early-access.

Programmatic dependent launch

Programmatic dependent launch can overlap the tail of a producer kernel with an independent preamble in the next kernel on the same stream. The producer signals when its dependents may start; the consumer is launched with dependent=true and waits before reading the producer's results:

function producer(a, producer_out)
    ct.grid_dependency_control_launch_dependents()
    tile = ct.load(a, 1, (32,))               # may overlap with the consumer
    ct.store(producer_out, 1, tile)
    return
endfunction consumer(b, producer_out, out)
    tile = ct.load(b, 1, (32,))               # independent preamble
    ct.grid_dependency_control_wait()
    ct.store(out, 1, tile + ct.load(producer_out, 1, (32,)))
    return
endstream = CUDA.stream()
@cuda backend=cuTile blocks=1 stream producer(a, producer_out)
@cuda backend=cuTile blocks=1 dependent=true stream consumer(b, producer_out, out)

The overlap is opportunistic, so correctness must never depend on the two kernels actually running concurrently. The feature requires Tile IR 13.4 and compute capability 9.0 or newer.

Other changes

  • Tile IR 13.4 is supported and emitted by default when tileiras accepts it. It brings ct.insert, the inverse of ct.extract, and check_bounds=false on ct.load/ct.store, an explicit promise that the whole tile is in bounds which drops the padding and selects Tile IR's unchecked encoding.

  • Explicit rounding modes on float-to-float conversions: Float32.(tile, RoundDown), with RoundNearest, RoundToZero, RoundDown, RoundUp and RoundNearestTiesAway. Supported modes and source/target pairs depend on the Tile IR version. Directed rounding generally requires 13.4; some conversions to Float8_E8M0FNU are available in 13.3.

  • 64-bit indexing. TileArray gained an index-type parameter, and arrays whose sizes or strides exceed the 32-bit range automatically use Int64. Smaller arrays continue to use Int32. ct.TileArray(a; index=Int64) selects wide indexing explicitly. Wide indexing requires Tile IR 13.3.

  • Array construction syntax works in kernels: [a, b, c], typed forms like Float32[a, b], bracket concatenation ([A; B], [a b; C], [A;;; B]) and cat(A, B...; dims). ct.cat was removed in favor of these.

  • Multi-dimensional reductions. dims now takes an integer, an iterable of integers, or :, for both tile-level and host-level ct.Tiled reductions.

  • Memory ordering on plain loads and stores, not just atomics: ct.load and ct.store accept memory_order and memory_scope.

  • Configuration moved to preferences. The JULIA_CUTILE_CACHE_DIR and JULIA_CUTILE_CACHE_SIZE environment variables were replaced by the disk_cache, cache_dir, and cache_size_bytes preferences. A new compiler_timeout_seconds preference bounds each tileiras invocation.

cuTile.jl 1.0 requires CUDA.jl 6.3. See NEWS.md for the user-facing release history and the release notes for the merged pull requests. Please file an issue if you run into a problem.

CUDA.jl 6.3: Compiler caching, a new cuDNN, and dependent launches

By: Tim Besard

Re-posted from: https://juliagpu.org/post/2026-08-19-cuda_6.3/index.html

CUDA.jl 6.3 features better integration with Julia's compiler caches, so that GPU-side inference done while a package precompiles survives across sessions. The cuDNN wrappers have been rebuilt on cuDNN 9's backend graph API, and there is also support for programmatic dependent launch.

Kernel compilation with CompilerCaching.jl

When you launch a kernel, CUDA.jl has to find the compiled code for it. Until now it kept that mapping itself: a dictionary per CUDA context, from (method instance, world age, compiler configuration) to a CuFunction. It worked, but it duplicated bookkeeping Julia already does for the same method instances, and the cached entries did not survive across Julia sessions.

CUDA.jl 6.3 adopts GPUCompiler 2, which builds on CompilerCaching.jl, and drops that dictionary. Compilation results are now stored in the CodeInstance that Julia caches anyway. This makes it possible to cache on disk, by saving into system or package images.

Right now, we only store inferred code. Work is underway to make the generated LLVM IR and machine code relocatable, which will enable caching those as well. However, just caching the inference results is already a big win. Let's demonstrate using a simple package:

module Blurusing CUDA
using PrecompileToolsfunction blur_kernel!(dst, src, ::Val{R}) where R
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    if i <= length(dst)
        acc = zero(eltype(src))
        for k in -R:R
            @inbounds acc += src[clamp(i + k, 1, length(src))] / (1 + abs(k))
        end
        @inbounds dst[i] = sqrt(abs(acc))
    end
    return
endfunction blur(src, ::Val{R} = Val(4)) where R
    dst = similar(src)
    @cuda threads=256 blocks=cld(length(dst), 256) blur_kernel!(dst, src, Val(R))
    return dst
end@setup_workload begin
    @compile_workload begin
        blur(CUDA.zeros(Float32, 1024))
    end
endend

Timing the first call in a fresh session, on an RTX 5080 with Julia 1.12 and CUDA 13.3:

julia> using Blur, CUDAjulia> src = CUDA.rand(Float32, 1024);julia> @time Blur.blur(src);
  0.129540 seconds (14.94 k allocations: 2.432 MiB, 71.87% compilation time: <1% of which was recompilation)

Delete the @setup_workload block, precompile again, and the same call in a fresh session costs this instead:

julia> @time Blur.blur(src);
  1.292274 seconds (5.43 M allocations: 262.556 MiB, 14.88% gc time, 95.62% compilation time: 10% of which was recompilation)

Note that this requires Julia 1.11 or later.

cuDNN, rebuilt on the graph API

cuDNN has two programming models: the legacy API, a fixed set of fixed-function operations and fusion patterns with a C entry point each, and the graph API, where you describe a computation as a graph of operations and let cuDNN pick an engine for the whole thing. The graph API can be reached two ways: directly through the C back-end API, or through NVIDIA's cudnn-frontend, whose C++ and Python layers provide a simplified programming model that covers most use cases.

cuDNN.jl was written against the fixed-function API. In version 6.3, it is rebuilt on the back-end API, with a front-end mimicking cudnn-frontend: a graph API and a set of operations implemented on top of it. The fixed-function wrappers are unchanged and remain available.

Graph front-end

Graph and Tensor describe a computation, build! lowers it, runs cuDNN's heuristics and selects an execution plan, and execute! binds arrays and runs it. Intermediate tensors are marked virtual, which is how the engine knows it may fuse instead of materializing them. As an example, a batched matrix multiply followed by a bias add and a ReLU, as one plan:

using CUDA, cuDNN
using cuDNN: Graph, tensor!, matmul!, pointwise!, build!, execute!A    = CUDA.rand(Float16, 256, 256, 8)
B    = CUDA.rand(Float16, 256, 256, 8)
bias = CUDA.rand(Float16, 256, 1, 8)
C    = CUDA.zeros(Float16, 256, 256, 8)g = Graph(io_dtype=Float16, intermediate_dtype=Float32, compute_dtype=Float32)
ta, tb = tensor!(g, A; name="A"), tensor!(g, B; name="B")
tbias  = tensor!(g, bias; name="Bias")
tc     = tensor!(g, C; name="C")tmm  = matmul!(g, ta, tb; name="MM")    # virtual
tsum = pointwise!(g, :add, tmm, tbias)  # virtual
pointwise!(g, :relu, tsum; y=tc)        # writes Cbuild!(g)
execute!(g, Dict(ta => A, tb => B, tbias => bias, tc => C))

Operations layer

On top of the frontend sits a higher-level API that's easier to use: attention! and attention_backward!, convolution! with its two gradients, maxpool!/meanpool! and their gradients, and the batchnorm_* family. These take CuArrays in Julia memory order and hide the graph entirely.

Both of these APIs are very new, and minor changes to the design or implementation are to be expected in future releases. Feedback is very welcome, so please report issues or missing features on the CUDA.jl bug tracker.

Programmatic dependent launch

Two kernels back-to-back in the same stream are fully serialized: the second one does not start until the last block of the first one retires. That is often more ordering than needed. If the consumer starts with work that does not touch the producer's output, like loading weights or zeroing an accumulator, that work could have been running while the producer's last few blocks were still draining.

CUDA calls the escape hatch programmatic dependent launch, and CUDA.jl 6.3 supports it. The producer signals when its dependents may start, the consumer is launched with dependent=true, and the consumer waits before it touches anything the producer wrote:

@inline function busy(x::Float32, n::Int)   # stand-in for real work
    for _ in 1:n
        x = fma(x, 1.0000001f0, 1f-7)
    end
    return x
endfunction producer!(out, n)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    trigger_programmatic_launch_completion()
    @inbounds out[i] = busy(Float32(i), n)   # the tail
    return
endfunction consumer!(out, in, n)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    pre = busy(Float32(i) * 0.5f0, n)        # independent preamble
    grid_dependency_synchronize()
    @inbounds out[i] = pre + in[i]
    return
end@cuda threads=256 blocks=32 producer!(a, 20_000)
@cuda threads=256 blocks=32 dependent=true consumer!(b, a, 20_000)

Each of these kernels takes about 36 µs on its own, and the grid is small enough that both fit on the device at once. Run back to back they cost 69 µs; with the trigger, the wait and dependent=true they cost 38 µs, so the consumer's preamble hides almost entirely behind the producer.

The trigger belongs at the point in the producer after which nothing else has to run before dependents may start, which is usually the top; a block that exits without calling it triggers completion implicitly. grid_dependency_synchronize is what makes the producer's writes visible, so the consumer needs it even when the trigger has already run. And the overlap is opportunistic: code whose correctness depends on the two kernels running concurrently can deadlock. Programmatic dependent launch requires compute capability 9.0 or higher.

Other changes

  • Support for CUDA 13.4. Since this version is still in early-access, it needs explicit opt-in by calling CUDA.set_runtime_version! or by configuring LocalPreferences.toml.

  • cuTENSOR.jl has been updated to cuTENSOR v2.7. Block-sparse contract! and plan_contraction take a reproducible keyword argument for bitwise reproducible contractions, and the compute-descriptor list gained the 16BF and FP-emulation descriptors that Hopper and Blackwell use.

  • There is a low-level API for conversion-free launches: KernelCall converts a kernel's function and arguments once, kernel_compile compiles the call, kernel_launch launches it without converting again, and rebind replaces a single argument. The KernelAbstractions back-end uses it when selecting a workgroup size, which removes a second conversion from operations such as broadcast.=

The full list is in NEWS.md and the release notes. If something in here breaks for you, please file an issue.

cuTile.jl 0.3: CUDA.jl integration, and even better performance &amp; latency

By: Tim Besard

Re-posted from: https://juliagpu.org/post/2026-05-05-cutile_0.3/index.html

cuTile.jl v0.3 integrates with CUDA.jl, making it even easier to write and run CUDA Tile kernels in Julia. Performance has also been greatly improved, closing the gap with cuTile Python on every benchmark we ship. Added features include a random number generator, and support for array slicing.

Performance: matching cuTile Python

Three months ago, several of our benchmarks lagged cuTile Python by 5–15%. Today, cuTile.jl matches or outperforms cuTile Python on every kernel we ship. The headline numbers (RTX 5080, tileiras 13.2.51):

Kernel Julia Python Δ
Vector Addition 845 GB/s 846 GB/s =
Matrix Transpose 812 GB/s 814 GB/s =
Layer Norm fwd 983 GB/s 716 GB/s +37%
Layer Norm bwd 248 GB/s 251 GB/s -1%
Matrix Multiplication 47.5 TFLOPS 43.5 TFLOPS +9%
Batch Matrix Multiply 34.0 TFLOPS 30.8 TFLOPS +10%
FFT (3-stage Cooley-Tukey) 529 μs 554 μs +5%
Mixture of Experts 27.0 TFLOPS 20.1 TFLOPS +34%
Attention (FMHA, causal) 103.6 TFLOPS 63.4 TFLOPS +63%
Softmax (TMA) 849 GB/s 857 GB/s -1%
Softmax (Chunked) 1684 GB/s 1640 GB/s +3%

Most of the gains come from extending the IR-level optimization pipeline introduced in v0.2 with a new dataflow framework that now powers several analyses and transformations.

CUDA.jl integration: @cuda backend=cuTile

Until v0.3, launching a cuTile kernel meant calling cuTile.launch(...) directly. cuTile.jl now plugs into CUDA.jl's existing @cuda macro as a first-class backend, making it much easier to launch cuTile.jl kernels:

using CUDA, cuTile
import cuTile as ctfunction vadd(a::ct.TileArray{Float32,1}, b::ct.TileArray{Float32,1},
              c::ct.TileArray{Float32,1})
    pid = ct.bid(1)
    ct.store(c; index=pid, tile=ct.load(a; index=pid, shape=(128,)) +
                                ct.load(b; index=pid, shape=(128,)))
    return
enda = CUDA.rand(Float32, 1024)
b = CUDA.rand(Float32, 1024)
c = CUDA.zeros(Float32, 1024)@cuda backend=cuTile blocks=8 vadd(a, b, c)

Time-to-first-launch

Compiling a cuTile kernel goes through several stages: Julia type inference, our IR rewriting passes, Tile IR bytecode emission, and finally tileiras-driven CUBIN generation. None of these are fast. Significant effort in v0.3 went into reducing the time-to-first-launch, and the latency is now comparable to a typical CUDA.jl kernel launch on the same hardware:

Benchmark 1: julia -e 'using CUDACore;
                       @cuda identity(nothing)'
  Time (mean ± σ):      1.882 s ±  0.012 s    [User: 2.554 s, System: 0.305 s]
  Range (min … max):    1.867 s …  1.906 s    10 runsBenchmark 2: julia -e 'using CUDACore, cuTile;
                       @cuda backend=cuTile identity(nothing)'
  Time (mean ± σ):      1.840 s ±  0.009 s    [User: 2.488 s, System: 0.329 s]
  Range (min … max):    1.827 s …  1.859 s    10 runs

Array slicing

view and @view now derive sub-range TileArrays from existing ones:

function copy_rows!(A::ct.TileArray{Float32,2}, B::ct.TileArray{Float32,2},
                    i::Int32, j::Int32)
    sub = @view A[i:j, :]                         # sub-range TileArray
    t = ct.load(sub; index=(1, 1), shape=(8, 8))
    ct.store(B; index=(1, 1), tile=t)
    return
end@cuda backend=cuTile copy_rows!(A, B, Int32(3), Int32(10))

Each index must be : or a UnitRange; other forms (StepRange, scalar indexes, CartesianIndex, …) are currently rejected at compile time. The result is itself a TileArray, and can be passed to ct.load / ct.store (or sliced again, for nested views). The new divisibility analysis sees through the slicing chain so contiguous-axis fast paths are preserved, while literal slice sizes fold to compile-time-constant shape operands.

Random number generation

cuTile.jl now ships a tile-vectorized Philox2x32-7 RNG, both as in-kernel intrinsics and as a host-side cuTile.RNG handle for filling CuArrays. The kernel API mirrors Base.Random:

function noise!(out::ct.TileArray{Float32,1})
    pid = ct.bid(1)
    t = randn(Float32, (128,))                 # in-kernel randn
    ct.store(out; index=pid, tile=t)
    return
end@cuda backend=cuTile blocks=cld(N, 128) noise!(A)

rand covers all of Int{8,16,32,64}, UInt{8,16,32,64}, Float16, BFloat16, Float32, and Float64; randn (via Box-Muller, sharing its uniforms with the existing rand path) and randexp (via -log(U)) cover the four floating-point types. ct.DeviceRNG() opens an independent stream inside a kernel; Random.seed! re-seeds.

The host-side cuTile.RNG integrates with Random.rand! / Random.randn! / Random.randexp! and auto-advances its counter, so consecutive fills produce disjoint streams:

A = CUDACore.zeros(Float32, 1 << 20)
rng = ct.RNG(42)
randn!(rng, A)                                 # fill via fused tile kernel
B = rand(rng, Float64, 16)                     # out-of-place

Performance of both the in-kernel and host-side APIs is excellent, matching or exceeding the performance of cuRAND and GPUArrays.jl' new generator.

What's next

If you've been watching cuTile.jl from a distance: now's a good time to try it out: add cuTile from the Julia REPL, or grab the examples to see how the moving parts fit together.

There is a webinar scheduled on May 12, 2026 at 1 PM ET, where Tim Besard (JuliaHub) and Andy Terrel (NVIDIA) will present cuTile.jl in a joint webinar, covering the design of CUDA Tile, how cuTile.jl is built, and several relevant examples. Click here to sign up.