Skip to content

linear-algebra-library.md

Latest commit

 

History

History
262 lines (199 loc) · 13.1 KB

linear-algebra-library.md

File metadata and controls

262 lines (199 loc) · 13.1 KB
Raw

Linear Algebra Library (<linalg>)

Executive Abstract

The C++26 Linear Algebra Library (<linalg>, std::linalg) is a standardized, high-performance linear algebra framework (proposal P1673R13 "A C++ Linear Algebra Library Based on BLAS", authored by Mark Hoemmen, Christian Trott, et al.). It standardizes Basic Linear Algebra Subprograms (BLAS Levels 1, 2, and 3) directly in ISO C++. Operating entirely over non-owning C++23 std::mdspan multidimensional array views, <linalg> replaces non-standard C/Fortran BLAS bindings (Intel MKL, Apple Accelerate, OpenBLAS) and hand-rolled matrix loops with a type-safe, generic C++ interface that enables hardware-accelerated SIMD and GPU vendor optimizations.

Motivation & Historical Context

Pre-C++26 Deficiencies: The Non-Standard Linear Algebra Landscape

Despite C++ being the primary language for high-performance computing (HPC), robotics, game physics, computer vision, and machine learning, C++ lacked a standard linear algebra library for over 40 years.

Prior to C++26, developers executing vector and matrix operations relied on third-party frameworks (Eigen, Armadillo, Boost.Numeric.Ublas) or raw C/Fortran BLAS C-APIs:

// Pre-C++26 Raw C BLAS Matrix-Matrix Multiply (cblas_dgemm) Anti-Pattern
#include <cblas.h> // Non-standard external library dependency

double A[6] = {1, 2, 3, 4, 5, 6}; // 2x3 matrix
double B[6] = {1, 2, 3, 4, 5, 6}; // 3x2 matrix
double C[4] = {0};                // 2x2 output

// Fragile C-API Call: Manually specifies dimensions, leading dimension strides (LDA, LDB),
// and raw memory pointers with ZERO C++ type safety!
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
            2, 2, 3, 1.0, A, 3, B, 2, 0.0, C, 2);

Deficiencies of Pre-C++26 Linear Algebra Solutions:

  1. Third-Party Library Lock-In: Codebases depended on non-standard libraries (Eigen, Armadillo), creating build-system friction and incompatible matrix type representations across different modules.
  2. Fragile C/Fortran BLAS Wrappers: Invoking hardware-accelerated vendor BLAS routines (Intel MKL, Apple Accelerate cblas_dgemm, cuBLAS) required passing raw pointers, explicit integer dimension counts, and manual stride leading dimensions (LDA, LDB). A single miscalculated dimension integer triggered memory corruption or invalid matrix results.
  3. Pessimized Hand-Rolled Loops: Writing naive triple-nested loops (for i... for j... for k...) failed to trigger optimal SIMD register vectorization, suffered from CPU cache line misses, and failed to exploit hardware tensor cores.

Solution: C++26 introduced <linalg>, standardizing BLAS Levels 1, 2, and 3 over C++23 std::mdspan matrix views.

Technical Mechanics & Library Architecture

1. The Three BLAS Operational Tiers (std::linalg)

                           std::linalg Framework (P1673R13)
                                          │
     ┌────────────────────────────────────┼────────────────────────────────────┐
     ▼                                    ▼                                    ▼
1. BLAS LEVEL 1                      2. BLAS LEVEL 2                      3. BLAS LEVEL 3
Vector-Vector Operations             Matrix-Vector Operations             Matrix-Matrix Operations
• std::linalg::dot(v1, v2)           • matrix_vector_product(A, x, y)     • matrix_product(A, B, C)
• std::linalg::vector_norm2(v)       • symmetric_matrix_vector_product    • GEMM (C = A * B + C)
• std::linalg::add(v1, v2, dest)     • triangular_matrix_vector_solve     • triangular_matrix_solve

Detailed BLAS API Functional Table

BLAS Level C++26 Standard Function API Mathematical Expression Time Complexity
Level 1 std::linalg::dot(v1, v2) Vector Dot Product: $r = \sum v_{1,i} v_{2,i}$ $\mathcal{O}(N)$
std::linalg::vector_norm2(v) Euclidean Norm: $|v|_2 = \sqrt{\sum v_i^2}$ $\mathcal{O}(N)$
std::linalg::add(v1, v2, out) Vector Addition: $\text{out} \gets v_1 + v_2$ $\mathcal{O}(N)$
Level 2 std::linalg::matrix_vector_product(A, x, y) Matrix-Vector Product: $y \gets A x$ or $y \gets A x + y$ $\mathcal{O}(M \times N)$
std::linalg::symmetric_matrix_vector_product Symmetric Matrix-Vector Multiply $\mathcal{O}(N^2)$
Level 3 std::linalg::matrix_product(A, B, C) General Matrix Multiply (GEMM): $C \gets A B$ $\mathcal{O}(N^3)$
std::linalg::triangular_matrix_product Triangular Matrix Multiply $\mathcal{O}(N^3)$

2. Integration with std::mdspan Views

std::linalg functions operate strictly over non-owning std::mdspan views:

$$\mathbf{\text{std::linalg::matrix_product( }} \underbrace{\text{std::mdspan } A}_{\text{Input Matrix } M \times K} , , \quad \underbrace{\text{std::mdspan } B}_{\text{Input Matrix } K \times N} , , \quad \underbrace{\text{std::mdspan } C}_{\text{Output Matrix } M \times N} \mathbf{\text{ )}}$$

  • Layout Invariance: std::linalg handles Row-Major (std::layout_right) and Column-Major (std::layout_left) layouts transparently without requiring explicit matrix transposition passes in user code!

Evolution Across Standards (C++26 to C++29)

┌──────────────────────────────────────────────────────────────────┐
│ C++23                                                            │
├──────────────────────────────────────────────────────────────────┤
│ • Introduced std::mdspan in <mdspan> (P0009R18)                  │
│ • Multidimensional subscripting operator[r, c]                   │
└──────────────────────────────────────────────────────────────────┘
                                 │
                                 ▼
┌──────────────────────────────────────────────────────────────────┐
│ C++26                                                            │
├──────────────────────────────────────────────────────────────────┤
│ • Introduced <linalg> header and std::linalg namespace (P1673R13)│
│ • Standardized BLAS Levels 1, 2, and 3                           │
│ • Full constexpr support for compile-time matrix operations      │
│ • Vendor hardware acceleration mapping (MKL, Accelerate, cuBLAS) │
└──────────────────────────────────────────────────────────────────┘
                                 │
                                 ▼
┌──────────────────────────────────────────────────────────────────┐
│ C++29 (Candidate)                                                │
├──────────────────────────────────────────────────────────────────┤
│ • High-level LAPACK matrix factorizations (LU, QR, Cholesky, SVD)│
└──────────────────────────────────────────────────────────────────┘

Comprehensive Code Examples

Example 1: Matrix-Vector Product (std::linalg::matrix_vector_product)

Demonstrates multiplying a $2 \times 2$ matrix by a 2-element vector using std::mdspan views and C++26 std::linalg.

#include <iostream>
#include <vector>
#include <array>

// C++26 <linalg> and C++23 <mdspan> Headers
#if defined(__cpp_lib_linalg) && __cpp_lib_linalg >= 202311L
#include <mdspan>
#include <linalg>

int main() {
    // 1. Storage buffers
    std::array<double, 4> A_data = {1.0, 2.0, 3.0, 4.0}; // 2x2 Matrix
    std::array<double, 2> x_data = {5.0, 6.0};           // 2x1 Vector
    std::array<double, 2> y_data = {0.0, 0.0};           // Output Vector

    // 2. Wrap buffers in non-owning std::mdspan views
    std::mdspan A(A_data.data(), 2, 2); // 2D Matrix view
    std::mdspan x(x_data.data(), 2);    // 1D Vector view
    std::mdspan y(y_data.data(), 2);    // 1D Output view

    // 3. C++26 Matrix-Vector Product: y = A * x
    // Row 0: 1.0*5.0 + 2.0*6.0 = 17.0
    // Row 1: 3.0*5.0 + 4.0*6.0 = 39.0
    std::linalg::matrix_vector_product(A, x, y);

    std::cout << "Matrix-Vector Product Result:\n";
    std::cout << "  y[0] = " << y[0] << " (Expected: 17.0)\n";
    std::cout << "  y[1] = " << y[1] << " (Expected: 39.0)\n";

    return 0;
}
#else
int main() {
    std::cout << "C++26 <linalg> framework demonstrated.\n";
    return 0;
}
#endif

Example 2: General Matrix Multiplication (GEMM: std::linalg::matrix_product)

Demonstrates multiplying two $2 \times 2$ matrices using std::linalg::matrix_product.

#include <iostream>
#include <vector>

#if defined(__cpp_lib_linalg) && __cpp_lib_linalg >= 202311L
#include <mdspan>
#include <linalg>

int main() {
    // Matrix A (2x2): [[1, 2], [3, 4]]
    std::vector<double> a_buffer = {1.0, 2.0, 3.0, 4.0};
    // Matrix B (2x2): [[5, 6], [7, 8]]
    std::vector<double> b_buffer = {5.0, 6.0, 7.0, 8.0};
    // Matrix C (2x2 Output): [[0, 0], [0, 0]]
    std::vector<double> c_buffer(4, 0.0);

    std::mdspan A(a_buffer.data(), 2, 2);
    std::mdspan B(b_buffer.data(), 2, 2);
    std::mdspan C(c_buffer.data(), 2, 2);

    // C = A * B
    // C[0,0] = 1*5 + 2*7 = 19
    // C[0,1] = 1*6 + 2*8 = 22
    // C[1,0] = 3*5 + 4*7 = 43
    // C[1,1] = 3*6 + 4*8 = 50
    std::linalg::matrix_product(A, B, C);

    std::cout << "Matrix Multiply C = A * B:\n";
    std::cout << "  [" << C[0, 0] << ", " << C[0, 1] << "]\n";
    std::cout << "  [" << C[1, 0] << ", " << C[1, 1] << "]\n";

    return 0;
}
#else
int main() { return 0; }
#endif

Example 3: Vector Dot Product & Euclidean Norm (std::linalg::dot and vector_norm2)

Demonstrates Level 1 BLAS vector operations.

#include <iostream>
#include <vector>

#if defined(__cpp_lib_linalg) && __cpp_lib_linalg >= 202311L
#include <mdspan>
#include <linalg>

int main() {
    std::vector<double> vec1 = {3.0, 4.0};
    std::vector<double> vec2 = {1.0, 2.0};

    std::mdspan v1(vec1.data(), 2);
    std::mdspan v2(vec2.data(), 2);

    // 1. Vector Dot Product
    double dot_product = std::linalg::dot(v1, v2); // 3*1 + 4*2 = 11.0

    // 2. Vector L2 Norm (Euclidean Length)
    double norm = std::linalg::vector_norm2(v1); // sqrt(3^2 + 4^2) = 5.0

    std::cout << "Dot Product:    " << dot_product << " (Expected: 11.0)\n";
    std::cout << "Euclidean Norm: " << norm << " (Expected: 5.0)\n";

    return 0;
}
#else
int main() { return 0; }
#endif

Engineering Best Practices & Common Pitfalls

  1. Verify Inner Dimension Match Invariants:

    • In matrix_product(A, B, C), matrix $A$ (size $M \times K$) and matrix $B$ (size $K \times N$) must have matching inner dimension $K$ ($A.\text{extent}(1) == B.\text{extent}(0)$). Mismatched inner dimensions cause contract precondition failures or out-of-bounds memory accesses.

  2. Pass Non-Owning std::mdspan Views By Value:

    • std::mdspan is a lightweight view handle. Pass std::mdspan instances by value into std::linalg functions.

  3. Compiler Vendor Hardware Acceleration:

    • When compiling production HPC binaries, link against vendor hardware acceleration libraries (e.g. Intel MKL, Apple Accelerate, or OpenBLAS). Compilers map standard std::linalg::matrix_product calls directly to vendor hardware assembly without code changes.

Relationships

  • Related to:
  • Superseded by: None.
  • Supersedes: Raw C/Fortran BLAS function wrappers (cblas_dgemm), non-standard third-party matrix libraries, and hand-rolled triple-nested multiplication loops.

Back to C++26 section.