linear-massiv

Pure Haskell numerical linear algebra implementing Golub & Van Loan.

linear-massiv provides type-safe matrix and vector operations backed by massiv arrays, with compile-time dimension checking via GHC's type-level naturals (KnownNat). It requires no foreign FFI, no C dependencies, and no Fortran libraries — every algorithm is implemented in Haskell, from the GEMM micro-kernel to the divide-and-conquer SVD.

Dual-Mode Execution Architecture

To balance portability across compiler versions and architectures with high-throughput execution on vector hardware, linear-massiv employs a dual-mode execution architecture isolated behind an unboxed kernel interface:

1. Portable Unboxed Scalar Mode (Default):

   • High-performance scalar kernels written in pure Haskell over unboxed ByteArray# and MutableByteArray# buffers.
   • Fully portable across CPU architectures (x86_64, AArch64 / ARM64, Apple Silicon, RISC-V).
   • Compatible across GHC versions (9.6, 9.8, 9.10, 9.12, 9.14) without requiring LLVM or specific instruction set extensions.
   • Zero compilation panics with GHC Native Code Generator (NCG).

2. AVX2 SIMD Mode (Optional, flag +simd):

   • Vectorised kernels leveraging GHC 9.14's native DoubleX4# SIMD primops.
   • Requires compilation via LLVM (-fllvm -pgmlo opt-17 -pgmlc llc-17 -mavx2 -mfma).
   • Hardware accelerated on AVX2 + FMA capable x86_64 processors (Intel Haswell / AMD Zen or later).

All low-level kernel primops are strictly encapsulated within Numeric.LinearAlgebra.Massiv.Internal.Kernel, ensuring zero contamination of calling modules.

Performance and Algorithmic Framing

linear-massiv is designed as a transparent, auditable Haskell implementation of the standard algorithms presented in Golub & Van Loan's Matrix Computations (4th ed.).

Rather than competing against heavily tuned vendor libraries (OpenBLAS, Intel MKL, Apple Accelerate) which employ architecture-specific assembly heuristics and dynamic thread management, linear-massiv focuses on:

• Eliminating Haskell runtime and array indexing overhead in inner loops via unboxed memory buffers (ByteArray#).
• Cache-conscious panel factorisations (LU, Cholesky, QR).
• Exact compile-time dimension safety via type-level naturals.
• Complete algorithmic transparency without leaving the Haskell language.

Quick Start

import Numeric.LinearAlgebra.Massiv
import qualified Data.Massiv.Array as M

-- Create a 3x3 matrix (Primitive representation, Double elements)
let a = makeMatrix @3 @3 @M.P $ \i j ->
          fromIntegral (i * 3 + j + 1) :: Double

-- QR factorisation
let (q, r) = qr a

-- Solve Ax = b via LU
let b = makeVector @3 @M.P $ \i -> fromIntegral (i + 1) :: Double
let x = luSolve a b

-- Symmetric eigenvalue decomposition
let (eigenvalues, eigenvectors) = symmetricEigen a

-- Singular value decomposition
let (u, sigma, vt) = svd a

Requirements & Portability

• Compilers: GHC 9.6.7, 9.8.4, 9.10.3, 9.12.4, 9.14.1
• Build tool: Cabal 3.10+
• Default build: Pure scalar kernels (works on any platform with standard GHC)

  cabal build packages/linear-massiv
  cabal test packages/linear-massiv
  

• Optional SIMD acceleration: Requires GHC 9.14+ with LLVM 17 and an AVX2-capable CPU

  cabal build packages/linear-massiv --flags=simd --ghc-options="-fllvm -pgmlo opt-17 -pgmlc llc-17 -mavx2 -mfma"
  

Operations

BLAS-equivalent (GVL4 Ch. 1)

• Level 1: dot, axpy, scal, nrm2, asum
• Level 2: matrix-vector multiply (matvec), outer product (ger)
• Level 3: matrix multiply (matMul), transpose, matrix add/sub/scale — with register-blocked kernels, cache-blocked tiling, and micro-panel packing

Direct Solvers (GVL4 Chs. 3-4)

• LU: outer product LU with partial pivoting, panel factorisation, luSolve
• Cholesky: symmetric positive definite factorisation, gaxpy and outer-product variants, choleskySolve
• Banded and Tridiagonal: bandLU, bandForwardSub, bandBackSub, bandCholesky, and tridiagSolve
• Triangular solvers: forward substitution (forwardSub, forwardSubUnit), back substitution (backSub)

Orthogonal Factorisations (GVL4 Ch. 5)

• QR: Householder reflections with blocked WY accumulation (qr)
• Givens rotations: single and multi-column Givens elimination
• Least squares: linear least squares solving via QR factorisation

Eigenvalue Problems and SVD (GVL4 Chs. 7-8)

• Symmetric eigendecomposition: panel tridiagonalisation, QR iteration with aggressive early deflation, divide-and-conquer with secular equation solver
• SVD:
  • svd / svdAtAP: \(A^T A\) eigendecomposition path with D&C eigensolver and pre-scaled recovery
  • svdGKP: Golub-Kahan bidiagonalisation with panel blocking and D&C bidiagonal SVD
• General eigenvalue: Hessenberg reduction, Schur decomposition
• Iterative methods: power method (powerMethod), inverse iteration (inverseIteration), Rayleigh quotient iteration (rayleighQuotient)

Norms (GVL4 Ch. 2)

• Vector norms: Euclidean (\(L_2\)), Manhattan (\(L_1\)), Chebyshev (\(L_\infty\))
• Matrix norms: Frobenius (normFrob), 1-norm (norm1), \(\infty\)-norm (normInf), spectral 2-norm

Integration with the linear Library

Conversion functions are provided for interoperability with Edward Kmett's linear library:

import Numeric.LinearAlgebra.Massiv.Linear

-- Convert from linear's V3
let v = fromV3 @M.P (V3 1.0 2.0 3.0)

-- Convert from linear's V n
let w = fromLinearV @5 @M.B someLinearVector

Algorithmic References

The implementations follow:

• Golub, G. H. & Van Loan, C. F. (2013). Matrix Computations (4th ed.). Johns Hopkins University Press. (GVL4)
• Higham, N. J. (2002). Accuracy and Stability of Numerical Algorithms (2nd ed.). SIAM.

License

EUPL-1.2 OR Apache-2.0 (see LICENSE).