L3 — Matrix Operations#

General matrix-matrix products (gemm and its tiled / strided / batched / indexed variants), matrix inversion (single, partial-pivoting inv_pivoted, and K-way fused multi-matrix), Cholesky factorization (single and K-way fused), triangular solves, symmetric rank-k / rank-2k updates (syrk / syr2k), symmetric-indefinite LDLᵀ (ldlt), and SPD solves (posv / potrs). The GEMM family follows the standard BLAS convention (C is M×N, contraction K) with TRANSPOSE_A / TRANSPOSE_B operand flags and a single ROW_MAJOR_C output flag — a row-major operand is just a transpose, so per-operand row-major flags were removed; see Backend Dispatch for how the dispatch GEMM chooses a path. Single-warp (glass::warp::) variants of gemm, potrf, trsm, and inv render inline beside their block-scoped sibling; see Warp-scoped operations (glass::warp::).

gemm#

Defines

GLASS_TILE4_HELPERS_DEFINED#

Functions

void tile4_load(const float *p, float &a0, float &a1, float &a2, float &a3)#

void tile4_load(const double *p, double &a0, double &a1, double &a2, double &a3)#

template<typename T> void tile4_load(const T *p, T &a0, T &a1, T &a2, T &a3)#

template<typename T> bool tile4_aligned(const T *p)#

constexpr bool tile4_profitable(uint32_t m)#

template<typename T, bool TRANSPOSE_B, bool HAS_BETA, bool VEC> void gemm_tile4_loop(uint32_t rank, uint32_t size, uint32_t m_, uint32_t n_, uint32_t k_, T alpha, const T *A, const T *B, T beta, T *C)#

template<typename T, bool TRANSPOSE_B, bool HAS_BETA> void gemm_tile4(uint32_t rank, uint32_t size, uint32_t m_, uint32_t n_, uint32_t k_, T alpha, const T *A, const T *B, T beta, T *C)#

template<typename T, uint32_t M, uint32_t N, uint32_t K, bool TRANSPOSE_B, bool HAS_BETA, bool VEC> void gemm_tile4_loop_ct(uint32_t rank, uint32_t size, T alpha, const T *A, const T *B, T beta, T *C)#

template<typename T, uint32_t M, uint32_t N, uint32_t K, bool TRANSPOSE_B, bool HAS_BETA> void gemm_tile4_ct(uint32_t rank, uint32_t size, T alpha, const T *A, const T *B, T beta, T *C)#

template<typename T, uint32_t M, uint32_t N, uint32_t K, bool TRANSPOSE_A, bool TRANSPOSE_B, bool ROW_MAJOR_C> void gemm_impl_ct(uint32_t rank, uint32_t size, T alpha, const T *A, const T *B, T beta, T *C)#

template<typename T, uint32_t M, uint32_t N, uint32_t K, bool TRANSPOSE_A, bool TRANSPOSE_B, bool ROW_MAJOR_C> void gemm_impl_ct(uint32_t rank, uint32_t size, T alpha, const T *A, const T *B, T *C)#

template<typename T, bool TRANSPOSE_A = false, bool TRANSPOSE_B = false, bool ROW_MAJOR_C = false, bool TRAILING_SYNC = true> void gemm(uint32_t m, uint32_t n, uint32_t k, T alpha, const T *A, const T *B, T beta, T *C)#

General matrix-matrix multiply: C = alpha * op(A) * op(B) + beta * C (GEMM).

Standard BLAS convention: C is m×n, contraction k. Runtime-size, single-block, flat-element parallelism: each thread owns output elements strided over the block. NumPy: C = alpha * opA(A) @ opB(B) + beta * C; Eigen: C.noalias() = alpha*(opA(A)*opB(B)) + beta*C;.

Template Parameters:

T – Scalar type.
TRANSPOSE_A – If true, A is k×m and op(A)=Aᵀ (else A is m×k).
TRANSPOSE_B – If true, B is n×k and op(B)=Bᵀ (else B is k×n).
ROW_MAJOR_C – Output storage order (false = column-major / Fortran, LDC=m).

Parameters:

m, n, k – Dimensions: C is m×n, contraction k.
alpha – Scalar multiplier on the product.
A, B – Input matrices (column-major; shapes per the transpose flags).
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out result matrix.

template<typename T, bool TRANSPOSE_A = false, bool TRANSPOSE_B = false, bool ROW_MAJOR_C = false, bool TRAILING_SYNC = true> void gemm(uint32_t m, uint32_t n, uint32_t k, T alpha, const T *A, const T *B, T *C)#

GEMM with implicit beta = 0: C = alpha * op(A) * op(B) (overwrite).

Runtime-size overload that overwrites C (the existing C is not read), avoiding the beta * C term. NumPy: C = alpha * opA(A) @ opB(B).

Template Parameters:

T – Scalar type.
TRANSPOSE_A – If true, A is k×m and op(A)=Aᵀ.
TRANSPOSE_B – If true, B is n×k and op(B)=Bᵀ.
ROW_MAJOR_C – Output storage order (false = column-major).

Parameters:

m, n, k – Dimensions: C is m×n, contraction k.
alpha – Scalar multiplier on the product.
A, B – Input matrices.
C – Output result matrix (overwritten).

template<typename T, uint32_t M, uint32_t N, uint32_t K, bool TRANSPOSE_A = false, bool TRANSPOSE_B = false, bool ROW_MAJOR_C = false, bool TRAILING_SYNC = true> void gemm(T alpha, const T *A, const T *B, T beta, T *C)#

Compile-time-size GEMM: C = alpha * op(A) * op(B) + beta * C (GEMM).

Dimensions are template parameters so the compiler unrolls the inner loop and replaces the el % M / el / M index math with magic-number multiplies. Standard BLAS convention: C is M×N, contraction K. NumPy: C = alpha * opA(A) @ opB(B) + beta * C; Eigen: C.noalias() = alpha*(opA(A)*opB(B)) + beta*C;.

Template Parameters:

T – Scalar type.
M, N, K – C is M×N, contraction K.
TRANSPOSE_A – If true, A is K×M and op(A)=Aᵀ (else A is M×K).
TRANSPOSE_B – If true, B is N×K and op(B)=Bᵀ (else B is K×N).
ROW_MAJOR_C – Output storage order (false = column-major / Fortran, LDC=M).

Parameters:

alpha – Scalar multiplier on the product.
A, B – Input matrices.
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out result matrix.

template<typename T, uint32_t M, uint32_t N, uint32_t K, bool TRANSPOSE_A = false, bool TRANSPOSE_B = false, bool ROW_MAJOR_C = false, bool TRAILING_SYNC = true> void gemm(T alpha, const T *A, const T *B, T *C)#

Compile-time-size GEMM with implicit beta = 0: C = alpha * op(A) * op(B).

Compile-time-size overload that overwrites C (the existing C is not read). NumPy: C = alpha * opA(A) @ opB(B).

Template Parameters:

T – Scalar type.
M, N, K – C is M×N, contraction K.
TRANSPOSE_A – If true, A is K×M and op(A)=Aᵀ.
TRANSPOSE_B – If true, B is N×K and op(B)=Bᵀ.
ROW_MAJOR_C – Output storage order (false = column-major).

Parameters:

alpha – Scalar multiplier on the product.
A, B – Input matrices.
C – Output result matrix (overwritten).

template<typename T, int TILE = 8> void gemm_tiled(uint32_t m, uint32_t n, uint32_t k, T alpha, const T *A, const T *B, T beta, T *C, T *s_A, T *s_B)#

Tiled GEMM with shared-memory staging: C = alpha * A * B + beta * C.

Standard convention, column-major, no transpose. C is m×n, contraction k. Stages TILE-wide column blocks of A (m×TILE) and the matching row blocks of B (TILE×n) into the caller-provided shared scratch, accumulating across tiles. Single-block; best when A/B values can be reused from shared memory. NumPy: C = alpha * A @ B + beta * C.

Template Parameters:

T – Scalar type.
TILE – Column-block width staged per pass.

Parameters:

m, n, k – Dimensions: A is m×k, B is k×n, C is m×n.
alpha – Scalar multiplier on the product.
A, B – Input matrices (column-major).
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out result matrix.
s_A – Shared scratch of m * TILE elements for the A tile.
s_B – Shared scratch of TILE * n elements for the B tile.

template<typename T, int TILE = 8> void gemm_dispatch(uint32_t m, uint32_t n, uint32_t k, T alpha, const T *A, const T *B, T beta, T *C, T *s_A = nullptr, T *s_B = nullptr)#

Auto-dispatching GEMM: C = alpha * A * B + beta * C (column-major).

Selects gemm_tiled when shared-memory scratch is provided and one output element fits per thread (m * n <= blockDim); otherwise falls back to the plain gemm. Standard convention: C is m×n, contraction k. Single-block. NumPy: C = alpha * A @ B + beta * C.

Template Parameters:

T – Scalar type.
TILE – Tile width passed through to gemm_tiled.

Parameters:

m, n, k – Dimensions: A is m×k, B is k×n, C is m×n.
alpha – Scalar multiplier on the product.
A, B – Input matrices (column-major).
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out result matrix.
s_A – Optional shared scratch for the A tile (nullptr selects the plain path).
s_B – Optional shared scratch for the B tile.

template<typename T, uint32_t M, uint32_t N, uint32_t K, bool TRANSPOSE_A = false, bool TRANSPOSE_B = false, bool ROW_MAJOR_C = false> void gemm(T alpha, const T *A, const T *B, T beta, T *C)#

Single-warp compile-time-size GEMM: C = alpha * op(A) * op(B) + beta * C.

Single-thread compile-time-size GEMM: C = alpha * op(A) * op(B) + beta * C.

One 32-lane warp computes the product with flat per-element parallelism (lanes stride over the M*N outputs, serial-K inner loop) — same semantics as the block-scoped compile-time gemm, but scoped to a single warp for warp-per-problem kernels (e.g. 4×4 homogeneous-transform multiplies). No inter-lane communication, no sync. C must not alias A/B.

ONE thread computes the whole product, walking the M*N outputs serially (serial-K inner loop) — same semantics as the block/warp compile-time gemm, reusing the same gemm_impl_ct body with (rank=0, size=1). Generally performs worse.

Warning

The shared body switches to the 4-row tile4 path when tile4_profitable(M) (M % 4 == 0 && M >= 12), which issues float4 / double2 vector loads through a reinterpret_cast. Unlikely to happen in the DOF where the single-thread is valuable

Template Parameters:

T – Scalar type.
M, N, K – C is M×N, contraction K.
TRANSPOSE_A – If true, A is K×M and op(A)=Aᵀ.
TRANSPOSE_B – If true, B is N×K and op(B)=Bᵀ.
ROW_MAJOR_C – Output storage order (false = column-major).
T – Scalar type.
M, N, K – C is M×N, contraction K.
TRANSPOSE_A – If true, A is K×M and op(A)=Aᵀ.
TRANSPOSE_B – If true, B is N×K and op(B)=Bᵀ.
ROW_MAJOR_C – Output storage order (false = column-major).

Parameters:

alpha – Scalar multiplier on the product.
A, B – Input matrices.
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out result matrix.
alpha – Scalar multiplier on the product.
A, B – Input matrices.
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out result matrix.

template<typename T, uint32_t M, uint32_t N, uint32_t K, bool TRANSPOSE_A = false, bool TRANSPOSE_B = false, bool ROW_MAJOR_C = false> void gemm(T alpha, const T *A, const T *B, T *C)#

Single-warp compile-time-size GEMM with implicit beta = 0: C = alpha * op(A) * op(B).

Single-thread compile-time-size GEMM with implicit beta = 0: C = alpha * op(A) * op(B).

Overwrites C (the existing C is not read). Otherwise identical to the beta overload above.

Overwrites C (the existing C is not read). Otherwise identical to the beta overload above, including the tile4 caveat.

Template Parameters:

T – Scalar type.
M, N, K – C is M×N, contraction K.
TRANSPOSE_A – If true, A is K×M and op(A)=Aᵀ.
TRANSPOSE_B – If true, B is N×K and op(B)=Bᵀ.
ROW_MAJOR_C – Output storage order (false = column-major).
T – Scalar type.
M, N, K – C is M×N, contraction K.
TRANSPOSE_A – If true, A is K×M and op(A)=Aᵀ.
TRANSPOSE_B – If true, B is N×K and op(B)=Bᵀ.
ROW_MAJOR_C – Output storage order (false = column-major).

Parameters:

alpha – Scalar multiplier on the product.
A, B – Input matrices.
C – Output result matrix (overwritten).
alpha – Scalar multiplier on the product.
A, B – Input matrices.
C – Output result matrix.

template<typename T> struct tile4_has_vec#

Public Static Attributes

static constexpr bool value = false#

template<> struct tile4_has_vec<float>#

Public Static Attributes

static constexpr bool value = true#

template<> struct tile4_has_vec<double>#

Public Static Attributes

static constexpr bool value = true#

namespace warp

namespace thread

Contraction-parallel gemm (`gemm_reduced`)#

One warp owns each output element and its lanes split the contraction — a thread-utilization variant for small outputs. See Contraction-parallel ops (the *_reduced family) for the honest win-condition.

Functions

template<uint32_t n_out, uint32_t K_contract, uint32_t blockDim> constexpr bool suggested_use_reduced()#

Should a contraction-parallel *_reduced op be preferred over the serial one?

Codegen / launch-time picker seeded by the measured crossover sweep (bench/RESULTS.md, reduced section). On sm_120 the measured answer is NO everywhere: the quiet-GPU resweep of 2026-07-08 found 0 of 48 configurations where *_reduced beats serial by more than the ±5% tie margin — the family pays a warp-shuffle latency per output and idles most lanes at short contractions, and even the former long-contraction corner (n_out <= blockDim/32 && K_contract >= 32) collapsed into the noise band. So this returns false unconditionally; it keeps its original signature as the seam where a retune on different hardware (e.g. Jetson Orin, whose shuffle/FMA balance differs) can reinstate a data-derived corner without touching call sites. The *_reduced ops stay in the library for expressiveness and fusion, not speed. Not a device function (the choice is a launch/codegen decision); constexpr so the if constexpr at call sites folds to the serial path with zero cost.

Template Parameters:

n_out – Output element count (e.g. M*K for gemm, M for gemv).
K_contract – Length of the contracted dimension.
blockDim – Launch thread count.

Returns:

true to use the *_reduced variant, false to use the serial op.

template<typename T, uint32_t M, uint32_t N, uint32_t K, bool TRANSPOSE_A, bool TRANSPOSE_B, bool ROW_MAJOR_C, bool HAS_BETA> void gemm_reduced_impl_ct(uint32_t rank, uint32_t size, T alpha, T *A, T *B, T beta, T *C)#

template<typename T, uint32_t M, uint32_t N, uint32_t K, bool TRANSPOSE_A = false, bool TRANSPOSE_B = false, bool ROW_MAJOR_C = false, bool TRAILING_SYNC = true> void gemm_reduced(T alpha, T *A, T *B, T beta, T *C)#

Contraction-parallel GEMM: C = alpha * A * op(B) + beta * C.

Single-warp contraction-parallel GEMM: C = alpha * A * op(B) + beta * C.

Same math and layout as the compile-time glass::gemm, but parallelizes the length-N contraction: one warp owns each output element and its 32 lanes split the inner sum (combined with a single warp-shuffle reduce) instead of one thread summing serially. A utilization win when the output count is smaller than the block — see :doc:../../user_guide/concepts/contraction_parallel and glass::suggested_use_reduced. Total MAC work is unchanged.

Thread-count invariant: bit-identical at any block size (a trailing partial warp idles; below 32 threads a register path reproduces the same rounding).

One 32-lane warp computes the full product, parallelizing the contraction across its lanes (warp-shuffle reduce per output). The warp-per-problem analogue of the block glass::gemm_reduced; the caller must run a full 32-lane warp. C must not alias A/B.

Template Parameters:

T – Scalar type.
M, N, K – C is M×N, contraction K. op(A) is M×K, op(B) is K×N (see gemm.cuh).
TRANSPOSE_A – If true, A is K×M and op(A)=Aᵀ.
TRANSPOSE_B – If true, B is N×K and op(B)=Bᵀ.
ROW_MAJOR_C – Output storage order (false = column-major / Fortran).
TRAILING_SYNC – Emit a trailing __syncthreads() (default true) so callers can read C safely.
T – Scalar type.
M, N, K – C is M×N, contraction K. op(A) is M×K, op(B) is K×N (see gemm.cuh).
TRANSPOSE_A – If true, A is K×M and op(A)=Aᵀ.
TRANSPOSE_B – If true, B is N×K and op(B)=Bᵀ.
ROW_MAJOR_C – Output storage order (false = column-major / Fortran).
TRAILING_SYNC – Emit a trailing __syncwarp() (default true) so lanes can read C safely.

Parameters:

alpha – Scalar multiplier on the product.
A, B – Input matrices.
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out result matrix.
alpha – Scalar multiplier on the product.
A, B – Input matrices.
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out result matrix.

template<typename T, uint32_t M, uint32_t N, uint32_t K, bool TRANSPOSE_A = false, bool TRANSPOSE_B = false, bool ROW_MAJOR_C = false, bool TRAILING_SYNC = true> void gemm_reduced(T alpha, T *A, T *B, T *C)#

Contraction-parallel GEMM with implicit beta = 0: C = alpha * A * op(B).

Single-warp contraction-parallel GEMM with implicit beta = 0: C = alpha * A * op(B).

Overwrites C (the existing C is not read), avoiding the beta * C term. Otherwise identical to the beta overload above.

Overwrites C (the existing C is not read). Otherwise identical to the beta overload above; the caller must run a full 32-lane warp.

Template Parameters:

T – Scalar type.
M, N, K – C is M×N, contraction K. op(A) is M×K, op(B) is K×N (see gemm.cuh).
TRANSPOSE_A – If true, A is K×M and op(A)=Aᵀ.
TRANSPOSE_B – If true, B is N×K and op(B)=Bᵀ.
ROW_MAJOR_C – Output storage order (false = column-major / Fortran).
TRAILING_SYNC – Emit a trailing __syncthreads() (default true).
T – Scalar type.
M, N, K – C is M×N, contraction K. op(A) is M×K, op(B) is K×N (see gemm.cuh).
TRANSPOSE_A – If true, A is K×M and op(A)=Aᵀ.
TRANSPOSE_B – If true, B is N×K and op(B)=Bᵀ.
ROW_MAJOR_C – Output storage order (false = column-major / Fortran).
TRAILING_SYNC – Emit a trailing __syncwarp() (default true).

Parameters:

alpha – Scalar multiplier on the product.
A, B – Input matrices.
C – Output result matrix (overwritten).
alpha – Scalar multiplier on the product.
A, B – Input matrices.
C – Output result matrix (overwritten).

namespace warp

Tensor ⊗ vector contractions (`tensor_vec_contract` / `vec_tensor_vec`)#

Enums

enum class TensorAxis#

Values:

enumerator K#

enumerator A#

enumerator B#

Functions

template<typename T, uint32_t K, uint32_t A, uint32_t B, TensorAxis CONTRACT = TensorAxis::K, bool SYMMETRIC = false, bool ACCUMULATE = true, bool TIN_ROW_MAJOR = false, bool TRAILING_SYNC = true> void tensor_vec_contract(const T *Tns, const T *v, T *Mout)#

Tensor ⊗ vector contraction: Mout (+)= Σ_c v[c] · T[..c..].

Single-warp tensor ⊗ vector contraction: Mout (+)= Σ_c v[c] · T[..c..].

Contracts a (K, A, B) tensor against a vector along one axis (default the leading K axis), producing a matrix. With CONTRACT = TensorAxis::K: Mout[a + b*A] (+)= Σ_k v[k] · Tns[k,a,b] — the second-order Hessian-fold Hxx += Σ_i Vx[i]·fxx[i]. Contracting A or B instead gives a (K,B) or (K,A) result. One warp owns each output and its lanes split the contracted axis (warp-shuffle reduce); thread-count invariant at any block size.

Warp-per-problem analogue of glass::tensor_vec_contract; one full 32-lane warp performs the whole contraction. See the block version for semantics.

Template Parameters:

T – Scalar type.
K, A, B – Tensor dimensions (slabs K, each A x B).
CONTRACT – Axis contracted away (default K). Output is the other two axes, column-major.
SYMMETRIC – When the K-slabs are symmetric in (a,b): compute the lower triangle and mirror (requires CONTRACT==K, A==B).
ACCUMULATE – Add into Mout (true, default) vs overwrite (false).
TIN_ROW_MAJOR – Each tensor slab is row-major a*B+b (true) vs column-major a+b*A (false, default).
TRAILING_SYNC – Emit a trailing __syncthreads() (default true).
T, K, A, B, CONTRACT, SYMMETRIC, ACCUMULATE, TIN_ROW_MAJOR – See glass::tensor_vec_contract.
TRAILING_SYNC – Emit a trailing __syncwarp() (default true).

Parameters:

Tns – Input tensor (K slabs of A x B).
v – Contraction vector (length = contracted-axis size).
Mout – In/out result matrix (column-major; read only when ACCUMULATE).
Tns, v, Mout – See glass::tensor_vec_contract.

template<typename T, uint32_t K, uint32_t A, uint32_t B, bool ACCUMULATE = false, bool TIN_ROW_MAJOR = false, bool TRAILING_SYNC = true> void vec_tensor_vec(const T *Tns, const T *u, const T *w, T *s)#

Vector–tensor–vector triple product: s[k] (+)= u^T · T_k · w.

Single-warp vector–tensor–vector triple product: s[k] (+)= u^T · T_k · w.

For each slab k of a (K, A, B) tensor, forms the bilinear form s[k] = Σ_{a,b} u[a] · Tns[k,a,b] · w[b] (second-order curvature along each mode). One warp owns each s[k] and its lanes split the flattened (a,b) contraction; thread-count invariant at any block size.

Warp-per-problem analogue of glass::vec_tensor_vec.

Template Parameters:

T – Scalar type.
K, A, B – Tensor dimensions (slabs K, each A x B).
ACCUMULATE – Add into s (false, default = overwrite).
TIN_ROW_MAJOR – Each tensor slab row-major a*B+b (true) vs column-major a+b*A (false, default).
TRAILING_SYNC – Emit a trailing __syncthreads() (default true).
T, K, A, B, ACCUMULATE, TIN_ROW_MAJOR – See glass::vec_tensor_vec.
TRAILING_SYNC – Emit a trailing __syncwarp() (default true).

Parameters:

Tns – Input tensor (K slabs of A x B).
u – Left vector (length A).
w – Right vector (length B).
s – In/out result vector (length K; read only when ACCUMULATE).
Tns, u, w, s – See glass::vec_tensor_vec.

template<typename T, uint32_t K, uint32_t A, uint32_t B, TensorAxis CONTRACT = TensorAxis::K, bool SYMMETRIC = false, bool ACCUMULATE = true, bool TIN_ROW_MAJOR = false> void tensor_vec_contract(const T *Tns, const T *v, T *Mout)#

Single-thread tensor ⊗ vector contraction: Mout (+)= Σ_c v[c] · T[..c..].

Thread-per-problem analogue of glass::tensor_vec_contract: ONE thread performs the whole contraction serially — for packing 32 independent low-DOF Hessian-folds into a warp. No barriers, no shuffles, no threadIdx read; operands may be thread-local register arrays (the implied T[K*A*B] tensor and T[A*B]-scale output stay register-resident only at very small dims — see the thread-tier N<=7 element-count ceiling in CLAUDE.md; larger operands still compute correctly but spill to local memory). Same algorithm and operand order as the glass:: twin, agreeing to a few ULP (cross-tier bit-identity is NOT guaranteed).

Template Parameters:: T, K, A, B, CONTRACT, SYMMETRIC, ACCUMULATE, TIN_ROW_MAJOR – See glass::tensor_vec_contract.
Parameters:: Tns, v, Mout – See glass::tensor_vec_contract (Mout read only when ACCUMULATE).

template<typename T, uint32_t K, uint32_t A, uint32_t B, bool ACCUMULATE = false, bool TIN_ROW_MAJOR = false> void vec_tensor_vec(const T *Tns, const T *u, const T *w, T *s)#

Single-thread vector–tensor–vector triple product: s[k] (+)= u^T · T_k · w.

Thread-per-problem analogue of glass::vec_tensor_vec: ONE thread forms every slab’s bilinear form serially. No barriers, no shuffles, no threadIdx read; operands may be thread-local register arrays (subject to the tier’s element-count ceiling — see CLAUDE.md). Same algorithm and operand order as the glass:: twin, agreeing to a few ULP (cross-tier bit-identity is NOT guaranteed).

Template Parameters:: T, K, A, B, ACCUMULATE, TIN_ROW_MAJOR – See glass::vec_tensor_vec.
Parameters:: Tns, u, w, s – See glass::vec_tensor_vec (s read only when ACCUMULATE).

namespace tensor_detail#

Functions

template<typename T, TensorAxis C, uint32_t K, uint32_t A, uint32_t B, bool TIN_ROW_MAJOR> T tvc_term(const T *Tns, const T *v, uint32_t o0, uint32_t o1, uint32_t c)#

template<TensorAxis C, uint32_t K, uint32_t A, uint32_t B> struct tvc_dims#

Public Static Attributes

static constexpr uint32_t OUT0 = (C == TensorAxis::K) ? A : K #

static constexpr uint32_t OUT1 = (C == TensorAxis::K) ? B : (C == TensorAxis::A ? B : A)#

static constexpr uint32_t CDIM = (C == TensorAxis::K) ? K : (C == TensorAxis::A ? A : B)#

namespace warp

namespace thread

Congruence / bilinear forms (`congruence_sym` / `bilinear`)#

Functions

template<typename T, uint32_t N, uint32_t Kdim> constexpr std::size_t congruence_scratch_bytes()#

Shared-memory bytes needed by congruence_sym / bilinear scratch (M·X).

The scratch holds the intermediate N x Kdim product. Host- and device-callable.

Template Parameters:

T – Scalar type.
N – Rows of X / dimension of M.
Kdim – Columns of X (= columns of the scratch).

Returns:

Bytes for the s_scratch buffer (these are small single-block sizes).

template<typename T, uint32_t N, uint32_t Kdim, bool ACCUMULATE = false, bool TRAILING_SYNC = true> void congruence_sym(T alpha, const T *X, const T *M, T beta, T *Q, T *s_scratch)#

Symmetric congruence: Q = alpha * Xᵀ·M·X + beta * Q (Q symmetric).

Single-warp symmetric congruence: Q = alpha * Xᵀ·M·X + beta * Q.

Forms MX = M·X into s_scratch, then contracts Q = Xᵀ·MX over the shared N dimension, computing only the lower triangle and mirroring it (the result is symmetric when M is). Replaces two gemms + a temp + a transpose with one fused call. Column-major; single block; thread-count invariant.

Warp-per-problem analogue of glass::congruence_sym; one full 32-lane warp forms M·X and the Xᵀ·MX contraction. See the block version for semantics.

Template Parameters:

T – Scalar type.
N – Dimension of M (N x N) and rows of X.
Kdim – Columns of X — Q is Kdim x Kdim.
ACCUMULATE – Add into Q (true) vs overwrite (false, default).
TRAILING_SYNC – Emit a trailing __syncthreads() (default true).
T, N, Kdim, ACCUMULATE – See glass::congruence_sym.
TRAILING_SYNC – Emit a trailing __syncwarp() (default true).

Parameters:

alpha – Scalar on the product.
X – Input N x Kdim matrix (column-major).
M – Input N x N matrix (column-major; symmetric for a symmetric Q).
beta – Scalar on the existing Q (read only when ACCUMULATE).
Q – In/out Kdim x Kdim result (column-major).
s_scratch – Shared scratch of congruence_scratch_bytes<T,N,Kdim>() bytes (holds M·X).
alpha, X, M, beta, Q, s_scratch – See glass::congruence_sym.

template<typename T, uint32_t N, uint32_t P, uint32_t Qd, bool ACCUMULATE = false, bool TRAILING_SYNC = true> void bilinear(T alpha, const T *X, const T *M, const T *Y, T beta, T *R, T *s_scratch)#

General bilinear form: R = alpha * Xᵀ·M·Y + beta * R.

Single-warp general bilinear form: R = alpha * Xᵀ·M·Y + beta * R.

Like congruence_sym but with a distinct right operand Y, so the result is not symmetric and the full P x Qd matrix is computed. Forms MY = M·Y into s_scratch, then contracts R = Xᵀ·MY. Column-major; single block; invariant.

Warp-per-problem analogue of glass::bilinear.

Template Parameters:

T – Scalar type.
N – Dimension of M (N x N) and rows of X and Y.
P – Columns of X — R has P rows.
Qd – Columns of Y — R has Qd columns.
ACCUMULATE – Add into R (true) vs overwrite (false, default).
TRAILING_SYNC – Emit a trailing __syncthreads() (default true).
T, N, P, Qd, ACCUMULATE – See glass::bilinear.
TRAILING_SYNC – Emit a trailing __syncwarp() (default true).

Parameters:

alpha – Scalar on the product.
X – Input N x P matrix (column-major).
M – Input N x N matrix (column-major).
Y – Input N x Qd matrix (column-major).
beta – Scalar on the existing R (read only when ACCUMULATE).
R – In/out P x Qd result (column-major).
s_scratch – Shared scratch of congruence_scratch_bytes<T,N,Qd>() bytes (holds M·Y).
alpha, X, M, Y, beta, R, s_scratch – See glass::bilinear.

template<typename T, uint32_t P, uint32_t Q> constexpr std::size_t congruence_accum_scratch_bytes()#

Scratch size in bytes for congruence_accum s_scratch.

Holds the Q×P transpose Gᵀ plus the congruence_sym scratch (M·Gᵀ, also Q×P). Total 2*P*Q elements of T.

template<typename T, uint32_t P, uint32_t Q, bool ACCUMULATE = false, bool TRAILING_SYNC = true> void congruence_accum(T alpha, const T *G, const T *M, T beta, T *C, T *s_scratch)#

Accumulating congruence with a rectangular left factor: C = alpha*G*M*Gᵀ + beta*C.

Single-warp accumulating congruence C = alpha*G*M*Gᵀ + beta*C (G is P×Q).

The “other orientation” of congruence_sym: here the rectangular factor G is P×Q (the natural storage in e.g. GATO’s Schur assembly, where G = B and M = R⁻¹, giving B·R⁻¹·Bᵀ), M is Q×Q symmetric, and the symmetric result C is P×P. Mathematically G·M·Gᵀ = XᵀMX with X = Gᵀ, so this transposes G into scratch and defers to congruence_sym<Q,P> — inheriting its exact triangle+mirror symmetry and its honest FMA-order note. ACCUMULATE adds into C (the += GATO wants); default overwrites. Single block, column-major, thread-count invariant. NumPy: C = alpha*(G @ M @ G.T) + beta*C.

Warp-per-problem analogue of glass::congruence_accum: one 32-lane warp transposes G into scratch and defers to warp::congruence_sym<Q,P>. See the block version for semantics.

Template Parameters:

T – Scalar type.
P – Rows of G — C is P×P.
Q – Columns of G and dimension of M (Q×Q).
ACCUMULATE – Add into C (true) vs overwrite (false, default).
TRAILING_SYNC – Emit a trailing __syncthreads() (default true).
T, P, Q, ACCUMULATE – See glass::congruence_accum.
TRAILING_SYNC – Emit a trailing __syncwarp() (default true).

Parameters:

alpha – Scalar on the product.
G – Input P×Q matrix (column-major).
M – Input Q×Q matrix (column-major, symmetric).
beta – Scalar on the existing C (read only when ACCUMULATE).
C – In/out P×P symmetric result (column-major).
s_scratch – Shared scratch of congruence_accum_scratch_bytes<T,P,Q>() elements.
alpha, G, M, beta, C, s_scratch – See glass::congruence_accum.

template<typename T, uint32_t N, uint32_t Kdim, bool ACCUMULATE = false> void congruence_sym(T alpha, const T *X, const T *M, T beta, T *Q, T *scratch)#

Single-thread symmetric congruence: Q = alpha * Xᵀ·M·X + beta * Q.

Thread-per-problem analogue of glass::congruence_sym: ONE thread forms MX = M·X into scratch then contracts Q = Xᵀ·MX (lower triangle + mirror). No barriers, no shuffles, no threadIdx read; operands and scratch may be thread-local register arrays (the implied T[N*N] M stays register-resident only under the tier’s N<=7 element-count ceiling — see CLAUDE.md). Same algorithm and operand order as the glass:: twin, agreeing to a few ULP (cross-tier bit-identity is NOT guaranteed).

Template Parameters:

T, N, Kdim, ACCUMULATE – See glass::congruence_sym.

Parameters:

alpha, X, M, beta, Q – See glass::congruence_sym (Q read only when ACCUMULATE).
scratch – Workspace of congruence_scratch_bytes<T,N,Kdim>() bytes (holds M·X); a thread-local T[N*Kdim] is the intended form.

template<typename T, uint32_t N, uint32_t P, uint32_t Qd, bool ACCUMULATE = false> void bilinear(T alpha, const T *X, const T *M, const T *Y, T beta, T *R, T *scratch)#

Single-thread general bilinear form: R = alpha * Xᵀ·M·Y + beta * R.

Thread-per-problem analogue of glass::bilinear: ONE thread forms MY = M·Y into scratch then contracts the full P x Qd result. No barriers, no shuffles, no threadIdx read; operands and scratch may be thread-local register arrays (subject to the tier’s element-count ceiling — see CLAUDE.md). Same algorithm and operand order as the glass:: twin, agreeing to a few ULP (cross-tier bit-identity is NOT guaranteed).

Template Parameters:

T, N, P, Qd, ACCUMULATE – See glass::bilinear.

Parameters:

alpha, X, M, Y, beta, R – See glass::bilinear (R read only when ACCUMULATE).
scratch – Workspace of congruence_scratch_bytes<T,N,Qd>() bytes (holds M·Y); a thread-local T[N*Qd] is the intended form.

template<typename T, uint32_t P, uint32_t Q, bool ACCUMULATE = false> void congruence_accum(T alpha, const T *G, const T *M, T beta, T *C, T *scratch)#

Single-thread accumulating congruence C = alpha*G*M*Gᵀ + beta*C (G is P×Q).

Thread-per-problem analogue of glass::congruence_accum: ONE thread transposes G into scratch and defers to thread::congruence_sym<Q,P> — the same construction as the block/warp twins, minus their barrier (sequential program order makes Gᵀ visible). No barriers, no shuffles, no threadIdx read; scratch is algorithmic workspace (holds Gᵀ then M·Gᵀ), intended as a thread-local T[2*P*Q] (see the tier’s element-count ceiling in CLAUDE.md). Same algorithm and operand order as the glass:: twin, agreeing to a few ULP (cross-tier bit-identity is NOT guaranteed).

Template Parameters:

T, P, Q, ACCUMULATE – See glass::congruence_accum.

Parameters:

alpha, G, M, beta, C – See glass::congruence_accum (C read only when ACCUMULATE).
scratch – Workspace of congruence_accum_scratch_bytes<T,P,Q>() bytes; a thread-local T[2*P*Q] is the intended form.

namespace congruence_detail#

namespace warp

namespace thread

Contraction-parallel syrk (`syrk_reduced`)#

Functions

template<typename T, uint32_t ROWS, uint32_t COLS, bool TRANSPOSE = false, bool TRAILING_SYNC = true> void syrk_reduced(T alpha, const T *A, T beta, T *C)#

Contraction-parallel symmetric rank-k update: C = alpha * A·op(A) + beta * C.

Single-warp contraction-parallel SYRK: C = alpha * A·op(A) + beta * C.

Compile-time-size SYRK that parallelizes the contracted dimension across a warp’s lanes (one warp per output) and computes the lower triangle, mirroring it to produce a full symmetric C — the SYRK analogue of glass::gemm_reduced. Column-major. Thread-count invariant at any block size.

Warp-per-problem analogue of glass::syrk_reduced; one full 32-lane warp.

Template Parameters:

T – Scalar type.
ROWS, COLS – A is ROWS x COLS (column-major).
TRANSPOSE – If true, C = AᵀA (COLS x COLS); else C = AAᵀ (ROWS x ROWS).
TRAILING_SYNC – Emit a trailing __syncthreads() (default true).
T, ROWS, COLS, TRANSPOSE – See glass::syrk_reduced.
TRAILING_SYNC – Emit a trailing __syncwarp() (default true).

Parameters:

alpha – Scalar on the product.
A – Input matrix (ROWS x COLS, column-major).
beta – Scalar on the existing C (read only when beta != 0).
C – In/out symmetric result (full storage; OUT x OUT, OUT = TRANSPOSE?COLS:ROWS).
alpha, A, beta, C – See glass::syrk_reduced.

template<typename T, uint32_t ROWS, uint32_t COLS, bool TRANSPOSE = false, bool TRAILING_SYNC = true> void syrk_reduced(T alpha, const T *A, T *C)#

Contraction-parallel SYRK with implicit beta = 0: C = alpha * A·op(A).

Single-warp contraction-parallel SYRK, implicit beta = 0: C = alpha * A·op(A).

Overwrites C (not read). Otherwise identical to the beta overload.

Template Parameters:

T, ROWS, COLS, TRANSPOSE, TRAILING_SYNC – See the beta overload.
T, ROWS, COLS, TRANSPOSE, TRAILING_SYNC – See the beta overload.

Parameters:

alpha, A – See the beta overload.
C – Output (overwritten; full symmetric).
alpha, A, C – See the beta overload.

namespace syrk_reduced_detail#

namespace warp

Riccati feedback gain (`riccati_gain`)#

Functions

template<typename T, uint32_t NX, uint32_t NU> constexpr std::size_t riccati_scratch_bytes()#

Scratch size in bytes for riccati_gain s_scratch.

Holds the NU×NU control-Hessian S = R + BᵀPB plus the larger of the two congruence/bilinear products (P·B is NX×NU, P·A is NX×NX).

Template Parameters:

T – Element type.
NX – State dimension.
NU – Control dimension.

Returns:

Bytes to allocate for riccati_gain’s s_scratch.

template<typename T, uint32_t NX, uint32_t NU, bool REGULARIZE = false, bool TRAILING_SYNC = true> void riccati_gain(const T *P, const T *A, const T *B, const T *R, T *Kgain, T *s_scratch, T rho = T(0), int *s_fail = nullptr)#

LQR/iLQR feedback gain: K = (R + BᵀPB)⁻¹ (BᵀPA).

Single-warp LQR/iLQR feedback gain K = (R + BᵀPB)⁻¹ (BᵀPA).

Forms the control Hessian S = R + BᵀPB (symmetric congruence), the coupling G = BᵀPA (bilinear), then solves S·K = G for the NU×NX gain by Cholesky (multi-RHS). With REGULARIZE, shifts S by rho·I before factoring (and always reports a non-PD S via s_fail) so an iLQR caller can escalate rho and retry. Single block, column-major; thread-count invariant within the surface. On return Kgain holds K (the inputs P,A,B,R are unchanged).

Warp-per-knot parity with the block glass::riccati_gain: one 32-lane warp forms S = R + BᵀPB (warp::congruence_sym), G = BᵀPA (warp::bilinear), then solves S·K = G for the NU×NX gain with the checked, optionally regularized warp::posv (NRHS=NX). Every sub-op is __syncwarp-scoped, so independent warps may run distinct knots of a batched backward pass concurrently in one block. On return Kgain holds K; P,A,B,R unchanged.

Template Parameters:

T – Scalar type (prefer double for ill-conditioned S).
NX – State dimension (P is NX×NX, A is NX×NX, B is NX×NU).
NU – Control dimension (R is NU×NU, K is NU×NX). Assumes NX >= NU.
REGULARIZE – If true, add rho·I to S before the solve (default false).
TRAILING_SYNC – Emit a trailing __syncthreads() (default true).
T – Scalar type (prefer double for ill-conditioned S).
NX – State dimension (P,A are NX×NX, B is NX×NU).
NU – Control dimension (R is NU×NU, K is NU×NX). Assumes NX >= NU.
REGULARIZE – If true, add rho·I to S before the solve (default false).
TRAILING_SYNC – Emit a trailing __syncwarp() (default true).

Parameters:

P – Cost-to-go Hessian (NX×NX, symmetric, column-major).
A – State Jacobian (NX×NX, column-major).
B – Control Jacobian (NX×NU, column-major).
R – Control cost (NU×NU, SPD, column-major).
Kgain – Out gain K (NU×NX, column-major).
s_scratch – Shared scratch of riccati_scratch_bytes<T,NX,NU>() bytes.
rho – Diagonal shift on S when REGULARIZE (ignored otherwise).
s_fail – Optional flag: set to 1 if S (after the shift) is not PD, else 0.
P, A, B, R – Inputs (column-major; see the block overload).
Kgain – Out gain K (NU×NX, column-major).
s_scratch – Shared scratch of riccati_scratch_bytes<T,NX,NU>() bytes (per warp).
rho – Diagonal shift on S when REGULARIZE (ignored otherwise).
s_fail – Optional flag: set to 1 if S (after the shift) is not PD, else 0.

template<typename T, uint32_t NX, uint32_t NU, bool REGULARIZE = false> void riccati_gain(const T *P, const T *A, const T *B, const T *R, T *Kgain, T *scratch, T rho = T(0), int *s_fail = nullptr)#

Single-thread LQR/iLQR feedback gain K = (R + BᵀPB)⁻¹ (BᵀPA).

Thread-per-knot analogue of glass::riccati_gain: ONE thread forms S = R + BᵀPB (thread::congruence_sym), G = BᵀPA (thread::bilinear), then solves S·K = G — composing the tier’s own pieces exactly as the block/warp twins compose theirs. Because the multi-RHS posv bodies are block/warp-scoped (BlockBarrier / warp trsm) and thread::posv is single-RHS, the solve leg is spelled out from the tier’s existing primitives with the identical algorithm and order the flagged posv runs: optional rho·I shift → checked thread::potrf → per-column forward/back substitution (thread::potrs, NX columns). No barriers, no shuffles, no threadIdx read, so it is safe in ragged-tail thread-per-problem launches; operands and scratch may be thread-local register arrays (the implied T[NX*NX] P stays register-resident only under the tier’s N<=7 element-count ceiling — see CLAUDE.md; larger dims still compute correctly but spill to local memory). Same algorithm and operand order as the glass:: twin, agreeing to a few ULP (cross-tier bit-identity is NOT guaranteed). On return Kgain holds K; P,A,B,R are unchanged.

scratch is ALGORITHMIC workspace (it holds S and the congruence / bilinear products the contraction re-reads), not cross-lane staging — the thread tier keeps the caller-provided pointer, intended as a thread-local array of riccati_scratch_bytes<T,NX,NU>() bytes. No TRAILING_SYNC parameter, matching the tier’s precedent.

Template Parameters:

T – Scalar type (prefer double for ill-conditioned S).
NX – State dimension (P,A are NX×NX, B is NX×NU).
NU – Control dimension (R is NU×NU, K is NU×NX). Assumes NX >= NU.
REGULARIZE – If true, add rho·I to S before the solve (default false).

Parameters:

P, A, B, R – Inputs (column-major; see the block overload).
Kgain – Out gain K (NU×NX, column-major).
scratch – Workspace of riccati_scratch_bytes<T,NX,NU>() bytes (per thread).
rho – Diagonal shift on S when REGULARIZE (ignored otherwise).
s_fail – Optional flag: set to 1 if S (after the shift) is not PD, else 0.

namespace warp

namespace thread

Strided gemm#

Functions

template<typename T, uint32_t M, uint32_t N, uint32_t K, uint32_t A_RS = M, uint32_t B_RS = K> void gemm_strided(T alpha, const T *A, const T *B, T beta, T *C)#

Strided compile-time GEMM: C = alpha * A * B + beta * C with custom leading dims.

Column-major GEMM (standard convention: C is M×N, contraction K) where A and B carry explicit leading dimensions (column strides): A[m][k] = A[m + k*A_RS], B[k][n] = B[k + n*B_RS]. Output C is standard column-major with LDC = M. When A_RS == M and B_RS == K this is identical to glass::gemm<T,M,N,K>. Single-block, flat-element parallelism; the inner K-loop is fully unrolled. NumPy: C = alpha * A @ B + beta * C on the strided sub-views.

Template Parameters:

T – Scalar type.
M, N, K – Compile-time dimensions: A is M×K, B is K×N, C is M×N (contraction K).
A_RS – Column stride (leading dimension) of A (default M).
B_RS – Column stride (leading dimension) of B (default K).

Parameters:

alpha – Scalar multiplier on the product.
A, B – Input matrices (column-major, strided).
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out result matrix (column-major, LDC = M).

template<typename T, uint32_t M, uint32_t N, uint32_t K, uint32_t A_RS = M, uint32_t B_RS = K> void gemm_strided(T alpha, const T *A, const T *B, T *C)#

Strided compile-time GEMM with implicit beta = 0: C = alpha * A * B.

Same as the beta overload but overwrites C (the existing C is not read). Column-major with explicit leading dims A_RS / B_RS; LDC = M. NumPy: C = alpha * A @ B.

Template Parameters:

T – Scalar type.
M, N, K – Compile-time dimensions: A is M×K, B is K×N, C is M×N (contraction K).
A_RS – Column stride (leading dimension) of A (default M).
B_RS – Column stride (leading dimension) of B (default K).

Parameters:

alpha – Scalar multiplier on the product.
A, B – Input matrices (column-major, strided).
C – Output result matrix (overwritten; column-major, LDC = M).

Indexed / batched gemm#

Functions

template<typename T, uint32_t DIM = 4, bool TRANSPOSE_A = false, bool TRANSPOSE_B = false, bool ATOMIC_C = false, typename IDX_T = int> void gemm_batched_indexed(uint32_t pairs, const IDX_T *a_idx, const IDX_T *b_idx, const IDX_T *c_idx, const T *A_base, const T *B_base, T *C_base)#

Indexed/gather batched square GEMM: C[c_idx[p]] = op(A[a_idx[p]]) * op(B[b_idx[p]]).

For each pair p in [0, pairs) multiplies two DIM x DIM column-major matrices selected by index into flat base buffers (a_idx[p] is a MATRIX index, so the matrix lives at offset a_idx[p] * DIM * DIM), computing all pairs concurrently in a single block. This is the indexed/gather analogue of gemm_strided, useful for assembling many independent small (e.g. 4x4 SE(3)) products from index lists. No alpha/beta — a pure overwrite (C = op(A) * op(B)) unless ATOMIC_C is set.

Layout flags read the factors transposed in place (matrices stay square DIM x DIM, so the output is always DIM x DIM): TRANSPOSE_A gives A_p^T * B_p, TRANSPOSE_B gives A_p * B_p^T, and both give A_p^T * B_p^T. Without ATOMIC_C, distinct pairs MUST target distinct c_idx slots (each output written once); a_idx / b_idx may alias freely.

Template Parameters:

T – Scalar type.
DIM – Compile-time matrix dimension (square; inner loop fully unrolled).
TRANSPOSE_A – If true, the left factor is read transposed (C_p = A_p^T * ...).
TRANSPOSE_B – If true, the right factor is read transposed (... * B_p^T).
ATOMIC_C – If true, accumulate via atomicAdd (C[c_idx[p]] += ...), allowing several pairs to share a c_idx slot; the caller must PRE-ZERO (or pre-load) the touched C slots, and no beta is applied.
IDX_T – Index type of the *_idx arrays.

Parameters:

pairs – Number of independent GEMMs.
a_idx – Per-pair matrix slot of the left factor in A_base (length pairs).
b_idx – Per-pair matrix slot of the right factor in B_base (length pairs).
c_idx – Per-pair matrix slot of the destination in C_base (length pairs).
A_base – Flat array of DIM x DIM left-factor matrices.
B_base – Flat array of DIM x DIM right-factor matrices.
C_base – Flat array of DIM x DIM destination matrices (written, or accumulated if ATOMIC_C).

Matrix inverse#

Functions

template<typename T, bool TRAILING_SYNC = true> void inv(uint32_t dimA, T *A, T *s_scratch)#

template<typename T, uint32_t N, bool TRAILING_SYNC = true> void inv(T *A, T *s_scratch)#

Compile-time-size in-place matrix inverse (augmented [A | I] Gauss-Jordan).

Same as the runtime inv but with the dimension as a template parameter. NumPy equivalent: Ainv = np.linalg.inv(A).

Template Parameters:

T – Scalar type.
N – Matrix dimension (A is N x N).

Parameters:

A – In/out augmented [A | I] buffer (column-major, N x 2*N); on return its right half holds A^-1.
s_scratch – Shared scratch of (2*N + 1) * sizeof(T) bytes.

template<typename T> constexpr std::size_t inv_scratch_bytes(uint32_t dimA)#

Scratch size in bytes for inv (augmented [A | I]).

The unpivoted Gauss-Jordan path saves the active pivot column + row window plus one slot: 2*dimA + 1 elements of T. Allocate inv_scratch_bytes<T>(dimA) for the s_scratch argument.

Template Parameters:: T – Scalar type.
Parameters:: dimA – Matrix dimension (A is dimA x dimA).
Returns:: Bytes to allocate for inv’s s_scratch.

template<typename T> constexpr std::size_t inv_pivoted_scratch_bytes(uint32_t dimA)#

Scratch size in bytes for inv_pivoted.

Row pivoting permutes the already-built inverse columns, so (unlike the unpivoted path) the elimination cannot use the reduced active-column window — it must save and update the full 2*dimA-wide pivot row. Layout: dimA slots for the pivot column, 2*dimA for the full pivot row, and one trailing slot to broadcast the chosen pivot-row index from the argmax.

Total = 3*dimA + 1 elements of T.

Template Parameters:: T – Scalar type.
Parameters:: dimA – Matrix dimension (A is dimA x dimA).
Returns:: Bytes to allocate for inv_pivoted’s s_scratch.

template<typename T, bool TRAILING_SYNC = true> void inv_pivoted(uint32_t dimA, T *A, T *s_scratch)#

template<typename T, uint32_t N, bool TRAILING_SYNC = true> void inv_pivoted(T *A, T *s_scratch)#

Compile-time-size ROBUST (partial-pivoting) matrix inverse ([A | I]).

Same as the runtime inv_pivoted but with the dimension as a template parameter; partial-pivoting Gauss-Jordan, tolerant of small leading pivots that the plain inv mishandles. NumPy equivalent: Ainv = np.linalg.inv(A).

Template Parameters:

T – Scalar type.
N – Matrix dimension (A is N x N).

Parameters:

A – In/out augmented [A | I] buffer (column-major, N x 2*N); on return its right half holds A^-1.
s_scratch – Shared scratch of (3*N + 1) * sizeof(T) bytes (= inv_pivoted_scratch_bytes<T>(N)).

template<typename T> constexpr std::size_t inv_fused_scratch_bytes(uint32_t K, const uint32_t *dims)#

Scratch size in bytes for the K-way fused inv (Σ_m (2*dims[m]+1) elements).

Template Parameters:

T – Scalar type.

Parameters:

K – Number of matrices.
dims – Per-matrix dimensions.

Returns:

Bytes to allocate for the fused inv’s s_scratch.

template<typename T, bool TRAILING_SYNC = true>
void inv(uint32_t K, const uint32_t *dims, uint32_t MAX_DIM, T **mats, T *s_scratch)#

Fused in-place inverse of K independent matrices (augmented [V | I]).

Inverts K matrices simultaneously in one block by interleaving their Gauss-Jordan sweeps over a single shared MAX_DIM = max(dims) pivot loop: matrix m participates while pivRC < dims[m] and sits idle thereafter. Every matrix keeps the same augmented [V | I] convention as the single-matrix inv — buffer mats[m] is column-major dims[m] x (2*dims[m]) and on return its right half holds inv(mats[m]). Fewer barriers than K separate calls (one save→update barrier pair per pivot step, shared by all K matrices). Used by GATO’s Schur kernel (Q_k, Q_kp1, R_k → K=3).

Scratch layout: matrix m owns the contiguous span [Σ_{j<m}(2*dims[j]+1), Σ_{j<=m}(2*dims[j]+1)) of s_scratch (each matrix needs 2*dims[m]+1 slots: dims[m] for its pivot column, dims[m]+1 for its pivot row plus the augmented column). The per-matrix base offset is the prefix sum Σ_{j<m}(2*dims[j]+1), recomputed locally per thread by scanning dims[] (no shared write, so race-free). Total scratch = Σ_m (2*dims[m]+1) elements.

NumPy equivalent (per matrix m): inv(m) = np.linalg.inv(mats[m]).

Template Parameters:

T – Scalar type.

Parameters:

K – Number of matrices.
dims – Per-matrix dimensions (dims[m] for matrix m).
MAX_DIM – max(dims[0..K-1]) — the shared pivot-loop length (precondition).
mats – Array of K in/out augmented [V | I] buffers (column-major, dims[m] x 2*dims[m]); on return each right half holds its inverse.
s_scratch – Shared scratch of inv_fused_scratch_bytes<T>(K, dims) bytes (= (Σ_m (2*dims[m]+1)) * sizeof(T)).

template<typename T, bool TRAILING_SYNC = true> void inv(uint32_t dimA, uint32_t dimB, uint32_t MAX_DIM, T *A, T *B, T *s_scratch)#

Fused in-place inverse of TWO independent matrices (augmented [V | I]).

Thin wrapper over the K-way inv (K=2). Inverts A (dimA x dimA) and B (dimB x dimB) simultaneously in one block; same augmented [V | I] convention and output as the single-matrix inv. NumPy: Ainv, Binv = inv(A), inv(B).

Template Parameters:

T – Scalar type.

Parameters:

dimA, dimB – Matrix dimensions.
MAX_DIM – max(dimA, dimB) — the shared pivot-loop length.
A, B – In/out augmented [V | I] buffers (column-major, dim x 2*dim).
s_scratch – Shared scratch of (2*dimA + 2*dimB + 2) * sizeof(T) bytes.

template<typename T, bool TRAILING_SYNC = true> void inv(uint32_t dimA, uint32_t dimB, uint32_t dimC, uint32_t MAX_DIM, T *A, T *B, T *C, T *s_scratch)#

Fused in-place inverse of THREE independent matrices (augmented [V | I]).

Thin wrapper over the K-way inv (K=3). Inverts A,B,C simultaneously in one block; same augmented [V | I] convention and output as the single-matrix inv. Used by GATO’s Schur kernel (Q_k, Q_kp1, R_k). NumPy: invert each independently.

Template Parameters:

T – Scalar type.

Parameters:

dimA, dimB, dimC – Matrix dimensions.
MAX_DIM – max(dimA, dimB, dimC) — the shared pivot-loop length.
A, B, C – In/out augmented [V | I] buffers (column-major, dim x 2*dim).
s_scratch – Shared scratch of (2*dimA + 2*dimB + 2*dimC + 3) * sizeof(T) bytes.

template<typename T> constexpr std::size_t inv_dense_scratch_bytes(uint32_t dimA)#

Scratch size in bytes for inv_dense.

The dual-buffer dense path saves one 3*dimA-element pivot working set. Allocate inv_dense_scratch_bytes<T>(dimA) for the s_scratch argument.

Template Parameters:: T – Scalar type.
Parameters:: dimA – Matrix dimension (A is dimA x dimA).
Returns:: Bytes to allocate for inv_dense’s s_scratch.

template<typename T, bool TRAILING_SYNC = true> void inv_dense(uint32_t dimA, T *A, T *Ainv, T *s_scratch)#

template<typename T, uint32_t N, bool TRAILING_SYNC = true> void inv_dense(T *A, T *Ainv, T *s_scratch)#

Compile-time-size dense in-place matrix inverse (dual-update Gauss-Jordan).

Same as the runtime inv_dense but with the dimension as a template parameter; on return both A and Ainv hold A^{-1}. NumPy equivalent: Ainv = np.linalg.inv(A).

Template Parameters:

T – Scalar type.
N – Matrix dimension (A is N x N).

Parameters:

A – In/out column-major N x N matrix; on return holds A^{-1}.
Ainv – Workspace column-major N x N; on return also holds A^{-1}.
s_scratch – Shared scratch of 3 * N * sizeof(T) bytes.

template<typename T, uint32_t N> void inv(T *A, T *s_scratch)#

Single-warp in-place matrix inverse (unpivoted Gauss-Jordan, augmented [A | I]), compile-time size.

Single-thread in-place matrix inverse (unpivoted Gauss-Jordan, augmented [A | I]), compile-time size.

One 32-lane warp reduces a column-major augmented N x 2*N [A | I] buffer so that on return columns N..2*N-1 hold A^-1 — the same layout, phases, and arithmetic as the block glass::inv, scoped to a warp for warp-per-problem kernels (e.g. packing many small Schur-block inversions from GATO/MPCGPU into one block, one warp each). Per pivot: a lane-strided SAVE of the pivot column + active pivot-row window into s_scratch, __syncwarp(), then a lane-strided Gauss-Jordan cell UPDATE over the N x (N+1) active window, __syncwarp() — mirroring inv_impl’s two-phase structure exactly. The pivot reciprocal is computed redundantly by every lane from the SAVED shared value (1 / s_scratch[pivRC], the same bits every lane) — deterministic, and never a lane-0-register broadcast, so there is nothing for the __restrict__ stale-shared-reread miscompile (guide §1g) to bite. No __syncthreads. Unpivoted: like the block inv, it divides by the leading pivots as-is (no row exchange) — use the block inv_pivoted when robustness to small/zero leading pivots is needed. Fused K-way and pivoted warp forms are deliberately not provided (future work).

NumPy equivalent: Ainv = np.linalg.inv(A).

ONE thread reduces a column-major augmented N x 2*N [A | I] buffer so that on return columns N..2*N-1 hold A^-1 — the sequential Gauss-Jordan sweep, for thread-per-problem solvers that pack 32 independent low-DOF problems into a warp. No shared scratch requirement, no barriers, no threadIdx read.

Delegates to the same inv_impl body the block surface uses, via ThreadBarrier (rank=0, size=1, no-op sync) — the same algorithm and operand order as glass::inv<T, N> on one thread, agreeing to a few ULP (FMA-contraction jitter; bit-identity across the two instantiations is NOT guaranteed — see test/test_thread.py).

s_scratch keeps the caller-provided-pointer signature for surface uniformity with glass::inv / warp::inv, but on this tier the intended use is a THREAD-LOCAL array — T scratch[2*N + 1]; declared in the caller — so both the operand and the scratch can stay register-resident (per-thread shared-memory scratch would defeat the tier’s packing).

Unpivoted: divides by the leading pivots as-is (no row exchange) — the data-dependent inv_pivoted is deliberately excluded from this tier (its argmax branches diverge across a warp of independent problems; use the block inv_pivoted when robustness to small/zero leading pivots is needed). NumPy equivalent: Ainv = np.linalg.inv(A).

Template Parameters:

T – Scalar type.
N – Matrix dimension (A is N x N).
T – Scalar type.
N – Matrix dimension (A is N x N). NOTE the augmented operand is T[2*N*N] — twice the T[N*N] the measured N<=7 register-residency ceiling refers to (see the thread-tier constraints in CLAUDE.md), so expect the local-memory cliff at a SMALLER N than for the factor/solve ops; larger sizes still compute correctly, they just forfeit the tier’s premise.

Parameters:

A – In/out augmented [A | I] buffer (column-major, N x 2*N); on return its right half holds A^-1.
s_scratch – Scratch of 2*N + 1 elements of T (= inv_scratch_bytes<T>(N) bytes), shared or global. Each warp needs its OWN 2*N + 1 span — when packing W warps into a block, give warp w s_scratch + w*(2*N+1).
A – In/out augmented [A | I] buffer (column-major, N x 2*N); on return its right half holds A^-1.
s_scratch – Scratch of 2*N + 1 elements of T (= inv_scratch_bytes<T>(N) bytes); a thread-local T scratch[2*N + 1] is the intended use here.

namespace warp

namespace thread

Cholesky#

Functions

template<typename T, bool CHECK = false, bool TRAILING_SYNC = true> void potrf(uint32_t n, T *s_A, int *s_fail = nullptr)#

template<typename T, bool TRAILING_SYNC = true>
void potrf(uint32_t K, const uint32_t *dims, uint32_t MAX_DIM, T **mats)#

Fused in-place Cholesky factorization of K independent SPD matrices (lower).

Factors K SPD matrices simultaneously in one block by interleaving their column sweeps over a single shared MAX_DIM = max(dims) row loop: matrix m participates while row < dims[m] and sits idle thereafter. Each matrix keeps the same column-major in-place A = L*L^T convention as the single-matrix potrf — on return the lower triangle of mats[m] holds its factor L (the upper triangle keeps its input values). Same two-barriers-per step structure as the single-matrix path; no shared scratch required.

The K diagonals of a given step are distributed across threads (for (m = rank; m < K; m += size), more parallel than rank-0 alone); then the trailing sub-diagonal column entries of each active matrix are updated in parallel. NumPy equivalent (per matrix m): cholesky(m) = np.linalg.cholesky(mats[m]) (lower).

Template Parameters:

T – Scalar type.

Parameters:

K – Number of matrices.
dims – Per-matrix dimensions (dims[m] for matrix m).
MAX_DIM – max(dims[0..K-1]) — the shared row-loop length (precondition).
mats – Array of K in/out column-major SPD buffers (dims[m] x dims[m]); on return each lower triangle holds its Cholesky factor L.

template<typename T, bool TRAILING_SYNC = true> void potrf(uint32_t dimA, uint32_t dimB, uint32_t MAX_DIM, T *A, T *B)#

Fused in-place Cholesky factorization of TWO SPD matrices (lower).

Thin wrapper over the K-way potrf (K=2). Same column-major in-place A = L*L^T convention and output as the single-matrix path. NumPy: La, Lb = cholesky(A), cholesky(B) (lower).

Template Parameters:

T – Scalar type.

Parameters:

dimA, dimB – Matrix dimensions.
MAX_DIM – max(dimA, dimB) — the shared row-loop length.
A, B – In/out column-major SPD buffers (dim x dim); lower triangles hold L.

template<typename T, bool TRAILING_SYNC = true> void potrf(uint32_t dimA, uint32_t dimB, uint32_t dimC, uint32_t MAX_DIM, T *A, T *B, T *C)#

Fused in-place Cholesky factorization of THREE SPD matrices (lower).

Thin wrapper over the K-way potrf (K=3). Same column-major in-place A = L*L^T convention and output as the single-matrix path. NumPy: invert-free factor each independently (lower).

Template Parameters:

T – Scalar type.

Parameters:

dimA, dimB, dimC – Matrix dimensions.
MAX_DIM – max(dimA, dimB, dimC) — the shared row-loop length.
A, B, C – In/out column-major SPD buffers (dim x dim); lower triangles hold L.

template<typename T, uint32_t N, bool CHECK = false, bool TRAILING_SYNC = true> void potrf(T *s_A, int *s_fail = nullptr)#

Compile-time-size in-place Cholesky factorization (LAPACK potrf, lower).

Same as the runtime overload but with the dimension as a template parameter, letting the compiler bake N in. Factors the SPD matrix A = L * L^T in place, writing only the lower triangle. NumPy equivalent: L = np.linalg.cholesky(A).

When CHECK is true and s_fail is non-null, reports a non-PD / NaN pivot via *s_fail (see the runtime overload). CHECK defaults false and compiles out, so the unchecked instantiation is byte-identical to the original.

Template Parameters:

T – Scalar type.
N – Matrix dimension (A is N x N).
CHECK – If true, detect a non-PD pivot and report it via s_fail (default false, compiles out).

Parameters:

s_A – In/out N x N matrix (column-major); on return its lower triangle holds L.
s_fail – Optional flag (CHECK only): set to 1 on a non-PD / NaN pivot, else 0. Ignored when null.

template<typename T, uint32_t N, bool CHECK = false> void potrf(T *s_A, int *s_fail = nullptr)#

Single-warp in-place Cholesky factorization (LAPACK potrf, lower), compile-time size.

Single-thread in-place Cholesky factorization (LAPACK potrf, lower), compile-time size.

One 32-lane warp factors the SPD matrix A = L * L^T in place, writing only the lower triangle (column-major). For warp-per-problem solvers on small systems (e.g. N≈7 normal equations). Lane 0 computes each diagonal; the remaining sub-diagonal entries of the column are filled by the warp’s lanes (stride 32), synchronized with __syncwarp. No shared scratch, no __syncthreads. A must be SPD. NumPy equivalent: L = np.linalg.cholesky(A).

When CHECK is true and s_fail is non-null, reports a non-PD / NaN pivot via *s_fail (lane 0 writes it, mirroring the block overload). CHECK defaults false and compiles out, so the unchecked instantiation is byte-identical to the original.

ONE thread factors the SPD matrix A = L * L^T in place, writing only the lower triangle (column-major) — the sequential algorithm, for thread-per-problem solvers that pack 32 independent low-DOF problems into a warp (e.g. N≈7 IK normal equations, one seed per lane). No shared scratch, no barriers, no threadIdx read; A may live in a thread-local array and stay register-resident. A must be SPD. NumPy: L = np.linalg.cholesky(A).

Delegates to the same potrf_impl body the block/warp surfaces use, via ThreadBarrier (rank=0, size=1, no-op sync) — the same algorithm and operand order as glass::potrf<T, N> on a single thread. NOT guaranteed bit-identical across the two instantiations: the no-op sync removes the optimization fences, so FMA contraction may differ by a last ULP (test/test_thread.py pins the bound).

COMPILE-TIME SIZE ONLY (no runtime-n overload): the tier’s value is an A that nvcc can keep in registers, which requires fully-unrolled, compile-time-resolvable indexing. A runtime-n form would silently spill to local memory and be strictly worse than glass::warp::potrf.

Template Parameters:

T – Scalar type (use double for stability on ill-conditioned A).
N – Matrix dimension (A is N x N).
CHECK – If true, detect a non-PD pivot and report it via s_fail (default false, compiles out).
T – Scalar type (use double for stability on ill-conditioned A).
N – Matrix dimension (A is N x N). N<=7 keeps A register-resident (measured ceiling, both dtypes — see the thread-tier constraints in CLAUDE.md); larger N still computes correctly but demotes A to local memory, forfeiting the tier’s premise.
CHECK – If true, detect a non-PD pivot and report it via s_fail (default false, compiles out).

Parameters:

s_A – In/out N x N matrix (column-major); on return its lower triangle holds L.
s_fail – Optional flag (CHECK only): set to 1 on a non-PD / NaN pivot, else 0. Ignored when null.
s_A – In/out N x N matrix (column-major); on return its lower triangle holds L.
s_fail – Optional flag (CHECK only): set to 1 on a non-PD / NaN pivot, else 0. Ignored when null.

namespace warp

namespace thread

LU with partial pivoting (getrf / getrs / gesv / laswp)#

The robustness path for general non-SPD systems: in-place pivoted LU (getrf, SciPy lu_factor), the matching solves (getrs / gesv, SciPy lu_solve / NumPy solve), and the LAPACK-style row-interchange helper laswp.

Functions

template<typename T, bool REVERSE = false, bool TRAILING_SYNC = true> void laswp(uint32_t n, uint32_t ncols, T *A, const uint32_t *piv, uint32_t k0, uint32_t k1)#

Apply LAPACK-style row interchanges to a rectangular matrix (LASWP).

Applies A[k,:] ↔ A[piv[k],:] sequentially for k in [k0, k1) to the n×ncols column-major matrix A (leading dimension n). This is how a getrf pivot vector is applied to right-hand sides before the triangular solves. REVERSE=true applies the swaps in the opposite order (k1−1 down to k0) — the INVERSE permutation, LAPACK laswp with INCX = −1. Each thread owns whole columns and applies the swap sequence serially within them (columns are independent), so the result is thread-count invariant with no interior barrier. LAPACK equivalent: ?laswp; NumPy: for k in range(k0, k1): A[[k, piv[k]], :] = A[[piv[k], k], :].

Template Parameters:

T – Scalar type.
REVERSE – Apply the swaps in reverse order (inverse permutation; default false).
TRAILING_SYNC – End on a barrier so the permuted A is block-visible (default true).

Parameters:

n – Leading dimension / number of rows of A.
ncols – Number of columns swaps are applied across.
A – In/out n×ncols matrix (column-major).
piv – Pivot indices (piv[k] = row swapped with row k; 0-based).
k0, k1 – Half-open range of swap steps to apply.

template<typename T, bool REVERSE = false, bool TRAILING_SYNC = true> void laswp(uint32_t n, T *A, const uint32_t *piv, uint32_t k0, uint32_t k1)#

Apply LAPACK-style row interchanges to a square n-column matrix (LASWP).

Square convenience form of the rectangular laswp: applies A[k,:] ↔ A[piv[k],:] sequentially for k in [k0, k1) to the n×n column-major matrix A. Swaps compose sequentially in k but are thread-strided across columns (each thread owns whole columns). NumPy: for k in range(k0, k1): A[[k, piv[k]], :] = A[[piv[k], k], :].

Template Parameters:

T – Scalar type.
REVERSE – Apply the swaps in reverse order (inverse permutation; default false).
TRAILING_SYNC – End on a barrier (default true).

Parameters:

n – Matrix dimension (A is n×n).
A – In/out n×n matrix (column-major).
piv – Pivot indices (piv[k] = row swapped with row k; 0-based).
k0, k1 – Half-open range of swap steps to apply.

template<typename T, bool REVERSE = false, bool TRAILING_SYNC = true> void laswp(const uint32_t *piv, uint32_t k0, uint32_t k1, T *x)#

Apply LAPACK-style row interchanges to a vector (LASWP, single column).

Vector form: applies x[k] ↔ x[piv[k]] sequentially for k in [k0, k1). The swaps compose (they must run in order), so one thread applies them serially; the routine ends on a barrier so the permuted x is block-visible. REVERSE=true undoes a forward application. NumPy: for k in range(k0, k1): x[[k, piv[k]]] = x[[piv[k], k]].

Template Parameters:

T – Scalar type.
REVERSE – Apply the swaps in reverse order (inverse permutation; default false).
TRAILING_SYNC – End on a barrier (default true).

Parameters:

piv – Pivot indices (piv[k] = element swapped with element k; 0-based).
k0, k1 – Half-open range of swap steps to apply.
x – In/out vector.

template<typename T, bool CHECK = false, bool TRAILING_SYNC = true> void getrf(uint32_t n, T *A, uint32_t *piv, int *s_fail = nullptr)#

In-place LU factorization with partial pivoting (LAPACK getrf).

Factors P·A = L·U, overwriting the n×n column-major A with L\\U (unit-lower L’s multipliers strictly below the diagonal, U on and above) and recording the row interchanges in piv (piv[k] = row swapped with row k at step k; LAPACK ipiv convention, 0-based, so piv[k] >= k). Partial pivoting (argmax |A[i,k]|, ties to the smallest index) makes this the robust factorization for GENERAL non-SPD matrices — it succeeds on any invertible A, including a zero leading pivot where no-pivot LU fails. SciPy equivalent: lu, piv = scipy.linalg.lu_factor(A) — the output pair is drop-in for scipy.linalg.lu_solve((lu, piv), b).

When CHECK is true and s_fail is non-null, a zero or non-finite pivot sets *s_fail = 1 and skips that column’s divide (A is singular to working precision; the factor is not usable), leaving *s_fail = 0 otherwise. CHECK defaults false and compiles out entirely (if constexpr), so the unchecked instantiation is byte-identical.

Ends on a barrier; deterministic pivot choice + barriers between every dependent phase make the output thread-count invariant.

Template Parameters:

T – Scalar type.
CHECK – If true, detect a zero/non-finite pivot and report it via s_fail (default false, compiles out).

Parameters:

n – Matrix dimension (A is n×n).
A – In/out n×n matrix (column-major); on return holds L\\U.
piv – Output pivot indices, length n (piv[k] = row swapped with row k; 0-based).
s_fail – Optional flag (CHECK only): set to 1 on a zero/non-finite pivot, else 0. Ignored when null.

template<typename T, uint32_t N, bool CHECK = false, bool TRAILING_SYNC = true> void getrf(T *A, uint32_t *piv, int *s_fail = nullptr)#

In-place LU factorization with partial pivoting, compile-time size (LAPACK getrf).

Compile-time-N overload of getrf, forwarding to the runtime form (same L\\U layout, 0-based LAPACK ipiv piv, and CHECK semantics). SciPy equivalent: lu, piv = scipy.linalg.lu_factor(A).

Template Parameters:

T – Scalar type.
N – Matrix dimension (A is N×N).
CHECK – If true, detect a zero/non-finite pivot and report it via s_fail (default false, compiles out).

Parameters:

A – In/out N×N matrix (column-major); on return holds L\\U.
piv – Output pivot indices, length N (0-based LAPACK ipiv).
s_fail – Optional flag (CHECK only): set to 1 on a zero/non-finite pivot, else 0. Ignored when null.

template<typename T, bool TRANSPOSE = false, bool TRAILING_SYNC = true> void getrs(uint32_t n, uint32_t nrhs, const T *LU, const uint32_t *piv, T *B)#

Solve op(A) X = B from a pivoted LU factorization, in place (LAPACK getrs).

Uses the (LU, piv) pair produced by getrf to solve for all nrhs columns of B (n×nrhs, column-major), overwriting B with X. TRANSPOSE=false solves A X = B: apply the row interchanges to B (laswp), forward-solve L Y = P·B (unit-lower trsm), back-solve U X = Y (upper trsm). TRANSPOSE=true solves Aᵀ X = B in the reverse order with transposed flags — Uᵀ Z = B, then Lᵀ W = Z, then the interchanges applied LAST in reverse order (the inverse permutation). Ends on a barrier. SciPy equivalent: X = scipy.linalg.lu_solve((lu, piv), B, trans=(1 if TRANSPOSE else 0)).

Template Parameters:

T – Scalar type.
TRANSPOSE – When true solve Aᵀ X = B (default false).

Parameters:

n – Dimension (LU is n×n; each column of B has length n).
nrhs – Number of right-hand sides (columns of B).
LU – Factorization from getrf (L\\U, column-major; read-only).
piv – Pivot indices from getrf (0-based LAPACK ipiv; read-only).
B – In/out right-hand sides (n×nrhs, column-major); on return holds X.

template<typename T, uint32_t N, uint32_t NRHS, bool TRANSPOSE = false, bool TRAILING_SYNC = true> void getrs(const T *LU, const uint32_t *piv, T *B)#

Solve op(A) X = B from a pivoted LU factorization, compile-time size (LAPACK getrs).

Compile-time-N/NRHS overload of getrs, forwarding to the runtime form. SciPy equivalent: X = scipy.linalg.lu_solve((lu, piv), B, trans=(1 if TRANSPOSE else 0)).

Template Parameters:

T – Scalar type.
N – Dimension (LU is N×N; each column of B has length N).
NRHS – Number of right-hand sides (columns of B).
TRANSPOSE – When true solve Aᵀ X = B (default false).

Parameters:

LU – Factorization from getrf (L\\U, column-major; read-only).
piv – Pivot indices from getrf (0-based LAPACK ipiv; read-only).
B – In/out right-hand sides (N×NRHS, column-major); on return holds X.

template<typename T, bool CHECK = false, bool TRAILING_SYNC = true> void gesv(uint32_t n, uint32_t nrhs, T *A, uint32_t *piv, T *B, int *s_fail = nullptr)#

Solve the general system A X = B via pivoted LU, in place (LAPACK gesv).

The composed general dense solve: getrf (in-place pivoted LU of A) then getrs (permute + two triangular solves on B). On return A holds L\\U, piv the interchanges, and B the solution X (n×nrhs, column-major). This is the robust path for GENERAL non-symmetric matrices — where posv/ldlt_solve require SPD/symmetry, gesv only requires invertibility (partial pivoting handles zero/small leading pivots). NumPy equivalent: X = np.linalg.solve(A, B).

When CHECK is true and s_fail is non-null, a zero/non-finite pivot in the factorization sets *s_fail = 1 (the “solution” is then meaningless — callers must test the flag), else 0. CHECK defaults false and compiles out. Ends on a barrier.

Template Parameters:

T – Scalar type.
CHECK – If true, report a zero/non-finite pivot via s_fail (default false, compiles out).

Parameters:

n – Dimension (A is n×n; each column of B has length n).
nrhs – Number of right-hand sides (columns of B).
A – In/out n×n matrix (column-major); on return holds L\\U.
piv – Output pivot indices, length n (0-based LAPACK ipiv).
B – In/out right-hand sides (n×nrhs, column-major); on return holds X.
s_fail – Optional flag (CHECK only): set to 1 on a zero/non-finite pivot, else 0. Ignored when null.

template<typename T, uint32_t N, uint32_t NRHS, bool CHECK = false, bool TRAILING_SYNC = true> void gesv(T *A, uint32_t *piv, T *B, int *s_fail = nullptr)#

Solve the general system A X = B via pivoted LU, compile-time size (LAPACK gesv).

Compile-time-N/NRHS overload of gesv, forwarding to the runtime form (same in-place L\\U/piv/X outputs and CHECK semantics). NumPy equivalent: X = np.linalg.solve(A, B).

Template Parameters:

T – Scalar type.
N – Dimension (A is N×N; each column of B has length N).
NRHS – Number of right-hand sides (columns of B).
CHECK – If true, report a zero/non-finite pivot via s_fail (default false, compiles out).

Parameters:

A – In/out N×N matrix (column-major); on return holds L\\U.
piv – Output pivot indices, length N (0-based LAPACK ipiv).
B – In/out right-hand sides (N×NRHS, column-major); on return holds X.
s_fail – Optional flag (CHECK only): set to 1 on a zero/non-finite pivot, else 0. Ignored when null.

Triangular solve (trsm)#

Functions

template<typename T, FillMode FILL = FillMode::Lower, Diag DIAG = Diag::NonUnit, bool TRANSPOSE = false, bool TRAILING_SYNC = true> void trsm(uint32_t n, uint32_t nrhs, const T *A, T *B)#

Triangular solve with multiple right-hand sides op(A) X = B, in place (TRSM).

Solves the triangular system for every column of B (n×nrhs, column-major), overwriting B with X. A is n×n column-major; only the triangle named by FILL is read. TRANSPOSE=true solves Aᵀ X = B against that same stored triangle; DIAG=Diag::Unit means an implicit unit diagonal. All right-hand sides share each elimination step’s two barriers, and the update is flat-strided over the (rows × nrhs) rectangle, so wide B keeps every thread busy even at small n. Ends on a barrier (composes cleanly). SciPy equivalent: X = scipy.linalg.solve_triangular(A, B, lower=(FILL==Lower), unit_diagonal=(DIAG==Unit), trans=(1 if TRANSPOSE else 0)).

Template Parameters:

T – Scalar type.
FILL – Which triangle of A holds the data (default FillMode::Lower).
DIAG – Diag::Unit for an implicit unit diagonal (default Diag::NonUnit).
TRANSPOSE – When true solve Aᵀ X = B (default false).

Parameters:

n – Dimension (A is n×n; each column of B has length n).
nrhs – Number of right-hand sides (columns of B).
A – Triangular matrix (column-major; read-only).
B – In/out right-hand sides (n×nrhs, column-major); on return holds X.

template<typename T, uint32_t N, uint32_t NRHS, FillMode FILL = FillMode::Lower, Diag DIAG = Diag::NonUnit, bool TRANSPOSE = false, bool TRAILING_SYNC = true> void trsm(const T *A, T *B)#

Triangular solve with multiple right-hand sides op(A) X = B, in place (TRSM), compile-time size.

Same as the runtime trsm but with the dimensions as template parameters. SciPy equivalent: X = scipy.linalg.solve_triangular(A, B, lower=(FILL==Lower), unit_diagonal=(DIAG==Unit), trans=(1 if TRANSPOSE else 0)).

Template Parameters:

T – Scalar type.
N – Dimension (A is N×N; each column of B has length N).
NRHS – Number of right-hand sides (columns of B).
FILL – Which triangle of A holds the data (default FillMode::Lower).
DIAG – Diag::Unit for an implicit unit diagonal (default Diag::NonUnit).
TRANSPOSE – When true solve Aᵀ X = B (default false).

Parameters:

A – Triangular matrix (column-major; read-only).
B – In/out right-hand sides (N×NRHS, column-major); on return holds X.

template<typename T, uint32_t N, FillMode FILL = FillMode::Lower, Diag DIAG = Diag::NonUnit, bool TRANSPOSE = false> void trsv(const T *A, T *b)

Single-warp triangular solve op(A) x = b in place (TRSV), compile-time size.

One 32-lane warp solves the triangular system for any {FILL, DIAG, TRANSPOSE} combination, overwriting b with x. A is column-major and only the triangle named by FILL is read; TRANSPOSE=true solves Aᵀx = b against that same stored triangle; DIAG=Diag::Unit skips the diagonal divide. Every pivot is broadcast from lane 0’s REGISTER via __shfl_sync (never a shared re-read of b[k]) — immune to the nvcc __restrict__ stale-cache miscompile (see warp::potrf). This is its OWN warp implementation (warp and block can’t share an impl: __shfl/__syncwarp vs __syncthreads). No shared scratch, no __syncthreads. SciPy: x = scipy.linalg.solve_triangular(A, b, lower=(FILL==Lower), unit_diagonal=(DIAG==Unit), trans=(1 if TRANSPOSE else 0)).

Template Parameters:

T – Scalar type.
N – Dimension (A is N×N, b has length N).
FILL – Which triangle of A holds the data (default FillMode::Lower).
DIAG – Diag::Unit for an implicit unit diagonal (default Diag::NonUnit).
TRANSPOSE – When true solve Aᵀx = b (default false).

Parameters:

A – Triangular matrix (column-major); only the FILL triangle read.
b – In/out right-hand side; on return holds the solution x.

template<typename T, uint32_t N, uint32_t NRHS, FillMode FILL = FillMode::Lower, Diag DIAG = Diag::NonUnit, bool TRANSPOSE = false> void trsm(const T *A, T *B)#

Single-warp triangular solve with multiple right-hand sides op(A) X = B (TRSM), compile-time size.

Single-thread triangular solve with multiple right-hand sides op(A) X = B, in place (TRSM), compile-time size.

Warp-per-problem parity with the block glass::trsm: one 32-lane warp solves all NRHS columns of B (N×NRHS, column-major) in place. Each elimination step resolves the pivot row across all columns (lane-strided) and flat-strides the update over the (rows × NRHS) rectangle, sharing the per-step __syncwarp() across every right-hand side. No shared scratch, no __syncthreads.

ONE thread solves the N×N triangular system for all NRHS columns of B (N×NRHS, column-major), overwriting B with X — for thread-per-problem solvers that pack 32 independent low-DOF problems into a warp. A is column-major and read-only; only the triangle named by FILL is read; TRANSPOSE=true solves Aᵀ X = B against that same stored triangle; DIAG=Diag::Unit skips the diagonal divide. No shared scratch, no barriers, no threadIdx read; operands may be thread-local register arrays. SciPy equivalent: X = scipy.linalg.solve_triangular(A, B, lower=(FILL==Lower), unit_diagonal=(DIAG==Unit), trans=(1 if TRANSPOSE else 0)).

Delegates to the same trsm_impl body the block surface uses, via ThreadBarrier (rank=0, size=1, no-op sync) — the same algorithm and operand order as glass::trsm<T, N, NRHS, …> on one thread, agreeing to a few ULP (FMA-contraction jitter; bit-identity across the two instantiations is NOT guaranteed — see test/test_thread.py).

Template Parameters:

T – Scalar type.
N – Dimension (A is N×N; each column of B has length N).
NRHS – Number of right-hand sides (columns of B).
FILL – Which triangle of A holds the data (default FillMode::Lower).
DIAG – Diag::Unit for an implicit unit diagonal (default Diag::NonUnit).
TRANSPOSE – When true solve Aᵀ X = B (default false).
T – Scalar type.
N – Dimension (A is N×N; each column of B has length N). N<=7 keeps a T[N*N] operand register-resident (measured ceiling, both dtypes — see the thread-tier constraints in CLAUDE.md); larger N — or a B wider than that element budget — still computes correctly but spills to local memory, forfeiting the tier’s premise.
NRHS – Number of right-hand sides (columns of B).
FILL – Which triangle of A holds the data (default FillMode::Lower).
DIAG – Diag::Unit for an implicit unit diagonal (default Diag::NonUnit).
TRANSPOSE – When true solve Aᵀ X = B (default false).

Parameters:

A – Triangular matrix (column-major; read-only).
B – In/out right-hand sides (N×NRHS, column-major); on return holds X.
A – Triangular matrix (column-major, N*N; read-only).
B – In/out right-hand sides (N×NRHS, column-major); on return holds X.

template<typename T, uint32_t N> void posv(T *A, T *b)#

Single-warp SPD solve A x = b via Cholesky (LAPACK posv), compile-time size.

One 32-lane warp solves the symmetric-positive-definite system A x = b in place: it factors A = L Lᵀ with warp::potrf (lower triangle overwrites A), then a forward solve L y = b and a back solve Lᵀ x = y (both warp::trsv). On return b holds x and the lower triangle of A holds L. This is the composed warp-per-problem solve — the proof that the warp L1/L2/L3 glue closes the gap. No shared scratch, no __syncthreads; every pivot broadcast from a register (§1g). A must be SPD (use double for ill-conditioned systems). NumPy equivalent: x = np.linalg.solve(A, b).

Template Parameters:

T – Scalar type.
N – Dimension (A is N x N, b has length N).

Parameters:

A – In/out SPD matrix (column-major); on return its lower triangle holds L.
b – In/out right-hand side; on return holds the solution x.

template<typename T, uint32_t N> void potrs(const T *L, T *b)#

Single-warp SPD solve from a precomputed Cholesky factor (LAPACK potrs), compile-time size.

Given the lower factor L (e.g. from warp::potrf), one 32-lane warp solves L Lᵀ x = b by forward then back substitution — the same two warp::trsv legs warp::posv composes, without the re-factor: the reusable-factor / multi-solve path. L is read-only; b is overwritten with x. No shared scratch, no __syncthreads; every pivot is broadcast from lane 0’s register via __shfl_sync (§1g — immune to the __restrict__ stale-shared-reread miscompile, see warp::trsv). SciPy equivalent: x = scipy.linalg.cho_solve((L, True), b).

Template Parameters:

T – Scalar type.
N – Dimension (L is N×N, b has length N).

Parameters:

L – Lower Cholesky factor (column-major, N*N; read-only).
b – In/out right-hand side; on return holds the solution x.

template<typename T, bool REG_DIAG = false, typename SizeT> void _posv_regularize(SizeT n, T *A, T rho)#

Add a diagonal regularization shift to A in place (single-warp helper).

Lane-strided over the n diagonal entries; REG_DIAG=false adds rho·I (Marquardt), REG_DIAG=true adds rho·diag(A) (Levenberg). Trailing __syncwarp() so the shifted A is warp-visible before factoring. Internal.

template<typename T, uint32_t N, uint32_t NRHS, bool REGULARIZE = false, bool CHECK = false, bool REG_DIAG = false> void posv(T *A, T *B, T rho = T(0), int *s_fail = nullptr)#

Single-warp regularized/checked multi-RHS SPD solve A X = B (LAPACK posv).

Warp-per-problem parity with the block multi-RHS glass::posv: one 32-lane warp optionally shifts A’s diagonal (REGULARIZE: rho·I, or rho·diag(A) when REG_DIAG), factors A = L Lᵀ via warp::potrf<…,CHECK> (reporting a non-PD pivot through s_fail), then forward/back-solves all NRHS columns of B at once with the multi-RHS warp::trsm. On return A holds L and B holds X. No shared scratch, no __syncthreads.

A flagged single-RHS solve is just NRHS=1 — the form HJCD’s LM step wants: warp::posv<T, DIM, 1, REGULARIZE=true, CHECK=true, REG_DIAG=true>(A, b, lambda, &s_fail) folds the A += lambda*diag(A) damping and the non-PD net into one call. (The unflagged 2-arg warp::posv<T,N>(A,b) stays; flags cannot live on it without colliding with this overload at NRHS in {0,1}.)

Template Parameters:

T – Scalar type (use double for ill-conditioned A).
N – Dimension (A is N x N, each column of B has length N).
NRHS – Number of right-hand sides (columns of B).
REGULARIZE – If true, shift A before factoring (default false, compiles out).
CHECK – If true, report a non-PD pivot via s_fail (default false, compiles out).
REG_DIAG – With REGULARIZE: shift by rho·diag(A) instead of rho·I (default false).

Parameters:

A – In/out SPD matrix (column-major); on return its lower triangle holds L.
B – In/out right-hand sides (N x NRHS, column-major); on return holds X.
rho – Diagonal shift applied when REGULARIZE (ignored otherwise).
s_fail – Optional non-PD flag when CHECK (set to 1 on a non-PD pivot, else 0).

namespace warp

namespace thread

Symmetric rank-k / rank-2k update (syrk / syr2k)#

Functions

bool syrk_in_canonical(FillMode fill, uint32_t row, uint32_t col)#

template<typename T, uint32_t N, uint32_t K, FillMode FILL, bool TRANSPOSE, bool ROW_MAJOR> void syrk_impl_ct(uint32_t rank, uint32_t size, T alpha, const T *A, T beta, T *C)#

template<typename T, uint32_t N, uint32_t K, FillMode FILL, bool TRANSPOSE, bool ROW_MAJOR> void syrk_impl_ct(uint32_t rank, uint32_t size, T alpha, const T *A, T *C)#

template<typename T, uint32_t N, uint32_t K, FillMode FILL, bool TRANSPOSE, bool ROW_MAJOR> void syr2k_impl_ct(uint32_t rank, uint32_t size, T alpha, const T *A, const T *B, T beta, T *C)#

template<typename T, uint32_t N, uint32_t K, FillMode FILL, bool TRANSPOSE, bool ROW_MAJOR> void syr2k_impl_ct(uint32_t rank, uint32_t size, T alpha, const T *A, const T *B, T *C)#

template<typename T, FillMode FILL = FillMode::Full, bool TRANSPOSE = false, bool ROW_MAJOR = false, bool TRAILING_SYNC = true> void syrk(uint32_t n, uint32_t k, T alpha, const T *A, T beta, T *C)#

Symmetric rank-k update: C = alpha * op(A) * op(A)^T + beta * C (SYRK).

Runtime-size, single-block, flat-element parallelism: each thread owns output cells of the n x n symmetric C strided over the block. The length-k dot is computed ONLY in the canonical triangle (the symmetry win, ~half the FLOPs of a GEMM); for Full the lower-cell-owning thread also writes the mirror C[col,row] (diagonal written once) so each cell is written exactly once and NO interior barrier is needed. Lower/Upper write only the named triangle and leave the other untouched.

NumPy equivalent: TRANSPOSE=false → alpha * A @ A.T + beta * C; TRANSPOSE=true → alpha * A.T @ A + beta * C.

Template Parameters:

T – Scalar type.
FILL – Which triangle of C to write (Lower / Upper / Full).
TRANSPOSE – If false, op(A)=A (A is n x k); if true, op(A)=A^T (A is k x n).
ROW_MAJOR – Storage order for A and C (false = column-major / Fortran).

Parameters:

n – Dimension of the symmetric result C (n x n).
k – Contraction length.
alpha – Scalar multiplier on the product.
A – Input matrix (n x k if TRANSPOSE=false, else k x n).
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out n x n symmetric result matrix.

template<typename T, FillMode FILL = FillMode::Full, bool TRANSPOSE = false, bool ROW_MAJOR = false, bool TRAILING_SYNC = true> void syrk(uint32_t n, uint32_t k, T alpha, const T *A, T *C)#

SYRK with implicit beta = 0: C = alpha * op(A) * op(A)^T (SYRK).

Runtime-size overload that overwrites C (C is overwritten, not read), avoiding the beta * C term — safe to write into uninitialized scratch. For Full, the full symmetric matrix is written; for Lower/Upper, only the named triangle is written and the other is left untouched. Single-block, flat-element parallelism; no interior barrier.

NumPy equivalent: TRANSPOSE=false → alpha * A @ A.T; TRANSPOSE=true → alpha * A.T @ A.

Template Parameters:

T – Scalar type.
FILL – Which triangle of C to write (Lower / Upper / Full).
TRANSPOSE – If false, op(A)=A (A is n x k); if true, op(A)=A^T (A is k x n).
ROW_MAJOR – Storage order for A and C (false = column-major / Fortran).

Parameters:

n – Dimension of the symmetric result C (n x n).
k – Contraction length.
alpha – Scalar multiplier on the product.
A – Input matrix (n x k if TRANSPOSE=false, else k x n).
C – Output n x n symmetric result matrix (overwritten, not read).

template<typename T, uint32_t N, uint32_t K, FillMode FILL = FillMode::Full, bool TRANSPOSE = false, bool ROW_MAJOR = false, bool TRAILING_SYNC = true> void syrk(T alpha, const T *A, T beta, T *C)#

Compile-time-size SYRK: C = alpha * op(A) * op(A)^T + beta * C (SYRK).

Dimensions are template parameters so the compiler unrolls the inner loop and replaces the el % N / el / N index math with magic-number multiplies. Single-block, flat-element parallelism, symmetry-exploiting (canonical triangle + mirror write), no interior barrier. C is read; caller must initialize it.

NumPy equivalent: TRANSPOSE=false → alpha * A @ A.T + beta * C; TRANSPOSE=true → alpha * A.T @ A + beta * C.

Template Parameters:

T – Scalar type.
N – Compile-time dimension of the symmetric result C (N x N).
K – Compile-time contraction length.
FILL – Which triangle of C to write (Lower / Upper / Full).
TRANSPOSE – If false, op(A)=A (A is N x K); if true, op(A)=A^T (A is K x N).
ROW_MAJOR – Storage order for A and C (false = column-major / Fortran).

Parameters:

alpha – Scalar multiplier on the product.
A – Input matrix (N x K if TRANSPOSE=false, else K x N).
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out N x N symmetric result matrix.

template<typename T, uint32_t N, uint32_t K, FillMode FILL = FillMode::Full, bool TRANSPOSE = false, bool ROW_MAJOR = false, bool TRAILING_SYNC = true> void syrk(T alpha, const T *A, T *C)#

Compile-time-size SYRK with implicit beta = 0: C = alpha * op(A) * op(A)^T.

Compile-time-size overload that overwrites C (C is overwritten, not read). Single-block, flat-element parallelism, symmetry-exploiting; no interior barrier. Safe to write into uninitialized scratch.

NumPy equivalent: TRANSPOSE=false → alpha * A @ A.T; TRANSPOSE=true → alpha * A.T @ A.

Template Parameters:

T – Scalar type.
N – Compile-time dimension of the symmetric result C (N x N).
K – Compile-time contraction length.
FILL – Which triangle of C to write (Lower / Upper / Full).
TRANSPOSE – If false, op(A)=A (A is N x K); if true, op(A)=A^T (A is K x N).
ROW_MAJOR – Storage order for A and C (false = column-major / Fortran).

Parameters:

alpha – Scalar multiplier on the product.
A – Input matrix (N x K if TRANSPOSE=false, else K x N).
C – Output N x N symmetric result matrix (overwritten, not read).

template<typename T, FillMode FILL = FillMode::Full, bool TRANSPOSE = false, bool ROW_MAJOR = false, bool TRAILING_SYNC = true> void syr2k(uint32_t n, uint32_t k, T alpha, const T *A, const T *B, T beta, T *C)#

Symmetric rank-2k update: C = alpha*(op(A)*op(B)^T + op(B)*op(A)^T) + beta*C (SYR2K).

Runtime-size, single-block, flat-element parallelism. The result is symmetric by construction; the length-k dot is computed only in the canonical triangle and (for Full) mirrored, so each cell is written once and NO interior barrier is needed. Lower/Upper write only the named triangle.

NumPy equivalent: TRANSPOSE=false → alpha*(A@B.T + B@A.T) + beta*C; TRANSPOSE=true → alpha*(A.T@B + B.T@A) + beta*C.

Template Parameters:

T – Scalar type.
FILL – Which triangle of C to write (Lower / Upper / Full).
TRANSPOSE – If false, op = identity (A,B are n x k); if true, op = transpose (A,B are k x n).
ROW_MAJOR – Storage order for A, B and C (false = column-major / Fortran).

Parameters:

n – Dimension of the symmetric result C (n x n).
k – Contraction length.
alpha – Scalar multiplier on the symmetrized product.
A, B – Input matrices (n x k if TRANSPOSE=false, else k x n).
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out n x n symmetric result matrix.

template<typename T, FillMode FILL = FillMode::Full, bool TRANSPOSE = false, bool ROW_MAJOR = false, bool TRAILING_SYNC = true> void syr2k(uint32_t n, uint32_t k, T alpha, const T *A, const T *B, T *C)#

SYR2K with implicit beta = 0: C = alpha*(op(A)*op(B)^T + op(B)*op(A)^T).

Runtime-size overload that overwrites C (C is overwritten, not read). Safe to write into uninitialized scratch. Single-block, flat-element parallelism, symmetry-exploiting; no interior barrier.

NumPy equivalent: TRANSPOSE=false → alpha*(A@B.T + B@A.T); TRANSPOSE=true → alpha*(A.T@B + B.T@A).

Template Parameters:

T – Scalar type.
FILL – Which triangle of C to write (Lower / Upper / Full).
TRANSPOSE – If false, op = identity (A,B are n x k); if true, op = transpose (A,B are k x n).
ROW_MAJOR – Storage order for A, B and C (false = column-major / Fortran).

Parameters:

n – Dimension of the symmetric result C (n x n).
k – Contraction length.
alpha – Scalar multiplier on the symmetrized product.
A, B – Input matrices (n x k if TRANSPOSE=false, else k x n).
C – Output n x n symmetric result matrix (overwritten, not read).

template<typename T, uint32_t N, uint32_t K, FillMode FILL = FillMode::Full, bool TRANSPOSE = false, bool ROW_MAJOR = false, bool TRAILING_SYNC = true> void syr2k(T alpha, const T *A, const T *B, T beta, T *C)#

Compile-time-size SYR2K: C = alpha*(op(A)*op(B)^T + op(B)*op(A)^T) + beta*C.

Dimensions are template parameters so the inner loop unrolls and el % N / el / N become magic-number multiplies. Single-block, flat-element parallelism, symmetry-exploiting; no interior barrier. C is read; caller must initialize it.

NumPy equivalent: TRANSPOSE=false → alpha*(A@B.T + B@A.T) + beta*C; TRANSPOSE=true → alpha*(A.T@B + B.T@A) + beta*C.

Template Parameters:

T – Scalar type.
N – Compile-time dimension of the symmetric result C (N x N).
K – Compile-time contraction length.
FILL – Which triangle of C to write (Lower / Upper / Full).
TRANSPOSE – If false, op = identity (A,B are N x K); if true, op = transpose (A,B are K x N).
ROW_MAJOR – Storage order for A, B and C (false = column-major / Fortran).

Parameters:

alpha – Scalar multiplier on the symmetrized product.
A, B – Input matrices (N x K if TRANSPOSE=false, else K x N).
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out N x N symmetric result matrix.

template<typename T, uint32_t N, uint32_t K, FillMode FILL = FillMode::Full, bool TRANSPOSE = false, bool ROW_MAJOR = false, bool TRAILING_SYNC = true> void syr2k(T alpha, const T *A, const T *B, T *C)#

Compile-time-size SYR2K with implicit beta = 0: C = alpha*(op(A)*op(B)^T + op(B)*op(A)^T).

Compile-time-size overload that overwrites C (C is overwritten, not read). Safe to write into uninitialized scratch. Single-block, flat-element parallelism, symmetry-exploiting; no interior barrier.

NumPy equivalent: TRANSPOSE=false → alpha*(A@B.T + B@A.T); TRANSPOSE=true → alpha*(A.T@B + B.T@A).

Template Parameters:

T – Scalar type.
N – Compile-time dimension of the symmetric result C (N x N).
K – Compile-time contraction length.
FILL – Which triangle of C to write (Lower / Upper / Full).
TRANSPOSE – If false, op = identity (A,B are N x K); if true, op = transpose (A,B are K x N).
ROW_MAJOR – Storage order for A, B and C (false = column-major / Fortran).

Parameters:

alpha – Scalar multiplier on the symmetrized product.
A, B – Input matrices (N x K if TRANSPOSE=false, else K x N).
C – Output N x N symmetric result matrix (overwritten, not read).

template<typename T, uint32_t N, uint32_t K, FillMode FILL = FillMode::Full, bool TRANSPOSE = false, bool ROW_MAJOR = false> void syrk(T alpha, const T *A, T beta, T *C)#

Single-warp SYRK C = alpha*op(A)*op(A)ᵀ + beta*C (compile-time size).

Single-thread compile-time-size SYRK: C = alpha * op(A) * op(A)^T + beta * C.

ONE thread computes the symmetric rank-k update serially (canonical triangle + mirror write for Full) — for thread-per-problem packing of low-DOF normal-equation builds (32 problems per warp). Reuses the same syrk_impl_ct body as the block/warp surfaces with (rank=0, size=1); no barriers, no shuffles, no threadIdx read, so operands may be thread-local register arrays (a T[N*N] C stays register-resident up to the tier’s measured N<=7 ceiling — see the thread-tier constraints in CLAUDE.md; larger N still computes correctly but spills to local memory). Same algorithm and operand order as the glass:: twin, agreeing to a few ULP (bit-identity across the two instantiations is NOT guaranteed).

See also

syrk (block form; identical math, (lane,32) element striping)

Template Parameters:

T – Scalar type.
N – Compile-time dimension of the symmetric result C (N x N).
K – Compile-time contraction length.
FILL – Which triangle of C to write (Lower / Upper / Full).
TRANSPOSE – If false, op(A)=A (A is N x K); if true, op(A)=A^T (A is K x N).
ROW_MAJOR – Storage order for A and C (false = column-major / Fortran).

Parameters:

alpha – Scalar multiplier on the product.
A – Input matrix (N x K if TRANSPOSE=false, else K x N).
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out N x N symmetric result matrix.

template<typename T, uint32_t N, uint32_t K, FillMode FILL = FillMode::Full, bool TRANSPOSE = false, bool ROW_MAJOR = false> void syrk(T alpha, const T *A, T *C)#

Single-warp SYRK with implicit beta = 0: C = alpha*op(A)*op(A)ᵀ (overwrite).

Single-thread SYRK with implicit beta = 0: C = alpha * op(A) * op(A)^T (overwrite).

Overwrites C (the existing C is not read — safe on uninitialized scratch). Otherwise identical to the beta overload above.

See also

syrk

Template Parameters:

T, N, K, FILL, TRANSPOSE, ROW_MAJOR – See the beta overload.

Parameters:

alpha – Scalar multiplier on the product.
A – Input matrix (N x K if TRANSPOSE=false, else K x N).
C – Output N x N symmetric result matrix (overwritten, not read).

template<typename T, uint32_t N, uint32_t K, FillMode FILL = FillMode::Full, bool TRANSPOSE = false, bool ROW_MAJOR = false> void syr2k(T alpha, const T *A, const T *B, T beta, T *C)#

Single-warp SYR2K C = alpha*(op(A)op(B)ᵀ + op(B)op(A)ᵀ) + beta*C (compile-time size).

Single-thread compile-time-size SYR2K: C = alpha*(op(A)*op(B)^T + op(B)*op(A)^T) + beta*C.

ONE thread computes the symmetric rank-2k update serially, reusing the shared syr2k_impl_ct body with (rank=0, size=1) — no barriers, no shuffles, no threadIdx read; operands may be thread-local register arrays (register-resident up to the tier’s N<=7 ceiling, CLAUDE.md). Same algorithm and operand order as the glass:: twin, agreeing to a few ULP (cross-tier bit-identity is NOT guaranteed).

See also

syr2k

Template Parameters:

T – Scalar type.
N – Compile-time dimension of the symmetric result C (N x N).
K – Compile-time contraction length.
FILL – Which triangle of C to write (Lower / Upper / Full).
TRANSPOSE – If false, op = identity (A,B are N x K); if true, op = transpose (A,B are K x N).
ROW_MAJOR – Storage order for A, B and C (false = column-major / Fortran).

Parameters:

alpha – Scalar multiplier on the symmetrized product.
A, B – Input matrices (N x K if TRANSPOSE=false, else K x N).
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out N x N symmetric result matrix.

template<typename T, uint32_t N, uint32_t K, FillMode FILL = FillMode::Full, bool TRANSPOSE = false, bool ROW_MAJOR = false> void syr2k(T alpha, const T *A, const T *B, T *C)#

Single-warp SYR2K with implicit beta = 0 (overwrite).

Single-thread SYR2K with implicit beta = 0 (overwrite).

Overwrites C (the existing C is not read — safe on uninitialized scratch). Otherwise identical to the beta overload above.

See also

syr2k

Template Parameters:

T, N, K, FILL, TRANSPOSE, ROW_MAJOR – See the beta overload.

Parameters:

alpha – Scalar multiplier on the symmetrized product.
A, B – Input matrices (N x K if TRANSPOSE=false, else K x N).
C – Output N x N symmetric result matrix (overwritten, not read).

namespace warp

namespace thread

Symmetric / triangular matrix-matrix multiply (symm / trmm)#

Left-side products against a triangle-stored matrix: symm (C = alpha*A_sym*B + beta*C, BLAS SYMM — only the FILL triangle of the symmetric A is stored, the other is read mirrored) and trmm (C = alpha*op(A_tri)*B, BLAS TRMM but deliberately OUT-of-place into a separate C so the flat one-output-per-thread loop stays race-free; FillMode / Diag / TRANSPOSE flags as in trsv/trsm).

SYMM (symmetric matrix-matrix multiply, left side) and TRMM (triangular matrix-matrix multiply, left side, out-of-place).

symm: C = alpha * A * B + beta * C where A is n x n SYMMETRIC with only the FILL triangle stored — the other triangle is never read, it is reconstructed by mirroring (A[i,k] reads A[k + i*n] when (i,k) falls outside the stored triangle). B and C are n x m column-major.

trmm: C = alpha * op(A) * B where A is n x n triangular (only the FILL triangle read; DIAG=Diag::Unit means an implicit unit diagonal that is never read) and op(A) = Aᵀ when TRANSPOSE. Out-of-place into C (deliberately NOT the BLAS in-place B := op(A) B — a separate output keeps the flat one-output-per-thread loop race-free with no interior barrier). C must not alias A or B.

Both use gemm’s plain-path parallelism: each thread owns disjoint output elements of the flat n*m space (el += size stride) with a serial ascending-k inner chain per element, so results are bit-identical at any thread count and NO interior barrier is needed (guide §1a counter-note). Block only (no warp:: variants yet — pack via warp::gemm on a materialized matrix in the meantime; future work).

Functions

template<typename T, FillMode FILL = FillMode::Lower, bool TRAILING_SYNC = true> void symm(uint32_t n, uint32_t m, T alpha, const T *A, const T *B, T beta, T *C)#

Symmetric matrix-matrix multiply (left side): C = alpha * A * B + beta * C (SYMM).

A is n x n symmetric with only the FILL triangle stored (column-major; the other triangle is never read — reconstructed by mirroring). B and C are n x m column-major. Single-block, flat one-output-per-thread parallelism (serial ascending-k chain per element ⇒ bit-identical at any thread count, no interior barrier). BLAS: SSYMM(‘L, uplo, …). NumPy (Lower):S = np.tril(A) + np.tril(A, -1).T; C = alpha * S @ B + beta * C`.

Template Parameters:

T – Scalar type.
FILL – Which triangle of A holds the data (default FillMode::Lower).
TRAILING_SYNC – Emit a trailing __syncthreads() (default true).

Parameters:

n – Dimension of A (n x n) and rows of B/C.
m – Columns of B/C.
alpha – Scalar multiplier on the product.
A – Symmetric matrix, FILL triangle stored (column-major; read-only).
B – Input matrix (n x m, column-major).
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out result matrix (n x m, column-major); must not alias A/B.

template<typename T, FillMode FILL = FillMode::Lower, bool TRAILING_SYNC = true> void symm(uint32_t n, uint32_t m, T alpha, const T *A, const T *B, T *C)#

SYMM with implicit beta = 0: C = alpha * A * B (overwrite).

Runtime-size overload that overwrites C (the existing C is never read) — safe to write into uninitialized scratch. NumPy (Lower): C = alpha * (np.tril(A) + np.tril(A, -1).T) @ B.

Template Parameters:

T – Scalar type.
FILL – Which triangle of A holds the data (default FillMode::Lower).
TRAILING_SYNC – Emit a trailing __syncthreads() (default true).

Parameters:

n, m – A is n x n; B/C are n x m.
alpha – Scalar multiplier on the product.
A – Symmetric matrix, FILL triangle stored (read-only).
B – Input matrix (n x m, column-major).
C – Output result matrix (overwritten); must not alias A/B.

template<typename T, uint32_t N, uint32_t M, FillMode FILL = FillMode::Lower, bool TRAILING_SYNC = true> void symm(T alpha, const T *A, const T *B, T beta, T *C)#

Compile-time-size SYMM: C = alpha * A * B + beta * C.

Same as the runtime symm but with the dimensions as template parameters (unrolled inner loop, magic-number %// indexing). BLAS: SSYMM(‘L, uplo, …). NumPy (Lower): C = alpha * (np.tril(A) + np.tril(A, -1).T) @ B + beta * C`.

Template Parameters:

T – Scalar type.
N – Dimension of A (N x N) and rows of B/C.
M – Columns of B/C.
FILL – Which triangle of A holds the data (default FillMode::Lower).
TRAILING_SYNC – Emit a trailing __syncthreads() (default true).

Parameters:

alpha – Scalar multiplier on the product.
A – Symmetric matrix, FILL triangle stored (read-only).
B – Input matrix (N x M, column-major).
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out result matrix (N x M, column-major); must not alias A/B.

template<typename T, uint32_t N, uint32_t M, FillMode FILL = FillMode::Lower, bool TRAILING_SYNC = true> void symm(T alpha, const T *A, const T *B, T *C)#

Compile-time-size SYMM with implicit beta = 0: C = alpha * A * B (overwrite).

Compile-time-size overload that overwrites C (never read) — safe into uninitialized scratch. NumPy (Lower): C = alpha * (np.tril(A) + np.tril(A, -1).T) @ B.

Template Parameters:

T – Scalar type.
N, M – A is N x N; B/C are N x M.
FILL – Which triangle of A holds the data (default FillMode::Lower).
TRAILING_SYNC – Emit a trailing __syncthreads() (default true).

Parameters:

alpha – Scalar multiplier on the product.
A – Symmetric matrix, FILL triangle stored (read-only).
B – Input matrix (N x M, column-major).
C – Output result matrix (overwritten); must not alias A/B.

template<typename T, FillMode FILL = FillMode::Lower, Diag DIAG = Diag::NonUnit, bool TRANSPOSE = false, bool TRAILING_SYNC = true> void trmm(uint32_t n, uint32_t m, T alpha, const T *A, const T *B, T *C)#

Triangular matrix-matrix multiply (left side): C = alpha * op(A) * B (TRMM, out-of-place).

A is n x n triangular, column-major; only the FILL triangle is read. TRANSPOSE=true multiplies by Aᵀ against that same stored triangle; DIAG=Diag::Unit means an implicit unit diagonal (A’s diagonal is never read). Unlike BLAS TRMM this writes a SEPARATE output C (n x m, column-major) instead of overwriting B in place — out-of-place keeps the flat one-output-per-thread loop race-free with no interior barrier. C must not alias A or B. The per-element k-chain covers only op(A)’s structural nonzeros, ascending ⇒ bit-identical at any thread count. BLAS: STRMM(‘L, uplo, transa, diag, …)(plus the copy). NumPy (Lower, NonUnit, no transpose):C = alpha * np.tril(A) @ B`.

Template Parameters:

T – Scalar type.
FILL – Which triangle of A holds the data (default FillMode::Lower).
DIAG – Diag::Unit for an implicit unit diagonal (default Diag::NonUnit).
TRANSPOSE – When true multiply by Aᵀ (default false).
TRAILING_SYNC – Emit a trailing __syncthreads() (default true).

Parameters:

n – Dimension of A (n x n) and rows of B/C.
m – Columns of B/C.
alpha – Scalar multiplier on the product.
A – Triangular matrix (column-major; only the FILL triangle read).
B – Input matrix (n x m, column-major; read-only).
C – Output result matrix (n x m, column-major, overwritten); no aliasing.

template<typename T, uint32_t N, uint32_t M, FillMode FILL = FillMode::Lower, Diag DIAG = Diag::NonUnit, bool TRANSPOSE = false, bool TRAILING_SYNC = true> void trmm(T alpha, const T *A, const T *B, T *C)#

Compile-time-size TRMM: C = alpha * op(A) * B (out-of-place).

Same as the runtime trmm but with the dimensions as template parameters (unrolled inner loop, magic-number %// indexing). NumPy (Lower, NonUnit, no transpose): C = alpha * np.tril(A) @ B.

Template Parameters:

T – Scalar type.
N – Dimension of A (N x N) and rows of B/C.
M – Columns of B/C.
FILL – Which triangle of A holds the data (default FillMode::Lower).
DIAG – Diag::Unit for an implicit unit diagonal (default Diag::NonUnit).
TRANSPOSE – When true multiply by Aᵀ (default false).
TRAILING_SYNC – Emit a trailing __syncthreads() (default true).

Parameters:

alpha – Scalar multiplier on the product.
A – Triangular matrix (column-major; only the FILL triangle read).
B – Input matrix (N x M, column-major; read-only).
C – Output result matrix (N x M, column-major, overwritten); no aliasing.

Diagonal-matrix multiply (dimm)#

Row/column scaling by a stored diagonal (cuBLAS dgmm analogue): C = alpha*diag(d)*B (default) or C = alpha*B*diag(d) (RIGHT). Contributed by Seyoung Yang.

Functions

template<typename T, bool RIGHT = false, bool TRAILING_SYNC = true> void dimm(uint32_t m, uint32_t n, T alpha, const T *d, const T *B, T *C)#

Diagonal-matrix scale of a dense matrix: C = alpha * diag(d) * B (RIGHT=false, scales ROWS) or C = alpha * B * diag(d) (RIGHT=true, scales COLUMNS).

cuBLAS analogue: cublasXdgmm.

B and C are m×n column-major; d holds the diagonal only — length m when RIGHT=false, length n when RIGHT=true. Pure elementwise: each output is owned by exactly one thread, so C may alias B (in-place scale). NumPy: C = alpha * np.diag(d) @ B / C = alpha * B @ np.diag(d).

Template Parameters:

T – Scalar type.
RIGHT – false: rows scaled by d[row] (diag on the left, default); true: columns scaled by d[col] (diag on the right).
TRAILING_SYNC – End on a barrier so C is valid for every thread on return (default true); callers owning the next barrier pass false to elide it.

Parameters:

m, n – Dimensions (B/C are m×n).
alpha – Scalar multiplier.
d – Diagonal entries (length m or n per RIGHT; read-only).
B – Input matrix (column-major; read-only, may alias C).
C – Output matrix (column-major).

template<typename T, uint32_t M, uint32_t N, bool RIGHT = false, bool TRAILING_SYNC = true> void dimm(T alpha, const T *d, const T *B, T *C)#

Compile-time-size diagonal-matrix scale (see the runtime dimm).

Template Parameters:

T – Scalar type.
M, N – Dimensions (B/C are M×N).
RIGHT – false: C = alpha*diag(d)*B; true: C = alpha*B*diag(d).
TRAILING_SYNC – End on a barrier (default true).

Parameters:

alpha – Scalar multiplier.
d – Diagonal entries (length M or N per RIGHT; read-only).
B – Input matrix (column-major; read-only, may alias C).
C – Output matrix (column-major).

LDLᵀ factorization (ldlt / ldlt_solve)#

Functions

template<typename T> constexpr std::size_t ldlt_scratch_bytes(uint32_t n)#

Scratch size in bytes for ldlt.

The pivot path uses n + 1 scratch elements (one broadcast slot for the argmax index + up to n working-row magnitudes fed to the Bunch–Kaufman rowmax scan); the non-pivoted path does not read it. Allocate ldlt_scratch_bytes<T>(n) bytes for the s_scratch argument so it is sized for both paths.

Template Parameters:: T – Scalar type.
Parameters:: n – Matrix dimension (A is n x n).
Returns:: Bytes to allocate for ldlt’s s_scratch.

template<typename T, bool CHECK = false, bool TRAILING_SYNC = true> void ldlt(uint32_t n, T *A, T *s_scratch, bool pivot = false, int32_t *piv = nullptr, int *s_fail = nullptr, int *s_inertia = nullptr)#

template<typename T, uint32_t N, bool CHECK = false, bool TRAILING_SYNC = true> void ldlt(T *A, T *s_scratch, bool pivot = false, int32_t *piv = nullptr, int *s_fail = nullptr, int *s_inertia = nullptr)#

Compile-time-size in-place LDLᵀ factorization (LAPACK sytrf, lower, optional Bunch–Kaufman pivoting).

Same as the runtime ldlt but with the dimension as a template parameter, letting the compiler bake N in. Factors a symmetric (possibly indefinite) A = L * D * Lᵀ in place. SciPy equivalence: lu, d, _ = scipy.linalg.ldl(A).

When CHECK is true, reports breakdowns via s_fail and the inertia via s_inertia (see the runtime overload). CHECK defaults false and compiles out, so the unchecked instantiation is byte-identical to the original.

Template Parameters:

T – Scalar type.
N – Matrix dimension (A is N x N).
CHECK – If true, report breakdowns and the inertia (default false, compiles out).

Parameters:

A – In/out N x N matrix (column-major); diagonal (+ 2×2 subdiagonal slots) holds D, strict lower holds L on return.
s_scratch – Shared scratch advertised as (N + 1) elements (used by the pivot path; non-pivoted path accepts nullptr).
pivot – If true, apply Bunch–Kaufman 1×1/2×2 pivoting (see the runtime overload).
piv – Out pivot array of N int32 entries (pivot path only); may be nullptr.
s_fail – Optional flag (CHECK only): 1 on a factorization breakdown, else 0. Ignored when null.
s_inertia – Optional 3 ints (CHECK only): {n_pos, n_neg, n_zero}. Ignored when null.

template<typename T, bool TRAILING_SYNC = true> void ldlt_solve(uint32_t n, const T *LD, T *b, const int32_t *piv = nullptr)#

template<typename T, uint32_t N, bool TRAILING_SYNC = true> void ldlt_solve(const T *LD, T *b, const int32_t *piv = nullptr)#

Compile-time-size LDLᵀ solve A x = b in place (LAPACK sytrs analogue).

Same as the runtime ldlt_solve but with the dimension as a template parameter. NumPy equivalence: x = np.linalg.solve(A, b). A non-null piv applies the recorded Bunch–Kaufman permutation and 2×2 diagonal blocks (factor made with pivot=true).

Template Parameters:

T – Scalar type.
N – Dimension (LD is N x N, b has length N).

Parameters:

LD – In LDLᵀ factor from ldlt (column-major; unit-L strict-lower, D diagonal).
b – In/out right-hand side; on return holds the solution x.
piv – Pivot array from the pivoted factorization, or nullptr (non-pivoted).

template<typename T, uint32_t N, bool CHECK = false> void ldlt(T *A, int *s_fail = nullptr, int *s_inertia = nullptr)#

Single-warp in-place LDLᵀ factorization (LAPACK sytrf, lower, NON-pivoted).

Single-thread in-place LDLᵀ factorization (LAPACK sytrf, lower, NON-pivoted), compile-time size.

Warp-per-problem parity with the block glass::ldlt: one 32-lane warp factors the symmetric (possibly INDEFINITE) A = L D Lᵀ in place — lane 0 runs the serial diagonal recurrence D_j = A_jj − Σ_{k<j} L_jk² D_k (broadcasting D_j from its register via __shfl_sync, never a shared re-read — immune to the __restrict__ stale-cache miscompile), lanes fill the trailing column L_ij = (A_ij − Σ_{k<j} L_ik D_k L_jk)/D_j strided by 32. On return the diagonal slots hold D, the strict lower triangle holds unit-L. No square root, so it factors KKT / saddle-point systems Cholesky cannot. No shared scratch, no __syncthreads. Non-pivoted (pivoting on the warp surface is deferred — it needs a warp::iamax over the working column; the block path covers the pivoted / Bunch–Kaufman case).

CHECK (compile-out) reports a zero/NaN pivot via s_fail and the inertia {n_pos, n_neg, n_zero} via s_inertia (lane 0 writes both). NumPy: lu, d, _ = scipy.linalg.ldl(A, lower=True) ⇒ A == lu @ np.diag(d) @ lu.T.

ONE thread factors the symmetric (possibly INDEFINITE) A = L D Lᵀ in place — the sequential recurrence D_j = A_jj − Σ_{k<j} L_jk² D_k, then L_ij = (A_ij − Σ_{k<j} L_ik D_k L_jk)/D_j down each column — for thread-per-problem solvers that pack 32 independent low-DOF problems into a warp. On return the diagonal slots hold D, the strict lower triangle holds unit-L (upper triangle untouched). No square root, so it factors KKT / saddle-point systems Cholesky cannot. No shared scratch, no barriers, no threadIdx read; A may live in a thread-local array and stay register-resident.

A FRESH serial body (not a ThreadBarrier instantiation of the block ldlt_impl): that impl’s runtime bool pivot branch would compile the whole Bunch–Kaufman path — shared-scratch broadcasts, iamax_lowmem, raw block-wide barriers — into the thread-tier function. Same algorithm and operand order as the NON-pivoted glass::ldlt<T, N> path (and warp::ldlt, itself a separate body), agreeing to a few ULP (FMA-contraction jitter; bit-identity across tiers is NOT guaranteed). Non-pivoted only — the tier is branch-free by contract (every lane owns a different problem, so the data-dependent Bunch–Kaufman branches would diverge across the warp; the block path covers the pivoted case). Every pivot D_j must be nonzero; a saddle like [[0,b],[b,0]] breaks down (use the block ldlt(..., pivot=true, piv) for those).

CHECK (compile-out, default false) reports a zero/NaN pivot via s_fail and the inertia {n_pos, n_neg, n_zero} via s_inertia. SciPy: lu, d, _ = scipy.linalg.ldl(A, lower=True) ⇒ A == lu @ np.diag(d) @ lu.T.

Template Parameters:

T – Scalar type (use double for ill-conditioned A).
N – Dimension (A is N x N).
CHECK – If true, report zero/NaN pivot + inertia (default false, compiles out).
T – Scalar type (use double for ill-conditioned A).
N – Dimension (A is N x N). N<=7 keeps A register-resident (measured ceiling, both dtypes — see the thread-tier constraints in CLAUDE.md); larger N still computes correctly but demotes A to local memory, forfeiting the tier’s premise.
CHECK – If true, report zero/NaN pivot + inertia (default false, compiles out).

Parameters:

A – In/out N x N symmetric matrix (column-major, lower); on return holds L (strict-lower, unit) and D (diagonal).
s_fail – Optional flag (CHECK only): set to 1 on a zero/NaN pivot, else 0.
s_inertia – Optional length-3 {n_pos, n_neg, n_zero} pivot-sign counts (CHECK only).
A – In/out N x N symmetric matrix (column-major, lower); on return holds L (strict-lower, unit) and D (diagonal).
s_fail – Optional flag (CHECK only): set to 1 on a zero/NaN pivot, else 0. Ignored when null.
s_inertia – Optional length-3 {n_pos, n_neg, n_zero} pivot-sign counts (CHECK only).

template<typename T, uint32_t N> void ldlt_solve(const T *LD, T *b)#

Single-warp LDLᵀ solve A x = b in place from an ldlt factor (NON-pivoted).

Single-thread LDLᵀ solve A x = b in place from an ldlt factor (LAPACK sytrs analogue, NON-pivoted), compile-time size.

Warp parity with block glass::ldlt_solve: one 32-lane warp runs the three sweeps — forward unit-L (L y = b), diagonal scale (z = y / D), back unit-Lᵀ (Lᵀ x = z) — over the factor LD from warp::ldlt. b is overwritten with x. No shared scratch; __syncwarp between dependent sweeps. Non-pivoted (matches the non-pivoted warp::ldlt). NumPy: x = np.linalg.solve(A, b).

ONE thread runs the three sweeps — forward unit-L (L y = b), diagonal scale (z = y / D), back unit-Lᵀ (Lᵀ x = z) — over the factor LD from thread::ldlt. LD is read-only; b is overwritten with x. No shared scratch, no barriers, no threadIdx read; operands may be thread-local register arrays. Non-pivoted only (matches thread::ldlt; see there for why the pivoted path is excluded from this tier). A fresh serial body for the same reason as thread::ldlt — the block ldlt_solve_impl interleaves runtime piv branches through every sweep. Same algorithm and operand order as the non-pivoted glass::ldlt_solve<T, N> path, agreeing to a few ULP (bit-identity across tiers is NOT guaranteed). NumPy: x = np.linalg.solve(A, b).

Template Parameters:

T – Scalar type.
N – Dimension (LD is N x N, b length N).
T – Scalar type.
N – Dimension (LD is N x N, b has length N). N<=7 keeps a T[N*N] factor register-resident (measured ceiling, both dtypes — see the thread-tier constraints in CLAUDE.md).

Parameters:

LD – LDLᵀ factor from warp::ldlt (column-major; unit-L strict-lower, D diagonal).
b – In/out right-hand side; on return holds the solution x.
LD – LDLᵀ factor from thread::ldlt (column-major; unit-L strict-lower, D diagonal).
b – In/out right-hand side; on return holds the solution x.

namespace warp

namespace thread

SPD solve (posv / potrs)#

SPD linear solve via Cholesky + two triangular solves (pure SIMT).

posv / potrs are thin single-block compositions of potrf (potrf.cuh) and trsv (trsv.cuh). Both callees end with a trailing __syncthreads(), so the factor and the two solves compose with NO inter-call barrier. Pure-SIMT companion to glass::nvidia::posv. Column-major throughout.

NOTE: glass::warp::posv is NOT in this file — it lives in trsm.cuh, after the warp::potrf/warp::trsm definitions it composes. glass::thread::posv IS here: it composes thread::potrf (potrf.cuh) with two single-RHS trsv_impl legs (trsv.cuh, L2), both of which glass.cuh already includes ahead of this file — so it needs nothing from trsm.cuh.

Functions

template<typename T, bool REG_DIAG = false, typename SizeT> void _posv_regularize(SizeT n, T *A, T rho)

Add a diagonal regularization shift to A in place (single-block helper).

REG_DIAG=false adds rho·I (Marquardt shift); REG_DIAG=true adds rho·diag(A), i.e. scales each diagonal by (1+rho) (Levenberg shift — scale-invariant across rows of very different magnitude, e.g. mixed prismatic/revolute Jacobians). Trailing __syncthreads() so the shifted A is block-visible before factoring. Internal; used by the flagged posv overloads.

template<typename T, bool TRAILING_SYNC = true> void posv(uint32_t n, T *A, T *b)#

template<typename T, uint32_t N> void posv(T *A, T *b)

Compile-time-size SPD solve A x = b (LAPACK posv).

Single-thread SPD solve A x = b via Cholesky (LAPACK posv), compile-time size.

Same as the runtime posv with the dimension as a template parameter. NumPy equivalent: x = np.linalg.solve(A, b) (A SPD). For the regularized / checked / Levenberg path use the multi-RHS overload with NRHS=1, e.g. posv<T, N, 1, true, true, true>(A, b, rho, s_fail).

ONE thread factors A = L Lᵀ in place, then forward-solves L y = b and back-solves Lᵀ x = y. On return A holds its lower Cholesky factor L and b holds the solution x. For thread-per-problem solvers packing 32 independent low-DOF systems into a warp (e.g. N≈7 IK normal equations, one seed per lane). No shared scratch, no barriers, no threadIdx read; operands may be thread-local register arrays. A must be SPD; behaviour on non-SPD input is undefined (the Cholesky step produces NaN, no info flag). NumPy equivalent: x = np.linalg.solve(A, b) (A SPD).

Unlike the block posv_impl — which hardcodes BlockBarrier and reads threadIdx for the trsv legs — this composes thread::potrf with two trsv_impl(0u, 1u, …) calls directly, so no barrier or threadIdx read survives. Same algorithm and operand order as glass::posv<T, N> on one thread, agreeing to within FMA-contraction jitter (a few ULP — see test/test_thread.py; bit-identity across instantiations is not guaranteed).

Template Parameters:

T – Scalar type.
N – Dimension (A is N×N, b has length N).
T – Scalar type (use double for stability on ill-conditioned A).
N – Dimension (A is N×N, b has length N). N<=7 keeps A register-resident (measured ceiling, both dtypes — see the thread-tier constraints in CLAUDE.md); larger N still computes correctly but demotes A to local memory, forfeiting the tier’s premise.

Parameters:

A – In/out SPD matrix (column-major); overwritten with its factor L.
b – In/out right-hand side; on return holds the solution x.
A – In/out SPD matrix (column-major); overwritten with its factor L.
b – In/out right-hand side; on return holds the solution x.

template<typename T, bool TRAILING_SYNC = true> void potrs(uint32_t n, const T *L, T *b)#

template<typename T, uint32_t N> void potrs(const T *L, T *b)

Compile-time-size SPD solve from a precomputed Cholesky factor (LAPACK potrs).

Single-thread SPD solve from a precomputed Cholesky factor (LAPACK potrs), compile-time size.

Given the lower factor L (e.g. from thread::potrf), solves L Lᵀ x = b by forward then back substitution — the reusable-factor path (no re-factor). L is read-only; b is overwritten with x. SciPy equivalent: x = scipy.linalg.cho_solve((L, True), b).

Template Parameters:

T – Scalar type.
N – Dimension.
T – Scalar type.
N – Dimension (L is N×N, b has length N).

Parameters:

L – Lower Cholesky factor (column-major, N*N; read-only).
b – In/out right-hand side; on return holds the solution x.
L – Lower Cholesky factor (column-major, N*N; read-only).
b – In/out right-hand side; on return holds the solution x.

template<typename T, bool REGULARIZE = false, bool CHECK = false, bool REG_DIAG = false, bool TRAILING_SYNC = true> void posv(uint32_t n, uint32_t nrhs, T *A, T *B, T rho = T(0), int *s_fail = nullptr)#

template<typename T, uint32_t N, uint32_t NRHS, bool REGULARIZE = false, bool CHECK = false, bool REG_DIAG = false, bool TRAILING_SYNC = true> void posv(T *A, T *B, T rho = T(0), int *s_fail = nullptr)#

Compile-time-size multi-RHS SPD solve A X = B (LAPACK posv).

Same as the runtime multi-RHS posv with the dimension and right-hand-side count as template parameters. B is N × NRHS column-major (column c at B + c*N). Factored once, solved per column. NumPy equivalent: X = np.linalg.solve(A, B) (A SPD).

The optional REGULARIZE / CHECK / REG_DIAG flags (default off, compile out) add a diagonal shift before factoring and report a non-PD pivot via s_fail — the fused regularize→factor→solve path posv<T, N, NRHS, true, true>(A, B, rho, s_fail). REG_DIAG (appended last so existing <…, true, true> callers are unaffected) switches the shift from rho·I to rho·diag(A) (Levenberg). A flagged single-RHS solve is just NRHS=1: posv<T, N, 1, true, true, true>(A, b, rho, s_fail).

Template Parameters:

T – Scalar type.
N – Dimension (A is N×N, each column of B has length N).
NRHS – Number of right-hand sides (columns of B).
REGULARIZE – If true, shift A before factoring (default false, compiles out).
CHECK – If true, report a non-PD pivot via s_fail (default false, compiles out).
REG_DIAG – With REGULARIZE: shift by rho·diag(A) instead of rho·I (default false).

Parameters:

A – In/out SPD matrix (column-major); overwritten with its factor L.
B – In/out right-hand sides (N×NRHS, column-major); on return holds X.
rho – Diagonal shift added to A when REGULARIZE (ignored otherwise).
s_fail – Optional non-PD flag when CHECK (set to 1 on a non-PD pivot, else 0).

template<typename T, bool TRAILING_SYNC = true> void potrs(uint32_t n, uint32_t nrhs, const T *L, T *B)#

template<typename T, uint32_t N, uint32_t NRHS, bool TRAILING_SYNC = true> void potrs(const T *L, T *B)#

Compile-time-size multi-RHS SPD solve from a precomputed Cholesky factor (LAPACK potrs).

B is N × NRHS column-major (column c at B + c*N). Solved per column, no re-factor. SciPy equivalent: X = scipy.linalg.cho_solve((L, True), B).

Template Parameters:

T – Scalar type.
N – Dimension.
NRHS – Number of right-hand sides (columns of B).

Parameters:

L – Lower Cholesky factor (column-major, N*N; read-only).
B – In/out right-hand sides (N×NRHS, column-major); on return holds X.

namespace thread

Symmetric eigenvalues (syev / eig_clamp)#

Small-matrix symmetric eigendecomposition via cyclic Jacobi (syev, NumPy np.linalg.eigh) and the eigenvalue-clamping consumer op (eig_clamp: decompose, floor the eigenvalues at eps, reconstruct V diag(max(W, eps)) Vᵀ in place) — the device-side replacement for a host Eigen::SelfAdjointEigenSolver round-trip when regularizing a Hessian mid-solve. Designed for n ≤ 32.

Symmetric eigendecomposition via cyclic Jacobi (syev) and the eigenvalue-clamping consumer op (eig_clamp).

Device-side small symmetric eigensolver (target sizes n ≤ 32 — robot state/control dimensions). Kills the host round-trip solvers like PDDP currently make to Eigen::SelfAdjointEigenSolver mid-solve just to clamp Hessian eigenvalues: eig_clamp does decompose → clamp → reconstruct entirely inside the block.

Algorithm (classical cyclic Jacobi): work on a copy B = A in scratch and accumulate V = I. For each pair (p, q), p < q, in a FIXED cyclic order, compute the Jacobi rotation (c, s) annihilating B[p,q] (standard stable formulas: theta = (B_qq - B_pp)/(2 B_pq), t = sign(theta)/(|theta| + sqrt(1 + theta^2)), c = 1/sqrt(1 + t^2), s = t*c) and apply it to rows/cols p, q of B and columns p, q of V. Sweeps repeat until off(B)_F <= eps_T * ||A||_F (checked once per sweep) or a fixed cap of 15 sweeps; typical convergence is 5–8 sweeps for n ≤ 32. Jacobi always converges for (finite) symmetric input, so there is no CHECK-style failure flag; NaN/Inf input propagates into W/V (garbage-in, NaN-out — the deterministic sort still terminates).

Thread-count invariance (bit-identical at any block size): the (p, q) pair loop is SERIAL in a fixed order; rank 0 computes (c, s) into shared slots followed by a barrier (one FMA chain, identical for every launch); the row/col rotation is thread-strided over the n affected indices, where index k owns exactly the entries {(k,p),(k,q),(p,k),(q,k)} of B (mirror writes keep B exactly symmetric) and row k of V — its new values read only slots owned by the same k (the pivot 2x2 block is owned by k == p), so the rotation phase has no cross-thread read-after-write hazard and needs no staging. The sweep-end off(B) probe and the final ascending selection sort run serially on rank 0 (deterministic), each published through a shared slot + barrier; the sorted copy-out permutes V’s columns in parallel through the no-longer-needed B scratch.

Functions

template<typename T> constexpr T syev_eps()#

Machine epsilon by scalar width (float / double).

Used to scale syev’s deterministic convergence / rotation-skip thresholds. Spelled locally (sizeof-keyed constants) so no <limits> lands inside namespace glass (this header is included inside the namespace).

Template Parameters:: T – Scalar type (4-byte -> FLT_EPSILON, 8-byte -> DBL_EPSILON).
Returns:: Machine epsilon of T.

template<typename T> constexpr std::size_t syev_scratch_bytes(uint32_t n)#

Scratch size in bytes for syev.

Exact layout (in T elements): n*n for the working copy B (reused at the end as the eigenvector permutation staging buffer) + n slots holding the ascending sort permutation (stored as uint32_t, one per T slot) + 4 control slots (c, s, the sweep-converged flag, one pad).

Template Parameters:: T – Scalar type.
Parameters:: n – Matrix dimension (A is n x n).
Returns:: Bytes to allocate for syev’s s_scratch.

template<typename T, bool TRAILING_SYNC = true> void syev(uint32_t n, const T *A, T *W, T *V, T *s_scratch)#

template<typename T, uint32_t N, bool TRAILING_SYNC = true> void syev(const T *A, T *W, T *V, T *s_scratch)#

Compile-time-size symmetric eigendecomposition (cyclic Jacobi).

Same as the runtime syev but with the dimension as a template parameter, letting the compiler bake N in (constant-folded trip counts / indexing). A is preserved; W ascending; V column i ↔ W[i]. NumPy equivalent: W, V = np.linalg.eigh(A).

Template Parameters:

T – Scalar type.
N – Matrix dimension (A is N x N; designed for N <= 32).

Parameters:

A – In: N x N symmetric matrix (column-major). NOT modified.
W – Out: N eigenvalues, ascending.
V – Out: N x N eigenvectors (column-major; column i ↔ W[i]).
s_scratch – Shared scratch of syev_scratch_bytes<T>(N) bytes.

template<typename T> constexpr std::size_t eig_clamp_scratch_bytes(uint32_t n)#

Scratch size in bytes for eig_clamp.

Exact layout (in T elements): n eigenvalues + n*n eigenvectors + syev’s own scratch (n*n + n + 4 — see syev_scratch_bytes), i.e. 2*n*n + 2*n + 4 elements total.

Template Parameters:: T – Scalar type.
Parameters:: n – Matrix dimension (A is n x n).
Returns:: Bytes to allocate for eig_clamp’s s_scratch.

template<typename T, bool TRAILING_SYNC = true> void eig_clamp(uint32_t n, T *A, T eps, T *s_scratch)#

template<typename T, uint32_t N, bool TRAILING_SYNC = true> void eig_clamp(T *A, T eps, T *s_scratch)#

Compile-time-size eigenvalue clamp A := V diag(max(W, eps)) Vᵀ.

Same as the runtime eig_clamp but with the dimension as a template parameter. NumPy equivalent: W, V = np.linalg.eigh(A); A = (V * np.maximum(W, eps)) @ V.T.

Template Parameters:

T – Scalar type.
N – Matrix dimension (A is N x N; designed for N <= 32).

Parameters:

A – In/out: N x N symmetric matrix (column-major); on return holds the eigenvalue-clamped reconstruction (SPD).
eps – Eigenvalue floor.
s_scratch – Shared scratch of eig_clamp_scratch_bytes<T>(N) bytes.

Cooperative-groups variants (`glass::cgrps::`)#

Functions

template<typename T, bool TRANSPOSE_A = false, bool TRANSPOSE_B = false, bool ROW_MAJOR_C = false, bool TRAILING_SYNC = true> void gemm(uint32_t m, uint32_t n, uint32_t k, T alpha, const T *A, const T *B, T beta, T *C, cgrps::thread_group g = cgrps::this_thread_block())#

GEMM: C = alpha * A * op(B) + beta * C (cooperative-groups variant).

Runtime-size, single-block; thread rank/size come from the cooperative group. Storage order is uniform across A, B, C (ROW_MAJOR; false = column-major). NumPy equivalent: C = alpha * A @ B + beta * C.

Template Parameters:

T – Scalar type.
TRANSPOSE_A – If true, A is k×m and op(A)=Aᵀ (else A is m×k).
TRANSPOSE_B – If true, B is n×k and op(B)=Bᵀ (else B is k×n).
ROW_MAJOR_C – Output storage order (false = column-major).

Parameters:

m, n, k – Dimensions: C is m×n, contraction k.
alpha – Scalar multiplier on the product.
A, B – Input matrices.
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out result matrix.
g – Cooperative thread group (defaults to the whole block).

template<typename T, bool TRANSPOSE_A = false, bool TRANSPOSE_B = false, bool ROW_MAJOR_C = false, bool TRAILING_SYNC = true> void gemm(uint32_t m, uint32_t n, uint32_t k, T alpha, const T *A, const T *B, T *C, cgrps::thread_group g = cgrps::this_thread_block())#

GEMM with implicit beta = 0: C = alpha * A * op(B) (cooperative-groups variant).

Runtime-size overload that overwrites C (the existing C is not read). NumPy equivalent: C = alpha * A @ B.

Template Parameters:

T – Scalar type.
TRANSPOSE_A – If true, A is k×m and op(A)=Aᵀ (else A is m×k).
TRANSPOSE_B – If true, B is n×k and op(B)=Bᵀ (else B is k×n).
ROW_MAJOR_C – Output storage order (false = column-major).

Parameters:

m, n, k – Dimensions: C is m×n, contraction k.
alpha – Scalar multiplier on the product.
A, B – Input matrices.
C – Output result matrix (overwritten).
g – Cooperative thread group (defaults to the whole block).

template<typename T, uint32_t M, uint32_t N, uint32_t K, bool TRANSPOSE_A = false, bool TRANSPOSE_B = false, bool ROW_MAJOR_C = false, bool TRAILING_SYNC = true> void gemm(T alpha, const T *A, const T *B, T beta, T *C, cgrps::thread_group g = cgrps::this_thread_block())#

Compile-time-size GEMM: C = alpha * op(A) * op(B) + beta * C (cooperative-groups variant).

Dimensions baked in as template parameters; standard convention (C is M×N, contraction K). NumPy: C = alpha * opA(A) @ opB(B) + beta * C.

Template Parameters:

T – Scalar type.
M, N, K – C is M×N, contraction K (see gemm.cuh).
TRANSPOSE_A – If true, A is K×M and op(A)=Aᵀ.
TRANSPOSE_B – If true, B is N×K and op(B)=Bᵀ.
ROW_MAJOR_C – Output storage order (false = column-major).

Parameters:

alpha – Scalar multiplier on the product.
A, B – Input matrices.
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out result matrix.
g – Cooperative thread group (defaults to the whole block).

template<typename T, uint32_t M, uint32_t N, uint32_t K, bool TRANSPOSE_A = false, bool TRANSPOSE_B = false, bool ROW_MAJOR_C = false, bool TRAILING_SYNC = true> void gemm(T alpha, const T *A, const T *B, T *C, cgrps::thread_group g = cgrps::this_thread_block())#

Compile-time-size GEMM with implicit beta = 0: C = alpha * op(A) * op(B) (cooperative-groups variant).

Overwrites C (the existing C is not read). NumPy: C = alpha * opA(A) @ opB(B).

Template Parameters:

T – Scalar type.
M, N, K – C is M×N, contraction K (see gemm.cuh).
TRANSPOSE_A – If true, A is K×M and op(A)=Aᵀ.
TRANSPOSE_B – If true, B is N×K and op(B)=Bᵀ.
ROW_MAJOR_C – Output storage order (false = column-major).

Parameters:

alpha – Scalar multiplier on the product.
A, B – Input matrices.
C – Output result matrix (overwritten).
g – Cooperative thread group (defaults to the whole block).

template<typename T> void inv(uint32_t dimA, T *A, T *s_scratch, cgrps::thread_group g = cgrps::this_thread_block())#

In-place matrix inverse via Gauss-Jordan on [A | I] (cooperative-groups variant).

Reduces a column-major augmented dimA x (2*dimA) matrix [A | I] so columns dimA..2*dimA-1 hold A^-1 on return. Serial pivot loop, block-parallel cell updates. NumPy equivalent: Ainv = np.linalg.inv(A).

Template Parameters:

T – Scalar type.

Parameters:

dimA – Matrix dimension (A is dimA x dimA).
A – In/out augmented [A | I] buffer (column-major, dimA x 2*dimA); on return its right half holds A^-1.
s_scratch – Shared scratch of (2*dimA + 1) * sizeof(T) bytes.
g – Cooperative thread group (defaults to the whole block).

template<typename T, bool CHECK = false> void potrf(uint32_t n, T *s_A, cgrps::thread_group g = cgrps::this_thread_block(), int *s_fail = nullptr)#

In-place Cholesky factorization of an SPD matrix (cooperative-groups variant).

Factors A = L * L^T and overwrites A with the lower-triangular factor L (only the lower triangle is written; the upper triangle keeps its input values). A must be symmetric positive-definite, column-major. NumPy equivalent: L = np.linalg.cholesky(A).

When CHECK is true and s_fail is non-null, rank 0 sets *s_fail = 1 on a non-PD / NaN pivot (else 0). CHECK defaults false and compiles out.

Template Parameters:

T – Scalar type.
CHECK – If true, detect a non-PD pivot and report it via s_fail (default false, compiles out).

Parameters:

n – Matrix dimension (A is n x n).
s_A – In/out n x n matrix (column-major); on return its lower triangle holds L.
g – Cooperative thread group (defaults to the whole block).
s_fail – Optional flag (CHECK only): set to 1 on a non-PD / NaN pivot, else 0. Ignored when null.

template<typename T, FillMode FILL = FillMode::Lower, Diag DIAG = Diag::NonUnit, bool TRANSPOSE = false> void trsm(uint32_t n, uint32_t nrhs, const T *A, T *B, cgrps::thread_group g = cgrps::this_thread_block())#

Triangular solve with multiple right-hand sides op(A) X = B, in place (cooperative-groups variant).

Solves the triangular system for every column of B (n×nrhs, column-major), overwriting B with X; flags match the block glass::trsm (FILL names the stored triangle, DIAG the implicit-unit choice, TRANSPOSE solves Aᵀ X = B). SciPy equivalent: X = scipy.linalg.solve_triangular(A, B, lower=(FILL==Lower), unit_diagonal=(DIAG==Unit), trans=(1 if TRANSPOSE else 0)).

Template Parameters:

T – Scalar type.
FILL – Which triangle of A holds the data (default FillMode::Lower).
DIAG – Diag::Unit for an implicit unit diagonal (default Diag::NonUnit).
TRANSPOSE – When true solve Aᵀ X = B (default false).

Parameters:

n – Dimension (A is n×n; each column of B has length n).
nrhs – Number of right-hand sides (columns of B).
A – Triangular matrix (column-major; read-only).
B – In/out right-hand sides (n×nrhs, column-major); on return holds X.
g – Cooperative thread group (defaults to the whole block).

template<typename T, uint32_t M, uint32_t N, uint32_t K, bool TRANSPOSE_A = false, bool TRANSPOSE_B = false, bool ROW_MAJOR_C = false, bool TRAILING_SYNC = true> void gemm_reduced(T alpha, T *A, T *B, T beta, T *C, cgrps::thread_group g = cgrps::this_thread_block())#

Contraction-parallel GEMM: C = alpha * A * op(B) + beta * C (cooperative-groups variant).

The cooperative-groups form of glass::gemm_reduced: one warp owns each output and its lanes split the contraction. Pass a warp-multiple group.

Template Parameters:

T – Scalar type.
M, N, K – C is M×N, contraction K (see gemm.cuh).
TRANSPOSE_A – If true, A is K×M and op(A)=Aᵀ.
TRANSPOSE_B – If true, B is N×K and op(B)=Bᵀ.
ROW_MAJOR_C – Output storage order (false = column-major).
TRAILING_SYNC – Emit a trailing g.sync() (default true).

Parameters:

alpha – Scalar multiplier on the product.
A, B – Input matrices.
beta – Scalar multiplier on the existing C (read only when beta != 0).
C – In/out result matrix.
g – Cooperative thread group (defaults to the whole block).

template<typename T, uint32_t M, uint32_t N, uint32_t K, bool TRANSPOSE_A = false, bool TRANSPOSE_B = false, bool ROW_MAJOR_C = false, bool TRAILING_SYNC = true> void gemm_reduced(T alpha, T *A, T *B, T *C, cgrps::thread_group g = cgrps::this_thread_block())#

Contraction-parallel GEMM with implicit beta = 0: C = alpha * A * op(B) (cooperative-groups variant).

Overwrites C (the existing C is not read). Otherwise identical to the beta overload above; pass a warp-multiple group.

Template Parameters:

T – Scalar type.
M, N, K – C is M×N, contraction K (see gemm.cuh).
TRANSPOSE_A – If true, A is K×M and op(A)=Aᵀ.
TRANSPOSE_B – If true, B is N×K and op(B)=Bᵀ.
ROW_MAJOR_C – Output storage order (false = column-major).
TRAILING_SYNC – Emit a trailing g.sync() (default true).

Parameters:

alpha – Scalar multiplier on the product.
A, B – Input matrices.
C – Output result matrix (overwritten).
g – Cooperative thread group (defaults to the whole block).

template<typename T, uint32_t K, uint32_t A, uint32_t B, TensorAxis CONTRACT = TensorAxis::K, bool SYMMETRIC = false, bool ACCUMULATE = true, bool TIN_ROW_MAJOR = false, bool TRAILING_SYNC = true> void tensor_vec_contract(const T *Tns, const T *v, T *Mout, cgrps::thread_group g = cgrps::this_thread_block())#

Tensor ⊗ vector contraction: Mout (+)= Σ_c v[c] · T[..c..] (cooperative-groups variant).

Cooperative-groups form of glass::tensor_vec_contract. See it for semantics.

Template Parameters:

T, K, A, B, CONTRACT, SYMMETRIC, ACCUMULATE, TIN_ROW_MAJOR – See glass::tensor_vec_contract.
TRAILING_SYNC – Emit a trailing g.sync() (default true).

Parameters:

Tns, v, Mout – See glass::tensor_vec_contract.
g – Cooperative thread group (defaults to the whole block; pass a warp-multiple group).

template<typename T, uint32_t K, uint32_t A, uint32_t B, bool ACCUMULATE = false, bool TIN_ROW_MAJOR = false, bool TRAILING_SYNC = true> void vec_tensor_vec(const T *Tns, const T *u, const T *w, T *s, cgrps::thread_group g = cgrps::this_thread_block())#

Vector–tensor–vector triple product: s[k] (+)= u^T · T_k · w (cooperative-groups variant).

Cooperative-groups form of glass::vec_tensor_vec. See it for semantics.

Template Parameters:

T, K, A, B, ACCUMULATE, TIN_ROW_MAJOR – See glass::vec_tensor_vec.
TRAILING_SYNC – Emit a trailing g.sync() (default true).

Parameters:

Tns, u, w, s – See glass::vec_tensor_vec.
g – Cooperative thread group (defaults to the whole block; pass a warp-multiple group).

template<typename T, uint32_t N, uint32_t Kdim, bool ACCUMULATE = false, bool TRAILING_SYNC = true> void congruence_sym(T alpha, const T *X, const T *M, T beta, T *Q, T *s_scratch, cgrps::thread_group g = cgrps::this_thread_block())#

Symmetric congruence: Q = alpha * Xᵀ·M·X + beta * Q (cooperative-groups variant).

Cooperative-groups form of glass::congruence_sym. See it for semantics.

Template Parameters:

T, N, Kdim, ACCUMULATE – See glass::congruence_sym.
TRAILING_SYNC – Emit a trailing g.sync() (default true).

Parameters:

alpha, X, M, beta, Q, s_scratch – See glass::congruence_sym.
g – Cooperative thread group (defaults to the whole block; pass a warp-multiple group).

template<typename T, uint32_t N, uint32_t P, uint32_t Qd, bool ACCUMULATE = false, bool TRAILING_SYNC = true> void bilinear(T alpha, const T *X, const T *M, const T *Y, T beta, T *R, T *s_scratch, cgrps::thread_group g = cgrps::this_thread_block())#

General bilinear form: R = alpha * Xᵀ·M·Y + beta * R (cooperative-groups variant).

Cooperative-groups form of glass::bilinear. See it for semantics.

Template Parameters:

T, N, P, Qd, ACCUMULATE – See glass::bilinear.
TRAILING_SYNC – Emit a trailing g.sync() (default true).

Parameters:

alpha, X, M, Y, beta, R, s_scratch – See glass::bilinear.
g – Cooperative thread group (defaults to the whole block; pass a warp-multiple group).

template<typename T, uint32_t ROWS, uint32_t COLS, bool TRANSPOSE = false, bool TRAILING_SYNC = true> void syrk_reduced(T alpha, const T *A, T beta, T *C, cgrps::thread_group g = cgrps::this_thread_block())#

Contraction-parallel symmetric rank-k update: C = alpha * A·op(A) + beta * C (cooperative-groups variant).

Cooperative-groups form of glass::syrk_reduced. Pass a warp-multiple group.

Template Parameters:

T, ROWS, COLS, TRANSPOSE – See glass::syrk_reduced.
TRAILING_SYNC – Emit a trailing g.sync() (default true).

Parameters:

alpha, A, beta, C – See glass::syrk_reduced.
g – Cooperative thread group (defaults to the whole block; pass a warp-multiple group).

template<typename T, uint32_t ROWS, uint32_t COLS, bool TRANSPOSE = false, bool TRAILING_SYNC = true> void syrk_reduced(T alpha, const T *A, T *C, cgrps::thread_group g = cgrps::this_thread_block())#

Contraction-parallel SYRK with implicit beta = 0: C = alpha * A·op(A) (cooperative-groups variant).

Template Parameters:: T, ROWS, COLS, TRANSPOSE, TRAILING_SYNC – See the beta overload.
Parameters:: alpha, A, C, g – See the beta overload.

L3 — Matrix Operations#

gemm#

Contraction-parallel gemm (gemm_reduced)#

Tensor ⊗ vector contractions (tensor_vec_contract / vec_tensor_vec)#

Congruence / bilinear forms (congruence_sym / bilinear)#

Contraction-parallel syrk (syrk_reduced)#

Riccati feedback gain (riccati_gain)#

Strided gemm#

Indexed / batched gemm#

Matrix inverse#

Cholesky#

LU with partial pivoting (getrf / getrs / gesv / laswp)#

Triangular solve (trsm)#

Symmetric rank-k / rank-2k update (syrk / syr2k)#

Symmetric / triangular matrix-matrix multiply (symm / trmm)#

Diagonal-matrix multiply (dimm)#

LDLᵀ factorization (ldlt / ldlt_solve)#

SPD solve (posv / potrs)#

Symmetric eigenvalues (syev / eig_clamp)#

Cooperative-groups variants (glass::cgrps::)#

Contraction-parallel gemm (`gemm_reduced`)#

Tensor ⊗ vector contractions (`tensor_vec_contract` / `vec_tensor_vec`)#

Congruence / bilinear forms (`congruence_sym` / `bilinear`)#

Contraction-parallel syrk (`syrk_reduced`)#

Riccati feedback gain (`riccati_gain`)#

Cooperative-groups variants (`glass::cgrps::`)#