Worked Examples#
Minimal, self-contained compile-and-run programs — one concept each. Every
file is a complete program: a __global__ kernel that calls a GLASS device
function, plus a main that allocates device memory, launches one block
(or a batch of blocks for the batched demos), copies the result back, and
verifies it — examples return non-zero on a numeric mismatch, and the test
suite compiles and runs every one of them on hardware. The sources live in
examples/ in the repository; examples/README.md carries the same table
plus build instructions (keep the two in sync).
All examples are pure SIMT (no external dependencies) except 05_nvidia_gemm,
which requires MathDx.
Example |
Shows |
Backend / deps |
|---|---|---|
|
L1 vector op |
pure SIMT |
|
THE GEMM example: standard-BLAS convention, both size overloads, all
four transpose combos, row-major-is-a-transpose (bit-identical), and
the |
pure SIMT |
|
|
pure SIMT |
|
|
pure SIMT |
|
the cuBLASDx-backed |
requires MathDx |
|
single-warp |
pure SIMT |
|
block-tridiagonal PCG solve |
pure SIMT |
|
|
pure SIMT |
|
GEMM on sub-blocks with explicit leading dims ( |
pure SIMT |
|
symmetric-indefinite |
pure SIMT |
|
the fused LQR gain |
pure SIMT |
|
augmented |
pure SIMT |
|
the |
pure SIMT |
|
Featherstone |
pure SIMT |
|
batched |
pure SIMT |
|
MPPI weight update = |
pure SIMT |
|
friction-cone AL: |
pure SIMT |
|
narrow-phase |
pure SIMT |
|
batched Wahba/Kabsch via |
pure SIMT |
Build and run#
cd examples
make -j ARCH=sm_120 # auto-detects ARCH if omitted
make run # runs every built example
# 05_nvidia_gemm builds only when MATHDX_ROOT is set
Sources#
01_axpy_simt#
// 01_axpy_simt.cu — basic L1 op (AXPY, runtime size), pure SIMT (no MathDx).
//
// Build (from this examples/ dir):
// nvcc -std=c++17 -arch=sm_75 -I.. 01_axpy_simt.cu -o axpy && ./axpy
//
// Computes y = alpha*x + y for a length-n vector inside a single block.
#include "glass.cuh"
#include <cstdio>
#include <cuda_runtime.h>
__global__ void axpy_kernel(float *x, float *y, int n) {
// Runtime size: every thread in the block strides over the n elements.
glass::block::axpy(static_cast<uint32_t>(n), 1.5f, x, y); // y = 1.5*x + y
}
int main() {
const int n = 8;
float hx[n], hy[n];
for (int i = 0; i < n; ++i) { hx[i] = static_cast<float>(i); hy[i] = 1.0f; }
float *dx, *dy;
cudaMalloc(&dx, n * sizeof(float));
cudaMalloc(&dy, n * sizeof(float));
cudaMemcpy(dx, hx, n * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(dy, hy, n * sizeof(float), cudaMemcpyHostToDevice);
// One block per independent data item (here: a single vector).
axpy_kernel<<<1, 256>>>(dx, dy, n);
cudaDeviceSynchronize();
cudaMemcpy(hy, dy, n * sizeof(float), cudaMemcpyDeviceToHost);
printf("y = 1.5*x + 1 ->");
for (int i = 0; i < n; ++i) printf(" %.1f", hy[i]); // expect 1.0 2.5 4.0 ...
printf("\n");
cudaFree(dx); cudaFree(dy);
return 0;
}
02_gemm_conventions#
// 02_gemm_conventions.cu — THE single-block GEMM example: the standard-BLAS
// convention, both size overloads, all four transpose combos, why row-major
// needs no flag, and the cooperative-groups spelling of the same call.
// (Merged from the former 02_gemm / 04_cgrps / 10_gemm_basics /
// 11_rowmajor_is_transpose examples, 2026-08-11.)
//
// Build (from this examples/ dir):
// nvcc -std=c++17 -arch=sm_75 -I.. 02_gemm_conventions.cu -o gemm && ./gemm
//
// GLASS gemm follows the standard BLAS / cuBLAS / NumPy / Eigen convention:
//
// C = alpha * op(A) * op(B) + beta * C (column-major)
// C is M×N, contraction K.
// op(A) is M×K: TRANSPOSE_A=false ⇒ A is M×K ; true ⇒ A is K×M (op(A)=Aᵀ).
// op(B) is K×N: TRANSPOSE_B=false ⇒ B is K×N ; true ⇒ B is N×K (op(B)=Bᵀ).
//
// NumPy: C = alpha * opA(A) @ opB(B) + beta * C
// Eigen: C.noalias() = alpha * (opA(A) * opB(B)) + beta * C; // col-major
//
// We deliberately use a NON-SQUARE shape (M=2, N=3, K=4): the dimension order
// matters, and a square example would hide a wrong mapping.
#include "glass-cgrps.cuh" // pulls glass.cuh; also enables the cgrps section
#include <cooperative_groups.h>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <cuda_runtime.h>
static constexpr int M = 2, N = 3, K = 4;
// ─── §1 the four transpose combos (compile-time-size overload) ──────────────
template <bool TA, bool TB>
__global__ void run(const float* A, const float* B, float* C) {
// beta = 0 overload: C is overwritten (never read).
glass::block::gemm<float, M, N, K, TA, TB>(1.0f, const_cast<float*>(A), const_cast<float*>(B), C);
}
// ─── §2 the runtime-size overload — sizes as arguments, same convention ─────
__global__ void run_rt(const float* A, const float* B, float* C, int m, int n, int k) {
glass::block::gemm(static_cast<uint32_t>(m), static_cast<uint32_t>(n),
static_cast<uint32_t>(k),
1.0f, const_cast<float*>(A), const_cast<float*>(B), 0.0f, C);
}
// ─── §3 row-major is just a transpose (why there is no ROW_MAJOR_A flag) ────
// A row-major M×K matrix occupies the SAME bytes as a column-major K×M matrix,
// so a row-major A is read with TRANSPOSE_A=true — bit-identically.
__global__ void run_ta(const float* A, const float* B, float* C) {
glass::block::gemm<float, M, N, K, /*TA=*/true, /*TB=*/false>(1.f, const_cast<float*>(A), const_cast<float*>(B), C);
}
// ─── §4 the cooperative-groups spelling — same numerics, group-driven ───────
__global__ void run_cgrps(const float* A, const float* B, float* C) {
// Whole block (default group). A sub-block tile also works:
// auto warp = cooperative_groups::tiled_partition<32>(
// cooperative_groups::this_thread_block());
// glass::cgrps::gemm<float, M, N, K>(1.f, ..., warp);
glass::cgrps::gemm<float, M, N, K>(1.f, const_cast<float*>(A), const_cast<float*>(B), 0.f, C);
}
// Host reference: logical op(A) is M×K, op(B) is K×N, C is M×N (all col-major).
static void ref(const float* opA, const float* opB, float* C) {
for (int m = 0; m < M; m++)
for (int n = 0; n < N; n++) {
float s = 0;
for (int k = 0; k < K; k++) s += opA[m + k*M] * opB[k + n*K];
C[m + n*M] = s;
}
}
int main() {
// opA (M×K) and opB (K×N) are the LOGICAL operands.
float opA[M*K], opB[K*N];
for (int i = 0; i < M*K; i++) opA[i] = 0.1f * (i + 1);
for (int i = 0; i < K*N; i++) opB[i] = 0.2f * (i + 1) - 0.5f;
// Physical storage per transpose flag (a transposed operand is op(_)ᵀ col-major).
float A_n[M*K], A_t[K*M], B_n[K*N], B_t[N*K];
for (int m = 0; m < M; m++) for (int k = 0; k < K; k++) { A_n[m + k*M] = opA[m + k*M]; A_t[k + m*K] = opA[m + k*M]; }
for (int k = 0; k < K; k++) for (int n = 0; n < N; n++) { B_n[k + n*K] = opB[k + n*K]; B_t[n + k*N] = opB[k + n*K]; }
float ref_C[M*N]; ref(opA, opB, ref_C);
float *dA, *dB, *dC; cudaMalloc(&dA, sizeof(float)*K*M); cudaMalloc(&dB, sizeof(float)*N*K); cudaMalloc(&dC, sizeof(float)*M*N);
const char* names[4] = {"C = A * B ", "C = AT * B ", "C = A * BT ", "C = AT * BT "};
int bad = 0;
float C[M*N];
// §1 — all four transpose combos vs the host reference.
for (int combo = 0; combo < 4; combo++) {
bool ta = combo & 2, tb = combo & 1;
cudaMemcpy(dA, ta ? A_t : A_n, sizeof(float)*(ta ? K*M : M*K), cudaMemcpyHostToDevice);
cudaMemcpy(dB, tb ? B_t : B_n, sizeof(float)*(tb ? N*K : K*N), cudaMemcpyHostToDevice);
if (!ta && !tb) run<false,false><<<1,64>>>(dA, dB, dC);
else if ( ta && !tb) run<true ,false><<<1,64>>>(dA, dB, dC);
else if (!ta && tb) run<false,true ><<<1,64>>>(dA, dB, dC);
else run<true ,true ><<<1,64>>>(dA, dB, dC);
cudaDeviceSynchronize();
cudaMemcpy(C, dC, sizeof(C), cudaMemcpyDeviceToHost);
float md = 0; for (int i = 0; i < M*N; i++) md = fmaxf(md, fabsf(C[i] - ref_C[i]));
printf(" %s max_err=%.2e %s\n", names[combo], md, md < 1e-5 ? "ok" : "FAIL");
bad += (md >= 1e-5);
}
// §2 — runtime-size overload matches the compile-time one.
cudaMemcpy(dA, A_n, sizeof(A_n), cudaMemcpyHostToDevice);
cudaMemcpy(dB, B_n, sizeof(B_n), cudaMemcpyHostToDevice);
run_rt<<<1,64>>>(dA, dB, dC, M, N, K); cudaDeviceSynchronize();
cudaMemcpy(C, dC, sizeof(C), cudaMemcpyDeviceToHost);
{ float md = 0; for (int i = 0; i < M*N; i++) md = fmaxf(md, fabsf(C[i] - ref_C[i]));
printf(" runtime-size overload max_err=%.2e %s\n", md, md < 1e-5 ? "ok" : "FAIL");
bad += (md >= 1e-5); }
// §3 — row-major A via TRANSPOSE_A is BIT-identical to col-major NN.
float C_nn[M*N], C_ta[M*N];
float A_rowmajor[M*K]; // A[m*K + k] == bytes of the K×M col-major transpose
for (int m = 0; m < M; m++) for (int k = 0; k < K; k++) A_rowmajor[m*K + k] = opA[m + k*M];
cudaMemcpy(dA, A_n, sizeof(A_n), cudaMemcpyHostToDevice);
run<false,false><<<1,64>>>(dA, dB, dC); cudaDeviceSynchronize();
cudaMemcpy(C_nn, dC, sizeof(C_nn), cudaMemcpyDeviceToHost);
cudaMemcpy(dA, A_rowmajor, sizeof(A_rowmajor), cudaMemcpyHostToDevice);
run_ta<<<1,64>>>(dA, dB, dC); cudaDeviceSynchronize();
cudaMemcpy(C_ta, dC, sizeof(C_ta), cudaMemcpyDeviceToHost);
{ bool identical = (memcmp(C_nn, C_ta, sizeof(C_nn)) == 0);
printf(" row-major-via-TRANSPOSE_A vs col-major NN: %s\n",
identical ? "BIT-IDENTICAL" : "DIFFER");
bad += !identical; }
// §4 — the cgrps spelling produces the same answer.
cudaMemcpy(dA, A_n, sizeof(A_n), cudaMemcpyHostToDevice);
run_cgrps<<<1,64>>>(dA, dB, dC); cudaDeviceSynchronize();
cudaMemcpy(C, dC, sizeof(C), cudaMemcpyDeviceToHost);
{ float md = 0; for (int i = 0; i < M*N; i++) md = fmaxf(md, fabsf(C[i] - ref_C[i]));
printf(" glass::cgrps::gemm max_err=%.2e %s\n", md, md < 1e-5 ? "ok" : "FAIL");
bad += (md >= 1e-5); }
cudaFree(dA); cudaFree(dB); cudaFree(dC);
printf(bad ? "FAIL\n" : "PASS\n");
return bad;
}
03_reductions_norms#
// 03_reductions_norms.cu — block reductions and norms: the in-place halving
// `reduce`, the warp-shuffle `reduce_fast` (with its scratch-sizing helper),
// and the nrm2 family across block + warp tiers. (Merged from the former
// 03_reduce / 12_nrm2 examples, 2026-08-11.)
//
// Build (from this examples/ dir):
// nvcc -std=c++17 -arch=sm_75 -I.. 03_reductions_norms.cu -o reductions && ./reductions
//
// reduce(x): x[0] = Σ xᵢ (in place, destructive)
// nrm2(x) = sqrt(Σ xᵢ²) NumPy: np.linalg.norm(x); Eigen: x.norm()
//
// The `_fast` suffix names the warp-shuffle REDUCTION STRATEGY (one scratch
// slot per warp — size with `reduce_fast_scratch_bytes<T>(blockDim)`);
// `glass::warp::nrm2` is the warp-tier form and returns its value in-register.
#include "glass.cuh"
#include <cstdio>
#include <cmath>
#include <cuda_runtime.h>
static constexpr int N = 8;
__global__ void k_reduce(float *x, int n) {
glass::block::reduce(static_cast<uint32_t>(n), x); // x[0] = sum(x)
}
__global__ void k_reduce_fast(float *x, int n) {
extern __shared__ float scratch[]; // one float per warp
glass::block::reduce_fast(static_cast<uint32_t>(n), x, scratch);
}
__global__ void k_nrm2_block(float* x, float* scratch) {
// nrm2_fast<T, N> — compile-time length, warp-reduced, DESTRUCTIVE (x[0] gets the result).
glass::block::nrm2_fast<float, N>(x, scratch);
}
__global__ void k_nrm2_warp(uint32_t n, const float* x, float* out) {
float r = glass::warp::nrm2<float>(n, x); // value-returning, non-destructive
if ((threadIdx.x & 31) == 0) *out = r;
}
int main() {
float hx[N];
for (int i = 0; i < N; ++i) hx[i] = static_cast<float>(i + 1); // sum = 36
float xn[N]; for (int i = 0; i < N; i++) xn[i] = 0.5f*i - 1.3f;
double s = 0; for (int i = 0; i < N; i++) s += (double)xn[i]*xn[i];
const float nrm_expected = (float)sqrt(s);
float *dx, *dout, *dscr;
cudaMalloc(&dx, N * sizeof(float));
cudaMalloc(&dout, sizeof(float));
cudaMalloc(&dscr, 8 * sizeof(float));
int bad = 0;
float out;
cudaMemcpy(dx, hx, sizeof(hx), cudaMemcpyHostToDevice);
k_reduce<<<1, 256>>>(dx, N); cudaDeviceSynchronize();
cudaMemcpy(&out, dx, sizeof(float), cudaMemcpyDeviceToHost);
printf(" reduce sum(1..8) = %.0f (expect 36)\n", out);
bad += (out != 36.f);
cudaMemcpy(dx, hx, sizeof(hx), cudaMemcpyHostToDevice);
const int threads = 256;
size_t smem = glass::reduce_fast_scratch_bytes<float>(threads);
k_reduce_fast<<<1, threads, smem>>>(dx, N); cudaDeviceSynchronize();
cudaMemcpy(&out, dx, sizeof(float), cudaMemcpyDeviceToHost);
printf(" reduce_fast sum(1..8) = %.0f (expect 36)\n", out);
bad += (out != 36.f);
float block_r, warp_r;
cudaMemcpy(dx, xn, sizeof(xn), cudaMemcpyHostToDevice);
k_nrm2_block<<<1, 64>>>(dx, dscr); cudaDeviceSynchronize();
cudaMemcpy(&block_r, dx, sizeof(float), cudaMemcpyDeviceToHost); // result in x[0]
cudaMemcpy(dx, xn, sizeof(xn), cudaMemcpyHostToDevice); // restore (destructive)
k_nrm2_warp<<<1, 32>>>(N, dx, dout); cudaDeviceSynchronize();
cudaMemcpy(&warp_r, dout, sizeof(float), cudaMemcpyDeviceToHost);
printf(" nrm2 block=%.6f warp=%.6f expected=%.6f\n", block_r, warp_r, nrm_expected);
bad += (fabsf(block_r - nrm_expected) >= 1e-5) + (fabsf(warp_r - nrm_expected) >= 1e-5);
cudaFree(dx); cudaFree(dout); cudaFree(dscr);
printf(bad ? "FAIL\n" : "PASS\n");
return bad;
}
04_gemm_dispatch#
// 04_gemm_dispatch.cu — glass::gemm_dispatch + dynamic shared memory (tiled path).
//
// Build (from this examples/ dir):
// nvcc -std=c++17 -arch=sm_75 -I.. 04_gemm_dispatch.cu -o dispatch && ./dispatch
//
// glass::gemm_dispatch auto-selects the shared-memory-tiled GEMM when scratch
// pointers are supplied (and m*n <= blockDim), else the plain path. The host
// helper glass_gemm_dispatch_smem() computes the bytes to launch with; it
// returns 0 when tiling is not warranted (then pass nullptr scratch).
#include "glass.cuh"
#include <cstdio>
#include <cuda_runtime.h>
// have_smem != 0 => the launch reserved dynamic shared memory for the tiles.
__global__ void dispatch_kernel(float *A, float *B, float *C,
int m, int n, int k, int have_smem) {
extern __shared__ float scratch[];
// TILE defaults to 8, so the B-tile starts m*8 floats past s_A.
float *s_A = have_smem ? scratch : nullptr;
float *s_B = have_smem ? scratch + m * 8 : nullptr;
glass::block::gemm_dispatch(static_cast<uint32_t>(m), static_cast<uint32_t>(n),
static_cast<uint32_t>(k), 1.f, A, B, 0.f, C, s_A, s_B);
}
int main() {
const int m = 4, n = 4, k = 4;
const int threads = 256;
float hA[m * n], hB[n * k], hC[m * k] = {0};
// Column-major: A = identity, B = ramp -> C = B.
for (int i = 0; i < m * n; ++i) hA[i] = 0.f;
for (int i = 0; i < m; ++i) hA[i + i * m] = 1.f; // A = I
for (int i = 0; i < n * k; ++i) hB[i] = static_cast<float>(i);
float *dA, *dB, *dC;
cudaMalloc(&dA, sizeof(hA));
cudaMalloc(&dB, sizeof(hB));
cudaMalloc(&dC, sizeof(hC));
cudaMemcpy(dA, hA, sizeof(hA), cudaMemcpyHostToDevice);
cudaMemcpy(dB, hB, sizeof(hB), cudaMemcpyHostToDevice);
// Host: how many shared bytes does the dispatched path need?
size_t smem = glass_gemm_dispatch_smem<float>(m, n, threads);
printf("dispatch smem = %zu bytes (%s)\n", smem,
smem ? "tiled path" : "plain path");
dispatch_kernel<<<1, threads, smem>>>(dA, dB, dC, m, n, k, smem > 0 ? 1 : 0);
cudaDeviceSynchronize();
cudaMemcpy(hC, dC, sizeof(hC), cudaMemcpyDeviceToHost);
printf("C = I*B -> C[0]=%.0f C[5]=%.0f C[15]=%.0f (expect 0 5 15)\n",
hC[0], hC[5], hC[15]);
cudaFree(dA); cudaFree(dB); cudaFree(dC);
return 0;
}
05_nvidia_gemm#
// 05_nvidia_gemm.cu — cuBLASDx-backed GEMM via glass::nvidia:: (REQUIRES MathDx).
//
// This is the ONLY example that needs NVIDIA MathDx (cuBLASDx). The pure-SIMT
// examples 01-05 build with plain nvcc; this one does not.
//
// Build (from this examples/ dir), with MATHDX_ROOT pointing at your MathDx
// install (see ../bench/INSTALL.md):
//
// nvcc -std=c++17 -arch=sm_86 -I.. \
// -DGLASS_BENCH_CUBLASDX -DSMS=860 \
// --expt-relaxed-constexpr -Xptxas -O1 \
// -I$MATHDX_ROOT/include \
// -I$MATHDX_ROOT/external/cutlass/include \
// 05_nvidia_gemm.cu -o nvidia_gemm && ./nvidia_gemm
//
// Notes:
// * -DGLASS_BENCH_CUBLASDX force-includes <cublasdx.hpp> from glass-nvidia.cuh.
// * -DSMS=XXX must match your -arch (860 for sm_86, 1200 for sm_120, ...);
// it selects the cuBLASDx-tuned config and the pre-instantiated GEMM table.
// * 16x16x16 is a pre-instantiated cuBLASDx shape (see glass-nvidia.cuh); the
// default form launches with EXACTLY gemm_threads<>() threads and
// gemm_scratch_bytes<>() bytes of shared memory — a mismatch deadlocks.
#include "glass-nvidia.cuh"
#include <cstdio>
#include <cuda_runtime.h>
constexpr int M = 16, N = 16, K = 16;
// Host-queryable, constexpr: the thread count + shared bytes cuBLASDx wants.
constexpr auto SMEM = glass::nvidia::block::gemm_scratch_bytes<float, M, N, K>();
constexpr auto THREADS = glass::nvidia::block::gemm_threads<float, M, N, K>();
__global__ void nvidia_gemm(float *A, float *B, float *C) {
extern __shared__ __align__(16) char smem_buf[];
glass::nvidia::block::gemm<float, M, N, K>(1.f, A, B, 0.f, C, smem_buf);
}
int main() {
float hA[M * N], hB[N * K], hC[M * K] = {0};
// Column-major: A = identity, B = ramp -> C = B.
for (int i = 0; i < M * N; ++i) hA[i] = 0.f;
for (int i = 0; i < M; ++i) hA[i + i * M] = 1.f; // A = I
for (int i = 0; i < N * K; ++i) hB[i] = static_cast<float>(i);
float *dA, *dB, *dC;
cudaMalloc(&dA, sizeof(hA));
cudaMalloc(&dB, sizeof(hB));
cudaMalloc(&dC, sizeof(hC));
cudaMemcpy(dA, hA, sizeof(hA), cudaMemcpyHostToDevice);
cudaMemcpy(dB, hB, sizeof(hB), cudaMemcpyHostToDevice);
// Launch with the EXACT thread count + smem cuBLASDx picked.
nvidia_gemm<<<1, THREADS, SMEM>>>(dA, dB, dC);
cudaDeviceSynchronize();
cudaMemcpy(hC, dC, sizeof(hC), cudaMemcpyDeviceToHost);
printf("glass::nvidia::gemm C = I*B (threads=%u smem=%zu)\n",
(unsigned)THREADS, (size_t)SMEM);
printf("C[0]=%.0f C[17]=%.0f C[255]=%.0f (expect 0 17 255)\n",
hC[0], hC[17], hC[255]);
cudaFree(dA); cudaFree(dB); cudaFree(dC);
return 0;
}
06_warp_ops#
// 06_warp_ops.cu — single-warp (glass::warp::) primitives, launched <<<1,32>>>.
//
// Build (from this examples/ dir):
// nvcc -std=c++17 -arch=sm_75 -I.. 06_warp_ops.cu -o warp_ops && ./warp_ops
//
// The glass::warp:: namespace holds warp-scoped SIMT variants (raw __shfl, one
// 32-lane warp, no shared scratch, no __syncthreads) for warp-per-problem
// kernels — e.g. a block that processes many independent problems, one per warp.
// No cooperative groups, no vendor (cuBLASDx/cuSOLVERDx) deps.
#include "glass.cuh"
#include <cstdio>
#include <cuda_runtime.h>
// Sum a vector within one warp; result in x[0].
__global__ void warp_reduce_kernel(float *x, int n) {
glass::warp::reduce(static_cast<uint32_t>(n), x);
}
// 4x4 column-major GEMM C = A*B within one warp (the mat4-multiply use case).
__global__ void warp_gemm_kernel(float *A, float *B, float *C) {
glass::warp::gemm<float, 4, 4, 4>(1.0f, A, B, 0.0f, C);
}
// SPD solve A x = b within one warp: factor A = L Lᵀ, then forward + transpose solve.
__global__ void warp_posv_kernel(float *A, float *b) {
glass::warp::potrf<float, 3>(A);
glass::warp::trsv<float, 3>(A, b); // forward: L y = b
glass::warp::trsv<float, 3, glass::block::FillMode::Lower, glass::block::Diag::NonUnit, true>(A, b); // back: Lᵀ x = y
}
int main() {
// ── warp::reduce ──
const int n = 8;
float hx[n];
for (int i = 0; i < n; ++i) hx[i] = static_cast<float>(i + 1); // sum = 36
float *dx; cudaMalloc(&dx, n * sizeof(float));
cudaMemcpy(dx, hx, n * sizeof(float), cudaMemcpyHostToDevice);
warp_reduce_kernel<<<1, 32>>>(dx, n);
cudaDeviceSynchronize();
float s = 0.f; cudaMemcpy(&s, dx, sizeof(float), cudaMemcpyDeviceToHost);
printf("glass::warp::reduce sum(1..8) = %.0f (expect 36)\n", s);
// ── warp::gemm (4x4, column-major; B = identity ⇒ C = A) ──
float hA[16], hB[16] = {0}, hC[16] = {0};
for (int i = 0; i < 16; ++i) hA[i] = static_cast<float>(i);
for (int i = 0; i < 4; ++i) hB[i*4 + i] = 1.0f; // identity
float *dA, *dB, *dC;
cudaMalloc(&dA, 16*sizeof(float)); cudaMalloc(&dB, 16*sizeof(float)); cudaMalloc(&dC, 16*sizeof(float));
cudaMemcpy(dA, hA, 16*sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(dB, hB, 16*sizeof(float), cudaMemcpyHostToDevice);
warp_gemm_kernel<<<1, 32>>>(dA, dB, dC);
cudaDeviceSynchronize();
cudaMemcpy(hC, dC, 16*sizeof(float), cudaMemcpyDeviceToHost);
printf("glass::warp::gemm C[0,5,10,15] = %.0f %.0f %.0f %.0f (expect 0 5 10 15)\n",
hC[0], hC[5], hC[10], hC[15]);
// ── warp:: SPD solve (3x3, column-major) ──
// A = [[4,1,0],[1,3,1],[0,1,2]] (SPD), b = [1,2,3]; solve A x = b.
float hAs[9] = {4,1,0, 1,3,1, 0,1,2}; // column-major == symmetric here
float hb[3] = {1,2,3};
float *dAs, *db;
cudaMalloc(&dAs, 9*sizeof(float)); cudaMalloc(&db, 3*sizeof(float));
cudaMemcpy(dAs, hAs, 9*sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(db, hb, 3*sizeof(float), cudaMemcpyHostToDevice);
warp_posv_kernel<<<1, 32>>>(dAs, db);
cudaDeviceSynchronize();
float hx3[3]; cudaMemcpy(hx3, db, 3*sizeof(float), cudaMemcpyDeviceToHost);
printf("glass::warp:: SPD solve x = %.3f %.3f %.3f\n", hx3[0], hx3[1], hx3[2]);
cudaFree(dx); cudaFree(dA); cudaFree(dB); cudaFree(dC); cudaFree(dAs); cudaFree(db);
return 0;
}
07_pcg_solve#
// 07_pcg_solve.cu — block-tridiagonal PCG solve (glass::pcg), pure SIMT.
//
// Build (from this examples/ dir):
// nvcc -std=c++17 -arch=sm_75 -I.. 07_pcg_solve.cu -o pcg && ./pcg
//
// Solves an SPD block-tridiagonal system S x = b in ONE CUDA block with
// preconditioned conjugate gradient, using a block-Jacobi preconditioner Pinv
// (inverse of each diagonal block). The matrix is stored as KP row-major
// [L | D | R] strips and the vectors use the padded (KP+2)*SS layout. See
// docs/source/user_guide/concepts/block_tridiagonal.rst.
#include "glass.cuh"
#include <cstdio>
#include <cuda_runtime.h>
constexpr int SS = 2; // state_size (block dimension)
constexpr int KP = 3; // knot_points (number of block-rows)
constexpr int BRL = 3 * SS; // columns per strip: [L|D|R]
constexpr int N = SS * KP; // unpadded system size
constexpr int VEC = (KP + 2) * SS; // padded vector length
__global__ void pcg_kernel(float *x, float *S, float *Pinv, float *b,
unsigned max_iters, float rel_tol, float abs_tol,
unsigned *iters) {
extern __shared__ float s_mem[];
glass::block::pcg<float, SS, KP>(x, S, Pinv, b, s_mem,
max_iters, rel_tol, abs_tol, iters);
}
// Write a 2x2 block into strip k at column offset `coloff` (0=L, SS=D, 2SS=R).
static void put2x2(float *M, int k, int coloff,
float a, float b, float c, float d) {
float *strip = M + k * BRL * SS; // SS rows x BRL cols, row-major
strip[0 * BRL + coloff + 0] = a; strip[0 * BRL + coloff + 1] = b;
strip[1 * BRL + coloff + 0] = c; strip[1 * BRL + coloff + 1] = d;
}
int main() {
// Diagonally-dominant SPD blocks: D = [[4,1],[1,4]], off-diagonals 0.1*I.
float S[KP * BRL * SS] = {0};
float Pinv[KP * BRL * SS] = {0};
for (int k = 0; k < KP; ++k) {
put2x2(S, k, SS, 4, 1, 1, 4); // D
if (k > 0) put2x2(S, k, 0, 0.1f, 0, 0, 0.1f); // L
if (k < KP - 1) put2x2(S, k, 2 * SS, 0.1f, 0, 0, 0.1f); // R
// Block-Jacobi preconditioner: inv(D) = (1/15)[[4,-1],[-1,4]].
put2x2(Pinv, k, SS, 4.f / 15, -1.f / 15, -1.f / 15, 4.f / 15);
}
// Padded RHS (b = 1 on the real entries, zeros in the pads) and zero guess.
float b[VEC] = {0}, x[VEC] = {0};
for (int i = 0; i < N; ++i) b[SS + i] = 1.0f;
float *dS, *dPinv, *db, *dx; unsigned *diters;
cudaMalloc(&dS, sizeof(S)); cudaMalloc(&dPinv, sizeof(Pinv));
cudaMalloc(&db, sizeof(b)); cudaMalloc(&dx, sizeof(x));
cudaMalloc(&diters, sizeof(unsigned));
cudaMemcpy(dS, S, sizeof(S), cudaMemcpyHostToDevice);
cudaMemcpy(dPinv, Pinv, sizeof(Pinv), cudaMemcpyHostToDevice);
cudaMemcpy(db, b, sizeof(b), cudaMemcpyHostToDevice);
cudaMemcpy(dx, x, sizeof(x), cudaMemcpyHostToDevice);
const int threads = 32; // must be a multiple of 32 (warp-dot)
size_t smem = glass::block::pcg_scratch_bytes<float, SS, KP>(threads);
pcg_kernel<<<1, threads, smem>>>(dx, dS, dPinv, db, 100, 1e-6f, 1e-12f, diters);
cudaDeviceSynchronize();
unsigned iters;
cudaMemcpy(x, dx, sizeof(x), cudaMemcpyDeviceToHost);
cudaMemcpy(&iters, diters, sizeof(unsigned), cudaMemcpyDeviceToHost);
printf("PCG converged in %u iters; x =", iters);
for (int i = 0; i < N; ++i) printf(" %.4f", x[SS + i]); // strip the pads
printf("\n");
cudaFree(dS); cudaFree(dPinv); cudaFree(db); cudaFree(dx); cudaFree(diters);
return 0;
}
08_backend_picker#
// 08_backend_picker.cu — choosing a backend + launch config with glass-defaults.cuh.
//
// Build (from this examples/ dir, pure SIMT — no MathDx needed):
// nvcc -std=c++17 -arch=sm_75 -I.. 08_backend_picker.cu -o picker && ./picker
// (to make the `nvidia` tier eligible, include glass-nvidia.cuh first + link MathDx.)
//
// glass-defaults.cuh exposes the measured thread/warp/block/nvidia ladder (bench/RESULTS.md)
// as constexpr helpers. The pick is host-/codegen-side because warp, block, and nvidia need
// DIFFERENT <<<grid,block>>> launches — so you query at compile time and branch the launch.
// With no MathDx linked (as here), the `nvidia` tier collapses to its warp/block runner-up.
#include "glass.cuh"
#include "glass-defaults.cuh"
#include <cstdio>
#include <cuda_runtime.h>
using glass::op;
using glass::backend;
static const char* name(backend b) {
return b == backend::warp ? "warp"
: b == backend::block ? "block"
: b == backend::thread ? "thread" : "nvidia";
}
// ── one SPD solve A x = b, dispatched to the picked backend ──────────────────
template <int N> __global__ void k_block_posv(float* A, float* b) { glass::block::posv<float, N>(A, b); }
template <int N> __global__ void k_warp_posv (float* A, float* b) {
int w = blockIdx.x * blockDim.y + threadIdx.y; // one warp per problem
glass::warp::posv<float, N>(A + (size_t)w*N*N, b + w*N);
}
template <int N> __global__ void k_thread_posv(float* A, float* b) {
int p = blockIdx.x * blockDim.x + threadIdx.x; // one problem per THREAD
float a[N*N], x[N]; // register-resident at N<=7
for (int i = 0; i < N*N; i++) a[i] = A[(size_t)p*N*N + i];
for (int i = 0; i < N; i++) x[i] = b[(size_t)p*N + i];
glass::thread::posv<float, N>(a, x);
for (int i = 0; i < N; i++) b[(size_t)p*N + i] = x[i];
}
template <int N>
static void solve_dispatch(float* dA, float* db) {
// Compile-time pick from the measured table (T=float, build's SM).
constexpr backend be = glass::suggested_backend<op::posv, N, float>();
printf(" posv N=%d -> backend=%s", N, name(be));
if constexpr (be == backend::thread) {
constexpr int TPB = glass::suggested_threads_per_block<op::posv, N, float>();
printf(" (TPB=%d)\n", TPB);
k_thread_posv<N><<<1, 1>>>(dA, db); // 1 problem here -> 1 thread
} else if constexpr (be == backend::warp) {
constexpr int WPB = glass::suggested_warps_per_block<op::posv>();
printf(" (WPB=%d)\n", WPB);
k_warp_posv<N><<<1, dim3(32, 1)>>>(dA, db); // 1 problem here -> 1 warp
} else { // block (or nvidia collapsed to block); a real nvidia tier would launch cuSOLVERDx
constexpr int TB = glass::suggested_block_threads<op::posv, N, float>();
printf(" (TB=%d)\n", TB);
k_block_posv<N><<<1, TB>>>(dA, db);
}
cudaDeviceSynchronize();
}
int main() {
// 1) Show what the picker chooses across ops/sizes (all compile-time constants).
printf("backend picks (T=float, this build's SM; no MathDx -> nvidia collapses):\n");
printf(" dot N=8 : %s\n", name(glass::suggested_backend<op::dot, 8, float>()));
printf(" dot N=64 : %s\n", name(glass::suggested_backend<op::dot, 64, float>()));
printf(" posv N=8 : %s\n", name(glass::suggested_backend<op::posv, 8, float>()));
printf(" gemv N=16 : %s\n", name(glass::suggested_backend<op::gemv, 16, float>()));
printf(" gemv N=64 : %s\n", name(glass::suggested_backend<op::gemv, 64, float>()));
printf(" gemm N=8 : %s\n", name(glass::suggested_backend<op::gemm, 8, float>()));
printf(" gemm N=32 : %s\n", name(glass::suggested_backend<op::gemm, 32, float>()));
printf(" chol N=8 : %s\n", name(glass::suggested_backend<op::chol, 8, float>()));
printf(" chol N=64 : %s\n", name(glass::suggested_backend<op::chol, 64, float>()));
// 2) Use the pick to dispatch a real solve. SPD A = M·Mᵀ + N·I (column-major), N=16.
const int N = 16;
float hA[N*N], hb[N];
for (int i=0;i<N;i++) for (int j=0;j<N;j++) {
float m=0; for (int k=0;k<N;k++) m += (((i+2*k)%5)*0.1f) * (((j+2*k)%5)*0.1f);
hA[i+j*N] = m + (i==j ? (float)N : 0.0f);
}
for (int i=0;i<N;i++) hb[i] = 1.0f + 0.1f*i;
float *dA,*db; cudaMalloc(&dA,N*N*4); cudaMalloc(&db,N*4);
cudaMemcpy(dA,hA,N*N*4,cudaMemcpyHostToDevice); cudaMemcpy(db,hb,N*4,cudaMemcpyHostToDevice);
printf("\ndispatch a real solve:\n");
solve_dispatch<N>(dA, db);
float hx[N]; cudaMemcpy(hx,db,N*4,cudaMemcpyDeviceToHost);
float res=0; for (int i=0;i<N;i++){ float Ax=0; for (int j=0;j<N;j++) Ax+=hA[i+j*N]*hx[j]; float r=Ax-hb[i]; res = r<0?(res<-r?res:-r):(res<r?r:res); }
printf(" residual ||A x - b||_inf = %.2e -> %s\n", res, res < 1e-3f ? "OK" : "FAIL");
cudaFree(dA); cudaFree(db);
return 0;
}
09_gemm_strided#
// 09_gemm_strided.cu — GEMM on column-major sub-blocks with explicit
// leading dimensions (the "strided" GEMM).
//
// Build (from this examples/ dir):
// nvcc -std=c++17 -arch=sm_75 -I.. 09_gemm_strided.cu -o rsgemm && ./rsgemm
//
// Same standard convention as glass::gemm (C is M×N, contraction K), but A and B
// live inside larger buffers with custom column strides (leading dimensions):
//
// A is M×K, leading dim A_RS ≥ M : A[m][k] = A_buf[m + k*A_RS]
// B is K×N, leading dim B_RS ≥ K : B[k][n] = B_buf[k + n*B_RS]
// C is M×N, standard column-major (LDC = M).
//
// alpha/beta are at the FRONT, matching every other GLASS op:
// glass::gemm_strided<T, M, N, K, A_RS, B_RS>(alpha, A, B, beta, C);
// NumPy: C = alpha * A @ B + beta * C (on the strided sub-views)
#include "glass.cuh"
#include <cstdio>
#include <cmath>
#include <cuda_runtime.h>
static constexpr int M = 5, N = 7, K = 3, A_RS = 8, B_RS = 6; // A_RS>M, B_RS>K
__global__ void run(const float* A, const float* B, float beta, float* C) {
glass::block::gemm_strided<float, M, N, K, A_RS, B_RS>(1.5f, const_cast<float*>(A), const_cast<float*>(B), beta, C);
}
int main() {
float A[A_RS*K], B[B_RS*N], C[M*N];
for (int i = 0; i < A_RS*K; i++) A[i] = 0.1f*i - 0.3f;
for (int i = 0; i < B_RS*N; i++) B[i] = 0.2f*i - 0.5f;
for (int i = 0; i < M*N; i++) C[i] = 1.0f;
const float alpha = 1.5f, beta = 0.25f;
float ref[M*N];
for (int m = 0; m < M; m++) for (int n = 0; n < N; n++) {
float s = 0; for (int k = 0; k < K; k++) s += A[m + k*A_RS] * B[k + n*B_RS];
ref[m + n*M] = alpha*s + beta*C[m + n*M];
}
float *dA, *dB, *dC; cudaMalloc(&dA, sizeof(A)); cudaMalloc(&dB, sizeof(B)); cudaMalloc(&dC, sizeof(C));
cudaMemcpy(dA, A, sizeof(A), cudaMemcpyHostToDevice);
cudaMemcpy(dB, B, sizeof(B), cudaMemcpyHostToDevice);
cudaMemcpy(dC, C, sizeof(C), cudaMemcpyHostToDevice);
run<<<1, 64>>>(dA, dB, beta, dC); cudaDeviceSynchronize();
float out[M*N]; cudaMemcpy(out, dC, sizeof(out), cudaMemcpyDeviceToHost);
float md = 0; for (int i = 0; i < M*N; i++) md = fmaxf(md, fabsf(out[i] - ref[i]));
printf(" gemm_strided %dx%dx%d A_RS=%d B_RS=%d max_err=%.2e\n", M, N, K, A_RS, B_RS, md);
cudaFree(dA); cudaFree(dB); cudaFree(dC);
printf(md < 1e-4 ? "PASS\n" : "FAIL\n");
return md < 1e-4 ? 0 : 1;
}
10_ldlt_solve#
// 10_ldlt_solve.cu — symmetric-INDEFINITE solve via LDLᵀ (no square root).
//
// Build (from this examples/ dir):
// nvcc -std=c++17 -arch=sm_75 -I.. 10_ldlt_solve.cu -o ldlt && ./ldlt
//
// Cholesky (potrf/posv) requires SPD. A symmetric matrix with NEGATIVE
// eigenvalues (a KKT / saddle-point system) has no Cholesky factor — but it
// does have A = L·D·Lᵀ with unit-lower L and a signed diagonal D. glass::ldlt
// factors it in place (SciPy: lu, d, _ = scipy.linalg.ldl(A, lower=True)) and
// glass::ldlt_solve runs the three sweeps L y = b, z = y/D, Lᵀ x = z.
//
// Also shown: the compile-out CHECK=true path. Factorizations default to
// CHECK=false and will silently produce NaN/Inf on a non-factorable input;
// with CHECK=true rank 0 reports a zero/NaN pivot via s_fail and the pivot
// sign counts {n_pos, n_neg, n_zero} (the matrix INERTIA) via s_inertia.
// Scratch is sized by glass::ldlt_scratch_bytes<T>(n) (covers the pivot path;
// the non-pivoted path would also accept nullptr).
#include "glass.cuh"
#include <cstdio>
#include <cmath>
#include <cuda_runtime.h>
static constexpr int N = 4; // indefinite 4x4 system
static constexpr int NZ = 2; // zero-pivot demo matrix
// Factor (checked, non-pivoted) then solve A x = b in place; x lands in b.
__global__ void k_factor_solve(float* A, float* b, int* fail, int* inertia) {
__shared__ float s_scratch[glass::block::ldlt_scratch_bytes<float>(N) / sizeof(float)];
glass::block::ldlt<float, N, /*CHECK=*/true>(A, s_scratch, /*pivot=*/false, nullptr, fail, inertia);
glass::block::ldlt_solve<float, N>(A, b);
}
// CHECK=true on a matrix whose FIRST pivot is exactly zero: D_0 = A_00 = 0.
// The factor itself is garbage (division by zero) — the point is fail == 1.
__global__ void k_factor_zero_pivot(float* A, int* fail) {
__shared__ float s_scratch[glass::block::ldlt_scratch_bytes<float>(NZ) / sizeof(float)];
glass::block::ldlt<float, NZ, /*CHECK=*/true>(A, s_scratch, /*pivot=*/false, nullptr, fail, nullptr);
}
int main() {
// Symmetric INDEFINITE A (column-major == row-major here, it's symmetric).
// Pivots come out {+, -, +, -} => inertia {2, 2, 0}: not SPD, potrf would NaN.
float hA[N*N] = { 2, 1, 0, 0,
1,-3, 1, 0,
0, 1, 4, 1,
0, 0, 1,-2 };
float x_true[N] = { 1.f, -2.f, 3.f, 0.5f };
float hb[N];
for (int i = 0; i < N; i++) { // b = A * x_true (col-major)
float s = 0; for (int j = 0; j < N; j++) s += hA[i + j*N] * x_true[j];
hb[i] = s;
}
float *dA, *db; int *dfail, *dinertia;
cudaMalloc(&dA, sizeof(hA)); cudaMalloc(&db, sizeof(hb));
cudaMalloc(&dfail, sizeof(int)); cudaMalloc(&dinertia, 3*sizeof(int));
cudaMemcpy(dA, hA, sizeof(hA), cudaMemcpyHostToDevice);
cudaMemcpy(db, hb, sizeof(hb), cudaMemcpyHostToDevice);
k_factor_solve<<<1, 64>>>(dA, db, dfail, dinertia);
cudaDeviceSynchronize();
float hx[N]; int fail, inertia[3];
cudaMemcpy(hx, db, sizeof(hx), cudaMemcpyDeviceToHost);
cudaMemcpy(&fail, dfail, sizeof(int), cudaMemcpyDeviceToHost);
cudaMemcpy(inertia, dinertia, sizeof(inertia), cudaMemcpyDeviceToHost);
float md = 0; for (int i = 0; i < N; i++) md = fmaxf(md, fabsf(hx[i] - x_true[i]));
printf(" ldlt + ldlt_solve x = %.3f %.3f %.3f %.3f max_err=%.2e %s\n",
hx[0], hx[1], hx[2], hx[3], md, md < 1e-4f ? "ok" : "FAIL");
printf(" CHECK (good A) fail=%d inertia={%d,%d,%d} (expect 0, {2,2,0})\n",
fail, inertia[0], inertia[1], inertia[2]);
bool ok = (md < 1e-4f) && (fail == 0)
&& inertia[0] == 2 && inertia[1] == 2 && inertia[2] == 0;
// Zero leading pivot: D_0 = 0 breaks the non-pivoted recurrence. Without
// CHECK this silently NaNs; with CHECK, s_fail flags it (a caller could
// then retry with ldlt(pivot=true) or escalate regularization).
float hZ[NZ*NZ] = { 0, 1,
1, 2 };
float *dZ; cudaMalloc(&dZ, sizeof(hZ));
cudaMemcpy(dZ, hZ, sizeof(hZ), cudaMemcpyHostToDevice);
k_factor_zero_pivot<<<1, 64>>>(dZ, dfail);
cudaDeviceSynchronize();
cudaMemcpy(&fail, dfail, sizeof(int), cudaMemcpyDeviceToHost);
printf(" CHECK (zero pivot) fail=%d (expect 1)\n", fail);
ok = ok && (fail == 1);
cudaFree(dA); cudaFree(db); cudaFree(dZ); cudaFree(dfail); cudaFree(dinertia);
printf(ok ? "PASS\n" : "FAIL\n");
return ok ? 0 : 1;
}
11_riccati_gain#
// 11_riccati_gain.cu — LQR feedback gain K = (R + BᵀPB)⁻¹ (BᵀPA) in one call.
//
// Build (from this examples/ dir):
// nvcc -std=c++17 -arch=sm_75 -I.. 11_riccati_gain.cu -o riccati && ./riccati
//
// The control-update solve at the heart of an LQR / iLQR backward pass.
// glass::riccati_gain composes three library primitives in one block:
// S = R + BᵀPB (congruence_sym, NU×NU control Hessian)
// G = BᵀPA (bilinear, NU×NX coupling)
// S·K = G (checked multi-RHS Cholesky posv; K overwrites G)
// Inputs P, A, B, R are unchanged; Kgain holds K (NU×NX, column-major).
// Scratch is dynamic shared memory sized by the host-callable
// glass::riccati_scratch_bytes<T, NX, NU>() (BYTES, pass as the launch smem).
// s_fail reports a non-PD S (an iLQR caller would escalate rho and retry via
// the REGULARIZE template flag; see 10_ldlt_solve.cu for the CHECK idea).
#include "glass.cuh"
#include <cstdio>
#include <cmath>
#include <cuda_runtime.h>
static constexpr int NX = 4, NU = 2;
__global__ void k_riccati(const float* P, const float* A, const float* B,
const float* R, float* Kgain, int* fail) {
extern __shared__ float s_scratch[];
glass::block::riccati_gain<float, NX, NU>(P, A, B, R, Kgain, s_scratch, 0.f, fail);
}
// CPU reference with plain loops (everything column-major, X[row + col*rows]).
static void ref_gain(const float* P, const float* A, const float* B,
const float* R, float* K) {
float PB[NX*NU], PA[NX*NX], S[NU*NU], G[NU*NX];
for (int c = 0; c < NU; c++) for (int r = 0; r < NX; r++) { // PB = P·B
float s = 0; for (int k = 0; k < NX; k++) s += P[r + k*NX] * B[k + c*NX];
PB[r + c*NX] = s;
}
for (int c = 0; c < NX; c++) for (int r = 0; r < NX; r++) { // PA = P·A
float s = 0; for (int k = 0; k < NX; k++) s += P[r + k*NX] * A[k + c*NX];
PA[r + c*NX] = s;
}
for (int c = 0; c < NU; c++) for (int r = 0; r < NU; r++) { // S = R + Bᵀ·PB
float s = 0; for (int k = 0; k < NX; k++) s += B[k + r*NX] * PB[k + c*NX];
S[r + c*NU] = R[r + c*NU] + s;
}
for (int c = 0; c < NX; c++) for (int r = 0; r < NU; r++) { // G = Bᵀ·PA
float s = 0; for (int k = 0; k < NX; k++) s += B[k + r*NX] * PA[k + c*NX];
G[r + c*NU] = s;
}
// Solve S·K = G (S is SPD NU×NU, NX right-hand sides): Gaussian elimination.
for (int p = 0; p < NU; p++) {
for (int r = p + 1; r < NU; r++) {
float m = S[r + p*NU] / S[p + p*NU];
for (int c = p; c < NU; c++) S[r + c*NU] -= m * S[p + c*NU];
for (int c = 0; c < NX; c++) G[r + c*NU] -= m * G[p + c*NU];
}
}
for (int c = 0; c < NX; c++)
for (int r = NU - 1; r >= 0; r--) {
float s = G[r + c*NU];
for (int k = r + 1; k < NU; k++) s -= S[r + k*NU] * K[k + c*NU];
K[r + c*NU] = s / S[r + r*NU];
}
}
int main() {
// P: symmetric PD cost-to-go; A: state Jacobian; B: control Jacobian; R: SPD.
float P[NX*NX], A[NX*NX], B[NX*NU];
for (int i = 0; i < NX; i++) for (int j = 0; j < NX; j++) {
P[i + j*NX] = 0.1f*(i + j) + (i == j ? 2.0f + i : 0.0f); // symmetric PD
A[i + j*NX] = 0.1f*i - 0.05f*j + (i == j ? 1.0f : 0.0f);
}
for (int i = 0; i < NX*NU; i++) B[i] = 0.2f*i - 0.3f;
float R[NU*NU] = { 1.0f, 0.2f,
0.2f, 0.8f };
float ref_K[NU*NX]; ref_gain(P, A, B, R, ref_K);
float *dP, *dA, *dB, *dR, *dK; int *dfail;
cudaMalloc(&dP, sizeof(P)); cudaMalloc(&dA, sizeof(A)); cudaMalloc(&dB, sizeof(B));
cudaMalloc(&dR, sizeof(R)); cudaMalloc(&dK, sizeof(float)*NU*NX); cudaMalloc(&dfail, sizeof(int));
cudaMemcpy(dP, P, sizeof(P), cudaMemcpyHostToDevice);
cudaMemcpy(dA, A, sizeof(A), cudaMemcpyHostToDevice);
cudaMemcpy(dB, B, sizeof(B), cudaMemcpyHostToDevice);
cudaMemcpy(dR, R, sizeof(R), cudaMemcpyHostToDevice);
const size_t smem = glass::block::riccati_scratch_bytes<float, NX, NU>(); // BYTES
k_riccati<<<1, 128, smem>>>(dP, dA, dB, dR, dK, dfail);
cudaDeviceSynchronize();
float K[NU*NX]; int fail;
cudaMemcpy(K, dK, sizeof(K), cudaMemcpyDeviceToHost);
cudaMemcpy(&fail, dfail, sizeof(int), cudaMemcpyDeviceToHost);
float md = 0; for (int i = 0; i < NU*NX; i++) md = fmaxf(md, fabsf(K[i] - ref_K[i]));
printf(" riccati_gain (NX=%d, NU=%d, smem=%zu B) fail=%d max_err vs CPU = %.2e\n",
NX, NU, smem, fail, md);
for (int r = 0; r < NU; r++)
printf(" K[%d,:] = %8.4f %8.4f %8.4f %8.4f\n",
r, K[r + 0*NU], K[r + 1*NU], K[r + 2*NU], K[r + 3*NU]);
bool ok = (md < 1e-4f) && (fail == 0);
cudaFree(dP); cudaFree(dA); cudaFree(dB); cudaFree(dR); cudaFree(dK); cudaFree(dfail);
printf(ok ? "PASS\n" : "FAIL\n");
return ok ? 0 : 1;
}
12_inv#
// 12_inv.cu — matrix inversion: the augmented [A | I] convention + the robust
// partial-pivoting variant.
//
// Build (from this examples/ dir):
// nvcc -std=c++17 -arch=sm_75 -I.. 12_inv.cu -o inv && ./inv
//
// glass::inv is Gauss-Jordan on an AUGMENTED buffer: you hand it a column-major
// N x 2N matrix laid out [A | I] (left half A, right half identity) and on
// return the RIGHT half (columns N..2N-1) holds A⁻¹. Scratch is sized by
// glass::inv_scratch_bytes<T>(N).
//
// The plain inv divides by each leading pivot AS-IS — a zero (or tiny) leading
// pivot produces Inf/NaN even when A is perfectly invertible. glass::inv_pivoted
// (scratch: inv_pivoted_scratch_bytes<T>(N)) row-pivots across the full
// augmented width, so the permutation is absorbed and the right half is still
// A⁻¹ directly. This example shows plain inv going non-finite on a zero leading
// pivot and inv_pivoted recovering the exact inverse.
#include "glass.cuh"
#include <cstdio>
#include <cmath>
#include <cuda_runtime.h>
static constexpr int N = 3; // well-behaved matrix
static constexpr int NP = 2; // zero-leading-pivot matrix
__global__ void k_inv(float* Aaug) {
__shared__ float s_scratch[glass::block::inv_scratch_bytes<float>(N) / sizeof(float)];
glass::block::inv<float, N>(Aaug, s_scratch);
}
__global__ void k_inv_plain_np(float* Aaug) { // mishandles pivot 0
__shared__ float s_scratch[glass::block::inv_scratch_bytes<float>(NP) / sizeof(float)];
glass::block::inv<float, NP>(Aaug, s_scratch);
}
__global__ void k_inv_pivoted_np(float* Aaug) { // robust
__shared__ float s_scratch[glass::block::inv_pivoted_scratch_bytes<float>(NP) / sizeof(float)];
glass::block::inv_pivoted<float, NP>(Aaug, s_scratch);
}
// Build the column-major dim x 2*dim augmented [A | I] buffer from A.
static void make_augmented(int dim, const float* A, float* aug) {
for (int c = 0; c < dim; c++)
for (int r = 0; r < dim; r++) {
aug[r + c*dim] = A[r + c*dim]; // left: A
aug[r + (dim + c)*dim] = (r == c) ? 1.0f : 0.0f; // right: I
}
}
// max |(A * Ainv - I)[i,j]| (all column-major dim x dim).
static float inverse_residual(int dim, const float* A, const float* Ainv) {
float md = 0;
for (int c = 0; c < dim; c++)
for (int r = 0; r < dim; r++) {
float s = 0; for (int k = 0; k < dim; k++) s += A[r + k*dim] * Ainv[k + c*dim];
md = fmaxf(md, fabsf(s - (r == c ? 1.0f : 0.0f)));
}
return md;
}
int main() {
// ── glass::inv on a well-behaved 3x3 (column-major) ──
float A[N*N] = { 4, 1, 2, // col 0
1, 3, 0, // col 1
2, 0, 5 }; // col 2
float aug[N*2*N]; make_augmented(N, A, aug);
float* dAug; cudaMalloc(&dAug, sizeof(aug));
cudaMemcpy(dAug, aug, sizeof(aug), cudaMemcpyHostToDevice);
k_inv<<<1, 64>>>(dAug);
cudaDeviceSynchronize();
cudaMemcpy(aug, dAug, sizeof(aug), cudaMemcpyDeviceToHost);
const float* Ainv = &aug[N*N]; // right half holds A⁻¹
float res = inverse_residual(N, A, Ainv);
printf(" inv (3x3) ||A·A⁻¹ - I||_max = %.2e %s\n",
res, res < 1e-5f ? "ok" : "FAIL");
bool ok = (res < 1e-5f);
// ── zero LEADING pivot: A invertible (det = -2) but A[0,0] = 0 ──
float Z[NP*NP] = { 0, 2, // col 0
1, 3 }; // col 1
float augZ[NP*2*NP];
// plain inv: divides by A[0,0] = 0 -> Inf/NaN contaminate the result.
make_augmented(NP, Z, augZ);
float* dAugZ; cudaMalloc(&dAugZ, sizeof(augZ));
cudaMemcpy(dAugZ, augZ, sizeof(augZ), cudaMemcpyHostToDevice);
k_inv_plain_np<<<1, 64>>>(dAugZ);
cudaDeviceSynchronize();
cudaMemcpy(augZ, dAugZ, sizeof(augZ), cudaMemcpyDeviceToHost);
bool nonfinite = false;
for (int i = 0; i < NP*NP; i++) nonfinite |= !isfinite(augZ[NP*NP + i]);
printf(" inv (zero pivot) non-finite result: %s (expected — plain inv mishandles it)\n",
nonfinite ? "yes" : "no");
ok = ok && nonfinite;
// inv_pivoted: row-swaps the largest |pivot| up first -> exact inverse.
make_augmented(NP, Z, augZ);
cudaMemcpy(dAugZ, augZ, sizeof(augZ), cudaMemcpyHostToDevice);
k_inv_pivoted_np<<<1, 64>>>(dAugZ);
cudaDeviceSynchronize();
cudaMemcpy(augZ, dAugZ, sizeof(augZ), cudaMemcpyDeviceToHost);
res = inverse_residual(NP, Z, &augZ[NP*NP]);
printf(" inv_pivoted (same A) ||A·A⁻¹ - I||_max = %.2e %s\n",
res, res < 1e-5f ? "ok" : "FAIL");
ok = ok && (res < 1e-5f);
cudaFree(dAug); cudaFree(dAugZ);
printf(ok ? "PASS\n" : "FAIL\n");
return ok ? 0 : 1;
}
13_thread_pack#
// 13_thread_pack.cu — the glass::thread:: tier: 32 low-DOF SPD solves packed per warp.
//
// Build (from this examples/ dir, pure SIMT — no MathDx needed):
// nvcc -std=c++17 -arch=sm_75 -I.. 13_thread_pack.cu -o thread_pack && ./thread_pack
//
// The low-DOF corner: at N=6 a warp-per-problem factor leaves ~26 of 32 lanes idle
// on the serial pivot steps. glass::thread:: flips the mapping — ONE problem per
// THREAD, sequential (no barriers, no shuffles, no threadIdx read inside the op),
// so a warp carries 32 independent problems at once. The price: compile-time sizes
// only, and the operands must stay register-resident (measured ceiling N<=7 —
// past it the thread-local arrays spill and the tier's premise is gone).
//
// The kernel is the caller's: it stages each problem global -> thread-local
// registers, runs the op on its own arrays, and writes back. Launch shape comes
// from glass::suggested_threads_per_block<>() (a seed heuristic, not a measured
// table entry — see glass-defaults.cuh).
#include "glass.cuh"
#include "glass-defaults.cuh"
#include <cstdio>
#include <cmath>
#include <cuda_runtime.h>
constexpr int N = 6; // per-problem size: a 6-DOF arm's normal equations
constexpr int P = 4096; // independent problems
__global__ void k_thread_posv(const float* A, const float* rhs, float* x, int np) {
int p = blockIdx.x * blockDim.x + threadIdx.x;
if (p >= np) return; // ragged tail: fine — no barriers inside
float a[N * N], b[N]; // register-resident at N<=7
for (int i = 0; i < N * N; i++) a[i] = A[(size_t)p * N * N + i];
for (int i = 0; i < N; i++) b[i] = rhs[(size_t)p * N + i];
glass::thread::posv<float, N>(a, b); // factor + fwd/back solve, one thread
for (int i = 0; i < N; i++) x[(size_t)p * N + i] = b[i];
}
int main() {
// Host: build P well-conditioned SPD systems (A = M Mᵀ + N·I, column-major).
static float hA[P * N * N], hb[P * N], hx[P * N];
for (int p = 0; p < P; p++) {
float M[N * N];
for (int i = 0; i < N * N; i++) M[i] = (float)(((p * 131 + i * 7919) % 200) - 100) / 100.f;
for (int c = 0; c < N; c++)
for (int r = 0; r < N; r++) {
float s = 0.f;
for (int k = 0; k < N; k++) s += M[r + k * N] * M[c + k * N];
hA[p * N * N + r + c * N] = s + (r == c ? (float)N : 0.f);
}
for (int i = 0; i < N; i++) hb[p * N + i] = (float)((p + i) % 5) - 2.f;
}
float *dA, *db, *dx;
cudaMalloc(&dA, sizeof(hA)); cudaMalloc(&db, sizeof(hb)); cudaMalloc(&dx, sizeof(hb));
cudaMemcpy(dA, hA, sizeof(hA), cudaMemcpyHostToDevice);
cudaMemcpy(db, hb, sizeof(hb), cudaMemcpyHostToDevice);
// One problem per thread; TPB from the defaults heuristic (N=6 -> 64).
constexpr uint32_t TPB = glass::suggested_threads_per_block<glass::op::posv, N, float>();
k_thread_posv<<<(P + TPB - 1) / TPB, TPB>>>(dA, db, dx, P);
cudaMemcpy(hx, dx, sizeof(hx), cudaMemcpyDeviceToHost);
// Verify every problem: ||A x - b||_inf against the untouched host copies.
float worst = 0.f;
for (int p = 0; p < P; p++) {
for (int r = 0; r < N; r++) {
float s = 0.f;
for (int c = 0; c < N; c++) s += hA[p * N * N + r + c * N] * hx[p * N + c];
worst = fmaxf(worst, fabsf(s - hb[p * N + r]));
}
}
printf("thread-packed posv: %d problems of N=%d at TPB=%u, worst |Ax-b| = %.2e -> %s\n",
P, N, TPB, worst, worst < 1e-3f ? "PASS" : "FAIL");
cudaFree(dA); cudaFree(db); cudaFree(dx);
return worst < 1e-3f ? 0 : 1;
}
14_spatial_dynamics#
// 14_spatial_dynamics.cu — Featherstone spatial cross products: the RNEA inner loop.
//
// Build: nvcc -std=c++17 -arch=sm_75 -I.. 14_spatial_dynamics.cu -o spatial_dynamics && ./spatial_dynamics
//
// USE CASE (rigid-body dynamics): every RNEA/ABA velocity/acceleration sweep is
// built from `v ×ₘ x` and `v ×* f` — a typical generated dynamics suite calls
// them at ~70 sites. This example runs one representative sweep step both ways:
//
// FUSED glass::motion_cross_mul / force_cross_mul — each output row is its
// 2-4 term formula; no 6x6 ever exists.
// COMPOSED glass::motion_cross → 6x6 in shared memory → glass::gemv.
//
// The two agree to rounding (the fused op IS the composed product, minus the
// 36-element temporary and the extra barrier) — that equivalence is also a
// pytest gate (test_robotics.py). A hand-rolled per-project copy of these row
// formulas computes the same thing at the same speed; what it can't give you is
// the pinned convention ([ω; v] angular-first, column-major), the sign-identity
// test suite (crf == −crmᵀ, the dual identity), and the three tiers for free.
//
// Convention: spatial vectors are ANGULAR-FIRST [ω(3); v(3)] (Featherstone) —
// see docs/source/user_guide/concepts/robotics_conventions.rst.
#include "glass.cuh"
#include <cstdio>
#include <cmath>
#include <cuda_runtime.h>
constexpr int P = 512; // links × timesteps in flight
// FUSED: y = crm(v)·x, f_out = crf(v)·f — no matrices materialized.
__global__ void k_fused(const float* v, const float* x, const float* f,
float* y, float* fo) {
int p = blockIdx.x;
glass::block::motion_cross_mul<float>(1.f, v + 6*p, x + 6*p, 0.f, y + 6*p);
glass::block::force_cross_mul<float>(1.f, v + 6*p, f + 6*p, 0.f, fo + 6*p);
}
// COMPOSED: materialize the 6x6s in shared memory, then gemv.
__global__ void k_composed(const float* v, const float* x, const float* f,
float* y, float* fo) {
__shared__ float M[36], F[36];
int p = blockIdx.x;
glass::block::motion_cross<float>(v + 6*p, M);
glass::block::force_cross<float>(v + 6*p, F);
glass::block::gemv<float, 6, 6>(1.f, M, x + 6*p, 0.f, y + 6*p);
glass::block::gemv<float, 6, 6>(1.f, F, f + 6*p, 0.f, fo + 6*p);
}
int main() {
static float hv[P*6], hx[P*6], hf[P*6];
for (int i = 0; i < P*6; i++) {
hv[i] = (float)((int)((i*2654435761u >> 8) % 2000) - 1000) / 500.f;
hx[i] = (float)((int)((i*40503u >> 4) % 2000) - 1000) / 500.f;
hf[i] = (float)((int)((i*9973u >> 2) % 2000) - 1000) / 500.f;
}
float *dv, *dx, *df, *dy1, *df1, *dy2, *df2;
cudaMalloc(&dv, sizeof(hv)); cudaMalloc(&dx, sizeof(hx)); cudaMalloc(&df, sizeof(hf));
cudaMalloc(&dy1, sizeof(hx)); cudaMalloc(&df1, sizeof(hf));
cudaMalloc(&dy2, sizeof(hx)); cudaMalloc(&df2, sizeof(hf));
cudaMemcpy(dv, hv, sizeof(hv), cudaMemcpyHostToDevice);
cudaMemcpy(dx, hx, sizeof(hx), cudaMemcpyHostToDevice);
cudaMemcpy(df, hf, sizeof(hf), cudaMemcpyHostToDevice);
k_fused<<<P, 32>>>(dv, dx, df, dy1, df1);
k_composed<<<P, 32>>>(dv, dx, df, dy2, df2);
cudaDeviceSynchronize();
static float y1[P*6], y2[P*6], f1[P*6], f2[P*6];
cudaMemcpy(y1, dy1, sizeof(y1), cudaMemcpyDeviceToHost);
cudaMemcpy(y2, dy2, sizeof(y2), cudaMemcpyDeviceToHost);
cudaMemcpy(f1, df1, sizeof(f1), cudaMemcpyDeviceToHost);
cudaMemcpy(f2, df2, sizeof(f2), cudaMemcpyDeviceToHost);
float maxerr = 0.f;
for (int i = 0; i < P*6; i++) {
maxerr = fmaxf(maxerr, fabsf(y1[i] - y2[i]));
maxerr = fmaxf(maxerr, fabsf(f1[i] - f2[i]));
}
// host identity check: crf(v)·f == −crm(v)ᵀ·f (sign structure of the duals)
printf("fused vs composed max |diff| = %.3g -> %s\n", maxerr,
maxerr < 1e-5f ? "PASS" : "FAIL");
return maxerr < 1e-5f ? 0 : 1;
}
15_floating_base_retract#
// 15_floating_base_retract.cu — batched SE(3) manifold integration, one state per THREAD.
//
// Build: nvcc -std=c++17 -arch=sm_75 -I.. 15_floating_base_retract.cu -o floating_base_retract && ./floating_base_retract
//
// USE CASE (floating-base dynamics / sampling control): a floating-base
// integrator step is a MANIFOLD update — position+quaternion pose ⊞ body twist
// — not a vector add. Naive `q += ω·dt` drifts off the unit sphere and off
// SO(3); the correct step is `glass::se3_retract` (exp on the group, matching
// Pinocchio's `integrate`). A batched rollout engine holds thousands of
// independent states, so the natural mapping is the THREAD tier: one state per
// thread, 32 states per warp, no barriers.
//
// This example integrates P=4096 rigid-body states through K steps of a
// constant body twist and checks the two invariants a hand-rolled integrator
// classically gets wrong:
// 1. the quaternion stays unit to rounding after K compounding steps
// (the retract renormalizes; additive updates diverge), and
// 2. integrating a constant twist for K steps equals ONE retract of K·(twist)
// for the rotation part (one-parameter-subgroup property of exp — a
// composed-vs-fused gate a Euler-style additive update fails badly).
#include "glass.cuh"
#include <cstdio>
#include <cmath>
#include <cuda_runtime.h>
constexpr int P = 4096;
constexpr int K = 200;
constexpr float DT = 0.01f;
__global__ void k_integrate(const float* pose0, const float* rho, const float* phi,
float* pose_out) {
int p = blockIdx.x * blockDim.x + threadIdx.x;
if (p >= P) return;
float pose[7], nxt[7], r[3], w[3];
for (int i = 0; i < 7; i++) pose[i] = pose0[7*p + i];
for (int i = 0; i < 3; i++) { r[i] = rho[3*p + i]*DT; w[i] = phi[3*p + i]*DT; }
for (int k = 0; k < K; k++) {
glass::thread::se3_retract<float>(pose, r, w, nxt);
for (int i = 0; i < 7; i++) pose[i] = nxt[i];
}
for (int i = 0; i < 7; i++) pose_out[7*p + i] = pose[i];
}
__global__ void k_one_shot(const float* pose0, const float* phi, float* q_out) {
int p = blockIdx.x * blockDim.x + threadIdx.x;
if (p >= P) return;
float w[3] = {phi[3*p]*DT*K, phi[3*p + 1]*DT*K, phi[3*p + 2]*DT*K};
float qn[4];
glass::thread::quat_retract<float>(pose0 + 7*p + 3, w, qn);
for (int i = 0; i < 4; i++) q_out[4*p + i] = qn[i];
}
int main() {
static float h0[P*7], hr[P*3], hw[P*3];
for (int p = 0; p < P; p++) {
// random-ish unit quaternion + position
float q[4], n = 0;
for (int i = 0; i < 4; i++) { q[i] = (float)((p*31 + i*97) % 200 - 100)/100.f + 0.01f; n += q[i]*q[i]; }
n = sqrtf(n);
h0[7*p + 0] = 0.1f*p; h0[7*p + 1] = -0.05f*p; h0[7*p + 2] = 0.5f;
for (int i = 0; i < 4; i++) h0[7*p + 3 + i] = q[i]/n;
for (int i = 0; i < 3; i++) {
hr[3*p + i] = (float)((p + i*7) % 100 - 50)/50.f;
hw[3*p + i] = (float)((p*3 + i*11) % 100 - 50)/100.f; // |ω| modest
}
}
float *d0, *dr, *dw, *dout, *dq1;
cudaMalloc(&d0, sizeof(h0)); cudaMalloc(&dr, sizeof(hr)); cudaMalloc(&dw, sizeof(hw));
cudaMalloc(&dout, sizeof(h0)); cudaMalloc(&dq1, P*4*sizeof(float));
cudaMemcpy(d0, h0, sizeof(h0), cudaMemcpyHostToDevice);
cudaMemcpy(dr, hr, sizeof(hr), cudaMemcpyHostToDevice);
cudaMemcpy(dw, hw, sizeof(hw), cudaMemcpyHostToDevice);
k_integrate<<<(P + 255)/256, 256>>>(d0, dr, dw, dout);
k_one_shot<<<(P + 255)/256, 256>>>(d0, dw, dq1);
cudaDeviceSynchronize();
static float out[P*7], q1[P*4];
cudaMemcpy(out, dout, sizeof(out), cudaMemcpyDeviceToHost);
cudaMemcpy(q1, dq1, sizeof(q1), cudaMemcpyDeviceToHost);
float max_norm_err = 0.f, max_sub_err = 0.f;
for (int p = 0; p < P; p++) {
float n = 0;
for (int i = 0; i < 4; i++) n += out[7*p + 3 + i]*out[7*p + 3 + i];
max_norm_err = fmaxf(max_norm_err, fabsf(sqrtf(n) - 1.f));
// one-parameter subgroup: K small steps == one K·(ω dt) retract (sign-free)
float dot = 0;
for (int i = 0; i < 4; i++) dot += out[7*p + 3 + i]*q1[4*p + i];
max_sub_err = fmaxf(max_sub_err, fabsf(fabsf(dot) - 1.f));
}
printf("unit-norm drift after %d steps: %.3g subgroup gap: %.3g -> %s\n",
K, max_norm_err, max_sub_err,
(max_norm_err < 1e-5f && max_sub_err < 1e-4f) ? "PASS" : "FAIL");
return (max_norm_err < 1e-5f && max_sub_err < 1e-4f) ? 0 : 1;
}
16_mppi_weights#
// 16_mppi_weights.cu — the path-integral (MPPI) weight update: softmax + argmin.
//
// Build: nvcc -std=c++17 -arch=sm_75 -I.. 16_mppi_weights.cu -o mppi_weights && ./mppi_weights
//
// USE CASE (sampling-based control): an MPPI/CEM controller rolls out N
// perturbed control sequences, scores each with a cost J_i, and blends them
// with the exponentially-weighted average
// w_i = exp(-λ(J_i - min J)) / Σ_j exp(-λ(J_j - min J)).
// That IS `glass::softmax(n, -λ, J, w, scratch)` — the baseline subtraction is
// the max shift, so the weights are overflow-safe and shift-invariant — plus
// `glass::argmin` for the best-rollout index (warm starts, elite selection).
// One BLOCK owns one controller instance here (the block tier); a batched
// multi-robot controller would drop the same calls one tier down.
//
// Every sampling planner hand-rolls this pair (baseline subtraction, the
// normalizer, the argmin tie rule). The GLASS ops pin the numerics: the
// reductions run fixed-order trees, so the weights are BIT-IDENTICAL at any
// block size — a property this example checks directly (64 vs 256 threads).
#include "glass.cuh"
#include <cstdio>
#include <cmath>
#include <cuda_runtime.h>
constexpr int NROLL = 480; // rollouts per controller
constexpr float LAMBDA = 3.0f; // temperature
__global__ void k_weights(const float* J, float* w, unsigned int* best) {
extern __shared__ float scr[];
glass::block::softmax<float>(NROLL, -LAMBDA, J, w, scr);
glass::block::argmin<float>(NROLL, J, best, scr);
}
int main() {
static float hJ[NROLL];
for (int i = 0; i < NROLL; i++)
hJ[i] = 5.f + 3.f*sinf(0.1f*i) + (float)((i*2654435761u >> 7) % 1000)/500.f;
float *dJ, *dw64, *dw256;
unsigned int* dbest;
cudaMalloc(&dJ, sizeof(hJ));
cudaMalloc(&dw64, sizeof(hJ));
cudaMalloc(&dw256, sizeof(hJ));
cudaMalloc(&dbest, sizeof(unsigned int));
cudaMemcpy(dJ, hJ, sizeof(hJ), cudaMemcpyHostToDevice);
size_t smem = glass::block::softmax_scratch_bytes<float>(NROLL);
size_t arg_smem = glass::block::argreduce_scratch_bytes<float>(256);
if (arg_smem > smem) smem = arg_smem; // the kernel reuses one buffer for both ops
k_weights<<<1, 64, smem>>>(dJ, dw64, dbest);
k_weights<<<1, 256, smem>>>(dJ, dw256, dbest);
cudaDeviceSynchronize();
static float w64[NROLL], w256[NROLL];
unsigned int best;
cudaMemcpy(w64, dw64, sizeof(w64), cudaMemcpyDeviceToHost);
cudaMemcpy(w256, dw256, sizeof(w256), cudaMemcpyDeviceToHost);
cudaMemcpy(&best, dbest, sizeof(best), cudaMemcpyDeviceToHost);
// host reference (double): the exact MPPI weights
double m = hJ[0];
for (int i = 1; i < NROLL; i++) m = fmin(m, (double)hJ[i]);
double Z = 0; int href = 0;
static double wref[NROLL];
for (int i = 0; i < NROLL; i++) {
wref[i] = exp(-LAMBDA*((double)hJ[i] - m));
Z += wref[i];
if (hJ[i] < hJ[href]) href = i;
}
float sum = 0, maxerr = 0;
bool bitinv = true;
for (int i = 0; i < NROLL; i++) {
sum += w64[i];
maxerr = fmaxf(maxerr, fabsf(w64[i] - (float)(wref[i]/Z)));
bitinv = bitinv && (w64[i] == w256[i]);
}
printf("Σw = %.6f max |w - ref| = %.3g best = %u (ref %d) "
"bit-identical 64 vs 256 threads: %s\n",
sum, maxerr, best, href, bitinv ? "yes" : "NO");
bool pass = fabsf(sum - 1.f) < 1e-5f && maxerr < 1e-6f
&& best == (unsigned int)href && bitinv;
printf("%s\n", pass ? "PASS" : "FAIL");
return pass ? 0 : 1;
}
17_cone_projection#
// 17_cone_projection.cu — friction-cone AL constraint step: soc_project + interval AL.
//
// Build: nvcc -std=c++17 -arch=sm_75 -I.. 17_cone_projection.cu -o cone_projection && ./cone_projection
//
// USE CASE (constrained trajectory optimization): a contact-rich solver carries
// friction-cone rows (‖f_tangential‖ <= μ·f_normal — a second-order cone) and
// state/control interval rows. The PHR augmented-Lagrangian machinery for both
// reduces to small per-row scalar ops:
// * `glass::soc_project` — the Euclidean cone projection (multiplier update
// λ ← Π_K(λ − ρg), and the AL gradient is −Π_K(λ − ρg));
// * `glass::al_soc_value` / `glass::al_interval_value` /
// `glass::al_interval_grad_hess` — the merit value and its GN derivatives.
// One (row-group, knot) pair per THREAD — exactly how GATO's bsqp solver runs
// this set (these ops are its cone wave, promoted).
//
// Checks: projection case split vs a host reference, idempotence
// (Π(Π(w)) == Π(w)), the convex-projection orthogonality <Π(w), Π(w) − w> == 0,
// and the m=1 degeneration of the conic AL value to the scalar hinge.
#include "glass.cuh"
#include <cstdio>
#include <cmath>
#include <cuda_runtime.h>
constexpr int P = 2048; // (row-group, knot) pairs
constexpr int M = 4; // cone dimension: normal + 3 tangential rows
constexpr float RHO = 1.7f;
__global__ void k_al_step(const float* g, const float* lam, float* proj,
float* val, float* hinge_gap) {
int p = blockIdx.x * blockDim.x + threadIdx.x;
if (p >= P) return;
// multiplier-update projection: Π_K(λ − ρg), then project AGAIN (idempotence
// is checked on host via this second image equaling the first)
float w[M], pr[M];
for (int i = 0; i < M; i++) w[i] = lam[M*p + i] - RHO*g[M*p + i];
glass::thread::soc_project<float>(w, pr, M);
for (int i = 0; i < M; i++) proj[M*p + i] = pr[i];
// AL merit value for the cone row
val[p] = glass::block::al_soc_value<float>(g + M*p, lam + M*p, RHO, M);
// m=1 cone == the g >= 0 hinge: the conic value must equal al_hinge_value
// on the sign convention bridge (hinge feasible c <= 0 ⇔ cone g >= 0).
float g1 = g[M*p], l1 = lam[M*p];
float conic = glass::block::al_soc_value<float>(&g1, &l1, RHO, 1);
float hinge = glass::block::al_hinge_value<float>(-g1, l1, RHO, 0.f);
hinge_gap[p] = fabsf(conic - hinge);
}
int main() {
static float hg[P*M], hl[P*M];
for (int i = 0; i < P*M; i++) {
hg[i] = (float)((int)((i*2654435761u >> 6) % 2000) - 1000)/700.f;
hl[i] = (float)((i*40503u >> 3) % 1000)/500.f;
}
float *dg, *dl, *dp, *dv, *dh;
cudaMalloc(&dg, sizeof(hg)); cudaMalloc(&dl, sizeof(hl));
cudaMalloc(&dp, sizeof(hg)); cudaMalloc(&dv, P*sizeof(float));
cudaMalloc(&dh, P*sizeof(float));
cudaMemcpy(dg, hg, sizeof(hg), cudaMemcpyHostToDevice);
cudaMemcpy(dl, hl, sizeof(hl), cudaMemcpyHostToDevice);
k_al_step<<<(P + 127)/128, 128>>>(dg, dl, dp, dv, dh);
cudaDeviceSynchronize();
static float proj[P*M], val[P], hinge_gap[P];
cudaMemcpy(proj, dp, sizeof(proj), cudaMemcpyDeviceToHost);
cudaMemcpy(val, dv, sizeof(val), cudaMemcpyDeviceToHost);
cudaMemcpy(hinge_gap, dh, sizeof(hinge_gap), cudaMemcpyDeviceToHost);
float max_orth = 0, max_hinge = 0, max_proj = 0;
for (int p = 0; p < P; p++) {
// host reference projection + orthogonality <Π(w), Π(w) − w>
double w[M], r = 0;
for (int i = 0; i < M; i++) w[i] = (double)hl[M*p + i] - RHO*(double)hg[M*p + i];
for (int i = 1; i < M; i++) r += w[i]*w[i];
r = sqrt(r);
double ref[M];
if (r <= w[0]) for (int i = 0; i < M; i++) ref[i] = w[i];
else if (r <= -w[0]) for (int i = 0; i < M; i++) ref[i] = 0;
else {
double a = 0.5*(w[0] + r);
ref[0] = a;
for (int i = 1; i < M; i++) ref[i] = (a/r)*w[i];
}
double orth = 0;
for (int i = 0; i < M; i++) {
max_proj = fmaxf(max_proj, fabsf(proj[M*p + i] - (float)ref[i]));
orth += (double)proj[M*p + i]*((double)proj[M*p + i] - w[i]);
}
max_orth = fmaxf(max_orth, fabsf((float)orth));
max_hinge = fmaxf(max_hinge, hinge_gap[p]);
}
printf("proj err %.3g orthogonality %.3g m=1 hinge gap %.3g\n",
max_proj, max_orth, max_hinge);
bool pass = max_proj < 1e-5f && max_orth < 1e-3f && max_hinge < 1e-6f;
printf("%s\n", pass ? "PASS" : "FAIL");
return pass ? 0 : 1;
}
18_collision_spheres#
// 18_collision_spheres.cu — sphere-based narrow phase: transform, distance, cost.
//
// Build: nvcc -std=c++17 -arch=sm_75 -I.. 18_collision_spheres.cu -o collision_spheres && ./collision_spheres
//
// USE CASE (motion generation / IK): the dominant GPU robot-collision
// representation decomposes the robot into spheres and scores each against
// world primitives — one (sphere, obstacle) pair per THREAD. A narrow-phase
// check is three GLASS calls:
// glass::transform_sphere FK pose applied to the link sphere (radius kept)
// glass::sphere_box_dist signed distance to a box, in the box frame
// (an OBB = rotate the center into the box frame
// with quat_rotate on the inverse pose)
// glass::smooth_hinge the C¹ activation turning distance into cost
// plus glass::sphere_sphere_dist for self-collision pairs. The gradient
// pipeline is the same three calls' grad outputs chained.
//
// Checks: the device signed distances against a double host reference over
// mixed inside/outside cases, and cost/gradient consistency of the hinge
// (cost decreasing in distance, zero beyond the activation band).
#include "glass.cuh"
#include <cstdio>
#include <cmath>
#include <cuda_runtime.h>
constexpr int P = 4096;
constexpr float ETA = 0.05f; // activation width
__global__ void k_narrow_phase(const float* q, const float* t, const float* sph,
const float* half, float* dist, float* cost, float* grad) {
int p = blockIdx.x * blockDim.x + threadIdx.x;
if (p >= P) return;
// link sphere -> world (FK pose), world == box frame here for brevity
float world[4];
glass::block::transform_sphere<float>(q + 4*p, t + 3*p, sph + 4*p, world);
float g[3];
float d = glass::block::sphere_box_dist<float>(world, world[3], half + 3*p, g);
dist[p] = d;
cost[p] = glass::block::smooth_hinge<float>(d, ETA);
float dcdd = glass::block::smooth_hinge_grad<float>(d, ETA);
for (int i = 0; i < 3; i++) grad[3*p + i] = dcdd*g[i]; // chain rule to the center
}
int main() {
static float hq[P*4], ht[P*3], hs[P*4], hh[P*3];
for (int p = 0; p < P; p++) {
float n = 0;
for (int i = 0; i < 4; i++) { hq[4*p + i] = (float)((p*13 + i*57) % 200 - 100)/100.f + 0.02f; n += hq[4*p + i]*hq[4*p + i]; }
for (int i = 0; i < 4; i++) hq[4*p + i] /= sqrtf(n);
for (int i = 0; i < 3; i++) {
ht[3*p + i] = (float)((p*7 + i*3) % 100 - 50)/40.f;
hs[4*p + i] = (float)((p*5 + i*11) % 60 - 30)/60.f;
hh[3*p + i] = 0.4f + (float)((p + i) % 40)/50.f;
}
hs[4*p + 3] = 0.03f + (float)(p % 20)/200.f; // radius
}
float *dq, *dt, *ds, *dh, *dd, *dc, *dg;
cudaMalloc(&dq, sizeof(hq)); cudaMalloc(&dt, sizeof(ht));
cudaMalloc(&ds, sizeof(hs)); cudaMalloc(&dh, sizeof(hh));
cudaMalloc(&dd, P*sizeof(float)); cudaMalloc(&dc, P*sizeof(float));
cudaMalloc(&dg, P*3*sizeof(float));
cudaMemcpy(dq, hq, sizeof(hq), cudaMemcpyHostToDevice);
cudaMemcpy(dt, ht, sizeof(ht), cudaMemcpyHostToDevice);
cudaMemcpy(ds, hs, sizeof(hs), cudaMemcpyHostToDevice);
cudaMemcpy(dh, hh, sizeof(hh), cudaMemcpyHostToDevice);
k_narrow_phase<<<(P + 255)/256, 256>>>(dq, dt, ds, dh, dd, dc, dg);
cudaDeviceSynchronize();
static float dist[P], cost[P];
cudaMemcpy(dist, dd, sizeof(dist), cudaMemcpyDeviceToHost);
cudaMemcpy(cost, dc, sizeof(cost), cudaMemcpyDeviceToHost);
// host reference in double
float maxerr = 0; int inside = 0; bool cost_ok = true;
for (int p = 0; p < P; p++) {
double x = hq[4*p], y = hq[4*p + 1], z = hq[4*p + 2], w = hq[4*p + 3];
double c[3] = {hs[4*p], hs[4*p + 1], hs[4*p + 2]}, rc[3];
// p' = p + 2w(v×p) + 2v×(v×p)
double cr1[3] = {y*c[2] - z*c[1], z*c[0] - x*c[2], x*c[1] - y*c[0]};
for (int i = 0; i < 3; i++) rc[i] = c[i] + 2*w*cr1[i];
double cr2[3] = {y*cr1[2] - z*cr1[1], z*cr1[0] - x*cr1[2], x*cr1[1] - y*cr1[0]};
for (int i = 0; i < 3; i++) rc[i] += 2*cr2[i] + ht[3*p + i];
double q[3], omax = -1e30, osum = 0;
for (int i = 0; i < 3; i++) {
q[i] = fabs(rc[i]) - hh[3*p + i];
omax = fmax(omax, q[i]);
osum += fmax(q[i], 0.0)*fmax(q[i], 0.0);
}
double ref = sqrt(osum) + fmin(omax, 0.0) - hs[4*p + 3];
maxerr = fmaxf(maxerr, fabsf(dist[p] - (float)ref));
if (ref < 0) inside++;
float want = dist[p] <= 0 ? -dist[p] + ETA/2
: (dist[p] >= ETA ? 0.f : (dist[p] - ETA)*(dist[p] - ETA)/(2*ETA));
cost_ok = cost_ok && fabsf(cost[p] - want) < 1e-6f;
}
printf("max |dist - ref| = %.3g over %d pairs (%d penetrating) hinge: %s\n",
maxerr, P, inside, cost_ok ? "ok" : "BROKEN");
bool pass = maxerr < 2e-5f && cost_ok && inside > 0;
printf("%s\n", pass ? "PASS" : "FAIL");
return pass ? 0 : 1;
}
19_best_fit_rotation#
// 19_best_fit_rotation.cu — batched point alignment + rotation cleanup (est kit).
//
// Build: nvcc -std=c++17 -arch=sm_75 -I.. 19_best_fit_rotation.cu -o best_fit_rotation && ./best_fit_rotation
//
// USE CASE (estimation / registration): the inner op of ICP and point-cloud
// alignment is Wahba's problem — given correspondences (a_i, b_i), find the
// rotation minimizing Σ‖b_i − R·a_i‖². The GPU-batched recipe, one problem
// per THREAD (thousands of independent alignments in one launch):
// 1. accumulate the 3x3 cross covariance M = Σ b_i·a_iᵀ (nine FMAs per
// correspondence — glass::ger shaped, done inline here),
// 2. glass::thread::closest_rotation(M, R) — the SVD-based Kabsch solution
// with the det fix (never returns a reflection).
// The SAME op is the rotation-matrix re-orthonormalizer: feed a drifted
// product of many incremental rotations and it returns the nearest proper
// rotation — the classic cleanup after long integrations.
//
// Checks (self-verifying): recovered rotations match the known ground truth
// to f32 tolerance for every problem; re-orthonormalized matrices are
// orthonormal with det +1 and stay near the drifted input.
#include "glass.cuh"
#include <cstdio>
#include <cmath>
#include <cuda_runtime.h>
constexpr int P = 4096; // independent alignment problems
constexpr int NC = 8; // correspondences per problem
__global__ void k_align(const float* a, const float* b, float* R) {
int p = blockIdx.x * blockDim.x + threadIdx.x;
if (p >= P) return;
// M = Σ b_i·a_iᵀ (column-major: M[c*3 + r] += b[r]·a[c])
float M[9];
for (int i = 0; i < 9; i++) M[i] = 0.f;
for (int i = 0; i < NC; i++) {
const float* ai = a + (size_t)p*NC*3 + 3*i;
const float* bi = b + (size_t)p*NC*3 + 3*i;
for (int c = 0; c < 3; c++)
for (int r = 0; r < 3; r++)
M[c*3 + r] += bi[r] * ai[c];
}
glass::thread::closest_rotation<float>(M, R + (size_t)p*9);
}
__global__ void k_cleanup(const float* A, float* R) {
int p = blockIdx.x * blockDim.x + threadIdx.x;
if (p >= P) return;
glass::thread::closest_rotation<float>(A + (size_t)p*9, R + (size_t)p*9);
}
// host-side rotation from an axis/angle (double, column-major)
static void rot_from_aa(const double* ax, double th, double* R) {
double n = sqrt(ax[0]*ax[0] + ax[1]*ax[1] + ax[2]*ax[2]);
double u[3] = {ax[0]/n, ax[1]/n, ax[2]/n};
double c = cos(th), s = sin(th), oc = 1.0 - c;
R[0] = c + u[0]*u[0]*oc; R[3] = u[0]*u[1]*oc - u[2]*s; R[6] = u[0]*u[2]*oc + u[1]*s;
R[1] = u[1]*u[0]*oc + u[2]*s; R[4] = c + u[1]*u[1]*oc; R[7] = u[1]*u[2]*oc - u[0]*s;
R[2] = u[2]*u[0]*oc - u[1]*s; R[5] = u[2]*u[1]*oc + u[0]*s; R[8] = c + u[2]*u[2]*oc;
}
int main() {
static double Rtrue[P][9];
float *ha = (float*)malloc((size_t)P*NC*3*sizeof(float));
float *hb = (float*)malloc((size_t)P*NC*3*sizeof(float));
float *hA = (float*)malloc((size_t)P*9*sizeof(float));
unsigned s = 12345u;
auto frand = [&]() { s = s*1664525u + 1013904223u; return ((s >> 8) / 8388608.0) * 2.0 - 1.0; };
for (int p = 0; p < P; p++) {
double ax[3] = {frand(), frand(), frand()};
rot_from_aa(ax, 0.1 + 2.8 * ((frand() + 1.0) / 2.0), Rtrue[p]);
for (int i = 0; i < NC; i++) {
double ai[3] = {frand(), frand(), frand()};
for (int k = 0; k < 3; k++) ha[(size_t)p*NC*3 + 3*i + k] = (float)ai[k];
for (int r = 0; r < 3; r++) // b = R_true·a (exact correspondences)
hb[(size_t)p*NC*3 + 3*i + r] =
(float)(Rtrue[p][r]*ai[0] + Rtrue[p][3+r]*ai[1] + Rtrue[p][6+r]*ai[2]);
}
for (int k = 0; k < 9; k++) // drifted rotation for the cleanup leg
hA[(size_t)p*9 + k] = (float)(Rtrue[p][k] + 5e-3 * frand());
}
float *da, *db, *dA, *dR1, *dR2;
cudaMalloc(&da, (size_t)P*NC*3*sizeof(float)); cudaMalloc(&db, (size_t)P*NC*3*sizeof(float));
cudaMalloc(&dA, (size_t)P*9*sizeof(float));
cudaMalloc(&dR1, (size_t)P*9*sizeof(float)); cudaMalloc(&dR2, (size_t)P*9*sizeof(float));
cudaMemcpy(da, ha, (size_t)P*NC*3*sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(db, hb, (size_t)P*NC*3*sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(dA, hA, (size_t)P*9*sizeof(float), cudaMemcpyHostToDevice);
k_align<<<(P + 127)/128, 128>>>(da, db, dR1);
k_cleanup<<<(P + 127)/128, 128>>>(dA, dR2);
if (cudaDeviceSynchronize() != cudaSuccess) { printf("FAIL: kernel error\n"); return 1; }
float *hR1 = (float*)malloc((size_t)P*9*sizeof(float));
float *hR2 = (float*)malloc((size_t)P*9*sizeof(float));
cudaMemcpy(hR1, dR1, (size_t)P*9*sizeof(float), cudaMemcpyDeviceToHost);
cudaMemcpy(hR2, dR2, (size_t)P*9*sizeof(float), cudaMemcpyDeviceToHost);
double worst_fit = 0.0, worst_orth = 0.0, worst_det = 0.0, worst_drift = 0.0;
for (int p = 0; p < P; p++) {
for (int k = 0; k < 9; k++) {
double d = fabs((double)hR1[(size_t)p*9 + k] - Rtrue[p][k]);
if (d > worst_fit) worst_fit = d;
double dd = fabs((double)hR2[(size_t)p*9 + k] - Rtrue[p][k]);
if (dd > worst_drift) worst_drift = dd;
}
const float* R = hR2 + (size_t)p*9; // orthonormality + det of the cleanup
double det = (double)R[0]*((double)R[4]*R[8] - (double)R[5]*R[7])
- (double)R[3]*((double)R[1]*R[8] - (double)R[2]*R[7])
+ (double)R[6]*((double)R[1]*R[5] - (double)R[2]*R[4]);
if (fabs(det - 1.0) > worst_det) worst_det = fabs(det - 1.0);
for (int i = 0; i < 3; i++)
for (int j = 0; j < 3; j++) {
double dot = 0.0;
for (int k = 0; k < 3; k++) dot += (double)R[i*3 + k]*R[j*3 + k];
double e = fabs(dot - (i == j ? 1.0 : 0.0));
if (e > worst_orth) worst_orth = e;
}
}
printf("align: worst |R - R_true| = %.3e (exact correspondences)\n", worst_fit);
printf("clean: worst orthonormality = %.3e, worst |det-1| = %.3e\n", worst_orth, worst_det);
printf("clean: worst |R - R_true| = %.3e (5e-3 drift input)\n", worst_drift);
bool ok = worst_fit < 2e-4 && worst_orth < 1e-5 && worst_det < 1e-5 && worst_drift < 2e-2;
printf("%s\n", ok ? "PASS" : "FAIL");
return ok ? 0 : 1;
}