Skip to content
BranchesTags

Repository files navigation

Legolas++: Building Blocks for Linear Algebra Solvers

Legolas++ Tensor SIMD Engine Logo

Header-Only Zero Dependencies C++ Standard Documentation Platforms Build & Test License: MIT

High-Performance Modern C++ Tensor Engine for Automatic SIMD Vectorization of Recurrences via Data Layout Interleaving (DLI).

🚀 1 Single Scalar Code $\longrightarrow$ A Multitude of Hardware Targets

Write your core numerical algorithm once using standard, readable scalar C++ math. Legolas++ automatically maps and executes it at peak hardware efficiency across SIMD vector execution units, multi-core CPU threads, and dedicated GPUs:

Write Once in Scalar C++, Target Any Hardware (CPU SIMD, Multi-Core, GPU Vulkan/Metal)

  • Single Natural Scalar Code: Written with standard loops and scalar arithmetic ($Y_i = A \cdot X_i + B \cdot Y_{i-1}$). Zero assembly intrinsics, zero compiler pragmas, zero vendor lock-in.
  • Legolas++ DLI Engine: Transforms your problem layout in memory via Data Layout Interleaving (DLI), turning impossible loop-carried recurrences into contiguous vector streams.
  • A Multitude of Hardware Targets:
    • ⚡ CPU SIMD Vectorized: Full vector width utilization on ARM NEON (4 floats / 128-bit), x86 AVX2 (8 floats / 256-bit), and x86 AVX-512 (16 floats / 512-bit) at 100% hardware line-rate.
    • 🚀 Multi-Core Parallel: Built-in, zero-dependency StaticThreadPool with dynamic chunk scheduling scaling across 4 to 128+ CPU cores, with vector SIMD active on every thread.
    • 🔥 GPU Acceleration: Optional header-only Vulkan Compute backend for Linux and Windows, and native Apple Metal backend for macOS with unified memory.

🎯 When to Use Legolas++: The Core Problem

The Universal Scenario

You need to apply an intrinsically sequential algorithm (recurrence relation, recursive filter, time-stepping scheme, tridiagonal solver $y_n = f(y_{n-1}, x_n)$) to a massive batch of independent problem instances of the same size (thousands of 1D grid lines, dozens of audio tracks, video streams, or neural network channels).

The Classic Dilemma

  • Standard Multi-Threading (OpenMP / Threads): Parallelizes across CPU cores, but within each core, compilers cannot vectorize across sequential dependencies. Hardware SIMD execution units (AVX2, AVX-512, NEON) sit idle—wasting 75% to 93% of the CPU's theoretical compute capacity.
  • Manual SIMD Intrinsics (_mm256_..., NEON): Attempting to vectorize manually across instances requires writing hundreds of lines of assembly-like intrinsics. The code becomes unreadable, non-portable, and a maintenance nightmare.

The Legolas++ Solution

  1. 100% Machine Utilization: Fully utilizes all CPU cores and 100% of SIMD vector register widths simultaneously.
  2. Natural Scalar Notation: You write your core algorithm once, as a simple sequential loop in standard scalar math.
  3. Hardware Vectorization by Construction: Through Data Layout Interleaving (DLI), vectorization is structural in memory and guaranteed—no reliance on fragile compiler heuristics.
  4. Zero-Overhead Portability: Pure header-only C++14 running with peak efficiency across Apple Silicon (NEON), Linux (x86_64 AVX2 / AVX-512), and Windows MSVC.

💡 Why SIMD is the Core Innovation of Legolas++

While multi-threading (multi-core thread scheduling) is a widely available and commoditized feature in modern computing, hardware SIMD vectorization across linear recurrences has historically remained an intractable barrier.

Standard optimizing compilers (Clang, GCC, MSVC, Intel oneAPI) systematically fail to auto-vectorize loops with loop-carried dependencies (e.g. tridiagonal solvers, recursive IIR filters, ADI sweeps).

Legolas++ breaks this recurrence barrier at the CPU register level. Using Data Layout Interleaving (DLI), Legolas++ reorganizes memory across problem instances so that the exact same generic scalar loop maps directly to full-width hardware vector registers (ARM NEON, x86 AVX2 / AVX-512) without writing a single line of intrinsics or inline assembly.

1. The Recurrence Barrier: Why Compilers Give Up

Many essential algorithms in scientific simulation, digital signal processing, edge AI, and deep learning feature strict loop-carried data dependencies where iteration $i$ requires the output of iteration $i-1$:

// Tridiagonal elimination (Thomas forward sweep), recursive IIR filters, Gauss-Seidel:
for (int i = 1; i < N; ++i) {
    X[i] = (B[i] - L[i] * X[i-1]) * invD[i]; // Strict recurrence: RAW hazard!
}

Because of this sequential dependency chain: $$X[0] \longrightarrow X[1] \longrightarrow X[2] \longrightarrow X[3] \longrightarrow \dots \longrightarrow X[N-1]$$

Every optimizing compiler falls back to scalar execution, leaving up to 90% of the CPU's vector compute capacity completely idle.

2. The Breakthrough: SIMD via Data Layout Interleaving (DLI)

In real-world applications, engineers rarely solve a single isolated recurrence. Instead, they process ensembles of independent problem instances:

  • Scientific Computing & PDEs: Thousands of 1D tridiagonal systems across 2D/3D ADI grids, heat diffusion, or fluid flow.
  • Real-Time Audio DSP: Filtering 32, 64, or 128 audio channels concurrently with recursive IIR/Biquad filters.
  • Edge AI & Computer Vision: Depthwise Separable Convolutions across channels (MobileNet, ConvNeXt).

Instead of struggling to vectorize sequentially along $i$, Legolas++ interleaves $P$ independent problem instances directly in memory:

Animated Data Layout Interleaving (DLI) Mechanism

Figure: Data Layout Interleaving (DLI) in Action. Elements at step $i$ across $P=4$ independent problem instances are mapped contiguously into physical memory, transforming strided access into single-instruction aligned SIMD streaming.

🌐 An Ubiquitous Pattern Across Science & Industry

This computing pattern appears everywhere across high-performance engineering. Legolas++ includes dedicated, self-contained examples and tutorials for each:

Domain Intrinsically Sequential Kernel Batch Dimension (Interleaved) Live Example & Tutorial
Scientific Computing & PDEs Tridiagonal Gaussian elimination (Thomas algorithm, ADI sweeps) Thousands of 1D spatial lines in 2D/3D grids 📄 Tridiagonal Thomas Tutorial (code)
Real-Time Audio Processing Recursive IIR Biquad filter (sample $t$ depends on $t-1$, $t-2$) 64+ concurrent audio channels / DAW mixer tracks 🎧 Audio Biquad Showcase (code)
Computer Vision & Video Temporal recursive motion differencing / spatial convolution 32 concurrent 720p HD live camera streams 🎥 Video Pipeline Showcase (code)
Deep Learning & Edge AI MobileNet depthwise 2D convolutions 32–512 feature map channels ⚡ Depthwise Conv Showcase (code)
Quantitative Finance 1D finite-difference PDE time-stepping (Black-Scholes / Dupire) Tens of thousands of independent option contracts 📈 Option Pricing Engine

⚡ Examples & Showcases

30-Second Quick Example

Write your numerical solver once using natural scalar notation; Legolas++ auto-vectorizes it in hardware SIMD and scales across CPU cores simultaneously:

#include <iostream>
#include "Legolas/Array/Array.hxx"
#include "Legolas/Array/Map.hxx"

// 1. Write numerical algorithm once using standard scalar math:
struct Scaler {
    template <class A2D>
    void operator()(int begin, int end, A2D X, A2D Y) const {
        using Scalar = typename A2D::RealType;
        Scalar factor(2.5f);
        for (int j = begin; j < end; ++j) {
            for (int i = 0; i < X[j].size(); ++i) {
                Y[j][i] = factor * X[j][i]; // Compiles to NEON / AVX2 FMA automatically!
            }
        }
    }
};

int main() {
    // 2. Interleave 1024 systems of size 256 with pack size P=4:
    using Tensor2D = Legolas::Array<float, 2, 4, 2>;
    Tensor2D X(1024, 256), Y(1024, 256);
    X.fill(1.0f);

    // 3. Parallelize across cores AND vectorize in hardware SIMD simultaneously:
    Legolas::parmap(Scaler(), X, Y);

    std::cout << "Computed 1024 systems at full hardware vector line-rate!" << std::endl;
    return 0;
}

🏆 Industrial Showcases & Energy Efficiency (Green Computing)

High-performance computing is fundamentally an energy efficiency challenge. On mobile platforms, drones, and data centers, unvectorized scalar loops waste energy by leaving SIMD hardware idle while drawing baseline power. Through Data Layout Interleaving (DLI), Legolas++ achieves peak arithmetic density, reducing energy consumption per operation by up to 87%:

Application Domain Scale Performance (Apple M1 Max) Speedup Energy & Power Footprint Showcase
Video Pipeline (CPU) Vision & Broadcast 32 HD 720p streams 8,442 FPS (7.78 GPix/s) 5.34× 3.55 µJ / frame (0.259 GPix/W) Tutorial
Video Pipeline (Metal GPU) Vision & Broadcast 32 HD 720p streams 17,264 FPS (15.91 GPix/s) 57.1× 2.43 µJ / frame (0.379 GPix/W) · 87.2% energy saved Tutorial
Audio IIR Biquad Audio DSP & Acoustics 64 parallel channels 5,832 MSamples/sec 14.7× 5.14 nJ / sample (< 0.05 W total draw) Tutorial
Depthwise 2D Conv Edge AI (MobileNet) 512 channels, $3\times3$ 212.3 GFlops (8 Cores) 5.83× 7.08 GFlops/Watt · 82.8% energy saved Tutorial
MultiThomas Solver Scientific Computing 262,144 systems, 134M unknowns 69.40 GFlops (8 Cores) 33.1× 2.31 GFlops/Watt (33× more math per Joule) Tutorial

Real-Time Multi-Stream Video Pipeline (CPU & Metal GPU)

Location: examples/VideoPipeline/VideoPipeline.cxx & VideoPipelineMetal.mm

  • 3×3 spatial Sobel edge detection fused with quadratic temporal motion differencing across 32 concurrent 720p HD feeds.
  • CPU Line-Rate: 8,442.1 FPS (7.78 GPixels/s), 3.79 ms latency per 32-frame batch (5.34× speedup vs scalar).
  • Metal GPU Line-Rate: 17,263.9 FPS (15.91 GPixels/s), 1.85 ms latency (57.15× speedup vs scalar, 7.57× vs CPU).
  • Energy Impact: Drops energy per HD frame from 19.0 µJ (scalar) down to 2.43 µJ on Metal GPU (87.2% reduction). A single 40 W laptop sustains 287 concurrent 60 FPS feeds without thermal throttling.

Digital Audio Processing: 64-Track IIR Biquad Filter

Location: examples/AudioBiquad/AudioBiquad.cxx

  • 2nd-order Direct Form II recursive IIR filter ($y[n] = b_0 x[n] + b_1 x[n-1] + b_2 x[n-2] - a_1 y[n-1] - a_2 y[n-2]$).
  • Temporal recurrence vectorized across 64 audio channels using $P=4$ NEON SIMD packing.
  • Result: 5,832 Megasamples/sec (14.66× speedup).
  • Energy Impact: Filters 64 channels of studio 192 kHz audio using less than 0.05 W of power (5.14 nJ per sample). Replaces multi-kilowatt dedicated hardware DSP racks with pure C++ software.

AI & Edge Computer Vision: Depthwise Separable 2D Convolution

Location: examples/DepthwiseConv/DepthwiseConv.cxx

  • Core compute primitive of lightweight vision neural networks (MobileNet, ConvNeXt, EfficientNet).
  • Evaluates $3\times3$ spatial filters across 512 channels simultaneously with NEON fmla.4s vector instructions.
  • Result: 212.3 GFlops (5.83× speedup over scalar).
  • Energy Impact: Delivers 7.08 GFlops/Watt on CPU with zero PCIe memory copy latency, cutting inference power consumption by 82.8% to extend battery life in robotics and drones.

Scientific Computing & PDEs: MultiThomas Recurrence

Location: tst/MultiThomas/MultiThomas.cxx

  • 262,144 tridiagonal systems (134.2 million unknowns) for Alternating Direction Implicit (ADI) heat diffusion and Navier-Stokes sweeps.
  • Result: 69.40 GFlops (33.05× speedup over scalar).
  • Energy Impact: Delivers 2.31 GFlops/Watt on tightly coupled recurrences where optimizing compilers drop to 0.07 GFlops/Watt, achieving 33× higher compute density per Joule.

📐 Roofline Performance Models: Proximity to Physical Hardware Ceilings

How close does Legolas++ operate to the physical limits of modern silicon? We model every workload using the Roofline Model ($P_{\text{attainable}} = \min(P_{\text{peak}}, I \times B_{\text{peak}})$):

Benchmark / Kernel Arithmetic Intensity ($I$) Physical Bottleneck Hardware Efficiency (Measured)
MultiThomas ($N_x=512$) $0.65 \text{ FLOP/Byte}$ Memory Bandwidth (DRAM) 81% of GPU VRAM limit / 45% of CPU DRAM
VideoPipeline (32×720p) $1.50 \text{ FLOP/Byte}$ Streaming Memory Bus 79.5% of Apple UMA Peak (191 GB/s)
AudioBiquad (64 Tracks) $1.00 \text{ FLOP/Byte}$ Recurrence Dependency / FMA Pipe 100% of FMA Vector Issue Capacity
DepthwiseConv (128ch) $2.25 \text{ FLOP/Byte}$ L2/L3 Cache Bandwidth 85% of L2/L3 Peak Bandwidth (628 GFlops)
OptionPricing (16k opts) $\approx 18 \text{ FLOP/Byte}$ Compute-Bound (L3 Cache) 88% of Multi-Core FMA Execution (1.5M opts/s)
Reductions (squaredNorm) $0.50 \text{ FLOP/Byte}$ VRAM Streaming 75% of Physical VRAM Peak (335 GB/s)

For complete mathematical derivations, memory traffic equations, and vector plots, see the Roofline Performance Model Guide.

📊 How Legolas++ Compares to Existing Solutions

Feature Standard Compilers (GCC/Clang/MSVC) Traditional Linear Algebra (BLAS, Armadillo, Blaze) Frameworks (PyTorch C++ ATen, oneDNN) Legolas++
Vectorize Recurrences ($X_i = f(X_{i-1})$) ❌ Fails (falls back to $1\times$ scalar) ❌ Scalar loops ($1\times$) ❌ Requires custom handwritten kernels ✅ Automatic SIMD ($4\times - 16\times$)
Write Once, Vectorize Everywhere ❌ Requires manual intrinsics / pragmas ⚠️ Limited to non-recursive ops ❌ Separate CPU/GPU/SIMD code paths ✅ Single generic scalar syntax
Data Layout Interleaving (DLI) ❌ Manual array-of-structs reshuffling ❌ Fixed row/column major layouts ⚠️ Heavy tensor transpose ops ✅ Native in the tensor type system
External Dependencies None BLAS / LAPACK Heavy (LibTorch, MKL, Python) ✅ Zero (Pure Standard C++14)
Build & Integration Model N/A Variable Multi-gigabyte binaries ✅ 100% Header-Only (INTERFACE)

Two-Level Decoupled Architecture

┌─────────────────────────────────────────────────────────────┐
│                 Level 1: Multi-Core Scaling                 │
│   Legolas::parmap(...) via persistent StaticThreadPool      │
│   Dynamic chunks claimed from an atomic counter (4/worker)  │
│   Hybrid spin-then-sleep idle policy (no hot-path locks)    │
└──────────────────────────────┬──────────────────────────────┘
                               │ (blocked ranges)
┌──────────────────────────────▼──────────────────────────────┐
│           Level 2: Hardware SIMD Vectorization (Core)       │
│   Legolas::Array<T, D, P, DP> (Data Layout Interleaving)    │
│   Packed views mapped directly to ARM NEON / x86 AVX2/512   │
│   Zero memory permute overhead during computation           │
└─────────────────────────────────────────────────────────────┘

🛠️ Installation & Integration

Legolas++ is designed with a strict 100% header-only and zero-dependency philosophy. No binary libraries (.a, .so, .dylib, .dll) need to be compiled or linked.

Prerequisites

  • Standard C++14 compliant compiler:
    • Apple Clang $\ge 12$
    • GCC $\ge 7$
    • LLVM Clang $\ge 8$
    • Microsoft Visual C++ (MSVC) $\ge 2017$
  • CMake $\ge 3.5$
  • Zero external dependencies required (native standard threads and SIMD wrappers)

Quick Integration (CMake FetchContent)

Add Legolas++ to your CMakeLists.txt in 4 lines:

include(FetchContent)
FetchContent_Declare(
  Legolas
  GIT_REPOSITORY https://github.com/LaurentPlagne/Legolas.git
  GIT_TAG        master
)
FetchContent_MakeAvailable(Legolas)

add_executable(my_solver main.cpp)
target_link_libraries(my_solver PRIVATE Legolas)

Alternatively, if using Git submodules or a local clone:

add_subdirectory(path/to/Legolas)
target_link_libraries(my_solver PRIVATE Legolas)

Or simply add the include directory to your compiler invocation:

c++ -O3 -std=c++14 -I/path/to/Legolas -I/path/to/Legolas/Legolas/include main.cpp -o my_solver

Standard Build & Testing

# Configure Release build
cmake -B build -DCMAKE_BUILD_TYPE=Release

# Compile tests and benchmarks
cmake --build build -j

# Run the automated test suite (100% passing)
ctest --test-dir build --output-on-failure

🔥 GPU Acceleration: Optional Vulkan Compute Backend (Linux & Windows)

Starting with version 2.1, Legolas++ includes an optional, fully self-contained Vulkan compute backend designed for Linux and Windows platforms (macOS utilizes the native Apple Metal backend).

Architectural Highlights

  • Zero Link-Time Dependencies: The Vulkan loader is dynamically opened at runtime via dlopen (Linux) or LoadLibrary (Windows) through vkGetInstanceProcAddr (VK_NO_PROTOTYPES). The backend links only ${CMAKE_DL_LIBS}. No Vulkan SDK or import libraries are needed at link time.
  • In-Tree Versioned SPIR-V: Compute shaders are precompiled and versioned directly inside the repository (Legolas/Vulkan/spv/*.hxx) as const uint32_t[] arrays. Neither glslc nor any shader compiler SDK is required to build or run.
  • Transparent CPU Fallback: If no Vulkan-capable GPU or driver is present on the system, Legolas::Vulkan::Context::available() evaluates to false, and execution falls back cleanly to the CPU path.
  • GPU Mapping of Data Layout Interleaving (DLI):
    • A GPU thread naturally replaces a CPU SIMD lane.
    • Transposing the problem ensemble to an [step][system] (Structure of Arrays / SoA) layout guarantees 100% memory coalescing across GPU warps/wavefronts.
    • For tridiagonal recurrences, each GPU thread processes 4 systems simultaneously via vec4, reducing address calculations, exposing instruction-level parallelism (ILP), and boosting compute throughput by 4.8× over scalar thread mappings.

Benchmark Results (NVIDIA GeForce RTX 2060 SUPER)

Workloads measured against single-core scalar, single-core AVX2 DLI, and a 12-thread CPU baseline (AMD Ryzen 5 3600):

Workload GPU Kernel Time Speedup vs Scalar Speedup vs 12-Thread CPU Verdict
VideoPipeline (32×720p Sobel + Temporal) 0.98 ms 57.3× 18.0× 🚀 Peak line-rate (32,659 FPS)
DepthwiseConv (128ch, 3×3 MobileNet) 0.044 ms 12.8× 4.9× 🚀 628 GFlops sustained
MultiThomas ($N_x=512$, 262,144 systems) 13.8 ms 58.0× 9.1× 🚀 127 GFlops sustained
AudioBiquad (64ch × 960k samples, blocked scan) 9.9 ms 15.7× 2.8× 🚀 6,176 MSamples/s
Reductions (squaredNorm, dot, 16.7M floats) 0.20–0.32 ms — 39–62× 🚀 Bandwidth saturated

⚠️ Key Takeaway: Resident Data vs. Discrete PCIe Transfers
On discrete GPUs, single-shot isolated executions are bounded by the PCIe bus transfer rate (~4.5 GB/s). The GPU backend delivers massive speedups when data stays resident in VRAM (iterative PDE solvers, multi-pass rendering/video pipelines, or audio streaming), or on Unified Memory Architectures (UMA) such as Apple Silicon or integrated APUs where zero-copy host-visible buffers are used (Context::isUnifiedMemory()).

How to Build & Run the Vulkan Benchmarks

Enabling the Vulkan backend is opt-in via CMake:

# 1. Configure with Vulkan enabled (requires vulkan headers on Linux/Windows)
cmake -B build -DCMAKE_BUILD_TYPE=Release -DLEGOLAS_ENABLE_VULKAN=ON

# 2. Build the benchmark suite and tests
cmake --build build -j

# 3. Run all tests (including Vulkan verification)
ctest --test-dir build --output-on-failure

# 4. Run the comprehensive Vulkan benchmark suite:
./build/examples/VulkanBench                # Full suite
./build/examples/VulkanBench --skip-heavy    # Fast smoke test

For complete implementation notes, architectural details, and regression logs, see vulkan.md.

🍎 Native Apple Metal Backend (macOS)

On macOS systems, Legolas++ leverages native Apple Metal compute (examples/VideoPipeline/VideoPipelineMetal.mm). Thanks to Apple Silicon's Unified Memory Architecture (UMA), the CPU and GPU share the same physical memory space:

  • Zero PCIe Transfer Bottleneck: Frames and tensors are written directly into shared buffers (MTLResourceStorageModeShared), eliminating upload and download overhead.
  • Peak End-to-End Speedup: Achieves 17,264 FPS on 32 concurrent 720p streams (57.15× speedup over scalar CPU, 7.57× faster than 8-core CPU), consuming only 2.43 µJ per frame.

📖 Documentation Website

Comprehensive tutorials, mathematical proofs, architecture guides, and API reference are available on the official documentation website:

🔗 https://laurentplagne.github.io/Legolas/

To preview the documentation site locally:

pip install -r requirements-docs.txt
mkdocs serve
# Open http://127.0.0.1:8000 in your browser

📜 Academic Background

Legolas++ is based on research presented at ACM SIGPLAN ARRAY:

Portable vectorization and parallelization of C++ multi-dimensional array computations
Laurent Plagne & Kavoos Bojnourdi
Proceedings of the 4th ACM SIGPLAN International Workshop on Libraries, Languages, and Compilers for Array Programming (ARRAY 2017), Pages 47–54.
DOI: 10.1145/3091966.3091973 · Available on ResearchGate

The foundational multi-target abstraction and data-structure concepts of Legolas++ originated in the PhD thesis of Wilfried Kirschenmann:

Vers des noyaux de calcul intensif pérennes (Towards sustainable high performance compute kernels)
Wilfried Kirschenmann (directed by Stéphane Vialle, advised by Laurent Plagne)
PhD Thesis in Computer Science, Université de Lorraine (2012), NNT: 2012LORR0203.
Available on HAL Theses · theses.fr

📄 License

This project is distributed under the terms of the MIT "Expat" License (Copyright (c) 2019-2026 EDF-R&D, TriScale innov). See License.md for details.

About

No description, website, or topics provided.

Resources

Readme
Activity

Stars

2 stars

Watchers

0 watching

Forks

2 forks
Report repository

Releases

Packages

Contributors

Languages