Linear Algebra Library

A numerical linear algebra library in modern C++. The project is designed to demonstrate ownership-aware C++ design, numerical reasoning, benchmarking, profiling, cache-aware optimization, SIMD, and concurrency through incremental, measured development.

The project currently has its owning Matrix and Vector foundations plus the first Phase 2 operations: arithmetic, broadcasting, transpose, trace, and vector/matrix norms. Performance experiments compare alternative traversal and access strategies before changing public implementations. More advanced numerical algorithms and optimized matrix multiplication kernels are intentionally being added incrementally.

Implemented functionality

Matrix

• contiguous row-major double storage;
• default, dimension, fill, flat-sequence, and nested initializer-list construction;
• copy and move ownership semantics;
• rows(), cols(), size(), empty(), and contiguous data() access;
• unchecked operator()(row, col) and checked at(row, col) access;
• in-place matrix addition and subtraction;
• in-place scalar addition, subtraction, and multiplication;
• value-returning addition, subtraction, and scalar multiplication;
• matrix-vector and reference matrix-matrix multiplication;
• adding a vector to every row or every column;
• constructing a matrix by repeating a vector as rows or columns;
• transpose and trace;
• one, Frobenius, and infinity norms selected with MatrixNorm.

Vector

• contiguous owning double storage;
• default, size, fill, initializer-list, and std::span construction;
• copy and move ownership semantics;
• size(), empty(), and contiguous data() access;
• unchecked operator[] and checked at() access;
• in-place vector addition and subtraction;
• in-place scalar addition, subtraction, and multiplication;
• value-returning addition, subtraction, and scalar multiplication;
• L2 normalization with validation for zero and non-finite vectors;
• L1, L2, and infinity norms selected with VectorNorm.

Basic usage

#include <linalg/linalg.hpp>

int main() {
    linalg::Matrix matrix{{1.0, 2.0, 3.0},
                          {4.0, 5.0, 6.0}};
    const linalg::Vector row_values{10.0, 20.0, 30.0};

    const linalg::Matrix adjusted =
        linalg::add_to_each_row(matrix, row_values);
    const linalg::Matrix transposed = linalg::transpose(matrix);

    const double matrix_size =
        linalg::norm(matrix, linalg::MatrixNorm::frobenius);

    linalg::Vector vector{1.0, 2.0, 3.0};
    const linalg::Vector other{4.0, 5.0, 6.0};

    vector += other;
    vector -= 1.0;
    vector *= 2.0;

    const double vector_length =
        linalg::norm(vector, linalg::VectorNorm::l2);

    return adjusted.at(0, 0) == 11.0 &&
                   transposed.at(0, 1) == 4.0 &&
                   matrix_size > 0.0 && vector_length > 0.0
               ? 0
               : 1;
}

For Vector, braces contain values while parentheses specify a size:

linalg::Vector one_value{3.0}; // one element containing 3.0
linalg::Vector three_zeroes(3); // three zero-initialized elements

Vector addition and subtraction require equal sizes. Matrix addition and subtraction require equal shapes. Broadcast operations require a vector length matching the relevant matrix dimension. Trace requires a square, non-empty matrix. Invalid dimensions throw std::invalid_argument; checked element access throws std::out_of_range.

Requirements

• CMake 3.20 or newer
• A C++20 compiler such as GCC or Clang
• Git and network access when GoogleTest or Google Benchmark are not installed locally

GoogleTest 1.17.0 and Google Benchmark 1.9.5 are development-only dependencies. CMake first looks for installed packages and otherwise fetches the pinned releases.

Build and test

cmake --preset debug
cmake --build --preset debug
ctest --preset debug --output-on-failure

Sanitizers

cmake --preset sanitized
cmake --build --preset sanitized
ctest --preset sanitized --output-on-failure

Release build and benchmarks

cmake --preset release
cmake --build --preset release
ctest --preset release --output-on-failure
./build/release/benchmarks/linalg_benchmarks

List registered benchmarks before choosing a filter:

./build/release/benchmarks/linalg_benchmarks --benchmark_list_tests

Run a quick comparison of one operation and size:

./build/release/benchmarks/linalg_benchmarks \
    --benchmark_filter='^BM_Transpose.*Kernel/512/512$' \
    --benchmark_min_time=0.1s \
    --benchmark_repetitions=7 \
    --benchmark_enable_random_interleaving=true \
    --benchmark_report_aggregates_only=true

Run all six matmul loop orders for one shape. Arguments are M/N/K, representing an M-by-K matrix multiplied by a K-by-N matrix:

./build/release/benchmarks/linalg_benchmarks \
    --benchmark_filter='^BM_MatMul.*Kernel/256/256/256$' \
    --benchmark_min_time=0.1s \
    --benchmark_repetitions=7 \
    --benchmark_enable_random_interleaving=true \
    --benchmark_report_aggregates_only=true

Plotting benchmark results

The optional Matplotlib tool under tools/plot_visualizer/ turns Google Benchmark JSON into SVG or PNG figures without adding a dependency to the C++ library. A focused matmul plot can be generated with:

python3 -m venv .venv
source .venv/bin/activate
python3 -m pip install -r tools/plot_visualizer/requirements.txt
python3 -m tools.plot_visualizer benchmark-results/matmul.json \
    --shape=512/512/512 \
    --metric='FLOP/s' \
    --output=docs/images/matmul-512.svg

See the plot visualizer guide for JSON capture, installation, CLI options, reusable Python API, plot selection, Markdown embedding, and troubleshooting.

Benchmarks are separated by operation family in benchmarks/, while one executable and Google Benchmark's regular-expression filter provide a common runner. Current experiments include broadcasting, transpose, trace, and all six scalar matmul loop orders. Kernel-only timings use a preallocated result, while end-to-end timings also include result allocation.

The main findings so far are:

• transpose performance is dominated by traversal order and contiguous destination writes;
• after inlining, indexed and pointer implementations with the same traversal compile and perform similarly;
• manually replacing indexed access with pointers did not produce a consistent improvement for add_to_each_row;
• direct-data and operator() trace implementations are currently indistinguishable within the observed measurement noise;
• row-major matmul loop order changes performance by nearly 3x at 64x64x64 and more than 12x at 512x512x512; kernels with contiguous inner-loop writes form the consistently fastest group.

Initial matmul loop-order results

The six scalar kernels perform the same arithmetic in different loop orders. For square matrices, the ikj and kij variants keep j innermost, making reads from B and updates to C contiguous. They remain near 5-7 GFLOP/s across the measured sizes. At 512x512x512, observed median CPU time was 50.4 ms for ikj versus 622.7 ms for the slowest jki variant, a 12.35x difference.

The exact ranking between ikj and kij is not yet decisive: their differences are often similar to the run-to-run variability. Rectangular matrices also demonstrate that the best outer-loop order depends on shape even when the inner loop is identical. See the full experiment for the second figure, measured tables, methodology, caveats, and cache-reuse hypothesis.

See docs/performance.md for commands, environment details, measured tables, limitations, and interpretation. Published performance claims should always be based on repeatable Release measurements on documented hardware.

Repository layout

include/linalg/   Public headers
src/              Library implementation
tests/            GoogleTest correctness tests
benchmarks/       Google Benchmark experiments
examples/         Small API usage programs
cmake/            Project CMake modules
docs/             Architecture and measured performance notes
tools/            Optional developer and documentation utilities