packages feed

keel-linalg (empty) → 0.1.0.0

raw patch · 9 files changed

+1732/−0 lines, 9 filesdep +basedep +keel-dyndep +keel-linalg

Dependencies added: base, keel-dyn, keel-linalg, process, vector

Files

+ CHANGELOG.md view
@@ -0,0 +1,5 @@+# Revision history for keel-linalg++## 0.1.0.0 -- 2026-08-18++* First release.
+ LICENSE view
@@ -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.
+ README.md view
@@ -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).
+ keel-linalg.cabal view
@@ -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
+ src/Keel/Linalg.hs view
@@ -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)
+ src/Keel/Linalg/Backend.hs view
@@ -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)
+ test/Oracle.hs view
@@ -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))
+ test/Smoke.hs view
@@ -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)
+ test/TestBackend.hs view
@@ -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 <> ")")