Block-tridiagonal Ops (glass::bdmv / glass::bdsv)

Single-block matrix-vector product for a block-tridiagonal matrix — the sparsity pattern of the KKT / Schur systems that arise in trajectory optimization and MPC. The matrix is stored as NumBlockRows contiguous BlockSize × (3·BlockSize) row-major strips laid out [L | D | R] (left / diagonal / right blocks), and the vectors use a padded layout (NumBlockRows + 2)·BlockSize (one BlockSize pad block on each end) so the edge block-rows need no special case — their absent L / R simply multiply the zero pad.

See Block-tridiagonal Solves for the layout in detail, and Preconditioned Conjugate Gradient (glass::pcg) for the conjugate-gradient solver built on top of this matvec.

Block-tridiagonal matrix-vector product (glass::bdmv).

Single-block, thread-count-invariant matvec for a block-tridiagonal matrix stored as NumBlockRows contiguous strips. Each block-row strip is a BlockSize x (3*BlockSize) row-major tile laid out [L | D | R] (left/diagonal/right BlockSize x BlockSize blocks). The block-row at index br starts at s_matrix + br * (3*BlockSize) * BlockSize.

The input/output vectors use the padded layout (NumBlockRows + 2) * BlockSize: one BlockSize pad block on each end. Block-row br reads the window s_vector[br*BlockSize : br*BlockSize + 3*BlockSize) (= previous / current / next state blocks) and writes its result at s_output[(br+1)*BlockSize + row]. Edge block-rows (0 and N-1) work with no special case because their absent L / R multiply the zero pad — the caller must pre-zero the leading/trailing pad blocks of s_vector.

No trailing __syncthreads() — the caller barriers before reusing the output (matches the rest of the GLASS surface).

Functions

template<typename T, uint32_t NumBlockRows, uint32_t BlockSize, bool TRAILING_SYNC = true>
void bdmv(T *s_output, const T *s_matrix, const T *s_vector)

Block-tridiagonal matvec: s_output = A_bd * s_vector.

Template Parameters:

• T – Scalar type (e.g. float, double).

• NumBlockRows – Number of block-rows.

• BlockSize – Block dimension (rows/cols of each L/D/R block).

Parameters:

• s_output – Padded output, length (NumBlockRows+2)*BlockSize; result for block-row br written at (br+1)*BlockSize + row.

• s_matrix – Block-tridiagonal strips, [L|D|R] row-major per block-row.

• s_vector – Padded input, length (NumBlockRows+2)*BlockSize.

template<typename T, uint32_t NumBlockRows, uint32_t BlockSize, bool TRAILING_SYNC = true>
void bdmv(T *s_output_1, T *s_output_2, const T *s_matrix, const T *s_vector)

Block-tridiagonal matvec writing the result to two buffers at once.

Computes A_bd * s_vector and stores the (identical) result into both s_output_1 and s_output_2 in a single pass — e.g. the PCG initialization z = p = Pinv * r.

Template Parameters:

• T – Scalar type (e.g. float, double).

• NumBlockRows – Number of block-rows.

• BlockSize – Block dimension.

Parameters:

• s_output_1 – First padded output, length (NumBlockRows+2)*BlockSize.

• s_output_2 – Second padded output, length (NumBlockRows+2)*BlockSize.

• s_matrix – Block-tridiagonal strips, [L|D|R] row-major per block-row.

• s_vector – Padded input, length (NumBlockRows+2)*BlockSize.

Direct solve (glass::bdsv)

Block-Cholesky (block-Thomas) direct factor + solve on the same strips — the exact one-sweep alternative to Preconditioned Conjugate Gradient (glass::pcg) for SPD block-tridiagonal systems: bdsv_factor (in place; MAIN ← Cholesky factor, LEFT ← coupling block), bdsv_solve (forward/backward block substitution, reusable per right-hand side), and the fused bdsv. Composed from potrf/trsm/syrk/gemv/trsv.

Direct block-tridiagonal SPD factor + solve (glass::bdsv).

Block-Cholesky (block-Thomas) sweep over a block-tridiagonal SPD matrix in the same [L|D|R] strip layout as glass::bdmv / glass::pcg: block-row br is a BlockSize × (3*BlockSize) row-major tile starting at s_matrix + br * 3*BlockSize*BlockSize, and vectors use the padded (NumBlockRows+2)*BlockSize layout (one zero pad block each end).

This is the DIRECT alternative to the iterative glass::pcg for the Riccati/KKT/Schur systems consumers assemble in this layout — exact in one sweep, no convergence tuning; serial over the NumBlockRows knots (the dependency chain is inherent) with full in-block parallelism inside each knot’s potrf/trsm/syrk/gemv/trsv.

Factorization (bdsv_factor, in place on the strips): F_0 = chol(D_0); for k>0: E_k = L_k F_{k-1}⁻ᵀ, F_k = chol(D_k − E_k E_kᵀ) After it returns, block-row k’s MAIN slot holds the Cholesky factor F_k (lower triangle; its upper triangle keeps stale D_k entries) and its LEFT slot holds E_k. RIGHT slots are never read or written (for a symmetric system they mirror LEFT of the next row). The strips no longer hold A — keep a copy if you still need bdmv products with A.

Solve (bdsv_solve, in place on the padded vector): forward sweep y_k = F_k⁻¹ (b_k − E_k y_{k-1}), then backward sweep x_k = F_k⁻ᵀ (y_k − E_{k+1}ᵀ x_{k+1}).

Scratch: bdsv_scratch_bytes<T, BlockSize>() bytes (two dense BlockSize×BlockSize staging buffers — strip blocks are strided inside the row-major tile, so each knot’s blocks are staged contiguous via load_block/store_block).

Thread-count invariant; composes entirely from existing primitives, so it ends on their trailing barriers.

Functions

template<typename T, uint32_t BlockSize>
constexpr std::size_t bdsv_scratch_bytes()

Scratch size in bytes for bdsv / bdsv_factor / bdsv_solve.

Template Parameters:

• T – Scalar type.

• BlockSize – Block dimension.

Returns:

Bytes to allocate for s_scratch (two dense BlockSize² blocks).

template<typename T, uint32_t NumBlockRows, uint32_t BlockSize, bool CHECK = false>
void bdsv_factor(T *s_matrix, T *s_scratch, int *s_fail = nullptr)

In-place block-Cholesky factorization of a block-tridiagonal SPD matrix.

See the file docs for the factor layout left in the strips. With CHECK, a non-positive-definite pivot in any knot’s Cholesky sets *s_fail = 1 (the sweep still runs to completion; treat the factor as invalid).

Template Parameters:

• T – Scalar type (prefer double for long/ill-conditioned chains).

• NumBlockRows – Number of block-rows (knots).

• BlockSize – Block dimension.

• CHECK – Report a non-PD pivot via s_fail (default false, compiles out).

Parameters:

• s_matrix – In/out [L|D|R] strips; on return MAIN←F_k (lower), LEFT←E_k.

• s_scratch – Shared scratch of bdsv_scratch_bytes<T, BlockSize>() bytes.

• s_fail – Optional non-PD flag when CHECK (set to 1 on failure, else untouched).

template<typename T, uint32_t NumBlockRows, uint32_t BlockSize>
void bdsv_solve(const T *s_matrix, T *s_vector, T *s_scratch)

Solve A x = b from a bdsv_factored block-tridiagonal system.

Forward then backward block substitution. s_vector uses the padded (NumBlockRows+2)*BlockSize layout with b in the interior blocks (block-row br at offset (br+1)*BlockSize); on return the interior holds x. Reusable: factor once, call this per right-hand side.

Template Parameters:

• T – Scalar type.

• NumBlockRows – Number of block-rows (knots).

• BlockSize – Block dimension.

Parameters:

• s_matrix – Factored strips from bdsv_factor (read-only).

• s_vector – In/out padded right-hand side; interior overwritten with x.

• s_scratch – Shared scratch of bdsv_scratch_bytes<T, BlockSize>() bytes.

template<typename T, uint32_t NumBlockRows, uint32_t BlockSize, bool CHECK = false>
void bdsv(T *s_matrix, T *s_vector, T *s_scratch, int *s_fail = nullptr)

Direct block-tridiagonal SPD solve: factor + solve in one call (in place).

bdsv_factor then bdsv_solve. On return the strips hold the factor (MAIN←F_k lower, LEFT←E_k) and the padded vector’s interior holds x. The exact one-sweep alternative to glass::pcg on the same layout. NumPy equivalent: x = np.linalg.solve(A_dense, b) (A SPD block-tridiag).

Template Parameters:

• T – Scalar type.

• NumBlockRows – Number of block-rows (knots).

• BlockSize – Block dimension.

• CHECK – Report a non-PD pivot via s_fail (default false, compiles out).

Parameters:

• s_matrix – In/out [L|D|R] strips (see bdsv_factor).

• s_vector – In/out padded right-hand side; interior overwritten with x.

• s_scratch – Shared scratch of bdsv_scratch_bytes<T, BlockSize>() bytes.

• s_fail – Optional non-PD flag when CHECK.

Block accessors (store_block / load_block)

Strided read/write of a single BlockSize × BlockSize block into / out of a [L | D | R] strip (with an optional scale), used to assemble and unpack the block-tridiagonal strips. Block- and warp-scoped forms.

Dense d×d block ↔ block-tridiagonal [L|D|R] strip movers (store_block / load_block).

GLASS owns the block-tridiagonal [L|D|R] strip layout used by bdmv / pcg: each block-row is a d × (3*d) row-major tile, so the (y,x) entry of the slot-s block lives at strip[y*band_width + s*d + x] (band_width = 3*d). Consumers (e.g. GATO’s Schur assembly) repeatedly hand-roll the transpose / negate remap between a dense d×d block and this strip; these two functions are that remap, once, correctly.

Each (y,x) element is written by exactly one thread/lane, so there is no race and the ops are trivially thread-count invariant. Block + warp::.

TRANSPOSE=false moves the dense block as src[y*d + x]

(the “store it as laid

out” case);

TRANSPOSE=true moves src[x*d + y] (store/load the transpose — GATO’s phi_k vs phi_kᵀ). scale folds in a multiplier (-1 negates, GATO’s common case). store_block writes the strip from a dense block; load_block is the exact inverse (strip → dense). With matching TRANSPOSE/scale=1, load_block ∘ store_block is the identity.

Enums

enum class BandSlot : uint32_t

Which sub-block of a [L|D|R] block-row strip a mover targets.

Column offset within the row-major strip: LEFT at 0, MAIN at d, RIGHT at 2*d (matching bdmv’s [L|D|R] order).

Values:

enumerator LEFT

enumerator MAIN

enumerator RIGHT

Functions

template<typename T, uint32_t d, uint32_t band_width = 3 * d, bool TRANSPOSE = false, bool TRAILING_SYNC = true>
void store_block(T *dst_strip, BandSlot slot, const T *src, T scale = T(1))

Store a dense d×d block into a [L|D|R] strip slot: strip[y*band_width + slot*d + x] = scale * src[(TRANSPOSE? x*d+y : y*d+x)].

Single-warp store_block — one 32-lane warp strides the d*d block.

Warp form of store_block; TRAILING_SYNC gates a closing __syncwarp(). Each element written once. Full 32 lanes required.

See also

store_block

Template Parameters:

• T – Scalar type (e.g. float, double).

• d – Block dimension (the L/D/R blocks are d×d).

• band_width – Row length of the strip (default 3*d).

• TRANSPOSE – Store the transpose of the dense block (default false).

• TRAILING_SYNC – Emit a trailing __syncthreads() (default true).

Parameters:

• dst_strip – Start of the block-row strip (row-major, row length band_width).

• slot – Which sub-block (LEFT/MAIN/RIGHT) to write.

• src – Dense d×d source block (read-only).

• scale – Multiplier folded into the store (e.g. -1 to negate).

template<typename T, uint32_t d, uint32_t band_width = 3 * d, bool TRANSPOSE = false, bool TRAILING_SYNC = true>
void load_block(T *dst, const T *src_strip, BandSlot slot, T scale = T(1))

Load a dense d×d block from a [L|D|R] strip slot (inverse of store_block): dst[(TRANSPOSE? x*d+y : y*d+x)] = scale * strip[y*band_width + slot*d + x].

Single-warp load_block — one 32-lane warp strides the d*d block.

Warp form of load_block; TRAILING_SYNC gates a closing __syncwarp(). Full 32 lanes required.

See also

load_block

Template Parameters:

• T – Scalar type (e.g. float, double).

• d – Block dimension.

• band_width – Row length of the strip (default 3*d).

• TRANSPOSE – Load the transpose (default false).

• TRAILING_SYNC – Emit a trailing __syncthreads() (default true).

Parameters:

• dst – Dense d×d destination block (overwritten).

• src_strip – Start of the block-row strip (row-major, row length band_width).

• slot – Which sub-block (LEFT/MAIN/RIGHT) to read.

• scale – Multiplier folded into the load.

namespace warp