Library

since C++26

Advanced

std::linalg

C++26 BLAS-level 1/2/3 free functions operating on std::mdspan views for portable dense linear algebra without external library dependencies.

std::linalgsince C++26

A namespace of BLAS-compatible free functions introduced in C++26 under <linalg> that perform dense vector, matrix–vector, and matrix–matrix operations directly on std::mdspan views, supporting row-major, column-major, and strided memory layouts without heap allocation inside the algorithm layer.

Overview

Before C++26, the standard library offered only scattered numerical primitives: std::inner_product, std::valarray, and the <numeric> algorithms. Production code needing true linear algebra reached for BLAS/LAPACK C bindings, Eigen, Blaze, or similar third-party libraries. C++26 addresses this gap with std::linalg: a standardised, namespace-scoped set of functions that mirrors the structure of BLAS Levels 1, 2, and 3.

The entire design rests on std::mdspan (C++23), which provides a non-owning, layout-parameterised view over contiguous or strided memory. Every std::linalg function takes mdspan arguments and dispatches work according to the view's layout and accessor policies. Ownership stays in caller-managed storage; layout policy lives in the mdspan type; algorithms live in linalg. This layering makes composition possible without heap allocation inside any algorithm.

The library ships four zero-copy view adaptors that wrap an existing mdspan to express mathematical structure:

Adaptor	Effect
`std::linalg::transposed(A)`	Swaps row and column extents/strides
`std::linalg::conjugated(A)`	Returns complex conjugates through the accessor
`std::linalg::conjugate_transposed(A)`	Hermitian transpose (†)
`std::linalg::scaled(α, A)`	Multiplies every element access by scalar `α`

Triangular and symmetric structure is communicated via tag-type arguments rather than a separate type hierarchy: std::linalg::upper_triangle / lower_triangle selects which triangle is meaningful; std::linalg::explicit_diagonal / implicit_unit_diagonal controls whether the diagonal is stored or assumed to be 1. This avoids allocating a dense copy just to use a triangular solver.

Syntax

All functions live in namespace std::linalg after #include <linalg>. Include <mdspan> (C++23) explicitly — the standard does not guarantee it is pulled in transitively.

BLAS Level 1 — vector/vector operations

cpp

// Dot product: result = Σ x[i] * y[i]  (real types)          // C++26
auto r = std::linalg::dot(x, y);

// Conjugate dot product: result = Σ conj(x[i]) * y[i]         // C++26
auto rc = std::linalg::dotc(x, y);

// Euclidean 2-norm, accumulated from init                      // C++26
auto n = std::linalg::vector_two_norm(x, 0.0);

// Overflow-safe sum of squares (returns {scale, ssq})         // C++26
auto [sc, ssq] = std::linalg::vector_sum_of_squares(x, {1.0, 0.0});

// In-place scale: x *= alpha                                   // C++26
std::linalg::scale(alpha, x);

// z = x + y                                                    // C++26
std::linalg::add(x, y, z);

// z = alpha*x + y  (compose scaled adaptor with add)          // C++26
std::linalg::add(std::linalg::scaled(alpha, x), y, z);

BLAS Level 2 — matrix/vector

cpp

std::linalg::matrix_vector_product(A, x, y);           // y = A*x     C++26
std::linalg::matrix_vector_product(A, x, z, y);        // y = A*x + z C++26

std::linalg::symmetric_matrix_vector_product(
    A, std::linalg::upper_triangle, x, y);              // y = A*x (symmetric, upper triangle) C++26

std::linalg::triangular_matrix_vector_product(
    A, std::linalg::upper_triangle,
    std::linalg::implicit_unit_diagonal, x, y);         // y = A*x (unit upper triangular) C++26

BLAS Level 3 — matrix/matrix

cpp

std::linalg::matrix_product(A, B, C);                  // C = A*B      C++26
std::linalg::matrix_product(A, B, E, C);               // C = A*B + E  C++26

Matrix norms

cpp

auto f = std::linalg::matrix_frob_norm(A, 0.0);        // Frobenius    C++26
auto o = std::linalg::matrix_one_norm(A, 0.0);         // max col-sum  C++26
auto i = std::linalg::matrix_inf_norm(A, 0.0);         // max row-sum  C++26

BLAS Level 2: Triangular Solve and Rank Updates

triangular_matrix_vector_solve solves the system Ax = b in-place for a triangular matrix A (BLAS TRSV). The triangle to use and whether the diagonal is stored are communicated via tag arguments, matching the same tag types used by the product functions.

cpp

// Solve A*x = b in place: x holds b on entry, solution on exit
// A is upper triangular with an explicitly stored diagonal        // C++26
std::linalg::triangular_matrix_vector_solve(
    A,
    std::linalg::upper_triangle,
    std::linalg::explicit_diagonal,
    b,          // right-hand side (in-vector)
    x);         // solution (out-vector)

// In-place overload: x holds b on entry, is overwritten in place  // C++26
std::linalg::triangular_matrix_vector_solve(
    A,
    std::linalg::lower_triangle,
    std::linalg::implicit_unit_diagonal,
    x);         // in-place: overwrites b with solution

The tag implicit_unit_diagonal tells the solver to treat diagonal entries as 1 without reading them — this is the convention used when storing an LU factorisation in-place where the unit diagonal of L is implicit.

symmetric_matrix_rank_k_update performs the BLAS SYRK operation C = C + α·A·A^T, updating only the triangle specified by the tag (BLAS 3, O(n²k) work):

cpp

// C = C + alpha * A * A^T  (update upper triangle of C)          // C++26
std::linalg::symmetric_matrix_rank_k_update(
    alpha,                           // Scalar
    A,                               // in-matrix  (m×k)
    C,                               // inout-matrix (m×m, symmetric)
    std::linalg::upper_triangle);    // which triangle to access

// Without an explicit alpha: C = C + A * A^T                     // C++26
std::linalg::symmetric_matrix_rank_k_update(
    A, C, std::linalg::upper_triangle);

Only the selected triangle of C is read and written; the opposite triangle is untouched. A Hermitian variant (hermitian_matrix_rank_k_update) exists for complex matrices, corresponding to BLAS HERK (C = C + α·A·A†).

Rotation and Transformation Utilities

Givens rotations are the building block of stable QR factorisation, SVD, and least-squares algorithms. std::linalg exposes both the setup step (computing the cosine/sine pair) and the application step (updating two vectors in place) as separate functions, matching BLAS DLARTG and DROT.

Step 1: compute the rotation parameters

cpp

#include <linalg>    // C++26
#include <complex>

// Real case: compute (c, s, r) such that [c  s]   [a]   [r]
//                                         [-s c] * [b] = [0]
auto res = std::linalg::setup_givens_rotation(a, b);
// res.c  — cosine (always real)
// res.s  — sine   (real for real inputs)
// res.r  — Euclidean norm of [a, b]

// Complex case: [c    s  ] [a]   [r]
//               [-conj(s) c] [b] = [0]
// c is always real even for complex a, b
std::complex<double> ca{1.0, 2.0}, cb{3.0, 4.0};
auto cres = std::linalg::setup_givens_rotation(ca, cb);
// cres.c is double, cres.s and cres.r are std::complex<double>

Step 2: apply the rotation to two vectors

cpp

// Apply rotation (c, s) to rows x and y in place               // C++26
// Each element pair (x[i], y[i]) is updated by the 2×2 rotation
std::linalg::apply_givens_rotation(x, y, res.c, res.s);

// ExecutionPolicy overload (may run in parallel)               // C++26
std::linalg::apply_givens_rotation(
    std::execution::par_unseq, x, y, res.c, res.s);

// Complex variant: s is std::complex<Real>
std::linalg::apply_givens_rotation(x, y, cres.c, cres.s);

setup_givens_rotation is noexcept and returns a givens_rotation_setup_result<Real> aggregate with members c, s, and r. The complex specialisation keeps c real while s and r are complex, matching the mathematical constraint c² + |s|² = 1.

Execution Policy Support

Most std::linalg functions that do not depend on sequential data access (such as apply_givens_rotation, scale, add, dot, matrix_vector_product, matrix_product, and symmetric_matrix_rank_k_update) provide ExecutionPolicy overloads. The policy is the first argument:

cpp

#include <linalg>
#include <execution>   // std::execution::par_unseq

// Parallel elementwise scale                                     // C++26
std::linalg::scale(std::execution::par_unseq, alpha, x);

// Parallel dot product (result type follows normal return rules)  // C++26
auto d = std::linalg::dot(std::execution::par_unseq, x, y, 0.0);

// Parallel general matrix-vector product                         // C++26
std::linalg::matrix_vector_product(
    std::execution::par_unseq, A, x, y);

// Parallel rank-k update                                         // C++26
std::linalg::symmetric_matrix_rank_k_update(
    std::execution::par_unseq, alpha, A, C, std::linalg::upper_triangle);

Not every function accepts an execution policy. In-place triangular solvers (triangular_matrix_vector_solve in-place overload) do have ExecutionPolicy overloads, but the solver is inherently sequential per column — providing a parallel policy permits an implementation to use parallelism internally (e.g., across independent right-hand sides) but does not guarantee speedup for a single solve. Always profile before assuming that par_unseq improves throughput on a given platform.

Examples

Dot product and 2-norm of a flat array

cpp

#include <linalg>    // C++26
#include <mdspan>    // C++23
#include <array>
#include <print>     // C++23

int main() {
    std::array<double, 4> data{1.0, 2.0, 3.0, 4.0};
    std::mdspan x(data.data(), 4);                             // C++23, rank-1

    auto d = std::linalg::dot(x, x);                          // 1+4+9+16 = 30
    auto n = std::linalg::vector_two_norm(x, 0.0);            // sqrt(30)

    std::println("dot={}, norm={:.6f}", d, n);
}

General matrix–vector product (row-major)

cpp

#include <linalg>
#include <mdspan>
#include <vector>

void gemv_row_major() {
    // Stored row-major: A = [[1,2,3],[4,5,6],[7,8,9]]
    std::vector<double> a_data{1,2,3, 4,5,6, 7,8,9};
    std::vector<double> x_data{1, 0, -1};
    std::vector<double> y_data(3, 0.0);

    // layout_right is the default — row-major, matching C array convention  // C++23
    std::mdspan<double, std::extents<std::size_t, 3, 3>> A(a_data.data());
    std::mdspan<double, std::extents<std::size_t, 3>>    x(x_data.data());
    std::mdspan<double, std::extents<std::size_t, 3>>    y(y_data.data());

    std::linalg::matrix_vector_product(A, x, y);               // C++26 — y = A*x
    // y == {-2.0, -2.0, -2.0}
}

Column-major layout for BLAS interoperability

Traditional BLAS and Fortran routines use column-major storage. Use std::layout_left to model this without copying:

cpp

#include <linalg>
#include <mdspan>
#include <vector>

void gemm_column_major() {
    // Column-major 2×2: a = {1,3,2,4} → [[1,2],[3,4]]
    std::vector<double> a{1,3,2,4};
    std::vector<double> b{5,7,6,8};
    std::vector<double> c(4, 0.0);

    using Ext = std::extents<std::size_t, 2, 2>;
    std::mdspan<double, Ext, std::layout_left> A(a.data()); // C++23
    std::mdspan<double, Ext, std::layout_left> B(b.data());
    std::mdspan<double, Ext, std::layout_left> C(c.data());

    std::linalg::matrix_product(A, B, C);                    // C++26 — C = A*B
}

Zero-copy A^T A using `transposed`

cpp

#include <linalg>
#include <mdspan>

template<class InMat, class OutMat>
void compute_gram_matrix(InMat A, OutMat ATA) {
    // A^T * A without allocating a transposed copy — C++26
    std::linalg::matrix_product(
        std::linalg::transposed(A),
        A,
        ATA
    );
}

Best Practices

Wrap data in mdspan at API boundaries and carry the type through. Building an mdspan from a raw pointer at every call site loses extent and layout information that enables bounds reasoning. Accept and return mdspan in internal interfaces.

Match storage layout to your access pattern. std::layout_right (row-major, the default) is cache-friendly for row-wise traversal. If downstream code expects Fortran or BLAS column-major conventions, create layout_left spans at the construction site rather than inserting a transpose step.

Prefer view adaptors over materialising temporaries. std::linalg::transposed(A) is a zero-cost type-level operation; allocating a transposed copy just to pass it to a product call wastes memory and bandwidth.

Use vector_sum_of_squares over a manual dot for norm computation on extreme-valued data. The two-argument vector_two_norm accumulates Σ x[i]² before taking the square root, which overflows for large elements. vector_sum_of_squares uses a scale/ssq pair that avoids intermediate overflow.

Initialise output spans before passing them to the three-argument accumulating forms. matrix_product(A, B, E, C) computes C = A*B + E, reading from E. A zero-initialised output fed as E is the correct pattern for a plain product that also accumulates into an existing matrix.

Common Pitfalls

Aliasing input and output is undefined behaviour. Unlike std::copy, which explicitly handles overlap, every std::linalg function requires that no output span alias any input span. There is no diagnostic. Passing the same buffer as both x and y in matrix_vector_product(A, x, y) silently corrupts results.

Using dot instead of dotc for complex inner products. For std::complex element types, dot(x, y) computes Σ x[i] * y[i] without conjugation. The standard sesquilinear inner product used in physics and signal processing is dotc(x, y) = Σ conj(x[i]) * y[i]. Using dot with complex vectors produces a bilinear result that is typically wrong in this context.

Misinterpreting extent order for column-major spans. std::mdspan<T, std::extents<size_t, M, N>, std::layout_left> always means M rows and N columns — the extent order is (rows, cols) regardless of storage order. BLAS C wrappers often use (N, M) conventions inherited from Fortran; transposing this mental model when constructing spans is a common source of silent shape errors.

Assuming all linalg functions are available in constant-evaluated contexts. The standard marks functions constexpr where the element type permits it, but compile-time linalg support in standard library implementations is still maturing in early C++26 toolchains. Verify before relying on it in constexpr or consteval contexts.

Matrix-vector operations

matrix_vector_product is the primary workhorse for BLAS Level 2 work. It computes y = A*x or the accumulating form y = A*x + z, and accepts any mdspan whose element type satisfies the usual arithmetic requirements. The triangular and symmetric variants narrow which elements of the matrix are read.

cpp

#include <linalg>   // C++26
#include <mdspan>
#include <vector>

void demo_matvec() {
    std::vector<double> a_data{2,1, 1,3};   // 2×2 stored row-major
    std::vector<double> x_data{1, 2};
    std::vector<double> y_data(2, 0.0);
    std::vector<double> z_data{10, 20};     // accumulation base

    using E2 = std::extents<std::size_t, 2, 2>;
    using E1 = std::extents<std::size_t, 2>;
    std::mdspan<double, E2> A(a_data.data());
    std::mdspan<double, E1> x(x_data.data());
    std::mdspan<double, E1> y(y_data.data());
    std::mdspan<double, E1> z(z_data.data());

    // y = A*x  → {4.0, 7.0}                                       // C++26
    std::linalg::matrix_vector_product(A, x, y);

    // y = A*x + z  → {14.0, 27.0}                                 // C++26
    std::linalg::matrix_vector_product(A, x, z, y);
}

The transposed adaptor lets you compute y = A^T * x without a copy: pass std::linalg::transposed(A) as the first argument. The resulting mdspan view swaps the extent and stride dimensions at zero runtime cost.

cpp

// y = A^T * x  — zero-copy transpose                              // C++26
std::linalg::matrix_vector_product(std::linalg::transposed(A), x, y);

BLAS-like interface

std::linalg is designed so that every function maps to a named BLAS routine. The table below shows the correspondence; knowing BLAS names makes it straightforward to translate legacy C-BLAS call sites.

`std::linalg` function	BLAS equivalent	Level
`dot` / `dotc`	DDOT / ZDOTC	1
`vector_two_norm`	DNRM2	1
`scale`	DSCAL	1
`add`	DAXPY (without α)	1
`matrix_vector_product`	DGEMV	2
`triangular_matrix_vector_product`	DTRMV	2
`triangular_matrix_vector_solve`	DTRSV	2
`symmetric_matrix_rank_k_update`	DSYRK	3
`matrix_product`	DGEMM	3

The naming drops the type-prefix convention (D/S/Z/C) because C++ templates already dispatch on element type. The scaled adaptor replaces the scalar alpha parameter found in BLAS signatures: instead of dgemv(alpha, A, x, beta, y), you write:

cpp

// Equivalent to BLAS: y = alpha * A * x + beta * y               // C++26
// Step 1: compute tmp = A * x
std::linalg::matrix_vector_product(A, x, tmp);
// Step 2: y = alpha * tmp + beta * y
std::linalg::add(
    std::linalg::scaled(alpha, tmp),
    std::linalg::scaled(beta,  y),
    y);

For the full-featured dgemv-style fused operation, implementations may provide it as a QoI extension; the standard expresses the same semantics via composition of adaptors and a three-argument matrix_vector_product.

In-place vs out-of-place operations

Most std::linalg functions have both an out-of-place form and an in-place overload. The in-place overload is identified by having one fewer argument for the output: the designated in/out span both reads the input value and is overwritten with the result.

cpp

// Out-of-place scale: reads x, writes result into y               // C++26
// (No such form exists — scale is always in-place)
std::linalg::scale(alpha, x);      // x *= alpha, no separate output

// Out-of-place triangular solve: reads b, writes solution to x    // C++26
std::linalg::triangular_matrix_vector_solve(
    A, std::linalg::upper_triangle, std::linalg::explicit_diagonal,
    b,    // in-vector (right-hand side, not modified)
    x);   // out-vector (solution)

// In-place triangular solve: x holds b on entry, solution on exit // C++26
std::linalg::triangular_matrix_vector_solve(
    A, std::linalg::upper_triangle, std::linalg::explicit_diagonal,
    x);   // in/out — b on entry, overwritten with solution

Aliasing is still prohibited even in the in-place forms. The in-place overload is syntactic shorthand, not a relaxation of the no-alias rule. The only operation permitted to read and write the same span is the in-place form where the standard explicitly defines that single-span argument as the inout operand.

Triangular and symmetric matrices

Triangular and symmetric structure is encoded via tag arguments, not separate matrix types. This avoids the need to pack a triangular matrix into a special class or to copy the full dense matrix into a dense temporary. The relevant tags live in namespace std::linalg.

cpp

#include <linalg>   // C++26

// Tags for triangle selection
constexpr auto upper = std::linalg::upper_triangle;
constexpr auto lower = std::linalg::lower_triangle;

// Tags for diagonal handling
constexpr auto expl = std::linalg::explicit_diagonal;      // read from storage
constexpr auto unit = std::linalg::implicit_unit_diagonal; // treat as 1, do not read

// Symmetric matrix-vector product (only reads the upper triangle)  // C++26
std::linalg::symmetric_matrix_vector_product(A, upper, x, y);

// Hermitian matrix-vector product (complex; reads upper triangle)  // C++26
std::linalg::hermitian_matrix_vector_product(A, upper, x, y);

// Triangular matrix-vector product: y = U*x, U upper unit triangular // C++26
std::linalg::triangular_matrix_vector_product(A, upper, unit, x, y);

// Rank-2 symmetric update: C = C + alpha*(x*y^T + y*x^T)          // C++26
std::linalg::symmetric_matrix_rank_2_update(alpha, x, y, C, upper);

Elements in the unselected triangle are never accessed. This means they can hold arbitrary values (e.g., the opposite factor of an LU decomposition) without affecting correctness. Pass the same tag consistently to all calls that operate on the same logical matrix.

Conjugate transpose

The conjugate_transposed adaptor composes conjugation and transposition in a single zero-copy wrapper. For a matrix A of std::complex<double> elements, std::linalg::conjugate_transposed(A) presents a view with swapped extents where each element access returns the complex conjugate. This is the Hermitian adjoint (†) required for inner products and unitary transformations in quantum mechanics, signal processing, and other complex-valued domains.

cpp

#include <linalg>
#include <mdspan>
#include <complex>
#include <vector>   // C++26

void hermitian_demo() {
    using cdbl = std::complex<double>;
    std::vector<cdbl> a_data{
        {1,0}, {2,1},
        {2,-1},{3,0}   // Hermitian 2×2: A† = A
    };
    std::vector<cdbl> x_data{{1,0},{0,1}};
    std::vector<cdbl> y_data(2, cdbl{});

    using E2 = std::extents<std::size_t, 2, 2>;
    using E1 = std::extents<std::size_t, 2>;
    std::mdspan<cdbl, E2> A(a_data.data());
    std::mdspan<cdbl, E1> x(x_data.data());
    std::mdspan<cdbl, E1> y(y_data.data());

    // y = A† * x  via conjugate_transposed — no allocation         // C++26
    std::linalg::matrix_vector_product(
        std::linalg::conjugate_transposed(A), x, y);
}

The conjugated adaptor (without transposition) is useful when you need element-wise conjugation of a vector passed to dotc-style computations manually. The two adaptors can be chained: std::linalg::conjugated(std::linalg::transposed(A)) is equivalent to std::linalg::conjugate_transposed(A).

Interop with Eigen and BLAS libraries

std::linalg and Eigen can share memory — no copy is required. Eigen's Map and Ref types address the same raw pointer as mdspan, so you can hand the same allocation to both APIs within the same translation unit. Use std::layout_left (column-major) for Eigen's default storage order.

cpp

#include <linalg>   // C++26
#include <mdspan>
#include <Eigen/Dense>

void eigen_interop_demo() {
    // Eigen allocates the storage
    Eigen::MatrixXd E = Eigen::MatrixXd::Random(4, 4);

    // Wrap with mdspan (column-major, matching Eigen::ColMajor)    // C++26
    using Ext = std::extents<std::size_t, std::dynamic_extent, std::dynamic_extent>;
    std::mdspan<double, Ext, std::layout_left> A(
        E.data(), E.rows(), E.cols());

    // Use std::linalg on the same memory — no copy
    auto frob = std::linalg::matrix_frob_norm(A, 0.0);             // C++26

    // Hand back to Eigen for operations not yet in std::linalg
    Eigen::MatrixXd result = E.llt().matrixL();  // Cholesky
}

For C BLAS (cblas) interoperability, construct an mdspan over the same pointer and use the same layout flag passed to the BLAS call. When calling cblas_dgemm with CblasColMajor, mirror that with std::layout_left; with CblasRowMajor use std::layout_right. Mixing layout flags between the BLAS call and the mdspan view produces silently wrong results.

cpp

#include <linalg>
#include <mdspan>
#include <cblas.h>   // platform BLAS header

void cblas_roundtrip(double* ptr, int m, int n) {
    // C BLAS column-major product
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans,
                m, n, n, 1.0, ptr, m, ptr, n, 0.0, ptr, m);

    // Inspect result via std::linalg (same pointer, same layout)   // C++26
    using Ext = std::extents<std::size_t, std::dynamic_extent, std::dynamic_extent>;
    std::mdspan<double, Ext, std::layout_left> A(ptr, m, n);
    auto inf = std::linalg::matrix_inf_norm(A, 0.0);
}

Note that std::linalg does not call an external BLAS library internally — it is a standalone implementation in the C++ standard library. An implementation is free to forward to an optimised BLAS as a quality-of-implementation decision, but that is not guaranteed by the standard. For maximum throughput on hardware-tuned BLAS (e.g., MKL, OpenBLAS), call the vendor BLAS directly and use std::linalg for composition, norms, and auxiliary operations.

See Also

std::mdspan — the rank-N non-owning view that all std::linalg functions consume
std::span — 1-D contiguous view; construct a rank-1 mdspan from it to use with linalg
std::inner_product — pre-C++26 dot-product and general binary-fold over two ranges in <numeric>
std::valarray — older numerical array type with element-wise operators; no mdspan integration