diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,5 @@
+# Revision history for keel-linalg
+
+## 0.1.0.0 -- 2026-08-18
+
+* First release.
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2026 Zhe Zhang
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/README.md b/README.md
new file mode 100644
--- /dev/null
+++ b/README.md
@@ -0,0 +1,31 @@
+# keel-linalg
+
+Dense linear algebra with zero build-time native dependencies: CBLAS
+level-3 and LAPACKE drivers over `Storable` vectors, resolved at run
+time from an OpenBLAS shared library through
+[keel-dyn](https://hackage.haskell.org/package/keel-dyn).
+
+```haskell
+import Keel.Linalg
+
+main :: IO ()
+main = do
+  Right be <- loadBackend defaultOpenBlasSpec
+  -- row-major dgemm: C := A(2x3) * B(3x2)
+  c <- dgemm be NoTrans NoTrans 2 2 3 1.0 a b 0.0
+  print c
+```
+
+The backend is located by the documented keel search policy
+(`KEEL_OPENBLAS` env override, the per-user keel data dir — populated
+by `keel setup openblas` from the umbrella package — then the system
+search path), probed for ILP64 misconfiguration via
+`openblas_get_config`, and pinned immutably in a `Backend` handle:
+no global state, no backend swapping under a pure API.
+
+Every driver is checked against numpy/scipy oracles in the test suite
+(GEMM, LU/solve/inverse, Cholesky, QR, SVD, symmetric and general
+eigenvalues, least squares — agreement within 1e-10).
+
+Part of the keel workspace — see the
+[project repository](https://github.com/skymanbp/keel).
diff --git a/keel-linalg.cabal b/keel-linalg.cabal
new file mode 100644
--- /dev/null
+++ b/keel-linalg.cabal
@@ -0,0 +1,71 @@
+cabal-version:      3.0
+name:               keel-linalg
+version:            0.1.0.0
+synopsis:           CBLAS/LAPACKE over a runtime-loaded OpenBLAS
+description:
+    Dense linear algebra for Haskell with zero build-time native
+    dependencies: CBLAS level-3 and LAPACKE drivers resolved at run time
+    from an OpenBLAS shared library through keel-dyn, over Storable
+    vectors.
+    .
+    The backend is located by the documented keel search policy
+    (@KEEL_OPENBLAS@ env override, then the per-user keel data dir, then
+    the system search path), probed for ILP64 misconfiguration via
+    @openblas_get_config@, and pinned immutably in a 'Backend' handle —
+    there is no global state and no backend swapping under a pure API.
+license:            MIT
+license-file:       LICENSE
+author:             Zhe Zhang
+maintainer:         Zhe Zhang
+category:           Math, Numeric
+build-type:         Simple
+homepage:           https://github.com/skymanbp/keel
+bug-reports:        https://github.com/skymanbp/keel/issues
+extra-doc-files:
+    CHANGELOG.md
+    README.md
+tested-with:        GHC ==9.10.3 || ==9.12.4 || ==9.14.1
+
+source-repository head
+    type:     git
+    location: https://github.com/skymanbp/keel
+
+common warnings
+    ghc-options:      -Wall
+    default-language: GHC2021
+
+library
+    import:           warnings
+    exposed-modules:
+        Keel.Linalg
+        Keel.Linalg.Backend
+    hs-source-dirs:   src
+    build-depends:
+        base >=4.20 && <4.23,
+        keel-dyn >=0.1 && <0.2,
+        vector >=0.13 && <0.14
+
+test-suite keel-linalg-smoke
+    import:           warnings
+    type:             exitcode-stdio-1.0
+    main-is:          Smoke.hs
+    other-modules:    TestBackend
+    hs-source-dirs:   test
+    build-depends:
+        base,
+        keel-dyn >=0.1 && <0.2,
+        keel-linalg >=0.1 && <0.2,
+        process >=1.6 && <1.7,
+        vector
+
+test-suite keel-linalg-oracle
+    import:           warnings
+    type:             exitcode-stdio-1.0
+    main-is:          Oracle.hs
+    other-modules:    TestBackend
+    hs-source-dirs:   test
+    build-depends:
+        base,
+        keel-linalg >=0.1 && <0.2,
+        process >=1.6 && <1.7,
+        vector
diff --git a/src/Keel/Linalg.hs b/src/Keel/Linalg.hs
new file mode 100644
--- /dev/null
+++ b/src/Keel/Linalg.hs
@@ -0,0 +1,739 @@
+-- | Dense linear algebra over 'VS.Vector' buffers, executed by a
+-- runtime-loaded OpenBLAS backend (see "Keel.Linalg.Backend").
+--
+-- Everything runs in 'IO' against an explicit 'Backend' handle — the
+-- immutable-pin discipline: no global state, no backend swapping under
+-- a pure API, results attributable to exactly one backend build.
+--
+-- Conventions, uniform across all drivers:
+--
+-- * matrices are dense, row-major, compact (leading dimension = column
+--   count); inputs are copied, never overwritten;
+-- * dimension mismatches throw 'DimensionMismatch' before any foreign
+--   call; a negative LAPACK @info@ throws 'LapackBadArgument' (a bug in
+--   these bindings, not in your data);
+-- * data-dependent failure (singular \/ not positive definite \/ rank
+--   deficient) is @Left info@ with LAPACK's 1-based index semantics —
+--   these are results, not exceptions;
+-- * every driver here requires the backend; none has a pure-Haskell
+--   fallback.
+module Keel.Linalg
+  ( -- * Backend
+    Backend
+  , BackendError (..)
+  , backendConfig
+  , openBackend
+  , openBackendWith
+  , closeBackend
+
+    -- * Errors and modes
+  , LinalgError (..)
+  , Transpose (..)
+  , Uplo (..)
+  , Diag (..)
+
+    -- * BLAS level 1
+  , ddot
+
+    -- * BLAS level 3
+  , dgemm
+
+    -- * Linear systems
+  , dgesv
+  , dposv
+  , dtrtrs
+
+    -- * Least squares
+  , dgels
+
+    -- * Factor \/ invert
+  , dgetrf
+  , dgetri
+  , dpotrf
+  , dpotri
+
+    -- * SVD
+  , dgesdd
+  , dgesvd
+
+    -- * Eigendecomposition
+  , dsyevd
+  , dgeev
+
+    -- * QR
+  , dgeqrf
+  , dorgqr
+  ) where
+
+import Control.Exception (Exception, throwIO)
+import Control.Monad (unless)
+import Data.Vector.Storable qualified as VS
+import Data.Vector.Storable.Mutable qualified as VSM
+import Foreign.C.String (castCharToCChar)
+import Foreign.C.Types (CChar (..), CInt (..))
+import Foreign.Ptr (FunPtr, Ptr, nullPtr)
+
+import Keel.Dyn (capOps)
+import Keel.Linalg.Backend
+
+-- | Programming errors — thrown, never returned. Data-dependent
+-- conditions (like a singular matrix) come back as @Either@ instead.
+data LinalgError
+  = DimensionMismatch String
+    -- ^ Operand sizes disagree; raised before any foreign call.
+  | LapackBadArgument String Int
+    -- ^ LAPACKE rejected argument /i/ (negative @info@) despite our
+    -- validation — indicates a bug in these bindings, please report.
+  deriving (Eq, Show)
+
+instance Exception LinalgError
+
+-- | Whether an operand is used as itself or its transpose.
+data Transpose = NoTrans | Trans
+  deriving (Eq, Show)
+
+-- | Which triangle of a symmetric\/triangular argument is stored\/read.
+data Uplo = Upper | Lower
+  deriving (Eq, Show)
+
+-- | Whether a triangular matrix has an implicit unit diagonal.
+data Diag = NonUnit | UnitDiag
+  deriving (Eq, Show)
+
+-- CBLAS enum values and LAPACKE char modes (frozen ABI: cblas.h, lapacke.h)
+cblasRowMajor, lapackRowMajor :: CInt
+cblasRowMajor = 101
+lapackRowMajor = 101
+
+transToC :: Transpose -> CInt
+transToC NoTrans = 111
+transToC Trans = 112
+
+transChar :: Transpose -> CChar
+transChar NoTrans = castCharToCChar 'N'
+transChar Trans = castCharToCChar 'T'
+
+uploChar :: Uplo -> CChar
+uploChar Upper = castCharToCChar 'U'
+uploChar Lower = castCharToCChar 'L'
+
+diagChar :: Diag -> CChar
+diagChar NonUnit = castCharToCChar 'N'
+diagChar UnitDiag = castCharToCChar 'U'
+
+-- Shared validation / info interpretation ------------------------------
+
+checkDim :: String -> String -> Int -> Int -> IO ()
+checkDim ctx what got want =
+  unless (got == want) . throwIO . DimensionMismatch $
+    ctx <> ": " <> what <> " has " <> show got <> " elements, want " <> show want
+
+-- Positive info = Left (data condition); zero = run the continuation;
+-- negative = bindings bug, thrown.
+interpretInfo :: String -> CInt -> IO a -> IO (Either Int a)
+interpretInfo ctx info onOk = case compare info 0 of
+  EQ -> Right <$> onOk
+  GT -> pure (Left (fromIntegral info))
+  LT -> throwIO (LapackBadArgument ctx (negate (fromIntegral info)))
+
+-- ---------------------------------------------------------------------
+-- BLAS
+
+foreign import ccall safe "dynamic"
+  callDdot
+    :: FunPtr (CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> IO Double)
+    -> CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> IO Double
+
+foreign import ccall safe "dynamic"
+  callDgemm
+    :: FunPtr
+         (  CInt -> CInt -> CInt
+         -> CInt -> CInt -> CInt
+         -> Double -> Ptr Double -> CInt
+         -> Ptr Double -> CInt
+         -> Double -> Ptr Double -> CInt
+         -> IO ()
+         )
+    -> CInt -> CInt -> CInt
+    -> CInt -> CInt -> CInt
+    -> Double -> Ptr Double -> CInt
+    -> Ptr Double -> CInt
+    -> Double -> Ptr Double -> CInt
+    -> IO ()
+
+-- | Dot product @x . y@ (@cblas_ddot@, unit strides).
+ddot :: Backend -> VS.Vector Double -> VS.Vector Double -> IO Double
+ddot be x y = do
+  checkDim "ddot" "y" (VS.length y) (VS.length x)
+  VS.unsafeWith x $ \px ->
+    VS.unsafeWith y $ \py ->
+      callDdot (opDdot (capOps be)) (fromIntegral (VS.length x)) px 1 py 1
+
+-- | @alpha * op(A) * op(B)@ (@cblas_dgemm@ with @beta = 0@ over compact
+-- row-major buffers): @op(A)@ is @m x k@, @op(B)@ is @k x n@, the
+-- result is @m x n@. @A@ itself is stored @m x k@ under 'NoTrans' and
+-- @k x m@ under 'Trans' (likewise @B@).
+dgemm
+  :: Backend
+  -> Transpose -- ^ op(A)
+  -> Transpose -- ^ op(B)
+  -> Int -- ^ m
+  -> Int -- ^ n
+  -> Int -- ^ k
+  -> Double -- ^ alpha
+  -> VS.Vector Double -- ^ A
+  -> VS.Vector Double -- ^ B
+  -> IO (VS.Vector Double)
+dgemm be ta tb m n k alpha a b = do
+  checkDim "dgemm" "A" (VS.length a) (m * k)
+  checkDim "dgemm" "B" (VS.length b) (k * n)
+  let lda = fromIntegral (if ta == NoTrans then k else m)
+      ldb = fromIntegral (if tb == NoTrans then n else k)
+      ldc = fromIntegral n
+  c <- VSM.new (m * n)
+  VS.unsafeWith a $ \pa ->
+    VS.unsafeWith b $ \pb ->
+      VSM.unsafeWith c $ \pc ->
+        callDgemm (opDgemm (capOps be))
+          cblasRowMajor (transToC ta) (transToC tb)
+          (fromIntegral m) (fromIntegral n) (fromIntegral k)
+          alpha pa lda pb ldb 0 pc ldc
+  VS.unsafeFreeze c
+
+-- ---------------------------------------------------------------------
+-- Linear systems
+
+foreign import ccall safe "dynamic"
+  callDgesv
+    :: FunPtr
+         (  CInt -> CInt -> CInt
+         -> Ptr Double -> CInt
+         -> Ptr CInt
+         -> Ptr Double -> CInt
+         -> IO CInt
+         )
+    -> CInt -> CInt -> CInt
+    -> Ptr Double -> CInt
+    -> Ptr CInt
+    -> Ptr Double -> CInt
+    -> IO CInt
+
+foreign import ccall safe "dynamic"
+  callDposv
+    :: FunPtr
+         (  CInt -> CChar -> CInt -> CInt
+         -> Ptr Double -> CInt
+         -> Ptr Double -> CInt
+         -> IO CInt
+         )
+    -> CInt -> CChar -> CInt -> CInt
+    -> Ptr Double -> CInt
+    -> Ptr Double -> CInt
+    -> IO CInt
+
+foreign import ccall safe "dynamic"
+  callDtrtrs
+    :: FunPtr
+         (  CInt -> CChar -> CChar -> CChar
+         -> CInt -> CInt
+         -> Ptr Double -> CInt
+         -> Ptr Double -> CInt
+         -> IO CInt
+         )
+    -> CInt -> CChar -> CChar -> CChar
+    -> CInt -> CInt
+    -> Ptr Double -> CInt
+    -> Ptr Double -> CInt
+    -> IO CInt
+
+-- | Solve @A X = B@ for square @A@ via LU with partial pivoting
+-- (@LAPACKE_dgesv@). @A@ is @n x n@, @B@ is @n x nrhs@. @Left i@:
+-- exactly singular, detected at @U(i,i)@.
+dgesv
+  :: Backend
+  -> Int -- ^ n
+  -> Int -- ^ nrhs
+  -> VS.Vector Double -- ^ A
+  -> VS.Vector Double -- ^ B
+  -> IO (Either Int (VS.Vector Double))
+dgesv be n nrhs a b = do
+  checkDim "dgesv" "A" (VS.length a) (n * n)
+  checkDim "dgesv" "B" (VS.length b) (n * nrhs)
+  aC <- VS.thaw a
+  bC <- VS.thaw b
+  ipiv <- VSM.new n :: IO (VSM.IOVector CInt)
+  info <-
+    VSM.unsafeWith aC $ \pa ->
+      VSM.unsafeWith bC $ \pb ->
+        VSM.unsafeWith ipiv $ \pp ->
+          callDgesv (opDgesv (capOps be))
+            lapackRowMajor (fromIntegral n) (fromIntegral nrhs)
+            pa (fromIntegral n) pp pb (fromIntegral nrhs)
+  interpretInfo "dgesv" info (VS.unsafeFreeze bC)
+
+-- | Solve @A X = B@ for symmetric positive-definite @A@ via Cholesky
+-- (@LAPACKE_dposv@); only the 'Uplo' triangle of @A@ is read. @Left i@:
+-- the leading minor of order @i@ is not positive definite.
+dposv
+  :: Backend
+  -> Uplo
+  -> Int -- ^ n
+  -> Int -- ^ nrhs
+  -> VS.Vector Double -- ^ A
+  -> VS.Vector Double -- ^ B
+  -> IO (Either Int (VS.Vector Double))
+dposv be uplo n nrhs a b = do
+  checkDim "dposv" "A" (VS.length a) (n * n)
+  checkDim "dposv" "B" (VS.length b) (n * nrhs)
+  aC <- VS.thaw a
+  bC <- VS.thaw b
+  info <-
+    VSM.unsafeWith aC $ \pa ->
+      VSM.unsafeWith bC $ \pb ->
+        callDposv (opDposv (capOps be))
+          lapackRowMajor (uploChar uplo) (fromIntegral n) (fromIntegral nrhs)
+          pa (fromIntegral n) pb (fromIntegral nrhs)
+  interpretInfo "dposv" info (VS.unsafeFreeze bC)
+
+-- | Solve @op(A) X = B@ for triangular @A@ (@LAPACKE_dtrtrs@); @A@ is
+-- read-only, only the 'Uplo' triangle is referenced ('UnitDiag' also
+-- skips the diagonal). @Left i@: @A(i,i)@ is exactly zero.
+dtrtrs
+  :: Backend
+  -> Uplo
+  -> Transpose -- ^ op(A)
+  -> Diag
+  -> Int -- ^ n
+  -> Int -- ^ nrhs
+  -> VS.Vector Double -- ^ A
+  -> VS.Vector Double -- ^ B
+  -> IO (Either Int (VS.Vector Double))
+dtrtrs be uplo ta diag n nrhs a b = do
+  checkDim "dtrtrs" "A" (VS.length a) (n * n)
+  checkDim "dtrtrs" "B" (VS.length b) (n * nrhs)
+  bC <- VS.thaw b
+  info <-
+    VS.unsafeWith a $ \pa ->
+      VSM.unsafeWith bC $ \pb ->
+        callDtrtrs (opDtrtrs (capOps be))
+          lapackRowMajor (uploChar uplo) (transChar ta) (diagChar diag)
+          (fromIntegral n) (fromIntegral nrhs)
+          pa (fromIntegral n) pb (fromIntegral nrhs)
+  interpretInfo "dtrtrs" info (VS.unsafeFreeze bC)
+
+-- ---------------------------------------------------------------------
+-- Least squares
+
+foreign import ccall safe "dynamic"
+  callDgels
+    :: FunPtr
+         (  CInt -> CChar
+         -> CInt -> CInt -> CInt
+         -> Ptr Double -> CInt
+         -> Ptr Double -> CInt
+         -> IO CInt
+         )
+    -> CInt -> CChar
+    -> CInt -> CInt -> CInt
+    -> Ptr Double -> CInt
+    -> Ptr Double -> CInt
+    -> IO CInt
+
+-- | Least squares \/ minimum norm via QR\/LQ (@LAPACKE_dgels@, no
+-- transpose): minimizes @||A X - B||@ when @m >= n@, finds the minimum
+-- norm solution when @m < n@. @A@ is @m x n@ full rank, @B@ is
+-- @m x nrhs@; the result is @n x nrhs@. @Left i@: @A@ is rank
+-- deficient (zero at diagonal element @i@ of the triangular factor).
+dgels
+  :: Backend
+  -> Int -- ^ m
+  -> Int -- ^ n
+  -> Int -- ^ nrhs
+  -> VS.Vector Double -- ^ A
+  -> VS.Vector Double -- ^ B
+  -> IO (Either Int (VS.Vector Double))
+dgels be m n nrhs a b = do
+  checkDim "dgels" "A" (VS.length a) (m * n)
+  checkDim "dgels" "B" (VS.length b) (m * nrhs)
+  let rows = max m n
+  aC <- VS.thaw a
+  -- LAPACKE wants B sized max(m,n) x nrhs: input in the first m rows,
+  -- solution comes back in the first n rows (row-major, ldb = nrhs,
+  -- so rows are contiguous and padding is a plain suffix)
+  bPad <- VSM.new (rows * nrhs)
+  VS.copy (VSM.slice 0 (m * nrhs) bPad) b
+  info <-
+    VSM.unsafeWith aC $ \pa ->
+      VSM.unsafeWith bPad $ \pb ->
+        callDgels (opDgels (capOps be))
+          lapackRowMajor (transChar NoTrans)
+          (fromIntegral m) (fromIntegral n) (fromIntegral nrhs)
+          pa (fromIntegral n) pb (fromIntegral nrhs)
+  interpretInfo "dgels" info (VS.take (n * nrhs) <$> VS.unsafeFreeze bPad)
+
+-- ---------------------------------------------------------------------
+-- Factor / invert
+
+foreign import ccall safe "dynamic"
+  callDgetrf
+    :: FunPtr (CInt -> CInt -> CInt -> Ptr Double -> CInt -> Ptr CInt -> IO CInt)
+    -> CInt -> CInt -> CInt -> Ptr Double -> CInt -> Ptr CInt -> IO CInt
+
+foreign import ccall safe "dynamic"
+  callDgetri
+    :: FunPtr (CInt -> CInt -> Ptr Double -> CInt -> Ptr CInt -> IO CInt)
+    -> CInt -> CInt -> Ptr Double -> CInt -> Ptr CInt -> IO CInt
+
+foreign import ccall safe "dynamic"
+  callDpotrf
+    :: FunPtr (CInt -> CChar -> CInt -> Ptr Double -> CInt -> IO CInt)
+    -> CInt -> CChar -> CInt -> Ptr Double -> CInt -> IO CInt
+
+foreign import ccall safe "dynamic"
+  callDpotri
+    :: FunPtr (CInt -> CChar -> CInt -> Ptr Double -> CInt -> IO CInt)
+    -> CInt -> CChar -> CInt -> Ptr Double -> CInt -> IO CInt
+
+-- | LU factorization with partial pivoting (@LAPACKE_dgetrf@) of an
+-- @m x n@ matrix: returns @(LU, ipiv)@ where @LU@ packs both factors
+-- and @ipiv@ (1-based) feeds 'dgetri'. @Left i@: @U(i,i)@ is exactly
+-- zero (the factor is dropped — downstream use would be invalid).
+dgetrf
+  :: Backend
+  -> Int -- ^ m
+  -> Int -- ^ n
+  -> VS.Vector Double -- ^ A
+  -> IO (Either Int (VS.Vector Double, VS.Vector CInt))
+dgetrf be m n a = do
+  checkDim "dgetrf" "A" (VS.length a) (m * n)
+  aC <- VS.thaw a
+  ipiv <- VSM.new (min m n) :: IO (VSM.IOVector CInt)
+  info <-
+    VSM.unsafeWith aC $ \pa ->
+      VSM.unsafeWith ipiv $ \pp ->
+        callDgetrf (opDgetrf (capOps be))
+          lapackRowMajor (fromIntegral m) (fromIntegral n)
+          pa (fromIntegral n) pp
+  interpretInfo "dgetrf" info
+    ((,) <$> VS.unsafeFreeze aC <*> VS.unsafeFreeze ipiv)
+
+-- | Matrix inverse from a 'dgetrf' factorization (@LAPACKE_dgetri@):
+-- pass the @n x n@ @LU@ and its @ipiv@. @Left i@: singular at
+-- @U(i,i)@.
+dgetri
+  :: Backend
+  -> Int -- ^ n
+  -> VS.Vector Double -- ^ LU from 'dgetrf'
+  -> VS.Vector CInt -- ^ ipiv from 'dgetrf'
+  -> IO (Either Int (VS.Vector Double))
+dgetri be n lu ipiv = do
+  checkDim "dgetri" "LU" (VS.length lu) (n * n)
+  checkDim "dgetri" "ipiv" (VS.length ipiv) n
+  luC <- VS.thaw lu
+  ipivC <- VS.thaw ipiv
+  info <-
+    VSM.unsafeWith luC $ \pa ->
+      VSM.unsafeWith ipivC $ \pp ->
+        callDgetri (opDgetri (capOps be))
+          lapackRowMajor (fromIntegral n) pa (fromIntegral n) pp
+  interpretInfo "dgetri" info (VS.unsafeFreeze luC)
+
+-- | Cholesky factorization (@LAPACKE_dpotrf@) of a symmetric
+-- positive-definite @n x n@ matrix; only the 'Uplo' triangle is read
+-- and only that triangle of the result is the factor (the other
+-- triangle keeps the input's bytes). @Left i@: leading minor of order
+-- @i@ not positive definite.
+dpotrf
+  :: Backend
+  -> Uplo
+  -> Int -- ^ n
+  -> VS.Vector Double -- ^ A
+  -> IO (Either Int (VS.Vector Double))
+dpotrf be uplo n a = do
+  checkDim "dpotrf" "A" (VS.length a) (n * n)
+  aC <- VS.thaw a
+  info <-
+    VSM.unsafeWith aC $ \pa ->
+      callDpotrf (opDpotrf (capOps be))
+        lapackRowMajor (uploChar uplo) (fromIntegral n) pa (fromIntegral n)
+  interpretInfo "dpotrf" info (VS.unsafeFreeze aC)
+
+-- | Inverse of a symmetric positive-definite matrix from its 'dpotrf'
+-- factor (@LAPACKE_dpotri@); only the 'Uplo' triangle of the result is
+-- the inverse. @Left i@: @factor(i,i)@ is exactly zero.
+dpotri
+  :: Backend
+  -> Uplo
+  -> Int -- ^ n
+  -> VS.Vector Double -- ^ Cholesky factor from 'dpotrf'
+  -> IO (Either Int (VS.Vector Double))
+dpotri be uplo n f = do
+  checkDim "dpotri" "factor" (VS.length f) (n * n)
+  fC <- VS.thaw f
+  info <-
+    VSM.unsafeWith fC $ \pa ->
+      callDpotri (opDpotri (capOps be))
+        lapackRowMajor (uploChar uplo) (fromIntegral n) pa (fromIntegral n)
+  interpretInfo "dpotri" info (VS.unsafeFreeze fC)
+
+-- ---------------------------------------------------------------------
+-- SVD
+
+foreign import ccall safe "dynamic"
+  callDgesdd
+    :: FunPtr
+         (  CInt -> CChar -> CInt -> CInt
+         -> Ptr Double -> CInt
+         -> Ptr Double
+         -> Ptr Double -> CInt
+         -> Ptr Double -> CInt
+         -> IO CInt
+         )
+    -> CInt -> CChar -> CInt -> CInt
+    -> Ptr Double -> CInt
+    -> Ptr Double
+    -> Ptr Double -> CInt
+    -> Ptr Double -> CInt
+    -> IO CInt
+
+foreign import ccall safe "dynamic"
+  callDgesvd
+    :: FunPtr
+         (  CInt -> CChar -> CChar -> CInt -> CInt
+         -> Ptr Double -> CInt
+         -> Ptr Double
+         -> Ptr Double -> CInt
+         -> Ptr Double -> CInt
+         -> Ptr Double
+         -> IO CInt
+         )
+    -> CInt -> CChar -> CChar -> CInt -> CInt
+    -> Ptr Double -> CInt
+    -> Ptr Double
+    -> Ptr Double -> CInt
+    -> Ptr Double -> CInt
+    -> Ptr Double
+    -> IO CInt
+
+-- Shared economy-size SVD wrapper: allocate s/U/VT, run the driver's
+-- foreign call, package the triple.
+svdWith
+  :: String
+  -> Int
+  -> Int
+  -> VS.Vector Double
+  -> (Ptr Double -> Ptr Double -> Ptr Double -> Ptr Double -> IO CInt)
+  -> IO (Either Int (VS.Vector Double, VS.Vector Double, VS.Vector Double))
+svdWith ctx m n a call = do
+  checkDim ctx "A" (VS.length a) (m * n)
+  let minmn = min m n
+  aC <- VS.thaw a
+  s <- VSM.new minmn
+  u <- VSM.new (m * minmn)
+  vt <- VSM.new (minmn * n)
+  info <-
+    VSM.unsafeWith aC $ \pa ->
+      VSM.unsafeWith s $ \ps ->
+        VSM.unsafeWith u $ \pu ->
+          VSM.unsafeWith vt $ \pvt ->
+            call pa ps pu pvt
+  interpretInfo ctx info
+    ((,,) <$> VS.unsafeFreeze s <*> VS.unsafeFreeze u <*> VS.unsafeFreeze vt)
+
+-- | Economy-size SVD by divide and conquer (@LAPACKE_dgesdd@,
+-- @jobz = \'S\'@): returns @(s, U, VT)@ with @s@ descending of length
+-- @min m n@, @U@ of @m x min m n@, @VT@ of @min m n x n@. @Left i@:
+-- the algorithm failed to converge.
+dgesdd
+  :: Backend
+  -> Int -- ^ m
+  -> Int -- ^ n
+  -> VS.Vector Double -- ^ A
+  -> IO (Either Int (VS.Vector Double, VS.Vector Double, VS.Vector Double))
+dgesdd be m n a =
+  svdWith "dgesdd" m n a $ \pa ps pu pvt ->
+    callDgesdd (opDgesdd (capOps be))
+      lapackRowMajor (castCharToCChar 'S')
+      (fromIntegral m) (fromIntegral n)
+      pa (fromIntegral n) ps
+      pu (fromIntegral (min m n))
+      pvt (fromIntegral n)
+
+-- | Economy-size SVD by QR iteration (@LAPACKE_dgesvd@,
+-- @jobu = jobvt = \'S\'@) — slower than 'dgesdd' but a different
+-- algorithm, useful as a cross-check. Same result shape as 'dgesdd'.
+-- @Left i@: @i@ superdiagonals failed to converge.
+dgesvd
+  :: Backend
+  -> Int -- ^ m
+  -> Int -- ^ n
+  -> VS.Vector Double -- ^ A
+  -> IO (Either Int (VS.Vector Double, VS.Vector Double, VS.Vector Double))
+dgesvd be m n a = do
+  superb <- VSM.new (max 1 (min m n - 1)) :: IO (VSM.IOVector Double)
+  svdWith "dgesvd" m n a $ \pa ps pu pvt ->
+    VSM.unsafeWith superb $ \psb ->
+      callDgesvd (opDgesvd (capOps be))
+        lapackRowMajor (castCharToCChar 'S') (castCharToCChar 'S')
+        (fromIntegral m) (fromIntegral n)
+        pa (fromIntegral n) ps
+        pu (fromIntegral (min m n))
+        pvt (fromIntegral n)
+        psb
+
+-- ---------------------------------------------------------------------
+-- Eigendecomposition
+
+foreign import ccall safe "dynamic"
+  callDsyevd
+    :: FunPtr
+         (  CInt -> CChar -> CChar -> CInt
+         -> Ptr Double -> CInt
+         -> Ptr Double
+         -> IO CInt
+         )
+    -> CInt -> CChar -> CChar -> CInt
+    -> Ptr Double -> CInt
+    -> Ptr Double
+    -> IO CInt
+
+foreign import ccall safe "dynamic"
+  callDgeev
+    :: FunPtr
+         (  CInt -> CChar -> CChar -> CInt
+         -> Ptr Double -> CInt
+         -> Ptr Double -> Ptr Double
+         -> Ptr Double -> CInt
+         -> Ptr Double -> CInt
+         -> IO CInt
+         )
+    -> CInt -> CChar -> CChar -> CInt
+    -> Ptr Double -> CInt
+    -> Ptr Double -> Ptr Double
+    -> Ptr Double -> CInt
+    -> Ptr Double -> CInt
+    -> IO CInt
+
+-- | Eigendecomposition of a symmetric matrix by divide and conquer
+-- (@LAPACKE_dsyevd@, @jobz = \'V\'@); only the 'Uplo' triangle is
+-- read. Returns @(w, V)@: eigenvalues ascending, and the eigenvector
+-- for @w[i]@ in /column/ @i@ of the row-major @n x n@ @V@ (i.e.
+-- @V[j*n + i]@), matching @numpy.linalg.eigh@. @Left i@: failed to
+-- converge.
+dsyevd
+  :: Backend
+  -> Uplo
+  -> Int -- ^ n
+  -> VS.Vector Double -- ^ A (symmetric)
+  -> IO (Either Int (VS.Vector Double, VS.Vector Double))
+dsyevd be uplo n a = do
+  checkDim "dsyevd" "A" (VS.length a) (n * n)
+  aC <- VS.thaw a
+  w <- VSM.new n
+  info <-
+    VSM.unsafeWith aC $ \pa ->
+      VSM.unsafeWith w $ \pw ->
+        callDsyevd (opDsyevd (capOps be))
+          lapackRowMajor (castCharToCChar 'V') (uploChar uplo)
+          (fromIntegral n) pa (fromIntegral n) pw
+  interpretInfo "dsyevd" info
+    ((,) <$> VS.unsafeFreeze w <*> VS.unsafeFreeze aC)
+
+-- | Eigendecomposition of a general matrix (@LAPACKE_dgeev@, right
+-- eigenvectors only). Returns @(wr, wi, VR)@ in LAPACK's packed real
+-- convention: eigenvalue @j@ is @wr[j] :+ wi[j]@; for a real
+-- eigenvalue, column @j@ of @VR@ is its eigenvector; a complex
+-- conjugate pair occupies columns @j@ (real part) and @j+1@ (imaginary
+-- part). @Left i@: the QR algorithm failed to converge.
+dgeev
+  :: Backend
+  -> Int -- ^ n
+  -> VS.Vector Double -- ^ A
+  -> IO (Either Int (VS.Vector Double, VS.Vector Double, VS.Vector Double))
+dgeev be n a = do
+  checkDim "dgeev" "A" (VS.length a) (n * n)
+  aC <- VS.thaw a
+  wr <- VSM.new n
+  wi <- VSM.new n
+  vr <- VSM.new (n * n)
+  info <-
+    VSM.unsafeWith aC $ \pa ->
+      VSM.unsafeWith wr $ \pwr ->
+        VSM.unsafeWith wi $ \pwi ->
+          VSM.unsafeWith vr $ \pvr ->
+            callDgeev (opDgeev (capOps be))
+              lapackRowMajor (castCharToCChar 'N') (castCharToCChar 'V')
+              (fromIntegral n) pa (fromIntegral n)
+              pwr pwi
+              nullPtr (fromIntegral n)
+              pvr (fromIntegral n)
+  interpretInfo "dgeev" info
+    ((,,) <$> VS.unsafeFreeze wr <*> VS.unsafeFreeze wi <*> VS.unsafeFreeze vr)
+
+-- ---------------------------------------------------------------------
+-- QR
+
+foreign import ccall safe "dynamic"
+  callDgeqrf
+    :: FunPtr (CInt -> CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt)
+    -> CInt -> CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt
+
+foreign import ccall safe "dynamic"
+  callDorgqr
+    :: FunPtr (CInt -> CInt -> CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt)
+    -> CInt -> CInt -> CInt -> CInt -> Ptr Double -> CInt -> Ptr Double -> IO CInt
+
+-- These two drivers define no positive info values, so success is the
+-- only non-throwing outcome.
+expectClean :: String -> CInt -> IO a -> IO a
+expectClean ctx info onOk
+  | info == 0 = onOk
+  | otherwise = throwIO (LapackBadArgument ctx (fromIntegral (abs info)))
+
+-- | QR factorization (@LAPACKE_dgeqrf@) of an @m x n@ matrix: returns
+-- @(QR, tau)@ where @QR@ packs @R@ in the upper triangle and the
+-- Householder reflectors below, and @tau@ has @min m n@ scalar
+-- factors. Feed both to 'dorgqr' to materialize @Q@.
+dgeqrf
+  :: Backend
+  -> Int -- ^ m
+  -> Int -- ^ n
+  -> VS.Vector Double -- ^ A
+  -> IO (VS.Vector Double, VS.Vector Double)
+dgeqrf be m n a = do
+  checkDim "dgeqrf" "A" (VS.length a) (m * n)
+  aC <- VS.thaw a
+  tau <- VSM.new (min m n)
+  info <-
+    VSM.unsafeWith aC $ \pa ->
+      VSM.unsafeWith tau $ \pt ->
+        callDgeqrf (opDgeqrf (capOps be))
+          lapackRowMajor (fromIntegral m) (fromIntegral n)
+          pa (fromIntegral n) pt
+  expectClean "dgeqrf" info
+    ((,) <$> VS.unsafeFreeze aC <*> VS.unsafeFreeze tau)
+
+-- | Materialize the first @n@ columns of @Q@ from a 'dgeqrf'
+-- factorization (@LAPACKE_dorgqr@): pass the packed @QR@ (@m x n@) and
+-- @tau@ (@k@ reflectors, @k = min m n@ from 'dgeqrf'); the result is
+-- @m x n@ with orthonormal columns.
+dorgqr
+  :: Backend
+  -> Int -- ^ m
+  -> Int -- ^ n
+  -> Int -- ^ k
+  -> VS.Vector Double -- ^ packed QR from 'dgeqrf'
+  -> VS.Vector Double -- ^ tau from 'dgeqrf'
+  -> IO (VS.Vector Double)
+dorgqr be m n k qr tau = do
+  checkDim "dorgqr" "QR" (VS.length qr) (m * n)
+  checkDim "dorgqr" "tau" (VS.length tau) k
+  qrC <- VS.thaw qr
+  tauC <- VS.thaw tau
+  info <-
+    VSM.unsafeWith qrC $ \pa ->
+      VSM.unsafeWith tauC $ \pt ->
+        callDorgqr (opDorgqr (capOps be))
+          lapackRowMajor (fromIntegral m) (fromIntegral n) (fromIntegral k)
+          pa (fromIntegral n) pt
+  expectClean "dorgqr" info (VS.unsafeFreeze qrC)
diff --git a/src/Keel/Linalg/Backend.hs b/src/Keel/Linalg/Backend.hs
new file mode 100644
--- /dev/null
+++ b/src/Keel/Linalg/Backend.hs
@@ -0,0 +1,333 @@
+-- | Locating, probing and pinning the BLAS\/LAPACK backend.
+--
+-- keel-linalg 0.1 targets OpenBLAS specifically: the library is loaded
+-- at run time (never linked), identified via @openblas_get_config@, and
+-- rejected loudly when the build is ILP64 (@USE64BITINT@) — a silently
+-- mis-matched integer width corrupts results above @2^31@ elements
+-- instead of failing. Symbol-renamed builds (e.g. the @scipy_@-prefixed
+-- or @64_@-suffixed wheels that numpy\/scipy bundle) do not resolve the
+-- standard names and are therefore rejected as 'BackendMissingSymbol' —
+-- point 'defaultBlasSpec' at a stock OpenBLAS instead.
+--
+-- All symbols are resolved eagerly at open time, so an OpenBLAS built
+-- without LAPACKE surfaces as one clear 'BackendMissingSymbol' up front
+-- instead of a crash mid-computation (the symbol-drift hazard).
+--
+-- The returned 'Backend' is an immutable pin: every operation runs
+-- against the handle you pass it, there is no global backend state and
+-- no swapping. Thread policy: unless the user has set
+-- @OPENBLAS_NUM_THREADS@ themselves, the backend is pinned to a single
+-- BLAS thread at open time — OpenBLAS's own pool fights the GHC RTS
+-- scheduler, and parallelism belongs to the caller.
+module Keel.Linalg.Backend
+  ( Backend
+  , Ops (..)
+  , BackendError (..)
+  , backendConfig
+  , defaultBlasSpec
+  , openBackend
+  , openBackendWith
+  , closeBackend
+
+    -- * Probes (exposed for tests and doctor)
+  , isILP64Config
+  ) where
+
+import Control.Exception (Exception)
+import Data.List (isInfixOf)
+import Foreign.C.String (CString, peekCString)
+import Foreign.C.Types (CChar, CInt (..))
+import Foreign.Ptr (FunPtr, Ptr)
+import System.Environment (lookupEnv)
+import System.Info (os)
+
+import Keel.Dyn
+import Keel.Dyn.Locate
+
+-- | C signatures of the resolved operations (CBLAS\/LAPACKE calling
+-- conventions, LP64 integers, @char@ mode arguments). Callers go
+-- through "Keel.Linalg"; the record is exposed so that layer can live
+-- in a separate module.
+data Ops = Ops
+  { opDdot :: FunPtr (CInt -> Ptr Double -> CInt -> Ptr Double -> CInt -> IO Double)
+  , opDgemm
+      :: FunPtr
+           (  CInt -> CInt -> CInt          -- order, transA, transB
+           -> CInt -> CInt -> CInt          -- m, n, k
+           -> Double -> Ptr Double -> CInt  -- alpha, A, lda
+           -> Ptr Double -> CInt            -- B, ldb
+           -> Double -> Ptr Double -> CInt  -- beta, C, ldc
+           -> IO ()
+           )
+  , opDgesv
+      :: FunPtr
+           (  CInt -> CInt -> CInt          -- layout, n, nrhs
+           -> Ptr Double -> CInt            -- A (overwritten with LU), lda
+           -> Ptr CInt                      -- ipiv
+           -> Ptr Double -> CInt            -- B (overwritten with X), ldb
+           -> IO CInt                       -- info
+           )
+  , opDposv
+      :: FunPtr
+           (  CInt -> CChar -> CInt -> CInt -- layout, uplo, n, nrhs
+           -> Ptr Double -> CInt            -- A (overwritten with factor), lda
+           -> Ptr Double -> CInt            -- B (overwritten with X), ldb
+           -> IO CInt
+           )
+  , opDgels
+      :: FunPtr
+           (  CInt -> CChar                 -- layout, trans
+           -> CInt -> CInt -> CInt          -- m, n, nrhs
+           -> Ptr Double -> CInt            -- A (overwritten with QR/LQ), lda
+           -> Ptr Double -> CInt            -- B (overwritten with X), ldb
+           -> IO CInt
+           )
+  , opDtrtrs
+      :: FunPtr
+           (  CInt -> CChar -> CChar -> CChar -- layout, uplo, trans, diag
+           -> CInt -> CInt                    -- n, nrhs
+           -> Ptr Double -> CInt              -- A (read-only), lda
+           -> Ptr Double -> CInt              -- B (overwritten with X), ldb
+           -> IO CInt
+           )
+  , opDgetrf
+      :: FunPtr
+           (  CInt -> CInt -> CInt          -- layout, m, n
+           -> Ptr Double -> CInt            -- A (overwritten with LU), lda
+           -> Ptr CInt                      -- ipiv
+           -> IO CInt
+           )
+  , opDgetri
+      :: FunPtr
+           (  CInt -> CInt                  -- layout, n
+           -> Ptr Double -> CInt            -- A (LU in, inverse out), lda
+           -> Ptr CInt                      -- ipiv from dgetrf
+           -> IO CInt
+           )
+  , opDpotrf
+      :: FunPtr
+           (  CInt -> CChar -> CInt         -- layout, uplo, n
+           -> Ptr Double -> CInt            -- A (overwritten with factor), lda
+           -> IO CInt
+           )
+  , opDpotri
+      :: FunPtr
+           (  CInt -> CChar -> CInt         -- layout, uplo, n
+           -> Ptr Double -> CInt            -- A (factor in, inverse out), lda
+           -> IO CInt
+           )
+  , opDgesdd
+      :: FunPtr
+           (  CInt -> CChar                 -- layout, jobz
+           -> CInt -> CInt                  -- m, n
+           -> Ptr Double -> CInt            -- A (destroyed), lda
+           -> Ptr Double                    -- s
+           -> Ptr Double -> CInt            -- U, ldu
+           -> Ptr Double -> CInt            -- VT, ldvt
+           -> IO CInt
+           )
+  , opDgesvd
+      :: FunPtr
+           (  CInt -> CChar -> CChar        -- layout, jobu, jobvt
+           -> CInt -> CInt                  -- m, n
+           -> Ptr Double -> CInt            -- A (destroyed), lda
+           -> Ptr Double                    -- s
+           -> Ptr Double -> CInt            -- U, ldu
+           -> Ptr Double -> CInt            -- VT, ldvt
+           -> Ptr Double                    -- superb workspace
+           -> IO CInt
+           )
+  , opDsyevd
+      :: FunPtr
+           (  CInt -> CChar -> CChar        -- layout, jobz, uplo
+           -> CInt                          -- n
+           -> Ptr Double -> CInt            -- A (in sym, out eigenvectors), lda
+           -> Ptr Double                    -- w (eigenvalues ascending)
+           -> IO CInt
+           )
+  , opDgeev
+      :: FunPtr
+           (  CInt -> CChar -> CChar        -- layout, jobvl, jobvr
+           -> CInt                          -- n
+           -> Ptr Double -> CInt            -- A (destroyed), lda
+           -> Ptr Double -> Ptr Double      -- wr, wi
+           -> Ptr Double -> CInt            -- VL, ldvl
+           -> Ptr Double -> CInt            -- VR, ldvr
+           -> IO CInt
+           )
+  , opDgeqrf
+      :: FunPtr
+           (  CInt -> CInt -> CInt          -- layout, m, n
+           -> Ptr Double -> CInt            -- A (out: packed QR), lda
+           -> Ptr Double                    -- tau
+           -> IO CInt
+           )
+  , opDorgqr
+      :: FunPtr
+           (  CInt -> CInt -> CInt -> CInt  -- layout, m, n, k
+           -> Ptr Double -> CInt            -- A (packed in, Q out), lda
+           -> Ptr Double                    -- tau
+           -> IO CInt
+           )
+  , opLapackeLibrary :: Maybe Library
+    -- ^ When LAPACKE lives in a separate shared library (Debian splits
+    -- OpenBLAS's CBLAS from @liblapacke@), the handle is kept here so
+    -- the resolved 'FunPtr's stay valid and 'closeBackend' can release
+    -- it; 'Nothing' when the main library carried LAPACKE itself.
+  , opSetNumThreads :: Maybe (FunPtr (CInt -> IO ()))
+    -- ^ @openblas_set_num_threads@ — optional so its absence degrades
+    -- only thread control, not the backend.
+  }
+
+-- | An immutably pinned OpenBLAS backend: the library handle, the
+-- @openblas_get_config@ string as the version tag, and the resolved
+-- operations.
+type Backend = Capability Ops
+
+-- | The backend's @openblas_get_config@ string, e.g.
+-- @\"OpenBLAS 0.3.30 DYNAMIC_ARCH NO_AFFINITY Cooperlake MAX_THREADS=64\"@.
+backendConfig :: Backend -> String
+backendConfig = capVersion
+
+-- | Why a backend could not be opened.
+data BackendError
+  = BackendNotFound DynError
+    -- ^ No library was found by the search policy.
+  | BackendNotOpenBLAS FilePath
+    -- ^ The library loaded but exports no @openblas_get_config@ —
+    -- keel-linalg 0.1 refuses to run un-probeable backends.
+  | BackendILP64 String
+    -- ^ The build is ILP64 (@USE64BITINT@ in the config string); these
+    -- bindings use 32-bit integers and would corrupt silently.
+  | BackendMissingSymbol DynError
+    -- ^ A required CBLAS\/LAPACKE symbol is absent (symbol-renamed or
+    -- LAPACKE-less builds land here).
+  deriving (Eq, Show)
+
+instance Exception BackendError
+
+-- | Where 'openBackend' looks: @KEEL_OPENBLAS@ override, the per-user
+-- keel data dir under the name @openblas@, then the system search path
+-- with the platform's stock library names.
+defaultBlasSpec :: LibrarySpec
+defaultBlasSpec =
+  LibrarySpec
+    { specName = "openblas"
+    , specEnvVar = "KEEL_OPENBLAS"
+    , specCandidates = case os of
+        "mingw32" -> ["libopenblas.dll", "openblas.dll"]
+        "darwin" -> ["libopenblas.dylib", "libopenblas.0.dylib"]
+        _ -> ["libopenblas.so.0", "libopenblas.so"]
+    }
+
+-- | @True@ when the config string names an ILP64 build.
+isILP64Config :: String -> Bool
+isILP64Config = ("USE64BITINT" `isInfixOf`)
+
+foreign import ccall unsafe "dynamic"
+  callGetConfig :: FunPtr (IO CString) -> IO CString
+
+foreign import ccall unsafe "dynamic"
+  callSetNumThreads :: FunPtr (CInt -> IO ()) -> CInt -> IO ()
+
+-- | 'openBackendWith' 'defaultBlasSpec'.
+openBackend :: IO (Either BackendError Backend)
+openBackend = openBackendWith defaultBlasSpec
+
+-- | Locate, probe and pin an OpenBLAS backend. See the module header
+-- for the probe and thread policy.
+openBackendWith :: LibrarySpec -> IO (Either BackendError Backend)
+openBackendWith spec = do
+  located <- locateLibrary spec
+  case located of
+    Left e -> pure (Left (BackendNotFound e))
+    Right loc -> do
+      let lib = locLibrary loc
+      cfgSym <- resolveSym lib "openblas_get_config"
+      case cfgSym of
+        Left _ -> do
+          closeLibrary lib
+          pure (Left (BackendNotOpenBLAS (libraryPath lib)))
+        Right cfgFp -> do
+          cfg <- peekCString =<< callGetConfig cfgFp
+          if isILP64Config cfg
+            then do
+              closeLibrary lib
+              pure (Left (BackendILP64 cfg))
+            else assemble lib cfg
+
+assemble :: Library -> String -> IO (Either BackendError Backend)
+assemble lib cfg = do
+  -- LAPACKE may live in the main library (stock OpenBLAS builds) or in
+  -- a separate liblapacke (Debian splits the packaging). Probe the main
+  -- library; on a miss, load liblapacke and resolve the drivers there.
+  -- The plain "liblapacke.so.3" name is the LP64 build by Debian's own
+  -- naming (the ILP64 variant is liblapacke64).
+  probe <- resolveSym lib "LAPACKE_dgesv" :: IO (Either DynError (FunPtr ()))
+  (lapackeLib, lapackeSource) <- case probe of
+    Right _ -> pure (Nothing, lib)
+    Left _ -> do
+      alt <- loadFirst ["liblapacke.so.3", "liblapacke.so"]
+      pure $ case alt of
+        Just l2 -> (Just l2, l2)
+        Nothing -> (Nothing, lib) -- resolution below reports the miss
+  let (<***>) :: IO (Either DynError (FunPtr a -> b)) -> String -> IO (Either DynError b)
+      (<***>) = resolveFrom lapackeSource
+  ops <-
+    Ops
+      <$$> "cblas_ddot"
+      <**> "cblas_dgemm"
+      <***> "LAPACKE_dgesv"
+      <***> "LAPACKE_dposv"
+      <***> "LAPACKE_dgels"
+      <***> "LAPACKE_dtrtrs"
+      <***> "LAPACKE_dgetrf"
+      <***> "LAPACKE_dgetri"
+      <***> "LAPACKE_dpotrf"
+      <***> "LAPACKE_dpotri"
+      <***> "LAPACKE_dgesdd"
+      <***> "LAPACKE_dgesvd"
+      <***> "LAPACKE_dsyevd"
+      <***> "LAPACKE_dgeev"
+      <***> "LAPACKE_dgeqrf"
+      <***> "LAPACKE_dorgqr"
+  threadsM <- resolveOptional lib "openblas_set_num_threads"
+  case (\f -> f lapackeLib threadsM) <$> ops of
+    Left e -> do
+      mapM_ closeLibrary lapackeLib
+      closeLibrary lib
+      pure (Left (BackendMissingSymbol e))
+    Right ops' -> do
+      userSet <- lookupEnv "OPENBLAS_NUM_THREADS"
+      case (userSet, opSetNumThreads ops') of
+        (Nothing, Just fp) -> callSetNumThreads fp 1
+        _ -> pure ()
+      pure (Right (Capability lib cfg ops'))
+  where
+    loadFirst :: [FilePath] -> IO (Maybe Library)
+    loadFirst [] = pure Nothing
+    loadFirst (n : ns) =
+      loadLibrary n >>= either (const (loadFirst ns)) (pure . Just)
+    -- applicative resolution over Either DynError, keeping the first
+    -- missing symbol's name in the error (explicit signatures: GHC2021's
+    -- MonoLocalBinds would otherwise monomorphise the FunPtr type).
+    -- <**> resolves from the main library; <***> (defined in the do
+    -- block, closing over lapackeSource) from wherever LAPACKE lives.
+    -- No fixity declarations: all three default to infixl 9, so the
+    -- chain associates left at one level.
+    (<$$>) :: (FunPtr a -> b) -> String -> IO (Either DynError b)
+    f <$$> name = fmap (fmap f) (resolveSym lib name)
+    resolveFrom :: Library -> IO (Either DynError (FunPtr a -> b)) -> String -> IO (Either DynError b)
+    resolveFrom src mf name = do
+      f <- mf
+      x <- resolveSym src name
+      pure (f <*> x)
+    (<**>) :: IO (Either DynError (FunPtr a -> b)) -> String -> IO (Either DynError b)
+    (<**>) = resolveFrom lib
+
+-- | Drop the pin (both libraries when LAPACKE was split out). All
+-- operations on this 'Backend' become invalid.
+closeBackend :: Backend -> IO ()
+closeBackend be = do
+  mapM_ closeLibrary (opLapackeLibrary (capOps be))
+  closeLibrary (capLibrary be)
diff --git a/test/Oracle.hs b/test/Oracle.hs
new file mode 100644
--- /dev/null
+++ b/test/Oracle.hs
@@ -0,0 +1,329 @@
+-- | Numerical oracle: cross-check keel-linalg results against numpy
+-- (LAPACK\/BLAS via a different build) out of process, on deterministic
+-- fixed-seed inputs.
+--
+-- Well-conditioned cases gate at 1e-10 relative error; the Hilbert 8x8
+-- system (cond ~1e10) gates on backward error instead — forward error
+-- is condition-limited across libraries, backward error is what a
+-- correct LU solve guarantees regardless of conditioning.
+--
+-- Needs python+numpy (reference values) and an OpenBLAS backend; either
+-- missing => SKIP unless @KEEL_LINALG_REQUIRE_ORACLE@ is set (CI does).
+module Main (main) where
+
+import Control.Exception (IOException, try)
+import Control.Monad (forM_, unless)
+import Data.Bits (shiftR)
+import Data.Vector.Storable qualified as VS
+import Data.Word (Word64)
+import System.Environment (lookupEnv)
+import System.Exit (ExitCode (..))
+import System.Process (readProcessWithExitCode)
+
+import Keel.Linalg
+import TestBackend (withTestBackend)
+
+expect :: Bool -> String -> IO ()
+expect ok msg = unless ok (fail msg)
+
+-- ---------------------------------------------------------------------
+-- Deterministic input generation
+
+lcg :: Word64 -> Word64
+lcg x = 6364136223846793005 * x + 1442695040888963407
+
+-- | @n@ doubles in [-1, 1).
+randDoubles :: Word64 -> Int -> [Double]
+randDoubles seed n = take n (map toD (drop 1 (iterate lcg seed)))
+  where
+    toD w = fromIntegral (w `shiftR` 11) / 4503599627370496 - 1 -- 2^52
+
+-- | Add @n@ to the diagonal of a row-major @n x n@ matrix.
+diagBoost :: Int -> VS.Vector Double -> VS.Vector Double
+diagBoost n =
+  VS.imap (\i v -> if i `div` n == i `mod` n then v + fromIntegral n else v)
+
+-- ---------------------------------------------------------------------
+-- numpy as reference implementation
+
+runNumpy :: String -> [Double] -> IO [Double]
+runNumpy script input = do
+  r <- try (readProcessWithExitCode "python" ["-c", script] (unwords (map show input)))
+        :: IO (Either IOException (ExitCode, String, String))
+  case r of
+    Right (ExitSuccess, out, _) -> pure (map read (words out))
+    _ -> fail "numpy reference run failed"
+
+numpyAvailable :: IO Bool
+numpyAvailable = do
+  r <- try (readProcessWithExitCode "python" ["-c", "import numpy"] "")
+        :: IO (Either IOException (ExitCode, String, String))
+  pure $ case r of
+    Right (ExitSuccess, _, _) -> True
+    _ -> False
+
+-- Script builders: stdin carries the flattened operands, sizes are
+-- baked into the source, output is one '%.17g' per line.
+pyHeader :: String
+pyHeader =
+  "import sys, numpy as np\n\
+  \d = np.array([float(t) for t in sys.stdin.read().split()])\n"
+
+pyEmit :: String -> String
+pyEmit expr = "print('\\n'.join('%.17g' % v for v in (" <> expr <> ").ravel()))\n"
+
+matmulScript :: Int -> Int -> Int -> String
+matmulScript m k n =
+  pyHeader
+    <> "m, k, n = " <> show m <> ", " <> show k <> ", " <> show n <> "\n"
+    <> "A = d[:m*k].reshape(m, k); B = d[m*k:].reshape(k, n)\n"
+    <> pyEmit "A @ B"
+
+solveScript :: Int -> Int -> String
+solveScript n nrhs =
+  pyHeader
+    <> "n, nrhs = " <> show n <> ", " <> show nrhs <> "\n"
+    <> "A = d[:n*n].reshape(n, n); B = d[n*n:].reshape(n, nrhs)\n"
+    <> pyEmit "np.linalg.solve(A, B)"
+
+triSolveScript :: Int -> Int -> String
+triSolveScript n nrhs =
+  pyHeader
+    <> "n, nrhs = " <> show n <> ", " <> show nrhs <> "\n"
+    <> "A = np.tril(d[:n*n].reshape(n, n)); B = d[n*n:].reshape(n, nrhs)\n"
+    <> pyEmit "np.linalg.solve(A, B)"
+
+lstsqScript :: Int -> Int -> Int -> String
+lstsqScript m n nrhs =
+  pyHeader
+    <> "m, n, nrhs = " <> show m <> ", " <> show n <> ", " <> show nrhs <> "\n"
+    <> "A = d[:m*n].reshape(m, n); B = d[m*n:].reshape(m, nrhs)\n"
+    <> pyEmit "np.linalg.lstsq(A, B, rcond=None)[0]"
+
+invScript :: Int -> String
+invScript n =
+  pyHeader
+    <> "n = " <> show n <> "\n"
+    <> pyEmit "np.linalg.inv(d.reshape(n, n))"
+
+cholScript :: Int -> String
+cholScript n =
+  pyHeader
+    <> "n = " <> show n <> "\n"
+    <> pyEmit "np.linalg.cholesky(d.reshape(n, n))"
+
+svdValsScript :: Int -> Int -> String
+svdValsScript m n =
+  pyHeader
+    <> "m, n = " <> show m <> ", " <> show n <> "\n"
+    <> pyEmit "np.linalg.svd(d.reshape(m, n), compute_uv=False)"
+
+eighValsScript :: Int -> String
+eighValsScript n =
+  pyHeader
+    <> "n = " <> show n <> "\n"
+    <> pyEmit "np.linalg.eigvalsh(d.reshape(n, n))"
+
+eigValsScript :: Int -> String
+eigValsScript n =
+  pyHeader
+    <> "n = " <> show n <> "\n"
+    <> "w = np.linalg.eigvals(d.reshape(n, n))\n"
+    <> "print('\\n'.join('%.17g\\n%.17g' % (v.real, v.imag) for v in w))\n"
+
+relErr :: Double -> Double -> Double
+relErr got ref = abs (got - ref) / max 1 (abs ref)
+
+checkAgainst :: String -> Double -> [Double] -> [Double] -> IO ()
+checkAgainst label tol got ref = do
+  expect (length got == length ref)
+    (label <> ": length " <> show (length got) <> " /= " <> show (length ref))
+  forM_ (zip3 [0 :: Int ..] got ref) $ \(i, g, r) ->
+    expect (relErr g r <= tol)
+      (label <> "[" <> show i <> "]: got " <> show g <> ", numpy " <> show r
+        <> ", rel " <> show (relErr g r))
+
+unwrap :: String -> Either Int a -> IO a
+unwrap ctx = either (\i -> fail (ctx <> ": unexpected info " <> show i)) pure
+
+-- ---------------------------------------------------------------------
+
+main :: IO ()
+main = do
+  np <- numpyAvailable
+  required <- lookupEnv "KEEL_LINALG_REQUIRE_ORACLE"
+  if not np
+    then case required of
+      Just v | v /= "" && v /= "0" ->
+        fail "KEEL_LINALG_REQUIRE_ORACLE set but python+numpy unavailable"
+      _ -> putStrLn "keel-linalg-oracle: SKIP - no python+numpy on this machine"
+    else withTestBackend "KEEL_LINALG_REQUIRE_ORACLE" run
+
+run :: Backend -> IO ()
+run be = do
+  -- 1. dgemm, random rectangular
+  let (m, k, n) = (40, 30, 20)
+      a = VS.fromList (randDoubles 42 (m * k))
+      b = VS.fromList (randDoubles 1337 (k * n))
+  ours <- dgemm be NoTrans NoTrans m n k 1 a b
+  ref <- runNumpy (matmulScript m k n) (VS.toList a <> VS.toList b)
+  checkAgainst "dgemm" 1e-10 (VS.toList ours) ref
+  putStrLn "oracle: dgemm 40x30.30x20 within 1e-10 of numpy"
+
+  -- 2. dgesv, well-conditioned random (diagonal boost)
+  let nn = 50
+      nrhs = 3
+      aSq = diagBoost nn (VS.fromList (randDoubles 7 (nn * nn)))
+      rhs = VS.fromList (randDoubles 99 (nn * nrhs))
+  solved <- unwrap "dgesv" =<< dgesv be nn nrhs aSq rhs
+  refX <- runNumpy (solveScript nn nrhs) (VS.toList aSq <> VS.toList rhs)
+  checkAgainst "dgesv" 1e-10 (VS.toList solved) refX
+  putStrLn "oracle: dgesv 50x50 within 1e-10 of numpy"
+
+  -- 3. Hilbert 8x8 (ill-conditioned): backward-error gate
+  let h = 8
+      hilbert = VS.fromList
+        [ 1 / fromIntegral (i + j + 1)
+        | i <- [0 .. h - 1], j <- [0 :: Int .. h - 1]
+        ]
+      bvec = VS.fromList
+        [ sum [1 / fromIntegral (i + j + 1) | j <- [0 :: Int .. h - 1]]
+        | i <- [0 .. h - 1]
+        ]
+  hx <- unwrap "dgesv(hilbert)" =<< dgesv be h 1 hilbert bvec
+  ax <- dgemm be NoTrans NoTrans h 1 h 1 hilbert hx
+  let residual = VS.maximum (VS.map abs (VS.zipWith (-) ax bvec))
+      backward = residual
+        / (VS.maximum (VS.map abs hilbert) * VS.maximum (VS.map abs hx)
+            + VS.maximum (VS.map abs bvec))
+  expect (backward <= 1e-12)
+    ("Hilbert backward error " <> show backward <> " > 1e-12")
+  putStrLn ("oracle: Hilbert 8x8 backward error " <> show backward <> " <= 1e-12")
+
+  -- 4. dposv on a real SPD matrix (M M^T + n I), vs numpy solve
+  let sn = 30
+      mMat = VS.fromList (randDoubles 11 (sn * sn))
+      srhs = VS.fromList (randDoubles 13 (sn * 2))
+  mmT <- dgemm be NoTrans Trans sn sn sn 1 mMat mMat
+  let spd = diagBoost sn mmT
+  psol <- unwrap "dposv" =<< dposv be Lower sn 2 spd srhs
+  pref <- runNumpy (solveScript sn 2) (VS.toList spd <> VS.toList srhs)
+  checkAgainst "dposv" 1e-10 (VS.toList psol) pref
+  putStrLn "oracle: dposv 30x30 SPD within 1e-10 of numpy"
+
+  -- 5. dgels, overdetermined 60x10, vs numpy lstsq
+  let (gm, gn, gr) = (60, 10, 2)
+      ga = VS.fromList (randDoubles 17 (gm * gn))
+      gb = VS.fromList (randDoubles 19 (gm * gr))
+  gsol <- unwrap "dgels" =<< dgels be gm gn gr ga gb
+  gref <- runNumpy (lstsqScript gm gn gr) (VS.toList ga <> VS.toList gb)
+  checkAgainst "dgels" 1e-10 (VS.toList gsol) gref
+  putStrLn "oracle: dgels 60x10 within 1e-10 of numpy lstsq"
+
+  -- 6. dgetrf + dgetri vs numpy inv, diag-boosted 40x40
+  let inn = 40
+      ia = diagBoost inn (VS.fromList (randDoubles 21 (inn * inn)))
+  (lu, piv) <- unwrap "dgetrf" =<< dgetrf be inn inn ia
+  inv <- unwrap "dgetri" =<< dgetri be inn lu piv
+  iref <- runNumpy (invScript inn) (VS.toList ia)
+  checkAgainst "dgetri" 1e-10 (VS.toList inv) iref
+  putStrLn "oracle: dgetrf+dgetri 40x40 within 1e-10 of numpy inv"
+
+  -- 7. dpotrf (lower) vs numpy cholesky — lower triangle only (the
+  -- upper triangle of our output keeps the input's bytes by contract)
+  chol <- unwrap "dpotrf" =<< dpotrf be Lower sn spd
+  cref <- runNumpy (cholScript sn) (VS.toList spd)
+  forM_ [(i, j) | i <- [0 .. sn - 1], j <- [0 .. i]] $ \(i, j) -> do
+    let g = chol VS.! (i * sn + j)
+        r = cref !! (i * sn + j)
+    expect (relErr g r <= 1e-10)
+      ("dpotrf[" <> show (i, j) <> "]: got " <> show g <> ", numpy " <> show r)
+  putStrLn "oracle: dpotrf 30x30 lower triangle within 1e-10 of numpy cholesky"
+
+  -- 8. dtrtrs on that Cholesky factor vs numpy solve over np.tril
+  let trhs = VS.fromList (randDoubles 23 sn)
+  tsol <- unwrap "dtrtrs" =<< dtrtrs be Lower NoTrans NonUnit sn 1 chol trhs
+  tref <- runNumpy (triSolveScript sn 1) (VS.toList chol <> VS.toList trhs)
+  checkAgainst "dtrtrs" 1e-10 (VS.toList tsol) tref
+  putStrLn "oracle: dtrtrs 30x30 within 1e-10 of numpy"
+
+  -- 9. SVD: singular values vs numpy (both algorithms), then a
+  -- reconstruction gate U diag(s) VT = A (sign-ambiguity-free check of
+  -- the factors themselves)
+  let (vm, vn) = (25, 15)
+      minmn = min vm vn
+      va = VS.fromList (randDoubles 29 (vm * vn))
+  (sv1, u1, vt1) <- unwrap "dgesdd" =<< dgesdd be vm vn va
+  svRef <- runNumpy (svdValsScript vm vn) (VS.toList va)
+  checkAgainst "dgesdd s" 1e-10 (VS.toList sv1) svRef
+  (sv2, _, _) <- unwrap "dgesvd" =<< dgesvd be vm vn va
+  checkAgainst "dgesvd s" 1e-10 (VS.toList sv2) svRef
+  let sVt = VS.fromList
+        [ (sv1 VS.! i) * (vt1 VS.! (i * vn + j))
+        | i <- [0 .. minmn - 1], j <- [0 .. vn - 1]
+        ]
+  recon <- dgemm be NoTrans NoTrans vm vn minmn 1 u1 sVt
+  checkAgainst "svd reconstruction" 1e-10 (VS.toList recon) (VS.toList va)
+  putStrLn "oracle: dgesdd/dgesvd 25x15 singular values + reconstruction within 1e-10"
+
+  -- 10. dsyevd: eigenvalues vs numpy eigvalsh, eigenvectors via the
+  -- residual A V = V diag(w) (signs are ambiguous, residuals are not)
+  let en = 20
+      base = VS.fromList (randDoubles 31 (en * en))
+      sym = VS.fromList
+        [ ((base VS.! (i * en + j)) + (base VS.! (j * en + i))) / 2
+        | i <- [0 .. en - 1], j <- [0 .. en - 1]
+        ]
+  (ew, ev) <- unwrap "dsyevd" =<< dsyevd be Lower en sym
+  ewRef <- runNumpy (eighValsScript en) (VS.toList sym)
+  checkAgainst "dsyevd w" 1e-10 (VS.toList ew) ewRef
+  av <- dgemm be NoTrans NoTrans en en en 1 sym ev
+  let vw = VS.imap (\idx x -> x * (ew VS.! (idx `mod` en))) ev
+  checkAgainst "dsyevd residual" 1e-10 (VS.toList av) (VS.toList vw)
+  putStrLn "oracle: dsyevd 20x20 eigenvalues + residual within 1e-10"
+
+  -- 11. dgeev: complex eigenvalues greedy-matched against numpy (order
+  -- differs between libraries; near-ties make positional compare wrong)
+  let gn2 = 12
+      gea = VS.fromList (randDoubles 37 (gn2 * gn2))
+  (wr, wi, _) <- unwrap "dgeev" =<< dgeev be gn2 gea
+  eigFlat <- runNumpy (eigValsScript gn2) (VS.toList gea)
+  let refPairs = pairUp eigFlat
+      gotPairs = zip (VS.toList wr) (VS.toList wi)
+  matchEigen refPairs gotPairs
+  putStrLn "oracle: dgeev 12x12 eigenvalues matched within 1e-10 of numpy"
+
+  -- 12. QR property gates: Q^T Q = I and Q R = A (both sign-free)
+  let (qm, qn) = (30, 12)
+      qa = VS.fromList (randDoubles 41 (qm * qn))
+  (packed, tau) <- dgeqrf be qm qn qa
+  q <- dorgqr be qm qn qn packed tau
+  qtq <- dgemm be Trans NoTrans qn qn qm 1 q q
+  let eye = [if i == j then 1 else 0 | i <- [0 .. qn - 1], j <- [0 :: Int .. qn - 1]]
+  checkAgainst "QtQ" 1e-12 (VS.toList qtq) eye
+  let r = VS.fromList
+        [ if i <= j then packed VS.! (i * qn + j) else 0
+        | i <- [0 .. qn - 1], j <- [0 .. qn - 1]
+        ]
+  qr <- dgemm be NoTrans NoTrans qm qn qn 1 q r
+  checkAgainst "QR=A" 1e-10 (VS.toList qr) (VS.toList qa)
+  putStrLn "oracle: dgeqrf/dorgqr 30x12 orthogonality + reconstruction gates passed"
+
+  putStrLn "keel-linalg-oracle: all oracle checks passed"
+
+pairUp :: [Double] -> [(Double, Double)]
+pairUp (x : y : rest) = (x, y) : pairUp rest
+pairUp _ = []
+
+-- Greedy nearest-neighbour matching of eigenvalue multisets.
+matchEigen :: [(Double, Double)] -> [(Double, Double)] -> IO ()
+matchEigen [] _ = pure ()
+matchEigen (r : rs) gs = do
+  let dists = [(dist r g, i) | (i, g) <- zip [0 :: Int ..] gs]
+      (dmin, imin) = minimum dists
+  expect (dmin <= 1e-10)
+    ("dgeev: ref eigenvalue " <> show r <> " nearest match at distance " <> show dmin)
+  let (before, after) = splitAt imin gs
+  matchEigen rs (before <> drop 1 after)
+  where
+    dist (a, b) (c, d) = sqrt ((a - c) ^ (2 :: Int) + (b - d) ^ (2 :: Int))
diff --git a/test/Smoke.hs b/test/Smoke.hs
new file mode 100644
--- /dev/null
+++ b/test/Smoke.hs
@@ -0,0 +1,149 @@
+-- | Smoke tests against a real OpenBLAS: probe machinery plus
+-- exact-value checks (small integer inputs stay exact in double
+-- precision, so equality is legitimate here). Backend discovery and the
+-- SKIP-vs-require policy live in "TestBackend".
+module Main (main) where
+
+import Control.Exception (try)
+import Control.Monad (unless)
+import Data.Vector.Storable qualified as VS
+import Foreign.C.Types (CInt (..))
+import Foreign.Ptr (FunPtr)
+import System.Environment (lookupEnv)
+
+import Keel.Dyn (capLibrary, resolveOptional)
+import Keel.Linalg
+import Keel.Linalg.Backend (isILP64Config)
+import TestBackend (withTestBackend)
+
+foreign import ccall unsafe "dynamic"
+  callGetNumThreads :: FunPtr (IO CInt) -> IO CInt
+
+expect :: Bool -> String -> IO ()
+expect ok msg = unless ok (fail msg)
+
+main :: IO ()
+main = do
+  -- unit-test the ILP64 classifier on both polarities first (no backend
+  -- needed; the local wheel is LP64 so only strings can cover this)
+  expect (isILP64Config "OpenBLAS 0.3.30 USE64BITINT DYNAMIC_ARCH NO_AFFINITY")
+    "ILP64 config not detected"
+  expect (not (isILP64Config "OpenBLAS 0.3.30 DYNAMIC_ARCH NO_AFFINITY Cooperlake"))
+    "LP64 config misclassified as ILP64"
+
+  withTestBackend "KEEL_LINALG_REQUIRE_OPENBLAS" run
+
+run :: Backend -> IO ()
+run be = do
+  putStrLn ("backend: " <> backendConfig be)
+  expect (take 8 (backendConfig be) == "OpenBLAS") "config string does not name OpenBLAS"
+
+  -- thread-pin hazard: unless the user chose a thread count themselves,
+  -- openBackend must have pinned the OpenBLAS pool to 1 (read back via
+  -- the optional openblas_get_num_threads)
+  userThreads <- lookupEnv "OPENBLAS_NUM_THREADS"
+  mGet <- resolveOptional (capLibrary be) "openblas_get_num_threads"
+  case (userThreads, mGet) of
+    (Nothing, Just fp) -> do
+      nthr <- callGetNumThreads fp
+      expect (nthr == 1) ("BLAS pool not pinned to 1 thread: " <> show nthr)
+    (Just _, _) -> putStrLn "note: OPENBLAS_NUM_THREADS set by user, pin check skipped"
+    (Nothing, Nothing) -> putStrLn "note: openblas_get_num_threads absent, pin unverifiable"
+
+  -- ddot: [1,2,3] . [4,5,6] = 32
+  d <- ddot be (VS.fromList [1, 2, 3]) (VS.fromList [4, 5, 6])
+  expect (d == 32) ("ddot: " <> show d)
+
+  -- ddot dimension mismatch must throw before calling BLAS
+  bad <- try (ddot be (VS.fromList [1, 2]) (VS.fromList [1, 2, 3]))
+        :: IO (Either LinalgError Double)
+  expect (either (const True) (const False) bad) "ddot length mismatch not rejected"
+
+  -- dgemm NoTrans/NoTrans: [[1,2],[3,4]] x [[5,6],[7,8]] = [[19,22],[43,50]]
+  c1 <- dgemm be NoTrans NoTrans 2 2 2 1 (VS.fromList [1, 2, 3, 4]) (VS.fromList [5, 6, 7, 8])
+  expect (VS.toList c1 == [19, 22, 43, 50]) ("dgemm NN: " <> show (VS.toList c1))
+
+  -- dgemm Trans on A: A stored [[1,3],[2,4]] (2x2), op(A) = [[1,2],[3,4]]
+  c2 <- dgemm be Trans NoTrans 2 2 2 1 (VS.fromList [1, 3, 2, 4]) (VS.fromList [5, 6, 7, 8])
+  expect (VS.toList c2 == [19, 22, 43, 50]) ("dgemm TN: " <> show (VS.toList c2))
+
+  -- rectangular: (1x3) x (3x1) = [[32]], alpha = 2 -> [[64]]
+  c3 <- dgemm be NoTrans NoTrans 1 1 3 2 (VS.fromList [1, 2, 3]) (VS.fromList [4, 5, 6])
+  expect (VS.toList c3 == [64]) ("dgemm alpha: " <> show (VS.toList c3))
+
+  -- dgesv: [[2,0],[0,4]] x = [3,8] -> x = [1.5, 2] (exact for powers of 2)
+  s1 <- dgesv be 2 1 (VS.fromList [2, 0, 0, 4]) (VS.fromList [3, 8])
+  case s1 of
+    Right x -> expect (VS.toList x == [1.5, 2]) ("dgesv: " <> show (VS.toList x))
+    Left i -> fail ("dgesv reported singular at " <> show i)
+
+  -- dgesv on an exactly singular matrix reports Left, not garbage
+  s2 <- dgesv be 2 1 (VS.fromList [1, 2, 2, 4]) (VS.fromList [1, 1])
+  expect (either (const True) (const False) s2) "singular matrix not reported"
+
+  -- dgesv dimension mismatch must throw before calling LAPACK
+  s3 <- try (dgesv be 2 1 (VS.fromList [1, 2, 3]) (VS.fromList [1, 1]))
+        :: IO (Either LinalgError (Either Int (VS.Vector Double)))
+  expect (either (const True) (const False) s3) "dgesv bad dims not rejected"
+
+  -- dposv on SPD [[4,2],[2,3]], b=[10,8]: x=[1.75,1.5] — near-equality,
+  -- the Cholesky pivot sqrt(2) is irrational so the result rounds
+  p1 <- unwrap "dposv" =<< dposv be Upper 2 1 (VS.fromList [4, 2, 2, 3]) (VS.fromList [10, 8])
+  expect (approxEq 1e-14 (VS.toList p1) [1.75, 1.5]) ("dposv: " <> show (VS.toList p1))
+
+  -- dposv rejects a non-positive-definite matrix
+  p2 <- dposv be Upper 2 1 (VS.fromList [-1, 0, 0, 1]) (VS.fromList [1, 1])
+  expect (either (const True) (const False) p2) "non-PD matrix not reported"
+
+  -- dtrtrs lower [[2,0],[1,1]] x = [2,3] -> x=[1,2] exact
+  t1 <- unwrap "dtrtrs" =<< dtrtrs be Lower NoTrans NonUnit 2 1
+          (VS.fromList [2, 0, 1, 1]) (VS.fromList [2, 3])
+  expect (VS.toList t1 == [1, 2]) ("dtrtrs: " <> show (VS.toList t1))
+
+  -- dgels overdetermined 3x2, consistent system -> x=[1,2]
+  g1 <- unwrap "dgels over" =<< dgels be 3 2 1
+          (VS.fromList [1, 0, 0, 1, 0, 0]) (VS.fromList [1, 2, 0])
+  expect (approxEq 1e-14 (VS.toList g1) [1, 2]) ("dgels over: " <> show (VS.toList g1))
+
+  -- dgels underdetermined 1x2 minimum norm: [[1,1]] x = [4] -> x=[2,2]
+  g2 <- unwrap "dgels under" =<< dgels be 1 2 1 (VS.fromList [1, 1]) (VS.fromList [4])
+  expect (approxEq 1e-14 (VS.toList g2) [2, 2]) ("dgels under: " <> show (VS.toList g2))
+
+  -- dgetrf + dgetri: inverse of diag(2,4) = diag(0.5,0.25) exact
+  (lu, piv) <- unwrap "dgetrf" =<< dgetrf be 2 2 (VS.fromList [2, 0, 0, 4])
+  inv <- unwrap "dgetri" =<< dgetri be 2 lu piv
+  expect (VS.toList inv == [0.5, 0, 0, 0.25]) ("dgetri: " <> show (VS.toList inv))
+
+  -- dgetrf reports exact singularity
+  f2 <- dgetrf be 2 2 (VS.fromList [1, 2, 2, 4])
+  expect (either (const True) (const False) f2) "dgetrf singularity not reported"
+
+  -- dpotrf lower of diag(4,9): factor diag(2,3), zeros preserved
+  ch <- unwrap "dpotrf" =<< dpotrf be Lower 2 (VS.fromList [4, 0, 0, 9])
+  expect (VS.toList ch == [2, 0, 0, 3]) ("dpotrf: " <> show (VS.toList ch))
+
+  -- dpotri from that factor: inverse diagonal [0.25, 1/9]
+  pin <- unwrap "dpotri" =<< dpotri be Lower 2 ch
+  expect (VS.head pin == 0.25 && abs (VS.last pin - 1 / 9) < 1e-15)
+    ("dpotri: " <> show (VS.toList pin))
+
+  -- dsyevd diag(3,1): w=[1,3] ascending; eigenvectors +/- unit basis
+  (w2, v2) <- unwrap "dsyevd" =<< dsyevd be Lower 2 (VS.fromList [3, 0, 0, 1])
+  expect (VS.toList w2 == [1, 3]) ("dsyevd w: " <> show (VS.toList w2))
+  expect (map abs (VS.toList v2) == [0, 1, 1, 0]) ("dsyevd v: " <> show (VS.toList v2))
+
+  -- dgesdd of diag(3,4): s = [4,3] descending exact
+  (sv, _, _) <- unwrap "dgesdd" =<< dgesdd be 2 2 (VS.fromList [3, 0, 0, 4])
+  expect (VS.toList sv == [4, 3]) ("dgesdd s: " <> show (VS.toList sv))
+
+  -- dgeev of the rotation [[0,-1],[1,0]]: eigenvalues +/- i
+  (wr2, wi2, _) <- unwrap "dgeev" =<< dgeev be 2 (VS.fromList [0, -1, 1, 0])
+  expect (approxEq 1e-14 (VS.toList wr2) [0, 0]) ("dgeev wr: " <> show (VS.toList wr2))
+  expect (approxEq 1e-14 (map abs (VS.toList wi2)) [1, 1]) ("dgeev wi: " <> show (VS.toList wi2))
+
+  putStrLn "keel-linalg-smoke: all checks passed against a real OpenBLAS"
+  where
+    unwrap :: String -> Either Int a -> IO a
+    unwrap ctx = either (\i -> fail (ctx <> ": unexpected info " <> show i)) pure
+    approxEq tol xs ys =
+      length xs == length ys && and (zipWith (\x y -> abs (x - y) <= tol) xs ys)
diff --git a/test/TestBackend.hs b/test/TestBackend.hs
new file mode 100644
--- /dev/null
+++ b/test/TestBackend.hs
@@ -0,0 +1,54 @@
+-- | Shared test harness: find a standard-symbol OpenBLAS and pin it.
+--
+-- Discovery order: 'openBackend' (env \/ data dir \/ system search); if
+-- that misses, a python site-packages sweep for a wheel-bundled stock
+-- OpenBLAS (faiss ships one on Windows; numpy\/scipy wheels are
+-- symbol-renamed and useless here) is tried through the @KEEL_OPENBLAS@
+-- override. No backend found => SKIP (exit 0) with the original
+-- failure printed, unless the given require-env-var is set —
+-- publish-stage CI sets it with a stock OpenBLAS installed.
+module TestBackend (withTestBackend) where
+
+import Control.Exception (IOException, finally, try)
+import System.Environment (lookupEnv, setEnv)
+import System.Exit (ExitCode (..))
+import System.Process (readProcessWithExitCode)
+
+import Keel.Linalg
+
+findWheelOpenblas :: IO (Maybe FilePath)
+findWheelOpenblas = do
+  r <- try (readProcessWithExitCode "python" ["-c", script] "")
+        :: IO (Either IOException (ExitCode, String, String))
+  pure $ case r of
+    Right (ExitSuccess, out, _) | not (null cleaned) -> Just cleaned
+      where cleaned = filter (`notElem` "\r\n") out
+    _ -> Nothing
+  where
+    script =
+      "import glob, os, sysconfig\n\
+      \sp = sysconfig.get_paths()['purelib']\n\
+      \g = glob.glob(os.path.join(sp, 'faiss_cpu.libs', 'libopenblas*.dll'))\n\
+      \print(g[0] if g else '')\n"
+
+withTestBackend :: String -> (Backend -> IO ()) -> IO ()
+withTestBackend requireVar action = do
+  first <- openBackend
+  case first of
+    Right be -> action be `finally` closeBackend be
+    Left firstErr -> do
+      wheel <- findWheelOpenblas
+      second <- case wheel of
+        Nothing -> pure Nothing
+        Just dll -> do
+          setEnv "KEEL_OPENBLAS" dll
+          either (const Nothing) Just <$> openBackend
+      required <- lookupEnv requireVar
+      case second of
+        Just be -> action be `finally` closeBackend be
+        Nothing -> case required of
+          Just v | v /= "" && v /= "0" ->
+            fail (requireVar <> " set but no usable OpenBLAS found: " <> show firstErr)
+          _ ->
+            putStrLn ("SKIP - no standard-symbol OpenBLAS on this machine ("
+              <> show firstErr <> ")")