ardea: introduction

Nicholas Cooley

2026-10-08

Introduction

ardea provides access to GPU devices and frameworks natively in R. It currently provides rudimentary access to OpenCL, CUDA, and Metal; letting users compile, test, and implement their own device specific code. ardea currently only requires existing OpenCL to install appropriately, but will selectively build out access to frameworks and devices that it can detect during package installation.

This package currently does not provide an interface for vendor supplied libraries. Functionality is currently focused on ingesting a user supplied framework specific function, compiling that function and preparing it for dispatch, executing function dispatch and buffer management, and returning that function’s result to R.

library(ardea)

opencl_is_available()
cuda_is_available()
metal_is_available()
## [1] TRUE
## [1] FALSE
## [1] TRUE

Using matrix multiplication as an example, a user can bring a naive mm kernel function conforming to the a given framework’s rules and requirements. In this case we can just write out our functions to temp files, and walk them forward through the frameworks and devices available on the current system. One of ardea’s current quirks is that it expects the first argument of a function to be the output. There is no keyword or tooling to enforce this, but the underlying buffer management currently only returns that first input from the device to R.

# write out a character vector to a tempfile for opencl, metal, and cuda
opencl_mm_naive <- '
__kernel void opencl_mm_naive(__global float* output,
                              const long M,
                              const long K,
                              const long N,
                              __global const float* A,
                              __global const float* B)
{
    long row = get_global_id(0);
    long col = get_global_id(1);

    if (row < M && col < N) {
        float sum = 0.0f;
        for (long k = 0; k < K; k++) {
            sum += A[row + k * M] * B[k + col * K];
        }
        output[row + col * M] = sum;
    }
}
'
tmp01 <- tempfile(fileext = ".cl")
writeLines(text = opencl_mm_naive,
           con = tmp01)

cuda_mm_naive <- '
extern "C" __global__ void cuda_mm_naive(float* output,
                                         const long long M,
                                         const long long K,
                                         const long long N,
                                         const float* A,
                                         const float* B)
{
    long long row = blockIdx.x * blockDim.x + threadIdx.x;
    long long col = blockIdx.y * blockDim.y + threadIdx.y;

    if (row < M && col < N) {
        float sum = 0.0f;
        for (long long k = 0; k < K; k++) {
            sum += A[row + k * M] * B[k + col * K];
        }
        output[row + col * M] = sum;
    }
}
'
tmp02 <- tempfile(fileext = ".cu")
tmp03 <- tempfile(fileext = ".ptx")
writeLines(text = cuda_mm_naive,
           con = tmp02)

metal_mm_naive <- '
kernel void metal_mm_naive(device float* output [[buffer(0)]],
                           constant uint& M [[buffer(1)]],
                           constant uint& K [[buffer(2)]],
                           constant uint& N [[buffer(3)]],
                           device const float* A [[buffer(4)]],
                           device const float* B [[buffer(5)]],
                           uint2 id [[thread_position_in_grid]])
{
    uint row = id.x;
    uint col = id.y;

    if (row < M && col < N) {
        float sum = 0.0f;
        for (uint k = 0; k < K; k++) {
            sum += A[row + k * M] * B[k + col * K];
        }
        output[row + col * M] = sum;
    }
}
'
tmp04 <- tempfile(fileext = ".metal")
tmp05 <- tempfile(fileext = ".metallib")
writeLines(text = metal_mm_naive,
           con = tmp04)

Matrix multiplication doesn’t need a lot of setup when using the builtin %*% function, which is calling on your R installation’s BLAS library and is highly optimized. However, because we’re reaching out to a device almost entirely ab initio with no helpers, we need to perform the work that those helpers would perform.

# generic inputs, regardless of framework
var00 <- vector(mode = "numeric",
               length = 250000)
var01 <- matrix(rnorm(250000),
               nrow = 500,
               ncol = 500)
var02 <- matrix(rnorm(250000),
               nrow = 500,
               ncol = 500)
dim01 <- nrow(var01)
dim02 <- ncol(var01)
# dim03 <- nrow(var02)
dim04 <- ncol(var02)

arg_list <- list(var00, # our expected output vector, we've just filled it with zeros
                 dim01,
                 dim02,
                 dim04,
                 as.vector(var01),
                 as.vector(var02))

Every framework has its own unique verbiage and intricacies, though they’re all kind of doing similar things. OpenCL is designed to be be deployable across a variety of device vendors and device types, not just GPUs. Metal is Apple’s proprietary GPU / accelerator framework, though it shares some conceptual choices with OpenCL as until recently Apple’s hardware stack was somewhat diverse. CUDA is NVIDIA’s well known framework for accessing their devices, and includes an intimidating array of libraries and functionalities.

This means that frameworks share concepts, but not necessarily the same keywords, type names, specific paths, or user interactions.

opencl_arg_types <- c("float",
                      "long",
                      "long",
                      "long",
                      "float",
                      "float")
metal_arg_types = c("float",
                    "uint",
                    "uint",
                    "uint",
                    "float",
                    "float")
cuda_arg_types <- c("float",
                    "long",
                    "long",
                    "long",
                    "float",
                    "float")

Though each framework handles program / module / library preparation somewhat uniquely, an effort has been made to unify how users interact with these steps. In this package OpenCL only supports API-side compilation and program generation, while Metal supports both API-side, and system-side compilation and program management, and CUDA only supports system-side compilation. Meaning that under the hood Metal (optionally) and CUDA are calling on system/vendor-specific compilers to build their libraries/modules (i.e. xcrun metal and nvcc respectively).

if (opencl_is_available() &
    opencl_devices_exist()) {
  print("OpenCL is available on this system!")
  cl_dvcs <- opencl_device_information()
  cl_ctx <- opencl_make_context(device = cl_dvcs[[1]])
  cl_program <- opencl_make_program(cl_file = tmp01,
                                    context = cl_ctx)
  cl_knl <- opencl_make_kernelptr(program = cl_program,
                                  kernel_names = "opencl_mm_naive")
}
if (cuda_is_available() &
    cuda_devices_exist()) {
  print("CUDA is available on this system!")
  cu_dvcs <- cuda_device_information()
  cu_ctx <- cuda_make_context(device = cu_dvcs[[1]])
  cu_program <- cuda_make_program(cuda_file = tmp02,
                                  context = cu_ctx,
                                  ptx_file = tmp03)
  cu_knl <- cuda_make_kernelptr(program = cu_program,
                                context = cu_ctx,
                                kernel_names = "cuda_mm_naive")
}
if (metal_is_available() &
    metal_devices_exist()) {
  print("Metal is available on this system!")
  mtl_dvcs <- metal_device_information()
  mtl_ctx <- metal_make_context(device = mtl_dvcs[[1]])
  if (metal_compiler_is_available()) {
    mtl_program <- metal_make_program(metal_file = tmp04,
                                      metallib_file = tmp05,
                                      context = mtl_ctx)
  } else {
    mtl_program <- metal_make_program(metal_file = tmp04,
                                      context = mtl_ctx)
  }
  
  mtl_knl <- metal_make_kernelptr(program = mtl_program,
                                  context = mtl_ctx,
                                  kernel_names = "metal_mm_naive")
}
## [1] "OpenCL is available on this system!"
## [1] "Metal is available on this system!"

A simple wrapper function is supplied with ardea, though it is mostly a thin wrapper for .Call() runner functions. It does perform some error checking, but for the most part users are being trusted with their judgements.

# our builtin optimized CPU implementation
system.time(res01 <- var01 %*% var02) 
if (opencl_is_available() &
    opencl_devices_exist()) {
  print("OpenCL implementation:")
  print(system.time(res02 <- simple_wrapper(framework = "opencl",
                                            context_ptr = cl_ctx,
                                            kernel_ptr = cl_knl$opencl_mm_naive,
                                            arg_types = opencl_arg_types,
                                            arg_list = arg_list,
                                            problem_dims = as.integer(c(dim01,
                                                                        dim04,
                                                                        1L)),
                                            group_dims = NULL,
                                            workers_per = NULL)))
  plot(as.vector(res01),
       res02,
       pch = 46,
       main = "opencl vs builtin values")
}

if (cuda_is_available() &
    cuda_devices_exist()) {
  print("CUDA implementation:")
  print(system.time(res03 <- simple_wrapper(framework = "cuda",
                                            context_ptr = cu_ctx,
                                            kernel_ptr = cu_knl$cuda_mm_naive,
                                            arg_types = cuda_arg_types,
                                            arg_list = arg_list,
                                            problem_dims = as.integer(c(dim01,
                                                                        dim04,
                                                                        1L)),
                                            group_dims = NULL,
                                            workers_per = NULL)))
  plot(as.vector(res01),
       res03,
       pch = 46,
       main = "cuda vs builtin values")
}
if (metal_is_available() &
    metal_devices_exist()) {
  print("Metal implementation:")
  print(system.time(res03 <- simple_wrapper(framework = "metal",
                                            context_ptr = mtl_ctx,
                                            kernel_ptr = mtl_knl$metal_mm_naive,
                                            arg_types = metal_arg_types,
                                            arg_list = arg_list,
                                            problem_dims = as.integer(c(dim01,
                                                                        dim04,
                                                                        1L)),
                                            group_dims = NULL,
                                            workers_per = NULL)))
  plot(as.vector(res01),
       res03,
       pch = 46,
       main = "metal vs builtin values")
}

##    user  system elapsed 
##   0.032   0.001   0.033 
## [1] "OpenCL implementation:"
##    user  system elapsed 
##   0.001   0.003   0.017 
## [1] "Metal implementation:"
##    user  system elapsed 
##   0.001   0.001   0.004

This is currently the limit of ardea’s functionality. This package began mostly as a curiosity project with Metal, and turning it into a CRAN acceptable package became a bit of a labor of love. If folks find this package useful I will consider expanding it to include other vendors, or more complicated dispatch routines. As it stands now, my primary concerns are conforming to CRAN submission guidelines, ensuring user interactions are clear and predictable, and working through a non-trivial TODO list.

sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Sequoia 15.7.9
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ardea_0.1.1
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.39   R6_2.6.1        fastmap_1.2.0   xfun_0.61      
##  [5] cachem_1.1.0    knitr_1.52      htmltools_0.5.9 rmarkdown_2.32 
##  [9] lifecycle_1.0.5 cli_3.6.6       sass_0.4.10     jquerylib_0.1.4
## [13] compiler_4.6.1  tools_4.6.1     evaluate_1.0.5  bslib_0.12.0   
## [17] yaml_2.3.12     rlang_1.3.0     jsonlite_2.0.0