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

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

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