Pure Rust BLAS/LAPACK implementation for the SciRS2 ecosystem

OxiBLAS is a production-grade, pure Rust implementation of BLAS (Basic Linear Algebra Subprograms) and LAPACK (Linear Algebra PACKage). Designed as the foundational linear algebra library for the SciRS2 scientific computing ecosystem.

• Pure Rust - No C dependencies, fully portable, safe by default
• SIMD Optimized - Custom SIMD layer using core::arch intrinsics
  • x86_64: AVX2/FMA (256-bit), AVX512F (512-bit)
  • AArch64: NEON (128-bit)
  • Fallback: Scalar operations
• Cache Aware - BLIS-style blocked algorithms for optimal cache usage
• Parallel - Optional rayon-based parallelization for multi-core systems
• Complete BLAS - Full Level 1, 2, 3 operations (including packed/banded)
• Extensive LAPACK - LU, Cholesky, QR, SVD, EVD, Schur, Hessenberg
• Extended Precision - f16, f128 (quad precision), Kahan/pairwise summation
• Tensor Operations - Einstein summation (24 patterns), batched matmul, batched BLAS
• Sparse Support - 9 formats (CSR, CSC, COO, ELL, DIA, BSR, BSC, HYB, SELL) with advanced solvers
• Advanced Preconditioners - AMG, SPAI, AINV, Schwarz, polynomial
• Recursive Factorizations - Cache-oblivious Cholesky, LU, QR (compute_recursive())
• Parallel Factorizations - Parallel blocked Cholesky and LU (compute_blocked_par())
• Batched BLAS - gemm_batched, gemm_strided_batched, axpy_batched, gemv_batched
• Mixed-Precision Refinement - f32 factorization + f64 residual for LU, Cholesky, QR
• Runtime Auto-Tuning - RuntimeAutoTuner selects optimal block sizes at runtime
• Multifrontal Sparse - MultifrontalCholesky, MultifrontalLU for large sparse systems
• no_std Support - oxiblas-core and oxiblas-matrix support #![no_std] with alloc
• Comprehensive Benchmarks - Direct performance comparison with OpenBLAS

Add to your Cargo.toml:

[dependencies]
oxiblas = "0.2"

# With parallelization
oxiblas = { version = "0.2", features = ["parallel"] }

# With extended precision
oxiblas = { version = "0.2", features = ["f16", "f128"] }

# All features
oxiblas = { version = "0.2", features = ["full"] }

use oxiblas_blas::level3::gemm;
use oxiblas_matrix::Mat;

// Create matrices
let a = Mat::from_rows(&[
    &[1.0, 2.0, 3.0],
    &[4.0, 5.0, 6.0],
]);
let b = Mat::from_rows(&[
    &[7.0, 8.0],
    &[9.0, 10.0],
    &[11.0, 12.0],
]);
let mut c = Mat::zeros(2, 2);

// GEMM: C = A * B
gemm(1.0, a.as_ref(), b.as_ref(), 0.0, c.as_mut());
// Result: [[58, 64], [139, 154]]

use oxiblas_lapack::lu::Lu;
use oxiblas_matrix::Mat;

let a = Mat::from_rows(&[
    &[2.0, 1.0],
    &[1.0, 3.0],
]);

let lu = Lu::compute(a.as_ref()).expect("Matrix is not singular");
let det = lu.determinant();
assert!((det - 5.0_f64).abs() < 1e-10); // det = 2*3 - 1*1 = 5

use oxiblas_lapack::svd::Svd;
use oxiblas_matrix::Mat;

let a = Mat::from_rows(&[
    &[3.0, 1.0],
    &[1.0, 3.0],
    &[1.0, 1.0],
]);

let svd = Svd::compute(a.as_ref()).expect("SVD failed");
let singular_values = svd.singular_values();  // [4.24, 2.00]
let u = svd.u();   // Left singular vectors
let vt = svd.vt(); // Right singular vectors (transposed)

OxiBLAS includes comprehensive examples demonstrating all major features:

# Basic BLAS operations (Level 1/2/3)
cargo run --example basic_blas

# LAPACK decompositions (LU, QR, Cholesky, SVD, EVD)
cargo run --example lapack_decompositions

# Extended precision (f128, Kahan, pairwise summation)
cargo run --example extended_precision --features f128

# Tensor operations & Einstein summation (24 patterns)
cargo run --example tensor_operations

# Sparse matrices (CSR, CSC, COO, iterative solvers, preconditioners)
cargo run --example sparse_matrices --features parallel

See crates/oxiblas/examples/ for complete source code.

OxiBLAS is organized as a workspace with specialized sub-crates:

OxiBLAS supports multiple precision levels:

use oxiblas::prelude::*;

// Quad-precision (f128) - ~31 decimal digits
#[cfg(feature = "f128")]
{
    use oxiblas_core::scalar::QuadFloat;
    let x = QuadFloat::from(2.0);
    let sqrt_x = x.sqrt();
    // Extremely high precision calculations
}

// Extended precision dot products
let x: Vec<f64> = vec![/* ... */];
let y: Vec<f64> = vec![/* ... */];

// Kahan summation (compensated)
let result = dot_kahan(&x, &y);

// Pairwise summation (divide-and-conquer)
let result = dot_pairwise(&x, &y);

// Mixed precision (f32 computation, f64 accumulation)
let x_f32: Vec<f32> = vec![/* ... */];
let y_f32: Vec<f32> = vec![/* ... */];
let result_f64 = dsdot(&x_f32, &y_f32);

OxiBLAS includes comprehensive tensor contraction support via einsum:

use oxiblas_blas::tensor::einsum;

// Matrix multiplication
let c = einsum("ij,jk->ik", &a, &[m, n], Some((&b, &[n, k])))?;

// Outer product
let c = einsum("i,j->ij", &x, &[m], Some((&y, &[n])))?;

// Transpose
let at = einsum("ij->ji", &a, &[m, n], None)?;

// Matrix trace
let trace = einsum("ii->", &a, &[n, n], None)?;

// Hadamard (element-wise) product
let c = einsum("ij,ij->ij", &a, &[m, n], Some((&b, &[m, n])))?;

// Tensor contractions (24 patterns supported)
let result = einsum("ijk,kl->ijl", &tensor3d, &[i, j, k], Some((&mat, &[k, l])))?;

Supported Patterns:

• Matrix operations: matmul, transpose, trace, diagonal
• Tensor transposes: 3 variants (ijk→ikj, ijk→jik, ijk→kji)
• Reductions: row/column sums, total sum, axis sums
• Products: outer, Hadamard, dot, Frobenius
• Advanced: tensor-matrix contraction, batched operations

use oxiblas_blas::tensor::{batched_matmul, Tensor3};

// Batch of matrix multiplications
let a = Tensor3::from_data(data_a, batch_size, m, k);
let b = Tensor3::from_data(data_b, batch_size, k, n);
let c = batched_matmul(&a, &b)?;
// c[i] = a[i] * b[i] for each batch i

• CSR - Compressed Sparse Row
• CSC - Compressed Sparse Column
• COO - Coordinate format
• ELL - ELLPACK format
• DIA - Diagonal format
• BSR - Block Sparse Row
• BSC - Block Sparse Column
• HYB - Hybrid ELL+COO
• SELL - Sliced ELLPACK (GPU-optimized)

• GMRES with restart (includes PGMRES, FGMRES)
• Conjugate Gradient (CG, PCG, Block-CG)
• BiCGStab
• MINRES (includes PMINRES)
• QMR (Quasi-Minimal Residual)
• TFQMR (Transpose-Free QMR)
• IDR(s) (Induced Dimension Reduction)
• Block-GMRES (multiple RHS)
• Lanczos iteration (symmetric)
• Arnoldi iteration (general)
• IRAM (Implicitly Restarted Arnoldi Method)

• ILU(0) / ILUT (Incomplete LU with threshold)
• ILUTP (ILUT with pivoting)
• Incomplete Cholesky (IC0, ICT)
• Jacobi / Block Jacobi
• Gauss-Seidel / SOR / SSOR
• AMG (Algebraic Multigrid - classical Ruge-Stüben)
• SA-AMG (Smoothed Aggregation AMG)
• SPAI (Sparse Approximate Inverse)
• AINV (Approximate Inverse)
• Additive Schwarz (domain decomposition)
• Polynomial preconditioners (Neumann, Chebyshev)

• Lanczos method (symmetric matrices)
• Block Lanczos
• Arnoldi iteration (general matrices)
• Block Arnoldi
• Shift-invert spectral transformation
• Polynomial filtering (Chebyshev)
• Interval eigenvalue computation (Sturm sequence)
• Truncated SVD
• Randomized SVD
• Incremental SVD (Brand algorithm)

• RCM (Reverse Cuthill-McKee)
• AMD (Approximate Minimum Degree)
• MMD (Multiple Minimum Degree)
• COLAMD (Column AMD for unsymmetric/rectangular)
• Nested Dissection (level-set based)

OxiBLAS implements BLIS-style blocked algorithms with architecture-specific SIMD kernels and platform-aware cache tuning.

Benchmarked: 2025-12-27 CPU: Intel Xeon E5-2623 v4 @ 2.60GHz (Broadwell-EP) Cache: L1D=32KB, L2=256KB, L3=10MB SIMD: AVX2, FMA, SSE4.1, SSE4.2

Peak Performance: 13.75 Gelem/s = 220 GFLOPS (256×256 f64) Key Achievement: OxiBLAS outperforms OpenBLAS by 2% on very large matrices (1024×1024)

Summary: OxiBLAS achieves 80-112% of OpenBLAS performance on Linux x86_64, with superior performance on small f32 and very large f64 matrices.

• 13-20% performance improvement after cache tuning for 256KB L2 systems
• Fine-tuned blocking parameters: KC=192, MC=128 (optimized for 256KB L2 cache)
• Increased prefetch distance: 12 iterations (optimized for Intel Xeon E5-2600 memory latency)
• Platform-aware cache detection: Linux sysfs, macOS sysctl, x86_64 CPUID fallback
• All 2,922 tests passing with zero warnings

Benchmarked: 2025-12-28 CPU: Apple M3 (ARM64 NEON) Cache: P-cores: L1D=128KB, L2=16MB SIMD: NEON (128-bit)

Peak Performance: 26.68 Gelem/s = 427 GFLOPS (1024×1024 f64) Key Achievement: OxiBLAS matches or exceeds OpenBLAS on Apple M3 for f64 GEMM

Breakthrough: OxiBLAS is 72% faster than OpenBLAS for large f32 matrices (1024×1024) on Apple M3!

Summary: OxiBLAS achieves 97-172% of OpenBLAS performance on Apple M3, with superior performance on f32 operations and competitive or better f64 performance.

• NEON micro-kernel excellence: 4×6 kernel optimized for f32, achieving 121-172% of OpenBLAS
• Large cache utilization: Optimal blocking for 16MB L2 cache on Apple M3 performance cores
• DOT product optimization: 65% faster than OpenBLAS (6.0 vs 3.6 Gelem/s)
• Production-ready: Competitive or superior performance across all major operations
• Platform-aware tuning: macOS sysctl cache detection working perfectly

Achievement: OxiBLAS matches or exceeds OpenBLAS on Apple M3 (97-172%), and is competitive on Linux x86_64 (80-112%). Pure Rust with NEON/AVX2 intrinsics.

Blocked algorithms using GEMM/TRSM as building blocks deliver dramatic speedups over the unblocked Level 2 baseline. All results measured on Apple M3 (f64).

Enable multi-core parallelization:

oxiblas = { version = "0.2", features = ["parallel"] }

Parallel operations automatically use all available cores via Rayon for large workloads.

OxiBLAS includes a comprehensive benchmarking suite with direct comparisons to OpenBLAS:

# Run all benchmarks
cargo bench --package oxiblas-benchmarks

# Compare with OpenBLAS (requires OpenBLAS installed)
cargo bench --package oxiblas-benchmarks --bench comparison --features compare-openblas

# View HTML reports
open target/criterion/report/index.html

See crates/oxiblas-benchmarks/README.md for detailed benchmarking guide.

OxiBLAS automatically detects and uses the best available SIMD:

• x86_64: AVX-512 > AVX2 > SSE4.2 > scalar
• aarch64: NEON (native 128-bit) + emulated 256-bit
• wasm32: SIMD128 (if enabled) + emulated 256-bit
• others: Scalar fallback

Runtime detection via SimdCapabilities::detect() with zero overhead after first call.

OxiBLAS ships a regress binary for tracking performance over time:

# Build the CLI
cargo build --release -p oxiblas-benchmarks

# Capture current performance as a baseline
./target/release/regress capture --output baseline.json

# On a subsequent commit/build, check for regressions (fail if >5% slower)
./target/release/regress check --baseline baseline.json --threshold 5.0

# Print a formatted performance report
./target/release/regress report --baseline baseline.json

# List all available benchmark names
./target/release/regress list

The check subcommand exits with a non-zero status when any operation regresses beyond the threshold, making it suitable for CI pipelines.

OxiBLAS provides multiple computation variants for dense factorizations, each suited to different use cases. The compute_auto() method is the recommended default — it selects between blocked and unblocked algorithms based on matrix size at zero cost.

use oxiblas_lapack::{cholesky::Cholesky, lu::Lu, qr::Qr};

// Recommended: automatic algorithm selection
let chol = Cholesky::compute_auto(a.as_ref())?;  // up to 9.75× faster than compute()
let lu   = Lu::compute_auto(a.as_ref())?;         // up to 23× faster than compute()
let qr   = Qr::compute_auto(a.as_ref())?;         // 3-7× faster for large matrices

// Explicit blocked with custom block size
let chol = Cholesky::compute_blocked(a.as_ref(), 64)?;

// Cache-oblivious recursive (no tuning required)
let lu = Lu::compute_recursive(a.as_ref())?;

// Parallel blocked (requires "parallel" feature)
let chol = Cholesky::compute_blocked_par(a.as_ref())?;

For applications requiring double-precision accuracy with single-precision speed:

use oxiblas_lapack::refine::{mixed_precision_solve, mixed_precision_solve_cholesky};

// LU-based: factorize in f32, refine residuals in f64
let x = mixed_precision_solve(a.as_ref(), b.as_ref())?;

// Cholesky-based (symmetric positive definite systems)
let x = mixed_precision_solve_cholesky(a.as_ref(), b.as_ref())?;

For batches of independent small matrix operations (e.g., neural network layers):

use oxiblas_blas::batched::{gemm_batched, gemm_strided_batched, axpy_batched};

// Batch of independent GEMMs: C[i] = A[i] * B[i]
gemm_batched(1.0, &a_slices, &b_slices, 0.0, &mut c_slices)?;

// Strided batched GEMM (contiguous memory layout)
gemm_strided_batched(1.0, &a, stride_a, &b, stride_b, 0.0, &mut c, stride_c, batch)?;

Benchmarks-specific features:

• Rust 1.85+ (Edition 2024)
• No external C dependencies

Version: 0.2.1 (2026-03-16)

• Lines of Code: ~223,935 Rust (371 files)
• Documentation: ~16,163 lines of comments, 12 comprehensive examples
• Tests: 2,922 passing lib tests + 287 doctests (100% success rate)
• Coverage:
  • ✅ Full BLAS Level 1/2/3 (including packed/banded variants)
  • ✅ Extensive LAPACK (LU, Cholesky, QR, SVD, EVD, Schur, Hessenberg)
  • ✅ Recursive cache-oblivious factorizations (Cholesky, LU, QR)
  • ✅ Parallel blocked factorizations (Cholesky, LU)
  • ✅ Fixed blocked QR with WY representation (3-7x speedup realized)
  • ✅ Complex bidiagonal reduction (ComplexBidiagFactors)
  • ✅ Batched BLAS operations (gemm_batched, axpy_batched, gemv_batched)
  • ✅ Mixed-precision iterative refinement (LU, Cholesky, symmetric, QR)
  • ✅ Runtime auto-tuning infrastructure (RuntimeAutoTuner)
  • ✅ Multifrontal sparse factorizations (MultifrontalCholesky, MultifrontalLU)
  • ✅ ndarray parallel GEMM and sparse integration
  • ✅ no_std support for oxiblas-core and oxiblas-matrix
  • ✅ Advanced sparse LU pivoting (threshold, static, Bunch-Kaufman LDL^T)
  • ✅ Sparse matrices (9 formats: CSR, CSC, COO, ELL, DIA, BSR, BSC, HYB, SELL)
  • ✅ Iterative solvers (GMRES, FGMRES, CG, BiCGStab, MINRES, QMR, TFQMR, IDR(s))
  • ✅ Advanced preconditioners (ILUT, IC, AMG, SA-AMG, SPAI, AINV, Schwarz)
  • ✅ Sparse eigenvalue/SVD (Lanczos, Arnoldi, IRAM, polynomial filtering, truncated/randomized SVD)
  • ✅ Tensor operations (einsum with 24 patterns, batched matmul)
  • ✅ Extended precision (f16, f128, Kahan/pairwise summation)
  • ✅ Comprehensive benchmarks with OpenBLAS comparison
  • ✅ Complete API documentation with examples

See TODO.md for the development roadmap.

OxiBLAS is part of the SciRS2 scientific computing ecosystem:

• SciRS2 - Scientific computing library
• NumRS2 - Numerical computing
• SkleaRS - Machine learning
• ToRSh - Tensor operations
• TrustformeRS - Transformers
• QuantRS2 - Quantum computing framework
• OxiRS - Semantic Web platform (SPARQL 1.2, GraphQL, Industrial Digital Twin, AI-augmented reasoning)

OxiBLAS is developed and maintained by COOLJAPAN OU (Team Kitasan).

If you find OxiBLAS useful, please consider sponsoring the project to support continued development of the Pure Rust ecosystem.

https://github.com/sponsors/cool-japan

Your sponsorship helps us:

• Maintain and improve the COOLJAPAN ecosystem
• Keep the entire ecosystem (OxiBLAS, OxiFFT, SciRS2, etc.) 100% Pure Rust
• Provide long-term support and security updates

Licensed under the Apache License, Version 2.0 (http://www.apache.org/licenses/LICENSE-2.0).

See LICENSE for details.

Contributions are welcome! Please feel free to submit issues and pull requests.