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 :doc:`../user_guide/concepts/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 :doc:`warp`.

gemm
----

.. doxygenfile:: src/base/L3/gemm.cuh

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
:doc:`../user_guide/concepts/contraction_parallel` for the honest win-condition.

.. doxygenfile:: src/base/L3/gemm_reduced.cuh

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

.. doxygenfile:: src/base/L3/tensor_contract.cuh

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

.. doxygenfile:: src/base/L3/congruence.cuh

Contraction-parallel syrk (``syrk_reduced``)
--------------------------------------------

.. doxygenfile:: src/base/L3/syrk_reduced.cuh

Riccati feedback gain (``riccati_gain``)
----------------------------------------

.. doxygenfile:: src/base/L3/riccati.cuh

Strided gemm
------------

.. doxygenfile:: src/base/L3/gemm_strided.cuh

Indexed / batched gemm
----------------------

.. doxygenfile:: src/base/L3/gemm_batched_indexed.cuh

Matrix inverse
--------------

.. doxygenfile:: src/base/L3/inv.cuh

Cholesky
--------

.. doxygenfile:: src/base/L3/potrf.cuh

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``.

.. doxygenfile:: src/base/L3/getrf.cuh

Triangular solve (trsm)
-----------------------

.. doxygenfile:: src/base/L3/trsm.cuh

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

.. doxygenfile:: src/base/L3/syrk.cuh

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``).

.. doxygenfile:: src/base/L3/symm.cuh

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.

.. doxygenfile:: src/base/L3/dimm.cuh

LDLᵀ factorization (ldlt / ldlt_solve)
--------------------------------------

.. doxygenfile:: src/base/L3/ldlt.cuh

SPD solve (posv / potrs)
------------------------

.. doxygenfile:: src/base/L3/posv.cuh

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.

.. doxygenfile:: src/base/L3/syev.cuh

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

.. doxygenfile:: src/cgrps/l3.cuh