linear-massiv-0.1.0.6: README.md

# linear-massiv

**Pure Haskell numerical linear algebra implementing Golub & Van Loan.**

`linear-massiv` provides type-safe matrix and vector operations backed by
[massiv](https://hackage.haskell.org/package/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

```haskell
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)
  ```bash
  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
  ```bash
  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](https://hackage.haskell.org/package/linear) library:

```haskell
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`).