Skip to content

[OpenBLAS] Port the LAPACK factorisations to a native managed C# backend #627

Description

@Nucs

Overview

NumSharp's matrix products are fully managed and value-gated, but the factorisations NumPy computes through OpenBLAS's bundled LAPACK (inv, solve, svd, qr, eig, det, cholesky, lstsq, …) currently throw NotSupportedException. This tracks implementing them as a pure-managed C# backend behind the existing IBlasBackend seam — no native dependency, correctness-not-byte-parity. NumSharp.Core stays 100% managed.

Problem

NumSharp.Core ships no managed LU/QR/SVD/eigensolver, so 18 np.linalg.* entry points + norm matrix orders {2, -2, 'nuc'} + matrix_power(n<0) have complete API / validation / error-parity scaffolding (gated by LinAlgSignatureParityTests and LinAlgErrorParityTests) but zero numerics — every real call raises.

Verified by three independent techniques:

  • Static throw-site classification: every factorisation file carries a NoLapack(...) numerics throw; the product files carry none.
  • Test contract: LinAlgEngineSeamTests asserts each throws NotSupportedException ("delete a case as each implementation lands").
  • Live execution (no backend installed): 17 products/shape ops compute; 18/18 factorisations throw NotSupported.

Proposal

Implement managed LAPACK, filling the 15 default-false Try* members in Backends/IBlasBackend.LinearAlgebra.cs (validation and verbatim errors already run before the seam, so only algorithms are needed). Phased so each lands independently:

  • Phase 1 — LU (getrf/getrs/getri/gesv) → solve, inv, det, slogdet, matrix_power(<0), tensorinv, tensorsolve — ~2.5k C#
  • Phase 2 — Cholesky + QR (potrf, geqrf/orgqr) → cholesky, qr — ~2.95k C#
  • Phase 3 — SVD (gesdd, gelsd) → svd, svdvals, pinv, matrix_rank, cond, norm{2,-2,'nuc'}, lstsq — ~9k C#
  • Phase 4 — Eigen (syevd/heevd, geev) → eigh, eigvalsh, eig, eigvals — ~10.2k C#
  • Kernels — add managed IL/SIMD trsm / syrk / trmm (the ~690 IL-hot lines; gemm/gemv/dot already exist)
  • Complex (c/z) paths for complex eig/eigh/svd — +10–13k C#
  • Gate — tolerance-based linalg differential-fuzz tier (factorisations are not byte-parity-gated)

Evidence — verified completeness audit

Surface Count Completeness Confidence
Products / CBLAS (dot, matmul, inner, tensordot, vdot, vecdot, matvec, vecmat, outer, multi_dot, matrix_power≥0, norm non-SVD, trace, diagonal, cross, …) 17 ~99% (managed + value-gated) high (static + live + products.jsonl 287 cases)
Partial (norm, matrix_power, einsum) 3 ~65% high (live-verified split)
Factorisations / LAPACK 18 ~30% (validation only) high (3 techniques agree)

LOC estimate (reference-LAPACK dependency closure, code-only Fortran × ~1.2 → C#): ~24.7k for float32/float64, +10–13k for complex. Only ~690 lines (~3%) are IL/SIMD-hot leaves (and gemm is already implemented in SimdMatMul); the remaining ~24k is scalar driver code where IL codegen buys nothing — its speed comes entirely from the Level-3 kernels underneath.

Reuse lever: MIT-licensed MathNet.Numerics ships managed LU/QR/SVD/Cholesky/Eigen. Adapting it behind the seam could cut net-new to glue + gaps rather than a full Fortran translation (its results aren't byte-identical to NumPy either — fine, since factorisations aren't parity-gated).

Scope / Non-goals

  • Byte-parity with NumPy's OpenBLAS float output — out of scope. That is NumSharp.Interop.OpenBLAS's job (a route-for-route port that calls the same binary). Factorisation results aren't parity-gated, and NumPy's own LAPACK bits aren't reproducible across builds/CPUs/thread-counts anyway.
  • Byte-perfect managed GEMM port — reproducing OpenBLAS's DYNAMIC_ARCH accumulation order is months of empirical reverse-engineering per CPU arch and inherently fragile to CPU/thread/version; not this issue.
  • Products — already managed-implemented and value-gated; not re-opened here.
  • np.einsum contraction — a summation kernel + path planner, NumSharp's own work, not LAPACK; tracked separately.
  • Keep NumSharp.Core 100% managed — this backend is pure C#, no native dependency; it plugs into the same TensorEngine.Blas property the OpenBLAS package uses.

Related issues

Metadata

Metadata

Assignees

No one assigned

    Labels

    NumPy 2.x ComplianceAligns behavior with NumPy 2.x (NEPs, breaking changes)architectureCross-cutting structural changes affecting multiple componentsenhancementNew feature or request

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions