Legolas++: Building Blocks for Linear Algebra Solvers
High-Performance Modern C++ Tensor Engine for Automatic SIMD Vectorization of Recurrences via Data Layout Interleaving (DLI).
⚡ Explore Industry Showcases 📐 Roofline Model (75–100% of Peak) 🚀 30-Second Quickstart 💻 GitHub Repository
🚀 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:
- Single Natural Scalar Code: Written with standard sequential math (\(Y_i = A \cdot X_i + B \cdot Y_{i-1}\)) with zero intrinsics and zero compiler pragmas.
- Legolas++ DLI Engine: Transforms your memory layout via Data Layout Interleaving (DLI), turning impossible loop-carried recurrences into contiguous vector streams.
- A Multitude of Hardware Targets:
- ⚡ CPU SIMD Vectorized: Full vector register width utilization on ARM NEON (4s), x86 AVX2 (8s), and x86 AVX-512 (16s) at 100% hardware line-rate.
- 🚀 Multi-Core Parallel: Built-in persistent thread pool with dynamic chunk scheduling scaling across 4 to 128+ CPU cores, with SIMD active on every thread.
- 🔥 GPU Acceleration: Optional header-only Vulkan Compute backend for Linux and Windows, and native Apple Metal backend for macOS.
- 📐 Physical Hardware Ceilings: Operates at 75% to 100% of theoretical peak limits across memory bandwidth and FMA compute pipelines (see the Roofline Performance Model).
🎯 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 nightmare to maintain.
The Legolas++ Solution:
- 100% Machine Utilization: Fully utilizes all CPU cores and 100% of SIMD vector register widths simultaneously.
- Natural Scalar Notation: You write your core algorithm once, as a simple sequential loop in standard scalar math.
- Hardware Vectorization by Construction: Through Data Layout Interleaving (DLI), vectorization is structural in memory and guaranteed—no reliance on fragile compiler heuristics.
- 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.
🌐 An Ubiquitous Pattern Across Science & Industry
This computing pattern appears everywhere across high-performance engineering. Explore our dedicated tutorials and showcases:
| 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 |
🏆 Industrial Showcases & Energy Efficiency
🌿 Green Computing & HPC Energy Efficiency
High-performance computing is no longer solely about raw clock speed—it is fundamentally an energy efficiency challenge. Whether in battery-constrained drones, automotive vision systems, or megawatt hyperscale data centers, unvectorized scalar loops waste energy by leaving silicon execution pipelines idle while burning power.Through Data Layout Interleaving (DLI), Legolas++ packs problem ensembles directly into hardware vector registers, slashing the energy consumed per frame, audio sample, and linear solver by up to 88%.
🎥 32-Stream Real-Time Video Pipeline
Spatial 3×3 Sobel edge extraction fused with quadratic temporal motion differencing across 32 concurrent 720p HD feeds.
• Metal GPU: Drops to 2.43 µJ per HD frame (0.379 GPixels/Watt).
• Green Impact: 87.2% energy reduction vs scalar. Sustains 287 concurrent 60 FPS feeds on a 40 W laptop.
🎧 64-Channel Studio Audio IIR Biquad
Recursive Direct Form II digital biquad filtering across 64 parallel audio tracks. Overcomes the recursive feedback barrier.
• Energy Density: Consumes just 5.14 nJ per sample (194.4 MSamples/Joule).
• Green Impact: Replaces dedicated 500 W hardware DSP racks with pure software on low-power CPUs.
🧠 Depthwise Separable 2D Convolution
Compute primitive of lightweight vision neural networks. Vectorizes per-channel spatial filters across large batches.
• Zero Transfer Overhead: Avoids high-power discrete GPU wake-up latency.
• Green Impact: Cuts inference power by 82.8%, maximizing battery life for robotics and drones.
🔬 MultiThomas Tridiagonal Solver
Solving 262,144 coupled tridiagonal linear systems (134.2 million unknowns) for Alternating Direction Implicit (ADI) heat diffusion, wave equations, and fluid dynamics.
• Scale: Solves 134 million unknowns with zero memory reordering overhead.
• Green Impact: 33× higher compute density per Joule, drastically shrinking cloud electricity bills.
📐 Theoretical Limits & Hardware Efficiency
How close does Legolas++ get to the physical limits of modern silicon? Across our 6 industrial and scientific benchmarks, Legolas++ operates at 75% to 100% of theoretical hardware ceilings (memory bandwidth saturation and vector FMA issue line-rate).
👉 Explore the Roofline Performance Model & Hardware Efficiency Matrix →
What is Legolas++?
Legolas++ solves one of the most stubborn performance bottlenecks in scientific and numeric computing: loop-carried recurrences (\(X_i = f(X_{i-1})\)).
Compilers cannot auto-vectorize sequential recurrences because each step strictly requires the output of the preceding step. Standard optimizing compilers systematically fall back to scalar execution, leaving up to 90% of the CPU's vector hardware completely idle.
Instead of fighting the compiler with hand-written assembly intrinsics, Legolas++ introduces Data Layout Interleaving (DLI): When computing ensembles of independent problem instances (video feeds, audio tracks, ADI mesh sweeps, neural channels), Legolas++ packs elements across instances contiguously into physical memory:
Figure: Data Layout Interleaving (DLI) in Action. Elements at step \(i\) across \(P=4\) independent systems are interleaved contiguously in memory, enabling direct single-instruction hardware SIMD vector streaming (ldr q / vmovaps).
The exact same generic algorithm code written using natural scalar math compiles into full hardware SIMD registers (ARM NEON, x86 AVX2 / AVX-512) with zero code modifications and zero overhead.
Key Technical Highlights
- 🚀 Break the Recurrence Barrier: Vectorize tridiagonal solvers (Thomas algorithm), IIR digital filters, and depthwise convolutions with 100% hardware SIMD register saturation.
- 📦 100% Header-Only & Zero Dependencies: Drop into any project with
#include <Legolas/Array/Array.hxx>. No precompiled libraries, no link flags, and zero external dependencies. Pure standard C++14. - ⚡ Write Once, Vectorize Everywhere: A single template implementation works for scalar floats and hardware SIMD vectors (
Legolas::NativeSimd<T, P>). - 🔀 Two-Level Decoupled Parallelism:
- Data-Level (SIMD): Automatic via Data Layout Interleaving.
- Thread-Level (Multi-Core): Built-in, header-only zero-allocation parallel loop scheduler.
- 🍏 Cross-Platform: Validated on Apple Silicon ARM64 (NEON), Linux x86_64 (AVX2 / AVX-512), and Windows MSVC.
# 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)
target_link_libraries(my_project PRIVATE Legolas)
How Legolas++ Compares
| Feature | Compilers (GCC/Clang) | Traditional Linear Algebra (BLAS, Armadillo) | Deep Learning Frameworks (PyTorch, oneDNN) | Legolas++ |
|---|---|---|---|---|
| Vectorize Recurrences (\(X_i = f(X_{i-1})\)) | ❌ Fails (scalar \(1\times\)) | ❌ Scalar loops (\(1\times\)) | ❌ Custom handwritten kernels | ✅ Automatic SIMD (\(4\times - 16\times\)) |
| Hardware Ceiling Efficiency | < 5% of peak | 40%–70% (non-recurrences only) | 60%–80% (dense GEMM only) | ✅ 75%–100% of Physical Peak |
| Write Once, Vectorize Everywhere | ❌ Manual intrinsics needed | ⚠️ Limited to non-recursive ops | ❌ Separate CPU/GPU implementations | ✅ Single generic scalar syntax |
| Data Layout Interleaving (DLI) | ❌ Manual memory transforms | ❌ Fixed row/col major | ⚠️ Heavy tensor reshape overhead | ✅ Native in the tensor type |
| External Dependencies | None | BLAS / LAPACK | Heavy (LibTorch, Python, MKL) | ✅ Zero (Pure C++14) |
| Integration Model | N/A | Variable | Multi-GB binaries | ✅ 100% Header-Only |
30-Second Quick Example
Write your numerical solver once using natural scalar math; Legolas++ executes it at full hardware vector line-rate and parallelizes across cores simultaneously:
#include <iostream>
#include "Legolas/Array/Array.hxx"
#include "Legolas/Array/Map.hxx"
// 1. Define your algorithm once using natural 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. Declare 2D array interleaved across dimension 2 with pack size P=4:
using Array2D = Legolas::Array<float, 2, 4, 2>; // 1024 systems of size 256
Array2D X(1024, 256), Y(1024, 256);
X.fill(1.0f);
// 3. Parallelize across cores and vectorize in SIMD simultaneously:
Legolas::parmap(Scaler(), X, Y);
std::cout << "Computed 1024 systems at full hardware vector line-rate!" << std::endl;
return 0;
}
Citation & Academic Background
Legolas++ is based on the research presented at ACM SIGPLAN ARRAY 2017:
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
Legolas++ is distributed under the terms of the permissive MIT License (Copyright (c) 2019-2026 EDF-R&D, TriScale innov).