GLASS: GPU Linear Algebra Simple Subroutines

GLASS is a comprehensive, header-only CUDA C++ __device__ template library for block-local linear algebra on GPUs. It is the foundational linear-algebra layer underneath GRiD, MPCGPU, GATO, HJCD-IK, and other A2R Lab GPU solvers.

Like Eigen on the CPU, GLASS aims to be comprehensive — it covers all block-local linear algebra in one consistent __device__ calling convention: BLAS (L1/L2/L3), LAPACK-style factorizations and triangular solves (Cholesky, LDLᵀ, and LU/QR via vendor backends), dense linear-system solvers (posv / ldlt / gesv), related algorithms — block-tridiagonal bdmv / pcg for trajectory optimization and MPC, plus a contraction-parallel + fused family — and robotics-specialized operators: Featherstone spatial 6-D cross products, coordinate transforms, and the 10-parameter inertia, the SO(3)/SE(3)/quaternion Lie family with its derivative chain and pose-error metrics, cone/augmented-Lagrangian projections, sphere-collision distance primitives, the 3x3 estimation kit (eig3/svd3/closest_rotation), and the sampling-planner softmax/argmin reductions (see Robotics operators & conventions). Everything runs inside one CUDA block — and you choose the granularity: the same operations exist as block-, warp-, or thread-scoped primitives (plus vendor-backed kernels), so a block can own one problem, pack one per warp, or pack 32 per warp with one problem per thread — whatever matches your problem size and batch count.

Interfaces

GLASS exposes four primary interfaces — pick the one that matches how your problem maps onto the GPU. They cover the same operations under one calling convention, so you switch between them by changing the namespace prefix. The ladder runs most→least problem packing: thread → warp → block → nvidia.

Block — glass::block::

The default. One block per problem; the block’s threads cooperate over shared/global data. Pure SIMT, no dependencies (#include "glass.cuh"). The contract tier: bit-exact, thread-count invariant, never re-dispatched. Choose this for a single moderate-to-large problem per block.

Library Overview

Warp — glass::warp::

One warp per problem (__shfl_*_sync, no __syncthreads), so warps run independently. Choose this to pack many small independent problems into one block for intra-block parallelism. Requires a full 32-lane warp.

Warp-scoped operations (glass::warp::)

Thread — glass::thread::

One problem per thread, 32 packed per warp — the low-DOF corner (compile-time sizes, register-resident up to N≤7). Sequential: no barriers, no shuffles. Choose this when a warp-per-problem factor at N≲7 would leave most lanes idle.

Thread-scoped operations (glass::thread::)

Nvidia — glass::nvidia::block::

CUB / cuBLASDx / cuSOLVERDx, auto-dispatched against SIMT by size — plus glass::nvidia::warp:: CUB WarpReduce L1 reductions (one full 32-lane warp per problem). Choose this when a vendor tensor-core kernel wins at your size (needs NVIDIA MathDx).

Backend Dispatch

Note

glass::cgrps:: is a convenience cooperative-groups alias of the Block interface — identical numerics (the same SIMT loop, indexed via a thread_group handle), for callers already in a cooperative-groups context or tiling arbitrary sub-block groups. It is not a separately-tuned backend. #include "glass-cgrps.cuh".

Note

Bare glass::op (and bare glass::nvidia::op) is the measured-default face: the same block-scope calling contract, with the implementation body chosen per (op, size, dtype) by glass::dispatch_body() (glass-dispatch.cuh, regenerated by the measured tune.py --legs body sweep). Cells with a robust measured win route to a warp- or thread-body inside the block; every other name is the same entity as glass::block::, and all pre-restructure spellings compile unchanged. Determinism-sensitive callers pin glass::block:: explicitly; see Namespaces, suffixes, and flags.

Performance

The four interfaces are numerically interchangeable, so GLASS can pick the fastest one per (operation, size, dtype) from a measured ladder: glass::suggested_backend<op, N, T>() returns the winning interface and a launch config for codegen and host-side dispatch. The shipped defaults are tuned on an RTX 5090 (sm_120); you can regenerate the table for your own GPU with the GLASS autotune workflow. The sm_120 tables now include the thread interface (2026-07-18 sweep): it wins the low-DOF corner of every operation except gemm — up to 7.5× on the small-N factor/solve chain in f64. See Tuning for Your Hardware for how the benchmarks drive the defaults, and the measured ladders at the bottom of this page.

The contraction-parallel + fused family (gemm_reduced / gemv_reduced / syrk_reduced, the tensor_vec_contract / vec_tensor_vec / congruence_sym / bilinear ops, riccati_gain, and compile-out robustness flags on potrf / ldlt / posv) is there for expressiveness and fusion only: on a quiet GPU the *_reduced decomposition measured slower than the plain serial in-thread contraction in 48 of 48 swept shapes (±5% margin, bench/RESULTS.md), and the suggested_use_reduced<> picker now declines it everywhere — prefer the plain ops for throughput. See Contraction-parallel ops (the *_reduced family) for the measurement and Namespaces, suffixes, and flags for the naming convention.

Get started

Header-only install, the single-block execution model, and an optional MathDx setup for the glass::nvidia:: backend.

Installation

API reference

The L1 / L2 / L3 and NVIDIA device functions, generated from the header doc-comments via Doxygen + Breathe.

API Reference

Quick start

#include "glass.cuh"

// One block solves one problem; threads stride over the data.
__global__ void saxpy_kernel(uint32_t n, float a, float *x, float *y) {
    glass::axpy(n, a, x, y);          // y = a*x + y
}

saxpy_kernel<<<1, 256>>>(n, 2.0f, d_x, d_y);

See Quickstart for a complete, compilable example, and Worked Examples for a worked program per concept.

Measured performance

The measured warp / block / nvidia ladder on an RTX 5090 (sm_120) — each op’s fastest interface across problem size, in ns/problem (the data behind glass::suggested_backend<>), shown here in the NPROB=8192 throughput regime:

GLASS also beats the standard host-batched recipe at robot sizes: against cublasGemmStridedBatched / cusolverDnPotrfBatched, gemm at N ≤ 24 and the factor-and-solve chain through N = 64 win at every batch size tested (up to 6.3× at saturation) — including with TF32 tensor cores permitted, which cuBLAS declines to engage below N = 24 anyway.

See Backend Sweep Results for the same ladder across the NPROB=64 / 1024 / 8192 batch regimes (the winner shifts with batch size), the host-batched cuBLAS/cuSOLVER and TF32 comparison, the fused riccati_gain case study, and the per-(op, N) winner table; see Tuning for Your Hardware to regenerate everything for your own GPU with bench/tune.py.