CUDACaml

About the project

CUDACaml is a compiler. You write ordinary OCaml — map, reduce, scan, matmul over arrays — and it produces CUDA C, compiles it at run time with NVRTC, and runs it on your GPU. You never write a kernel, choose a launch geometry, or decide by hand which operations to fuse. Five chained array operations become one kernel with no intermediate memory, because a fusion analysis decided that, not you.

That part is ordinary. Here is the part that isn't.

A wrong GPU kernel doesn't crash. It returns a number. Every optimisation a compiler performs — fusing two loops, keeping a value in a register, reassociating arithmetic — is a chance to change the answer, and the code that changed it is precisely the code that was supposed to be equivalent. Normally nothing checks you: there is no second implementation of your program for the first one to disagree with. A float is a float. It comes back either way.

So CUDACaml has one graph and two executors. Every program is a typed, immutable DAG, and that DAG can be run two completely different ways: through the optimising CUDA path, or through a pure-OCaml reference interpreter that runs no passes at all and shares no code with the path it is judging. Differential.check runs both and compares them elementwise. It works on any graph you write, not just our examples.

This caught a real bug, during the hackathon. Signed integer overflow is undefined behaviour in C++, and nvcc exploits it. Our emitter was printing integer arithmetic as plain signed C, so given i * 0xD2511F53, nvcc assumed the product could not go negative and folded a later < 0 test to a constant. Our random number generator — Philox-4x32-10, written in the DSL itself — overflows deliberately, every round. So every random number drawn on the GPU was wrong: uniforms outside [0,1), normals returning NaN, and a Black-Scholes Monte Carlo price reading 6.0461 against a correct 9.64478. Nothing crashed. Nothing warned. The only thing that noticed was that the interpreter disagreed with the GPU.

What's in it

• Fusion as an analysis, not a rewrite. A node owns a buffer only if it is a param, a reduction, a graph output, or read more than once. Everything else is inlined into its consumer.
• Reverse-mode automatic differentiation as a graph-to-graph transform. Adjoints are built out of ordinary DSL calls, so gradients fuse and lower exactly like forward code — no tape, no second runtime. A Black-Scholes Monte Carlo yields price, delta, vega and rho from one backward pass instead of one full repricing per Greek: 0.7 ms on an A100, about 4× the cost of the price alone.
• A Philox RNG written in the DSL, so a random tensor is a fusable subgraph whose only parameter is the seed — reproducible, and independent of how work is blocked across threads.
• Streams and events, pinned host memory, resident device values, per-device contexts and a data-parallel multi-GPU driver, and a Longstaff-Schwartz American option pricer built on top.

Why OCaml

The graph is an immutable value, not a sequence of side effects — so fusion, AD and lowering are three transformations over one value, and they compose because none of them mutates it. Algebraic data types with exhaustive matching mean adding a node kind turns every unhandled site into a compile error rather than a wrong number at 2am. And purity is what makes the oracle free: same graph, two interpretations.

Honest counterpoint: OCaml is not the obvious choice for kernel-level performance work, and CUDACaml still emits CUDA C and hands it to nvcc. The claim is about the layer that decides what to run, not about beating hand-tuned CUDA at running it.

Where it stands

Verified on two cards — an RTX 5080 Laptop (sm_120) and an A100-SXM4-40GB (sm_80, two-GPU box):

• 38 end-to-end tests, 0 failures, on both cards. 231 unit tests on the A100 with zero skips; the laptop skips only the devices >= 2 multi-GPU half it cannot run.
• All 13 bundled examples differential-tested green against the interpreter.
• Greeks agree with Black-Scholes closed form; the LSM American put lands within 0.071 of a 500-step binomial; Rng.u32 is bit-exact against the interpreter over 2²⁰ draws.

And an honest pair of throughput numbers, because one of them is unflattering: on saxpy we are 0.36× — slower than a plain OCaml loop, because at three flops per element it is all PCIe traffic. On a degree-64 polynomial with the same memory traffic and real arithmetic, 34×. The value is in the compiler layer, and low arithmetic intensity is where that layer cannot help you.

Not production-ready: no nested parallelism, dynamic shapes, bool tensors or autotuning. In production you would run CUDA alone — the interpreter is a deploy-time gate, not a runtime one, and because shapes are static per graph the set of graphs you ship is finite and can all be checked before it goes out.

One graph, two executors, and a reason to believe the number.

Built With

• compiler
• cuda
• ocaml