Nicholas Cooley 2026-06-02
Infrastructure for calling NVIDIA’s CUDA framework from R
This package is no longer under development. The capabilities presented here are now supported in ardea, which has recently been submitted to CRAN.
A brief introduction to CUDA, and why access to CUDA in R shouldn’t be as hard as it is now…
There is a split between the runtime and the driver APIs, and that distinction is both confusing and important.
ACFcuda itself is relatively bare on R functions, with an emphasis of
providing C level functions that can be used to divert intensive compute
steps to non-CPU hardware. Library and function compilation, device
discovery, and context construction are exposed in R; while general
buffer management and function dispatch are only exposed in C. A simple
R wrapper is present, but it is almost solely overhead checking to
ensure argument compliance.
Users should expect to interact with CUDA devices through R’s .Call
with an emphasis being placed on maintability and reliability for the
time being.
A brief code example for package usage is included below.
# write out a kernel function
cuda_mm_naive <- '
#include <stdint.h>
extern "C" __global__ void cuda_mm_naive(float* output,
const uint32_t M,
const uint32_t K,
const uint32_t N,
const float* A,
const float* B)
{
uint32_t row = blockIdx.x * blockDim.x + threadIdx.x;
uint32_t col = blockIdx.y * blockDim.y + threadIdx.y;
if (row < M && col < N) {
float sum = 0.0f;
for (uint32_t k = 0; k < K; k++) {
sum += A[row + k * M] * B[k + col * K];
}
output[row + col * M] = sum;
}
}
'
# write out our kernel function to a temporary file
# and compile it to a library
tmp01 <- tempfile(fileext = ".cu")
writeLines(text = cuda_mm_naive,
con = tmp01)
tmp02 <- tempfile(fileext = ".ptx")library(ACFcuda)## ACFcuda 0.1.1 - CUDA GPU acceleration is available
library(microbenchmark)
# make device context
dvc_ctx <- cuda_make_context()
cuda_source_to_ptx(cuda_file = tmp01,
ptx_file = tmp02,
verbose = TRUE)## character(0)
fun_ptr <- cuda_kernel_from_ptx(cuda_context = dvc_ctx,
ptx_file = tmp02,
kernel_name = "cuda_mm_naive")set.seed(1986)
dim_set <- c(100, 250, 500, 1000, 2500)
res_vals <- vector(mode = "list",
length = length(dim_set))
fun_vec <- c("internal",
"fun01")
for (a1 in seq_along(dim_set)) {
dim_size <- dim_set[a1]
var1 <- matrix(data = runif(dim_size^2),
nrow = dim_size,
ncol = dim_size)
var2 <- matrix(data = runif(dim_size^2),
nrow = dim_size,
ncol = dim_size)
dummy_vec <- matrix(data = numeric(),
nrow = dim_size,
ncol = dim_size)
arg_types <- c("float",
"uint",
"uint",
"uint",
"float",
"float")
# we can cheat here because these are symetrical matrices
arg_list <- list(dummy_vec,
dim_size,
dim_size,
dim_size,
var1,
var2)
# suppress the nanosecond timings warning...
timing_res <- suppressWarnings(microbenchmark("internal" = var1 %*% var2,
"fun01" = simple_cuda_wrapper(context_pointer = dvc_ctx,
kernel_pointer = fun_ptr,
arg_types = arg_types,
arg_list = arg_list),
times = 5))
res_vals[[a1]] <- vapply(X = c("internal",
"fun01"),
FUN = function(x) {
mean(timing_res$time[timing_res$expr == x])
},
FUN.VALUE = vector(mode = "numeric",
length = 1))
}
res_vals <- do.call(rbind,
res_vals)
colset <- c("black",
"blue")
par(mar = c(3,2.5,2,1),
mgp = c(1.45, .55, 0))
plot(x = 1e-100,
y = 1e-100,
xlim = range(dim_set),
ylim = range(unlist(res_vals)),
xlab = "dim size (log)",
ylab = "nanoseconds (log)",
type = "n",
log = "xy",
main = "Matrix Multiplication")
for (a1 in seq(ncol(res_vals))) {
points(x = dim_set,
y = res_vals[, a1],
pch = 20,
col = colset[a1],
cex = 2)
}
legend("topleft",
legend = c("internal",
"cuda"),
col = colset,
pch = 20,
cex = 1.5)
As should be expected (because compute is taking place with different precisions), values returned by compute on a CUDA devices are not going to be exactly the same as those generated on the CPU, but typically deviations here are within generally acceptable tolerances. Checking them for personal preferences is rarely a poor choice though. Unique implementations of the same operations should return matching values however.
dim_size <- 100
var1 <- matrix(data = runif(dim_size^2),
nrow = dim_size,
ncol = dim_size)
var2 <- matrix(data = runif(dim_size^2),
nrow = dim_size,
ncol = dim_size)
dummy_vec <- matrix(data = numeric(),
nrow = dim_size,
ncol = dim_size)
arg_types <- c("float",
"uint",
"uint",
"uint",
"float",
"float")
# we can cheat here because these are symetrical matrices
arg_list <- list(dummy_vec,
dim_size,
dim_size,
dim_size,
var1,
var2)
res01 <- var1 %*% var2
res02 <- simple_cuda_wrapper(context_pointer = dvc_ctx,
kernel_pointer = fun_ptr,
arg_types = arg_types,
arg_list = arg_list)
layout(mat = matrix(data = 1:2,
nrow = 1))
plot(as.vector(res01),
as.vector(res02),
pch = 46,
xlab = "CPU result",
ylab = "GPU result",
main = "CPU vs GPU")
plot(as.vector(res01) - as.vector(res02),
pch = 46,
xlab = "index",
ylab = "delta between implementations",
main = "implementation 1 vs 2")