{-# LANGUAGE ForeignFunctionInterface #-}
{-# LANGUAGE ScopedTypeVariables #-}

-- SPDX-License-Identifier: MPL-2.0

{- | Checked dense CPU kernels with explicit, garbage-collected provider handles.
Inputs and outputs use private contiguous copies. No pointer escapes to the
caller and the immutable matrices are never mutated. Auto selection tries
OpenBLAS LP64, then the pure implementation; explicit selection never silently
falls back. Provider paths are executable native libraries and must be trusted.

OpenBLAS reports its integer ABI via openblas_get_config; BLIS reports both
query and BLAS integer widths. oneMKL's process-global interface layer is set
on first explicit selection: applications sharing MKL must select this handle
before other MKL use and must not change that layer afterwards. Accelerate's
system CBLAS uses LP64. Thread setters affect the provider globally; callers
serialize reconfiguration with other uses of that provider.

LAPACKE solve/spectral dispatch is independently configured by the environment
documented in "HMath.Matrix.LAPACK". LU factors, QR, Cholesky and least squares expose
the checked pure reference algorithms; the provider affects BLAS kernels only.
Batched operations loop over GEMM calls (no vendor batch-extension requirement).
-}
module HMath.Backend (
    Provider (..),
    IntegerABI (..),
    BackendRequest (..),
    Backend,
    BackendInfo (..),
    selectBackend,
    backendInfo,
    pureBackend,
    setBackendThreads,
    lapackAvailable,
    BackendLayout (..),
    Transpose (..),
    matrixFromLayout,
    matrixToLayout,
    axpy,
    dot,
    nrm2,
    gemv,
    gemm,
    complexAxpy,
    complexDot,
    complexNrm2,
    complexGemv,
    complexGemm,
    batchedGemm,
    stridedBatchedGemm,
    backendLU,
    backendQR,
    backendSolve,
    backendCholesky,
    backendLeastSquares,
    backendSVD,
    backendSymmetricEigen,
    pureSVD,
    pureSymmetricEigen,
    lapackStatus,
) where

import Control.Exception (mask_)
import Control.Monad (foldM)
import Data.Complex (Complex ((:+)), conjugate, imagPart, realPart)
import Data.List (sortOn, transpose)
import Data.Ord (Down (..))
import Foreign hiding (rotate)
import Foreign.C
import HMath.Error
import qualified HMath.Matrix as M
import qualified HMath.Matrix.Decomposition as D
import qualified HMath.Matrix.Internal as I
import qualified HMath.Matrix.LAPACK as L
import qualified HMath.Matrix.LU as LU
import qualified HMath.Matrix.QR as QR
import System.Environment (lookupEnv)
import System.IO.Unsafe (unsafePerformIO)

-- | Pure reference, OpenBLAS, BLIS, Intel oneMKL SDL, or Apple's system framework.
data Provider = Pure | OpenBLAS | BLIS | OneMKL | Accelerate deriving (Provider -> Provider -> Bool
(Provider -> Provider -> Bool)
-> (Provider -> Provider -> Bool) -> Eq Provider
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Provider -> Provider -> Bool
== :: Provider -> Provider -> Bool
$c/= :: Provider -> Provider -> Bool
/= :: Provider -> Provider -> Bool
Eq, Eq Provider
Eq Provider =>
(Provider -> Provider -> Ordering)
-> (Provider -> Provider -> Bool)
-> (Provider -> Provider -> Bool)
-> (Provider -> Provider -> Bool)
-> (Provider -> Provider -> Bool)
-> (Provider -> Provider -> Provider)
-> (Provider -> Provider -> Provider)
-> Ord Provider
Provider -> Provider -> Bool
Provider -> Provider -> Ordering
Provider -> Provider -> Provider
forall a.
Eq a =>
(a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
$ccompare :: Provider -> Provider -> Ordering
compare :: Provider -> Provider -> Ordering
$c< :: Provider -> Provider -> Bool
< :: Provider -> Provider -> Bool
$c<= :: Provider -> Provider -> Bool
<= :: Provider -> Provider -> Bool
$c> :: Provider -> Provider -> Bool
> :: Provider -> Provider -> Bool
$c>= :: Provider -> Provider -> Bool
>= :: Provider -> Provider -> Bool
$cmax :: Provider -> Provider -> Provider
max :: Provider -> Provider -> Provider
$cmin :: Provider -> Provider -> Provider
min :: Provider -> Provider -> Provider
Ord, Int -> Provider -> ShowS
[Provider] -> ShowS
Provider -> String
(Int -> Provider -> ShowS)
-> (Provider -> String) -> ([Provider] -> ShowS) -> Show Provider
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Provider -> ShowS
showsPrec :: Int -> Provider -> ShowS
$cshow :: Provider -> String
show :: Provider -> String
$cshowList :: [Provider] -> ShowS
showList :: [Provider] -> ShowS
Show, Int -> Provider
Provider -> Int
Provider -> [Provider]
Provider -> Provider
Provider -> Provider -> [Provider]
Provider -> Provider -> Provider -> [Provider]
(Provider -> Provider)
-> (Provider -> Provider)
-> (Int -> Provider)
-> (Provider -> Int)
-> (Provider -> [Provider])
-> (Provider -> Provider -> [Provider])
-> (Provider -> Provider -> [Provider])
-> (Provider -> Provider -> Provider -> [Provider])
-> Enum Provider
forall a.
(a -> a)
-> (a -> a)
-> (Int -> a)
-> (a -> Int)
-> (a -> [a])
-> (a -> a -> [a])
-> (a -> a -> [a])
-> (a -> a -> a -> [a])
-> Enum a
$csucc :: Provider -> Provider
succ :: Provider -> Provider
$cpred :: Provider -> Provider
pred :: Provider -> Provider
$ctoEnum :: Int -> Provider
toEnum :: Int -> Provider
$cfromEnum :: Provider -> Int
fromEnum :: Provider -> Int
$cenumFrom :: Provider -> [Provider]
enumFrom :: Provider -> [Provider]
$cenumFromThen :: Provider -> Provider -> [Provider]
enumFromThen :: Provider -> Provider -> [Provider]
$cenumFromTo :: Provider -> Provider -> [Provider]
enumFromTo :: Provider -> Provider -> [Provider]
$cenumFromThenTo :: Provider -> Provider -> Provider -> [Provider]
enumFromThenTo :: Provider -> Provider -> Provider -> [Provider]
Enum)

-- | BLAS/LAPACK integer width; Haskell dimensions must additionally fit Int.
data IntegerABI = LP64 | ILP64 deriving (IntegerABI -> IntegerABI -> Bool
(IntegerABI -> IntegerABI -> Bool)
-> (IntegerABI -> IntegerABI -> Bool) -> Eq IntegerABI
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: IntegerABI -> IntegerABI -> Bool
== :: IntegerABI -> IntegerABI -> Bool
$c/= :: IntegerABI -> IntegerABI -> Bool
/= :: IntegerABI -> IntegerABI -> Bool
Eq, Eq IntegerABI
Eq IntegerABI =>
(IntegerABI -> IntegerABI -> Ordering)
-> (IntegerABI -> IntegerABI -> Bool)
-> (IntegerABI -> IntegerABI -> Bool)
-> (IntegerABI -> IntegerABI -> Bool)
-> (IntegerABI -> IntegerABI -> Bool)
-> (IntegerABI -> IntegerABI -> IntegerABI)
-> (IntegerABI -> IntegerABI -> IntegerABI)
-> Ord IntegerABI
IntegerABI -> IntegerABI -> Bool
IntegerABI -> IntegerABI -> Ordering
IntegerABI -> IntegerABI -> IntegerABI
forall a.
Eq a =>
(a -> a -> Ordering)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> Bool)
-> (a -> a -> a)
-> (a -> a -> a)
-> Ord a
$ccompare :: IntegerABI -> IntegerABI -> Ordering
compare :: IntegerABI -> IntegerABI -> Ordering
$c< :: IntegerABI -> IntegerABI -> Bool
< :: IntegerABI -> IntegerABI -> Bool
$c<= :: IntegerABI -> IntegerABI -> Bool
<= :: IntegerABI -> IntegerABI -> Bool
$c> :: IntegerABI -> IntegerABI -> Bool
> :: IntegerABI -> IntegerABI -> Bool
$c>= :: IntegerABI -> IntegerABI -> Bool
>= :: IntegerABI -> IntegerABI -> Bool
$cmax :: IntegerABI -> IntegerABI -> IntegerABI
max :: IntegerABI -> IntegerABI -> IntegerABI
$cmin :: IntegerABI -> IntegerABI -> IntegerABI
min :: IntegerABI -> IntegerABI -> IntegerABI
Ord, Int -> IntegerABI -> ShowS
[IntegerABI] -> ShowS
IntegerABI -> String
(Int -> IntegerABI -> ShowS)
-> (IntegerABI -> String)
-> ([IntegerABI] -> ShowS)
-> Show IntegerABI
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> IntegerABI -> ShowS
showsPrec :: Int -> IntegerABI -> ShowS
$cshow :: IntegerABI -> String
show :: IntegerABI -> String
$cshowList :: [IntegerABI] -> ShowS
showList :: [IntegerABI] -> ShowS
Show)

{- | Optional path overrides the provider default; Accelerate only accepts its
system framework. OpenBLAS custom symbol prefixes are unsupported; the usual
unsuffixed and 64_ suffixed builds are recognized.
-}
data BackendRequest = AutoBackend | RequireBackend Provider IntegerABI (Maybe FilePath)
    deriving (BackendRequest -> BackendRequest -> Bool
(BackendRequest -> BackendRequest -> Bool)
-> (BackendRequest -> BackendRequest -> Bool) -> Eq BackendRequest
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: BackendRequest -> BackendRequest -> Bool
== :: BackendRequest -> BackendRequest -> Bool
$c/= :: BackendRequest -> BackendRequest -> Bool
/= :: BackendRequest -> BackendRequest -> Bool
Eq, Int -> BackendRequest -> ShowS
[BackendRequest] -> ShowS
BackendRequest -> String
(Int -> BackendRequest -> ShowS)
-> (BackendRequest -> String)
-> ([BackendRequest] -> ShowS)
-> Show BackendRequest
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> BackendRequest -> ShowS
showsPrec :: Int -> BackendRequest -> ShowS
$cshow :: BackendRequest -> String
show :: BackendRequest -> String
$cshowList :: [BackendRequest] -> ShowS
showList :: [BackendRequest] -> ShowS
Show)

-- | Actual selected BLAS provider and integer-width ABI, available through backendInfo.
data BackendInfo = BackendInfo
    { BackendInfo -> Provider
selectedProvider :: Provider
    -- ^ Provider actually selected, including Pure when reference fallback is used.
    , BackendInfo -> IntegerABI
selectedIntegerABI :: IntegerABI
    -- ^ Selected BLAS integer width: LP64 uses 32-bit indices and ILP64 uses 64-bit indices.
    }
    deriving (BackendInfo -> BackendInfo -> Bool
(BackendInfo -> BackendInfo -> Bool)
-> (BackendInfo -> BackendInfo -> Bool) -> Eq BackendInfo
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: BackendInfo -> BackendInfo -> Bool
== :: BackendInfo -> BackendInfo -> Bool
$c/= :: BackendInfo -> BackendInfo -> Bool
/= :: BackendInfo -> BackendInfo -> Bool
Eq, Int -> BackendInfo -> ShowS
[BackendInfo] -> ShowS
BackendInfo -> String
(Int -> BackendInfo -> ShowS)
-> (BackendInfo -> String)
-> ([BackendInfo] -> ShowS)
-> Show BackendInfo
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> BackendInfo -> ShowS
showsPrec :: Int -> BackendInfo -> ShowS
$cshow :: BackendInfo -> String
show :: BackendInfo -> String
$cshowList :: [BackendInfo] -> ShowS
showList :: [BackendInfo] -> ShowS
Show)

data NativeBackend

{- | Opaque selected reference/native backend owner. Retain the value while using its
dynamically loaded provider.
-}
data Backend = Backend !BackendInfo !(Maybe (ForeignPtr NativeBackend))

-- | Provider identity is immutable for the lifetime of a handle.
backendInfo :: Backend -> BackendInfo
backendInfo :: Backend -> BackendInfo
backendInfo (Backend BackendInfo
info Maybe (ForeignPtr NativeBackend)
_) = BackendInfo
info

{- | Construct the deterministic reference backend without dynamically loading BLAS. Its
lifetime is garbage-collected; selected numeric reference decomposition limits still
apply.
-}
pureBackend :: Backend
pureBackend :: Backend
pureBackend = BackendInfo -> Maybe (ForeignPtr NativeBackend) -> Backend
Backend (Provider -> IntegerABI -> BackendInfo
BackendInfo Provider
Pure IntegerABI
ILP64) Maybe (ForeignPtr NativeBackend)
forall a. Maybe a
Nothing

foreign import ccall safe "hmath_backend_open"
    cOpen :: CInt -> CInt -> CString -> Ptr CInt -> IO (Ptr NativeBackend)
foreign import ccall unsafe "&hmath_backend_close" cClose :: FinalizerPtr NativeBackend
foreign import ccall safe "hmath_backend_threads" cThreads :: Ptr NativeBackend -> CInt -> IO CInt
foreign import ccall safe "hmath_backend_vector"
    cVector ::
        Ptr NativeBackend -> CInt -> Int64 -> Double -> Ptr Double -> Ptr Double -> Ptr Double -> IO CInt
foreign import ccall safe "hmath_backend_gemv"
    cGemv ::
        Ptr NativeBackend ->
        CInt ->
        Int64 ->
        Int64 ->
        Double ->
        Ptr Double ->
        Ptr Double ->
        Double ->
        Ptr Double ->
        IO CInt
foreign import ccall safe "hmath_backend_gemm"
    cGemm ::
        Ptr NativeBackend ->
        CInt ->
        CInt ->
        CInt ->
        Int64 ->
        Int64 ->
        Int64 ->
        Ptr a ->
        Ptr a ->
        Int64 ->
        Ptr a ->
        Int64 ->
        Ptr a ->
        Ptr a ->
        IO CInt
foreign import ccall safe "hmath_lapack_solve" cSolve :: CInt -> Ptr Double -> Ptr Double -> IO CInt
foreign import ccall safe "hmath_lapack_available" cLapackAvailable :: IO CInt

-- | Availability of the once-initialized independent LAPACKE loader.
lapackAvailable :: IO Bool
lapackAvailable :: IO Bool
lapackAvailable = (CInt -> CInt -> Bool
forall a. Eq a => a -> a -> Bool
/= CInt
0) (CInt -> Bool) -> IO CInt -> IO Bool
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> IO CInt
cLapackAvailable

{- | Given AutoBackend or explicit provider/integer-ABI/library request, select a reference
or native BLAS backend in IO. Auto honors HMATH_BACKEND, HMATH_BLAS_PROVIDER,
HMATH_BLAS_INTEGER_BITS and HMATH_BLAS_LIBRARY; unknown settings or unavailable
explicit providers fail. Keep the Backend alive while calling it.
-}
selectBackend :: BackendRequest -> IO (Either MathError Backend)
selectBackend :: BackendRequest -> IO (Either MathError Backend)
selectBackend BackendRequest
AutoBackend = do
    mode <- String -> IO (Maybe String)
lookupEnv String
"HMATH_BACKEND"
    provider <- lookupEnv "HMATH_BLAS_PROVIDER"
    width <- lookupEnv "HMATH_BLAS_INTEGER_BITS"
    path <- lookupEnv "HMATH_BLAS_LIBRARY"
    case (mode, width) of
        (Maybe String, Maybe String)
_
            | Bool -> (String -> Bool) -> Maybe String -> Bool
forall b a. b -> (a -> b) -> Maybe a -> b
maybe Bool
False (String -> [String] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [String
"auto", String
"pure", String
"openblas", String
"blis", String
"mkl", String
"accelerate"]) Maybe String
provider ->
                Either MathError Backend -> IO (Either MathError Backend)
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (MathError -> Either MathError Backend
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"unknown HMATH_BLAS_PROVIDER"))
        (Just String
unknown, Maybe String
_)
            | String
unknown String -> [String] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [String
"auto", String
"pure"] ->
                Either MathError Backend -> IO (Either MathError Backend)
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (MathError -> Either MathError Backend
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"HMATH_BACKEND must be auto or pure"))
        (Maybe String
_, Just String
unknown)
            | String
unknown String -> [String] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`notElem` [String
"32", String
"64"] ->
                Either MathError Backend -> IO (Either MathError Backend)
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (MathError -> Either MathError Backend
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"HMATH_BLAS_INTEGER_BITS must be 32 or 64"))
        (Maybe String, Maybe String)
_ ->
            let abi :: IntegerABI
abi = if Maybe String
width Maybe String -> Maybe String -> Bool
forall a. Eq a => a -> a -> Bool
== String -> Maybe String
forall a. a -> Maybe a
Just String
"64" then IntegerABI
ILP64 else IntegerABI
LP64
                requested :: Maybe String
requested = if Maybe String
mode Maybe String -> Maybe String -> Bool
forall a. Eq a => a -> a -> Bool
== String -> Maybe String
forall a. a -> Maybe a
Just String
"pure" then String -> Maybe String
forall a. a -> Maybe a
Just String
"pure" else Maybe String
provider
             in case Maybe String
requested of
                    Just String
"pure" -> BackendRequest -> IO (Either MathError Backend)
selectBackend (Provider -> IntegerABI -> Maybe String -> BackendRequest
RequireBackend Provider
Pure IntegerABI
abi Maybe String
forall a. Maybe a
Nothing)
                    Just String
"openblas" -> BackendRequest -> IO (Either MathError Backend)
selectBackend (Provider -> IntegerABI -> Maybe String -> BackendRequest
RequireBackend Provider
OpenBLAS IntegerABI
abi Maybe String
path)
                    Just String
"blis" -> BackendRequest -> IO (Either MathError Backend)
selectBackend (Provider -> IntegerABI -> Maybe String -> BackendRequest
RequireBackend Provider
BLIS IntegerABI
abi Maybe String
path)
                    Just String
"mkl" -> BackendRequest -> IO (Either MathError Backend)
selectBackend (Provider -> IntegerABI -> Maybe String -> BackendRequest
RequireBackend Provider
OneMKL IntegerABI
abi Maybe String
path)
                    Just String
"accelerate" -> BackendRequest -> IO (Either MathError Backend)
selectBackend (Provider -> IntegerABI -> Maybe String -> BackendRequest
RequireBackend Provider
Accelerate IntegerABI
abi Maybe String
path)
                    Just String
unknown | String
unknown String -> String -> Bool
forall a. Eq a => a -> a -> Bool
/= String
"auto" -> Either MathError Backend -> IO (Either MathError Backend)
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (MathError -> Either MathError Backend
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"unknown HMATH_BLAS_PROVIDER"))
                    Maybe String
_ -> do
                        native <- BackendRequest -> IO (Either MathError Backend)
selectBackend (Provider -> IntegerABI -> Maybe String -> BackendRequest
RequireBackend Provider
OpenBLAS IntegerABI
abi Maybe String
path)
                        pure (Right (either (const pureBackend) id native))
selectBackend (RequireBackend Provider
Pure IntegerABI
abi Maybe String
Nothing) = Either MathError Backend -> IO (Either MathError Backend)
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Backend -> Either MathError Backend
forall a b. b -> Either a b
Right (BackendInfo -> Maybe (ForeignPtr NativeBackend) -> Backend
Backend (Provider -> IntegerABI -> BackendInfo
BackendInfo Provider
Pure IntegerABI
abi) Maybe (ForeignPtr NativeBackend)
forall a. Maybe a
Nothing))
selectBackend (RequireBackend Provider
Pure IntegerABI
_ (Just String
_)) = Either MathError Backend -> IO (Either MathError Backend)
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (MathError -> Either MathError Backend
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"pure backend does not load a library path"))
selectBackend (RequireBackend Provider
_ IntegerABI
_ (Just String
path)) | Char
'\0' Char -> String -> Bool
forall a. Eq a => a -> [a] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` String
path = Either MathError Backend -> IO (Either MathError Backend)
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (MathError -> Either MathError Backend
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"backend path contains NUL"))
selectBackend (RequireBackend Provider
provider IntegerABI
abi Maybe String
path) = IO (Either MathError Backend) -> IO (Either MathError Backend)
forall a. IO a -> IO a
mask_ (IO (Either MathError Backend) -> IO (Either MathError Backend))
-> IO (Either MathError Backend) -> IO (Either MathError Backend)
forall a b. (a -> b) -> a -> b
$ String
-> (CString -> IO (Either MathError Backend))
-> IO (Either MathError Backend)
forall a. String -> (CString -> IO a) -> IO a
withCString (String -> ShowS -> Maybe String -> String
forall b a. b -> (a -> b) -> Maybe a -> b
maybe String
"" ShowS
forall a. a -> a
id Maybe String
path) ((CString -> IO (Either MathError Backend))
 -> IO (Either MathError Backend))
-> (CString -> IO (Either MathError Backend))
-> IO (Either MathError Backend)
forall a b. (a -> b) -> a -> b
$ \CString
p -> (Ptr CInt -> IO (Either MathError Backend))
-> IO (Either MathError Backend)
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr CInt -> IO (Either MathError Backend))
 -> IO (Either MathError Backend))
-> (Ptr CInt -> IO (Either MathError Backend))
-> IO (Either MathError Backend)
forall a b. (a -> b) -> a -> b
$ \Ptr CInt
status -> do
    pointer <- CInt -> CInt -> CString -> Ptr CInt -> IO (Ptr NativeBackend)
cOpen (Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Provider -> Int
forall a. Enum a => a -> Int
fromEnum Provider
provider)) (if IntegerABI
abi IntegerABI -> IntegerABI -> Bool
forall a. Eq a => a -> a -> Bool
== IntegerABI
LP64 then CInt
32 else CInt
64) CString
p Ptr CInt
status
    code <- peek status
    if pointer == nullPtr
        then
            pure
                ( Left
                    ( DomainError
                        ("backend unavailable: " ++ show provider ++ " " ++ show abi ++ " (" ++ reason code ++ ")")
                    )
                )
        else Right . Backend (BackendInfo provider abi) . Just <$> newForeignPtr cClose pointer
  where
    reason :: a -> String
reason a
2 = String
"integer ABI incompatible or cannot be verified"
    reason a
3 = String
"allocation failed"
    reason a
4 = String
"required CBLAS symbols absent"
    reason a
_ = String
"library not found or native dispatch disabled"

{- | Positive thread count; Pure supports one thread. Accelerate has no supported
runtime setter and returns an explicit error. No implicit global restoration.
-}
setBackendThreads :: Backend -> Int -> IO (Either MathError ())
setBackendThreads :: Backend -> Int -> IO (Either MathError ())
setBackendThreads Backend
b Int
n
    | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
1 Bool -> Bool -> Bool
|| Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> CInt -> Integer
forall a. Integral a => a -> Integer
toInteger (CInt
forall a. Bounded a => a
maxBound :: CInt) =
        Either MathError () -> IO (Either MathError ())
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"backend thread count must fit positive C int"))
    | Bool
otherwise = case Backend
b of
        Backend BackendInfo
_ Maybe (ForeignPtr NativeBackend)
Nothing -> Either MathError () -> IO (Either MathError ())
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 then () -> Either MathError ()
forall a b. b -> Either a b
Right () else MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"pure backend is single-threaded"))
        Backend BackendInfo
_ (Just ForeignPtr NativeBackend
fp) -> ForeignPtr NativeBackend
-> (Ptr NativeBackend -> IO (Either MathError ()))
-> IO (Either MathError ())
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr NativeBackend
fp ((Ptr NativeBackend -> IO (Either MathError ()))
 -> IO (Either MathError ()))
-> (Ptr NativeBackend -> IO (Either MathError ()))
-> IO (Either MathError ())
forall a b. (a -> b) -> a -> b
$ \Ptr NativeBackend
p -> CInt -> Either MathError ()
adapterStatus (CInt -> Either MathError ())
-> IO CInt -> IO (Either MathError ())
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Ptr NativeBackend -> CInt -> IO CInt
cThreads Ptr NativeBackend
p (Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)

adapterStatus :: CInt -> Either MathError ()
adapterStatus :: CInt -> Either MathError ()
adapterStatus CInt
0 = () -> Either MathError ()
forall a b. b -> Either a b
Right ()
adapterStatus CInt
2 = MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"backend integer extent or option invalid")
adapterStatus CInt
_ = MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"backend operation unavailable")

{- | Convert the LAPACK convention without losing the distinction between an
invalid argument and failed numerical convergence. For factorization pivots
callers may instead map a positive code to SingularMatrix (info - 1).
-}
lapackStatus :: String -> Int -> Either MathError ()
lapackStatus :: String -> Int -> Either MathError ()
lapackStatus String
_ Int
0 = () -> Either MathError ()
forall a b. b -> Either a b
Right ()
lapackStatus String
name Int
n
    | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument (String
name String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
": argument " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show (Integer -> Integer
forall a. Num a => a -> a
negate (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n))))
    | Bool
otherwise = MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
name Int
n)

{- | Physical matrix storage order: RowMajor or ColumnMajor, independent of mathematical
transpose.
-}
data BackendLayout = RowMajor | ColumnMajor deriving (BackendLayout -> BackendLayout -> Bool
(BackendLayout -> BackendLayout -> Bool)
-> (BackendLayout -> BackendLayout -> Bool) -> Eq BackendLayout
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: BackendLayout -> BackendLayout -> Bool
== :: BackendLayout -> BackendLayout -> Bool
$c/= :: BackendLayout -> BackendLayout -> Bool
/= :: BackendLayout -> BackendLayout -> Bool
Eq, Int -> BackendLayout -> ShowS
[BackendLayout] -> ShowS
BackendLayout -> String
(Int -> BackendLayout -> ShowS)
-> (BackendLayout -> String)
-> ([BackendLayout] -> ShowS)
-> Show BackendLayout
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> BackendLayout -> ShowS
showsPrec :: Int -> BackendLayout -> ShowS
$cshow :: BackendLayout -> String
show :: BackendLayout -> String
$cshowList :: [BackendLayout] -> ShowS
showList :: [BackendLayout] -> ShowS
Show)

-- | BLAS operand transformation: unchanged,transposed,or conjugate-transposed.
data Transpose = NoTranspose | Transpose | ConjugateTranspose deriving (Transpose -> Transpose -> Bool
(Transpose -> Transpose -> Bool)
-> (Transpose -> Transpose -> Bool) -> Eq Transpose
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Transpose -> Transpose -> Bool
== :: Transpose -> Transpose -> Bool
$c/= :: Transpose -> Transpose -> Bool
/= :: Transpose -> Transpose -> Bool
Eq, Int -> Transpose -> ShowS
[Transpose] -> ShowS
Transpose -> String
(Int -> Transpose -> ShowS)
-> (Transpose -> String)
-> ([Transpose] -> ShowS)
-> Show Transpose
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Transpose -> ShowS
showsPrec :: Int -> Transpose -> ShowS
$cshow :: Transpose -> String
show :: Transpose -> String
$cshowList :: [Transpose] -> ShowS
showList :: [Transpose] -> ShowS
Show)

-- | Layout conversion is a checked, explicit copy; zero shapes retain rank.
matrixFromLayout :: BackendLayout -> Int -> Int -> [a] -> Either MathError (M.Matrix a)
matrixFromLayout :: forall a.
BackendLayout -> Int -> Int -> [a] -> Either MathError (Matrix a)
matrixFromLayout BackendLayout
RowMajor Int
m Int
n = Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked Int
m Int
n
matrixFromLayout BackendLayout
ColumnMajor Int
m Int
n = (Matrix a -> Matrix a)
-> Either MathError (Matrix a) -> Either MathError (Matrix a)
forall a b. (a -> b) -> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Matrix a -> Matrix a
forall a. Matrix a -> Matrix a
M.transposeMatrix (Either MathError (Matrix a) -> Either MathError (Matrix a))
-> ([a] -> Either MathError (Matrix a))
-> [a]
-> Either MathError (Matrix a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked Int
n Int
m

{- | Given RowMajor/ColumnMajor layout then Matrix, return values in that physical order
without changing the logical matrix.
-}
matrixToLayout :: BackendLayout -> M.Matrix a -> [a]
matrixToLayout :: forall a. BackendLayout -> Matrix a -> [a]
matrixToLayout BackendLayout
RowMajor = Matrix a -> [a]
forall a. Matrix a -> [a]
M.toList
matrixToLayout BackendLayout
ColumnMajor = Matrix a -> [a]
forall a. Matrix a -> [a]
M.toList (Matrix a -> [a]) -> (Matrix a -> Matrix a) -> Matrix a -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix a -> Matrix a
forall a. Matrix a -> Matrix a
M.transposeMatrix

transCode :: Transpose -> CInt
transCode :: Transpose -> CInt
transCode Transpose
NoTranspose = CInt
111
transCode Transpose
Transpose = CInt
112
transCode Transpose
ConjugateTranspose = CInt
113
op :: (Num a) => (a -> a) -> Transpose -> M.Matrix a -> Either MathError (M.Matrix a)
op :: forall a.
Num a =>
(a -> a) -> Transpose -> Matrix a -> Either MathError (Matrix a)
op a -> a
_ Transpose
NoTranspose Matrix a
a = Matrix a -> Either MathError (Matrix a)
forall a b. b -> Either a b
Right Matrix a
a
op a -> a
_ Transpose
Transpose Matrix a
a = Matrix a -> Either MathError (Matrix a)
forall a b. b -> Either a b
Right (Matrix a -> Matrix a
forall a. Matrix a -> Matrix a
M.transposeMatrix Matrix a
a)
op a -> a
conj Transpose
ConjugateTranspose Matrix a
a = let t :: Matrix a
t = Matrix a -> Matrix a
forall a. Matrix a -> Matrix a
M.transposeMatrix Matrix a
a in (Int -> Int -> [a] -> Either MathError (Matrix a))
-> (Int, Int) -> [a] -> Either MathError (Matrix a)
forall a b c. (a -> b -> c) -> (a, b) -> c
uncurry Int -> Int -> [a] -> Either MathError (Matrix a)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
M.shape Matrix a
t) ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map a -> a
conj (Matrix a -> [a]
forall a. Matrix a -> [a]
M.toList Matrix a
t))
finite :: [Double] -> Either MathError ()
finite :: [Double] -> Either MathError ()
finite = [Double] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
I.finiteValues
finiteC :: [Complex Double] -> Either MathError ()
finiteC :: [Complex Double] -> Either MathError ()
finiteC = [Double] -> Either MathError ()
finite ([Double] -> Either MathError ())
-> ([Complex Double] -> [Double])
-> [Complex Double]
-> Either MathError ()
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Complex Double -> [Double]) -> [Complex Double] -> [Double]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (\Complex Double
z -> [Complex Double -> Double
forall a. Complex a -> a
realPart Complex Double
z, Complex Double -> Double
forall a. Complex a -> a
imagPart Complex Double
z])
extents :: Backend -> [Int] -> Either MathError ()
extents :: Backend -> [Int] -> Either MathError ()
extents Backend
b [Int]
dims
    | (Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0) [Int]
dims Bool -> Bool -> Bool
|| (Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any ((Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
limit) (Integer -> Bool) -> (Int -> Integer) -> Int -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Integer
forall a. Integral a => a -> Integer
toInteger) [Int]
dims =
        MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"backend extent exceeds selected integer ABI")
    | Bool
otherwise = () -> Either MathError ()
forall a b. b -> Either a b
Right ()
  where
    limit :: Integer
limit = if BackendInfo -> IntegerABI
selectedIntegerABI (Backend -> BackendInfo
backendInfo Backend
b) IntegerABI -> IntegerABI -> Bool
forall a. Eq a => a -> a -> Bool
== IntegerABI
LP64 then Integer
2147483647 else Int64 -> Integer
forall a. Integral a => a -> Integer
toInteger (Int64
forall a. Bounded a => a
maxBound :: Int64)
allocation :: Int -> [Int] -> Either MathError ()
allocation :: Int -> [Int] -> Either MathError ()
allocation Int
bytes [Int]
counts
    | (Int -> Bool) -> [Int] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any (\Int
n -> Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
1 Int
n) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
bytes Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int
forall a. Bounded a => a
maxBound :: Int)) [Int]
counts =
        MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"backend byte allocation overflow")
    | Bool
otherwise = () -> Either MathError ()
forall a b. b -> Either a b
Right ()
checkedList :: [Double] -> Either MathError [Double]
checkedList :: [Double] -> Either MathError [Double]
checkedList [Double]
xs = [Double] -> Either MathError ()
finite [Double]
xs Either MathError ()
-> Either MathError [Double] -> Either MathError [Double]
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right [Double]
xs
sameLength :: [a] -> [b] -> Either MathError Int
sameLength :: forall a b. [a] -> [b] -> Either MathError Int
sameLength [a]
a [b]
b
    | [a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= [b] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [b]
b = MathError -> Either MathError Int
forall a b. a -> Either a b
Left (Int -> Int -> MathError
VectorLengthMismatch ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
a) ([b] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [b]
b))
    | Bool
otherwise = Int -> Either MathError Int
forall a b. b -> Either a b
Right ([a] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [a]
a)
vectorOperation ::
    Backend -> Int -> Double -> [Double] -> [Double] -> IO (Either MathError [Double])
vectorOperation :: Backend
-> Int
-> Double
-> [Double]
-> [Double]
-> IO (Either MathError [Double])
vectorOperation Backend
b Int
kind Double
alpha [Double]
xs [Double]
ys = case Either MathError Int
validate of
    Left MathError
e -> Either MathError [Double] -> IO (Either MathError [Double])
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (MathError -> Either MathError [Double]
forall a b. a -> Either a b
Left MathError
e)
    Right Int
n -> case Backend
b of
        Backend BackendInfo
_ Maybe (ForeignPtr NativeBackend)
Nothing -> Either MathError [Double] -> IO (Either MathError [Double])
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Int -> Either MathError [Double]
forall {a}. (Eq a, Num a) => a -> Either MathError [Double]
reference Int
kind)
        Backend BackendInfo
_ (Just ForeignPtr NativeBackend
fp) | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 -> ForeignPtr NativeBackend
-> (Ptr NativeBackend -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr NativeBackend
fp ((Ptr NativeBackend -> IO (Either MathError [Double]))
 -> IO (Either MathError [Double]))
-> (Ptr NativeBackend -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. (a -> b) -> a -> b
$ \Ptr NativeBackend
p -> [Double]
-> (Ptr Double -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. Storable a => [a] -> (Ptr a -> IO b) -> IO b
withArray [Double]
xs ((Ptr Double -> IO (Either MathError [Double]))
 -> IO (Either MathError [Double]))
-> (Ptr Double -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. (a -> b) -> a -> b
$ \Ptr Double
px -> [Double]
-> (Ptr Double -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. Storable a => [a] -> (Ptr a -> IO b) -> IO b
withArray [Double]
ys ((Ptr Double -> IO (Either MathError [Double]))
 -> IO (Either MathError [Double]))
-> (Ptr Double -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. (a -> b) -> a -> b
$ \Ptr Double
py -> (Ptr Double -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr Double -> IO (Either MathError [Double]))
 -> IO (Either MathError [Double]))
-> (Ptr Double -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. (a -> b) -> a -> b
$ \Ptr Double
out -> do
            Ptr Double -> Double -> IO ()
forall a. Storable a => Ptr a -> a -> IO ()
poke Ptr Double
out Double
0
            status <- Ptr NativeBackend
-> CInt
-> Int64
-> Double
-> Ptr Double
-> Ptr Double
-> Ptr Double
-> IO CInt
cVector Ptr NativeBackend
p (Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
kind) (Int -> Int64
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) Double
alpha Ptr Double
px Ptr Double
py Ptr Double
out
            values <- if kind == 0 then peekArray n py else (: []) <$> peek out
            pure (adapterStatus status >> checkedList values)
        Backend
_ -> Either MathError [Double] -> IO (Either MathError [Double])
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Int -> Either MathError [Double]
forall {a}. (Eq a, Num a) => a -> Either MathError [Double]
reference Int
kind)
  where
    validate :: Either MathError Int
validate = do
        n <- [Double] -> [Double] -> Either MathError Int
forall a b. [a] -> [b] -> Either MathError Int
sameLength [Double]
xs [Double]
ys
        extents b [n]
        allocation 8 [n]
        finite (alpha : xs ++ ys)
        Right n
    reference :: a -> Either MathError [Double]
reference a
0 = [Double] -> Either MathError [Double]
checkedList ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
x Double
y -> Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y) [Double]
xs [Double]
ys)
    reference a
1 = [Double] -> Either MathError [Double]
checkedList [[Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
xs [Double]
ys)]
    reference a
_ = [Double] -> Either MathError [Double]
checkedList [[Double] -> Double
forall a. RealFloat a => [a] -> a
I.stableNorm [Double]
xs]

{- | Given backend,alpha,x,y, return alpha*x+y as a new Double list. Require equal lengths,
finite inputs/result and provider-supported integer extents; validation/provider
failures return Left MathError in IO.
-}
axpy :: Backend -> Double -> [Double] -> [Double] -> IO (Either MathError [Double])
axpy :: Backend
-> Double -> [Double] -> [Double] -> IO (Either MathError [Double])
axpy Backend
b = Backend
-> Int
-> Double
-> [Double]
-> [Double]
-> IO (Either MathError [Double])
vectorOperation Backend
b Int
0

{- | Real inner product. Finite input and output; accumulation can overflow even
when exact cancellation would produce a finite mathematical result.
-}
dot :: Backend -> [Double] -> [Double] -> IO (Either MathError Double)
dot :: Backend -> [Double] -> [Double] -> IO (Either MathError Double)
dot Backend
b [Double]
x [Double]
y = (Either MathError [Double] -> Either MathError Double)
-> IO (Either MathError [Double]) -> IO (Either MathError Double)
forall a b. (a -> b) -> IO a -> IO b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (([Double] -> Double)
-> Either MathError [Double] -> Either MathError Double
forall a b. (a -> b) -> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ([Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0)) (Backend
-> Int
-> Double
-> [Double]
-> [Double]
-> IO (Either MathError [Double])
vectorOperation Backend
b Int
1 Double
0 [Double]
x [Double]
y)

{- | Given backend then finite Double vector, return Euclidean norm through the selected
provider or reference path. Empty input gives zero; invalid/non-finite output or
provider bounds fail in IO.
-}
nrm2 :: Backend -> [Double] -> IO (Either MathError Double)
nrm2 :: Backend -> [Double] -> IO (Either MathError Double)
nrm2 Backend
b [Double]
x = (Either MathError [Double] -> Either MathError Double)
-> IO (Either MathError [Double]) -> IO (Either MathError Double)
forall a b. (a -> b) -> IO a -> IO b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (([Double] -> Double)
-> Either MathError [Double] -> Either MathError Double
forall a b. (a -> b) -> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ([Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0)) (Backend
-> Int
-> Double
-> [Double]
-> [Double]
-> IO (Either MathError [Double])
vectorOperation Backend
b Int
2 Double
0 [Double]
x (Int -> Double -> [Double]
forall a. Int -> a -> [a]
replicate ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
x) Double
0))

-- | y := alpha op(A) x + beta y. Finite inputs are checked even for zero alpha.
gemv ::
    Backend ->
    Transpose ->
    Double ->
    M.Matrix Double ->
    [Double] ->
    Double ->
    [Double] ->
    IO (Either MathError [Double])
gemv :: Backend
-> Transpose
-> Double
-> Matrix Double
-> [Double]
-> Double
-> [Double]
-> IO (Either MathError [Double])
gemv Backend
b Transpose
tr Double
alpha Matrix Double
a [Double]
x Double
beta [Double]
y = case Either MathError (Matrix Double)
validate of
    Left MathError
e -> Either MathError [Double] -> IO (Either MathError [Double])
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (MathError -> Either MathError [Double]
forall a b. a -> Either a b
Left MathError
e)
    Right Matrix Double
aa -> case Backend
b of
        Backend BackendInfo
_ (Just ForeignPtr NativeBackend
fp) | Int
m Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 -> ForeignPtr NativeBackend
-> (Ptr NativeBackend -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr NativeBackend
fp ((Ptr NativeBackend -> IO (Either MathError [Double]))
 -> IO (Either MathError [Double]))
-> (Ptr NativeBackend -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. (a -> b) -> a -> b
$ \Ptr NativeBackend
p -> [Double]
-> (Ptr Double -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. Storable a => [a] -> (Ptr a -> IO b) -> IO b
withArray (Matrix Double -> [Double]
forall a. Matrix a -> [a]
M.toList Matrix Double
a) ((Ptr Double -> IO (Either MathError [Double]))
 -> IO (Either MathError [Double]))
-> (Ptr Double -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. (a -> b) -> a -> b
$ \Ptr Double
pa -> [Double]
-> (Ptr Double -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. Storable a => [a] -> (Ptr a -> IO b) -> IO b
withArray [Double]
x ((Ptr Double -> IO (Either MathError [Double]))
 -> IO (Either MathError [Double]))
-> (Ptr Double -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. (a -> b) -> a -> b
$ \Ptr Double
px -> [Double]
-> (Ptr Double -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. Storable a => [a] -> (Ptr a -> IO b) -> IO b
withArray [Double]
y ((Ptr Double -> IO (Either MathError [Double]))
 -> IO (Either MathError [Double]))
-> (Ptr Double -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. (a -> b) -> a -> b
$ \Ptr Double
py -> do
            code <- Ptr NativeBackend
-> CInt
-> Int64
-> Int64
-> Double
-> Ptr Double
-> Ptr Double
-> Double
-> Ptr Double
-> IO CInt
cGemv Ptr NativeBackend
p (Transpose -> CInt
transCode Transpose
tr) (Int -> Int64
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
m) (Int -> Int64
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) Double
alpha Ptr Double
pa Ptr Double
px Double
beta Ptr Double
py
            output <- peekArray (length y) py
            pure (adapterStatus code >> checkedList output)
        Backend
_ ->
            Either MathError [Double] -> IO (Either MathError [Double])
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure ([Double] -> Either MathError [Double]
checkedList [Double
alpha Double -> Double -> Double
forall a. Num a => a -> a -> a
* [Double] -> Double
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(*) [Double]
row [Double]
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
beta Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v | ([Double]
row, Double
v) <- [[Double]] -> [Double] -> [([Double], Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip (Matrix Double -> [[Double]]
forall a. Matrix a -> [[a]]
I.entries Matrix Double
aa) [Double]
y])
  where
    (Int
m, Int
n) = Matrix Double -> (Int, Int)
forall a. Matrix a -> (Int, Int)
M.shape Matrix Double
a
    validate :: Either MathError (Matrix Double)
validate = do
        aa <- (Double -> Double)
-> Transpose -> Matrix Double -> Either MathError (Matrix Double)
forall a.
Num a =>
(a -> a) -> Transpose -> Matrix a -> Either MathError (Matrix a)
op Double -> Double
forall a. a -> a
id Transpose
tr Matrix Double
a
        let (r, c) = M.shape aa
        if length x /= c then Left (VectorLengthMismatch c (length x)) else Right ()
        if length y /= r then Left (VectorLengthMismatch r (length y)) else Right ()
        extents b [m, n]
        allocation 8 [length (M.toList a), length x, length y]
        finite (alpha : beta : M.toList a ++ x ++ y)
        Right aa

{- | C := alpha op(A) op(B) + beta C. Matrices and output must have compatible
shapes; all zero-dimensional cases use the reference arithmetic.
-}
gemm ::
    Backend ->
    Transpose ->
    Transpose ->
    Double ->
    M.Matrix Double ->
    M.Matrix Double ->
    Double ->
    M.Matrix Double ->
    IO (Either MathError (M.Matrix Double))
gemm :: Backend
-> Transpose
-> Transpose
-> Double
-> Matrix Double
-> Matrix Double
-> Double
-> Matrix Double
-> IO (Either MathError (Matrix Double))
gemm = CInt
-> (Double -> Double)
-> ([Double] -> Either MathError ())
-> Backend
-> Transpose
-> Transpose
-> Double
-> Matrix Double
-> Matrix Double
-> Double
-> Matrix Double
-> IO (Either MathError (Matrix Double))
forall a.
(Num a, Storable a) =>
CInt
-> (a -> a)
-> ([a] -> Either MathError ())
-> Backend
-> Transpose
-> Transpose
-> a
-> Matrix a
-> Matrix a
-> a
-> Matrix a
-> IO (Either MathError (Matrix a))
gemmWith CInt
0 Double -> Double
forall a. a -> a
id [Double] -> Either MathError ()
finite

aComplex :: Complex Double -> [Double]
aComplex :: Complex Double -> [Double]
aComplex Complex Double
z = [Complex Double -> Double
forall a. Complex a -> a
realPart Complex Double
z, Complex Double -> Double
forall a. Complex a -> a
imagPart Complex Double
z]

{- | Given backend,transpose A,transpose B,alpha,A,B,beta,C, return alpha*op(A)*op(B)+beta*C
as a new complex matrix. Validate compatible dimensions, finite components and provider
integer/allocation bounds.
-}
complexGemm ::
    Backend ->
    Transpose ->
    Transpose ->
    Complex Double ->
    M.Matrix (Complex Double) ->
    M.Matrix (Complex Double) ->
    Complex Double ->
    M.Matrix (Complex Double) ->
    IO (Either MathError (M.Matrix (Complex Double)))
complexGemm :: Backend
-> Transpose
-> Transpose
-> Complex Double
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Complex Double
-> Matrix (Complex Double)
-> IO (Either MathError (Matrix (Complex Double)))
complexGemm = CInt
-> (Complex Double -> Complex Double)
-> ([Complex Double] -> Either MathError ())
-> Backend
-> Transpose
-> Transpose
-> Complex Double
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Complex Double
-> Matrix (Complex Double)
-> IO (Either MathError (Matrix (Complex Double)))
forall a.
(Num a, Storable a) =>
CInt
-> (a -> a)
-> ([a] -> Either MathError ())
-> Backend
-> Transpose
-> Transpose
-> a
-> Matrix a
-> Matrix a
-> a
-> Matrix a
-> IO (Either MathError (Matrix a))
gemmWith CInt
1 Complex Double -> Complex Double
forall a. Num a => Complex a -> Complex a
conjugate [Complex Double] -> Either MathError ()
finiteC

gemmWith ::
    forall a.
    (Num a, Storable a) =>
    CInt ->
    (a -> a) ->
    ([a] -> Either MathError ()) ->
    Backend ->
    Transpose ->
    Transpose ->
    a ->
    M.Matrix a ->
    M.Matrix a ->
    a ->
    M.Matrix a ->
    IO (Either MathError (M.Matrix a))
gemmWith :: forall a.
(Num a, Storable a) =>
CInt
-> (a -> a)
-> ([a] -> Either MathError ())
-> Backend
-> Transpose
-> Transpose
-> a
-> Matrix a
-> Matrix a
-> a
-> Matrix a
-> IO (Either MathError (Matrix a))
gemmWith CInt
complexValues a -> a
conj [a] -> Either MathError ()
check Backend
b Transpose
ta Transpose
tb a
alpha Matrix a
a Matrix a
bb a
beta Matrix a
c = case Either MathError (Matrix a, Matrix a, Int, Int, Int)
validate of
    Left MathError
e -> Either MathError (Matrix a) -> IO (Either MathError (Matrix a))
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (MathError -> Either MathError (Matrix a)
forall a b. a -> Either a b
Left MathError
e)
    Right (Matrix a
aa, Matrix a
bbb, Int
m, Int
n, Int
k) -> case Backend
b of
        Backend BackendInfo
_ (Just ForeignPtr NativeBackend
fp) | Int
m Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 Bool -> Bool -> Bool
&& Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 -> ForeignPtr NativeBackend
-> (Ptr NativeBackend -> IO (Either MathError (Matrix a)))
-> IO (Either MathError (Matrix a))
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr NativeBackend
fp ((Ptr NativeBackend -> IO (Either MathError (Matrix a)))
 -> IO (Either MathError (Matrix a)))
-> (Ptr NativeBackend -> IO (Either MathError (Matrix a)))
-> IO (Either MathError (Matrix a))
forall a b. (a -> b) -> a -> b
$ \Ptr NativeBackend
p -> a
-> (Ptr a -> IO (Either MathError (Matrix a)))
-> IO (Either MathError (Matrix a))
forall a b. Storable a => a -> (Ptr a -> IO b) -> IO b
with a
alpha ((Ptr a -> IO (Either MathError (Matrix a)))
 -> IO (Either MathError (Matrix a)))
-> (Ptr a -> IO (Either MathError (Matrix a)))
-> IO (Either MathError (Matrix a))
forall a b. (a -> b) -> a -> b
$ \Ptr a
al -> a
-> (Ptr a -> IO (Either MathError (Matrix a)))
-> IO (Either MathError (Matrix a))
forall a b. Storable a => a -> (Ptr a -> IO b) -> IO b
with a
beta ((Ptr a -> IO (Either MathError (Matrix a)))
 -> IO (Either MathError (Matrix a)))
-> (Ptr a -> IO (Either MathError (Matrix a)))
-> IO (Either MathError (Matrix a))
forall a b. (a -> b) -> a -> b
$ \Ptr a
be -> [a]
-> (Ptr a -> IO (Either MathError (Matrix a)))
-> IO (Either MathError (Matrix a))
forall a b. Storable a => [a] -> (Ptr a -> IO b) -> IO b
withArray (Matrix a -> [a]
forall a. Matrix a -> [a]
M.toList Matrix a
a) ((Ptr a -> IO (Either MathError (Matrix a)))
 -> IO (Either MathError (Matrix a)))
-> (Ptr a -> IO (Either MathError (Matrix a)))
-> IO (Either MathError (Matrix a))
forall a b. (a -> b) -> a -> b
$ \Ptr a
pa -> [a]
-> (Ptr a -> IO (Either MathError (Matrix a)))
-> IO (Either MathError (Matrix a))
forall a b. Storable a => [a] -> (Ptr a -> IO b) -> IO b
withArray (Matrix a -> [a]
forall a. Matrix a -> [a]
M.toList Matrix a
bb) ((Ptr a -> IO (Either MathError (Matrix a)))
 -> IO (Either MathError (Matrix a)))
-> (Ptr a -> IO (Either MathError (Matrix a)))
-> IO (Either MathError (Matrix a))
forall a b. (a -> b) -> a -> b
$ \Ptr a
pb -> [a]
-> (Ptr a -> IO (Either MathError (Matrix a)))
-> IO (Either MathError (Matrix a))
forall a b. Storable a => [a] -> (Ptr a -> IO b) -> IO b
withArray (Matrix a -> [a]
forall a. Matrix a -> [a]
M.toList Matrix a
c) ((Ptr a -> IO (Either MathError (Matrix a)))
 -> IO (Either MathError (Matrix a)))
-> (Ptr a -> IO (Either MathError (Matrix a)))
-> IO (Either MathError (Matrix a))
forall a b. (a -> b) -> a -> b
$ \Ptr a
pc -> do
            code <-
                Ptr NativeBackend
-> CInt
-> CInt
-> CInt
-> Int64
-> Int64
-> Int64
-> Ptr a
-> Ptr a
-> Int64
-> Ptr a
-> Int64
-> Ptr a
-> Ptr a
-> IO CInt
forall a.
Ptr NativeBackend
-> CInt
-> CInt
-> CInt
-> Int64
-> Int64
-> Int64
-> Ptr a
-> Ptr a
-> Int64
-> Ptr a
-> Int64
-> Ptr a
-> Ptr a
-> IO CInt
cGemm
                    Ptr NativeBackend
p
                    CInt
complexValues
                    (Transpose -> CInt
transCode Transpose
ta)
                    (Transpose -> CInt
transCode Transpose
tb)
                    (Int -> Int64
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
m)
                    (Int -> Int64
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n)
                    (Int -> Int64
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k)
                    Ptr a
al
                    Ptr a
pa
                    (Int -> Int64
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Int, Int) -> Int
forall a b. (a, b) -> b
snd (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
M.shape Matrix a
a)))
                    Ptr a
pb
                    (Int -> Int64
forall a b. (Integral a, Num b) => a -> b
fromIntegral ((Int, Int) -> Int
forall a b. (a, b) -> b
snd (Matrix a -> (Int, Int)
forall a. Matrix a -> (Int, Int)
M.shape Matrix a
bb)))
                    Ptr a
be
                    Ptr a
pc
            output <- peekArray (m * n) pc
            pure (adapterStatus code >> check output >> M.fromListChecked m n output)
        Backend
_ -> Either MathError (Matrix a) -> IO (Either MathError (Matrix a))
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Either MathError (Matrix a) -> IO (Either MathError (Matrix a)))
-> Either MathError (Matrix a) -> IO (Either MathError (Matrix a))
forall a b. (a -> b) -> a -> b
$ do
            productM <- Matrix a -> Matrix a -> Either MathError (Matrix a)
forall a.
Num a =>
Matrix a -> Matrix a -> Either MathError (Matrix a)
M.multiplyChecked Matrix a
aa Matrix a
bbb
            output <- M.addMatrixChecked (M.scaleMatrix alpha productM) (M.scaleMatrix beta c)
            check (M.toList output)
            Right output
  where
    validate :: Either MathError (Matrix a, Matrix a, Int, Int, Int)
validate = do
        aa <- (a -> a) -> Transpose -> Matrix a -> Either MathError (Matrix a)
forall a.
Num a =>
(a -> a) -> Transpose -> Matrix a -> Either MathError (Matrix a)
op a -> a
conj Transpose
ta Matrix a
a
        bbb <- op conj tb bb
        let (m, k) = M.shape aa; (kk, n) = M.shape bbb
        if k /= kk then Left (InnerDimensionMismatch k kk) else Right ()
        if M.shape c /= (m, n) then Left (ShapeMismatch (m, n) (M.shape c)) else Right ()
        extents b [m, n, k, snd (M.shape a), snd (M.shape bb)]
        allocation
            (sizeOf (undefined :: a))
            [length (M.toList a), length (M.toList bb), length (M.toList c)]
        check (alpha : beta : M.toList a ++ M.toList bb ++ M.toList c)
        Right (aa, bbb, m, n, k)

{- | Given backend,complex alpha,x,y, compute alpha*x+y, requiring equal lengths and finite
real/imaginary components. Return a new list or checked shape/provider/numeric error in
IO.
-}
complexAxpy ::
    Backend ->
    Complex Double ->
    [Complex Double] ->
    [Complex Double] ->
    IO (Either MathError [Complex Double])
complexAxpy :: Backend
-> Complex Double
-> [Complex Double]
-> [Complex Double]
-> IO (Either MathError [Complex Double])
complexAxpy Backend
b Complex Double
alpha [Complex Double]
x [Complex Double]
y = case [Complex Double] -> [Complex Double] -> Either MathError Int
forall a b. [a] -> [b] -> Either MathError Int
sameLength [Complex Double]
x [Complex Double]
y Either MathError Int
-> (Int
    -> Either
         MathError
         (Matrix (Complex Double), Matrix (Complex Double),
          Matrix (Complex Double)))
-> Either
     MathError
     (Matrix (Complex Double), Matrix (Complex Double),
      Matrix (Complex Double))
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= \Int
n -> (,,) (Matrix (Complex Double)
 -> Matrix (Complex Double)
 -> Matrix (Complex Double)
 -> (Matrix (Complex Double), Matrix (Complex Double),
     Matrix (Complex Double)))
-> Either MathError (Matrix (Complex Double))
-> Either
     MathError
     (Matrix (Complex Double)
      -> Matrix (Complex Double)
      -> (Matrix (Complex Double), Matrix (Complex Double),
          Matrix (Complex Double)))
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int
-> Int
-> [Complex Double]
-> Either MathError (Matrix (Complex Double))
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked Int
n Int
1 [Complex Double]
x Either
  MathError
  (Matrix (Complex Double)
   -> Matrix (Complex Double)
   -> (Matrix (Complex Double), Matrix (Complex Double),
       Matrix (Complex Double)))
-> Either MathError (Matrix (Complex Double))
-> Either
     MathError
     (Matrix (Complex Double)
      -> (Matrix (Complex Double), Matrix (Complex Double),
          Matrix (Complex Double)))
forall a b.
Either MathError (a -> b)
-> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Int
-> Int
-> [Complex Double]
-> Either MathError (Matrix (Complex Double))
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked Int
1 Int
1 [Complex Double
1] Either
  MathError
  (Matrix (Complex Double)
   -> (Matrix (Complex Double), Matrix (Complex Double),
       Matrix (Complex Double)))
-> Either MathError (Matrix (Complex Double))
-> Either
     MathError
     (Matrix (Complex Double), Matrix (Complex Double),
      Matrix (Complex Double))
forall a b.
Either MathError (a -> b)
-> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Int
-> Int
-> [Complex Double]
-> Either MathError (Matrix (Complex Double))
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked Int
n Int
1 [Complex Double]
y of
    Left MathError
e -> Either MathError [Complex Double]
-> IO (Either MathError [Complex Double])
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (MathError -> Either MathError [Complex Double]
forall a b. a -> Either a b
Left MathError
e)
    Right (Matrix (Complex Double)
a, Matrix (Complex Double)
one, Matrix (Complex Double)
c) -> (Either MathError (Matrix (Complex Double))
 -> Either MathError [Complex Double])
-> IO (Either MathError (Matrix (Complex Double)))
-> IO (Either MathError [Complex Double])
forall a b. (a -> b) -> IO a -> IO b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ((Matrix (Complex Double) -> [Complex Double])
-> Either MathError (Matrix (Complex Double))
-> Either MathError [Complex Double]
forall a b. (a -> b) -> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Matrix (Complex Double) -> [Complex Double]
forall a. Matrix a -> [a]
M.toList) (Backend
-> Transpose
-> Transpose
-> Complex Double
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Complex Double
-> Matrix (Complex Double)
-> IO (Either MathError (Matrix (Complex Double)))
complexGemm Backend
b Transpose
NoTranspose Transpose
NoTranspose Complex Double
alpha Matrix (Complex Double)
a Matrix (Complex Double)
one Complex Double
1 Matrix (Complex Double)
c)

-- | Set the Boolean to True for the Hermitian dot product (conjugate x).
complexDot ::
    Backend -> Bool -> [Complex Double] -> [Complex Double] -> IO (Either MathError (Complex Double))
complexDot :: Backend
-> Bool
-> [Complex Double]
-> [Complex Double]
-> IO (Either MathError (Complex Double))
complexDot Backend
b Bool
hermitian [Complex Double]
x [Complex Double]
y = case [Complex Double] -> [Complex Double] -> Either MathError Int
forall a b. [a] -> [b] -> Either MathError Int
sameLength [Complex Double]
x [Complex Double]
y Either MathError Int
-> (Int
    -> Either
         MathError
         (Matrix (Complex Double), Matrix (Complex Double),
          Matrix (Complex Double)))
-> Either
     MathError
     (Matrix (Complex Double), Matrix (Complex Double),
      Matrix (Complex Double))
forall a b.
Either MathError a
-> (a -> Either MathError b) -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= \Int
n ->
    (,,)
        (Matrix (Complex Double)
 -> Matrix (Complex Double)
 -> Matrix (Complex Double)
 -> (Matrix (Complex Double), Matrix (Complex Double),
     Matrix (Complex Double)))
-> Either MathError (Matrix (Complex Double))
-> Either
     MathError
     (Matrix (Complex Double)
      -> Matrix (Complex Double)
      -> (Matrix (Complex Double), Matrix (Complex Double),
          Matrix (Complex Double)))
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int
-> Int
-> [Complex Double]
-> Either MathError (Matrix (Complex Double))
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked Int
1 Int
n (if Bool
hermitian then (Complex Double -> Complex Double)
-> [Complex Double] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map Complex Double -> Complex Double
forall a. Num a => Complex a -> Complex a
conjugate [Complex Double]
x else [Complex Double]
x)
        Either
  MathError
  (Matrix (Complex Double)
   -> Matrix (Complex Double)
   -> (Matrix (Complex Double), Matrix (Complex Double),
       Matrix (Complex Double)))
-> Either MathError (Matrix (Complex Double))
-> Either
     MathError
     (Matrix (Complex Double)
      -> (Matrix (Complex Double), Matrix (Complex Double),
          Matrix (Complex Double)))
forall a b.
Either MathError (a -> b)
-> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Int
-> Int
-> [Complex Double]
-> Either MathError (Matrix (Complex Double))
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked Int
n Int
1 [Complex Double]
y
        Either
  MathError
  (Matrix (Complex Double)
   -> (Matrix (Complex Double), Matrix (Complex Double),
       Matrix (Complex Double)))
-> Either MathError (Matrix (Complex Double))
-> Either
     MathError
     (Matrix (Complex Double), Matrix (Complex Double),
      Matrix (Complex Double))
forall a b.
Either MathError (a -> b)
-> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Int
-> Int
-> [Complex Double]
-> Either MathError (Matrix (Complex Double))
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked Int
1 Int
1 [Complex Double
0] of
    Left MathError
e -> Either MathError (Complex Double)
-> IO (Either MathError (Complex Double))
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (MathError -> Either MathError (Complex Double)
forall a b. a -> Either a b
Left MathError
e)
    Right (Matrix (Complex Double)
a, Matrix (Complex Double)
bb, Matrix (Complex Double)
c) -> (Either MathError (Matrix (Complex Double))
 -> Either MathError (Complex Double))
-> IO (Either MathError (Matrix (Complex Double)))
-> IO (Either MathError (Complex Double))
forall a b. (a -> b) -> IO a -> IO b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ((Matrix (Complex Double) -> Complex Double)
-> Either MathError (Matrix (Complex Double))
-> Either MathError (Complex Double)
forall a b. (a -> b) -> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (([Complex Double] -> Int -> Complex Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
0) ([Complex Double] -> Complex Double)
-> (Matrix (Complex Double) -> [Complex Double])
-> Matrix (Complex Double)
-> Complex Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix (Complex Double) -> [Complex Double]
forall a. Matrix a -> [a]
M.toList)) (Backend
-> Transpose
-> Transpose
-> Complex Double
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Complex Double
-> Matrix (Complex Double)
-> IO (Either MathError (Matrix (Complex Double)))
complexGemm Backend
b Transpose
NoTranspose Transpose
NoTranspose Complex Double
1 Matrix (Complex Double)
a Matrix (Complex Double)
bb Complex Double
0 Matrix (Complex Double)
c)

{- | Given backend then complex vector, return real Euclidean norm of all real/imaginary
components. Non-finite input/output or provider extent errors return Left MathError in
IO.
-}
complexNrm2 :: Backend -> [Complex Double] -> IO (Either MathError Double)
complexNrm2 :: Backend -> [Complex Double] -> IO (Either MathError Double)
complexNrm2 Backend
b = Backend -> [Double] -> IO (Either MathError Double)
nrm2 Backend
b ([Double] -> IO (Either MathError Double))
-> ([Complex Double] -> [Double])
-> [Complex Double]
-> IO (Either MathError Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Complex Double -> [Double]) -> [Complex Double] -> [Double]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap Complex Double -> [Double]
aComplex

{- | Given backend,transpose selector,alpha,A,x,beta,y, return alpha*op(A)*x+beta*y as a new
complex list. ConjugateTranspose conjugates matrix values; dimensions, finite values
and provider bounds are checked.
-}
complexGemv ::
    Backend ->
    Transpose ->
    Complex Double ->
    M.Matrix (Complex Double) ->
    [Complex Double] ->
    Complex Double ->
    [Complex Double] ->
    IO (Either MathError [Complex Double])
complexGemv :: Backend
-> Transpose
-> Complex Double
-> Matrix (Complex Double)
-> [Complex Double]
-> Complex Double
-> [Complex Double]
-> IO (Either MathError [Complex Double])
complexGemv Backend
b Transpose
tr Complex Double
alpha Matrix (Complex Double)
a [Complex Double]
x Complex Double
beta [Complex Double]
y = case (,) (Matrix (Complex Double)
 -> Matrix (Complex Double)
 -> (Matrix (Complex Double), Matrix (Complex Double)))
-> Either MathError (Matrix (Complex Double))
-> Either
     MathError
     (Matrix (Complex Double)
      -> (Matrix (Complex Double), Matrix (Complex Double)))
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int
-> Int
-> [Complex Double]
-> Either MathError (Matrix (Complex Double))
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked ([Complex Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Complex Double]
x) Int
1 [Complex Double]
x Either
  MathError
  (Matrix (Complex Double)
   -> (Matrix (Complex Double), Matrix (Complex Double)))
-> Either MathError (Matrix (Complex Double))
-> Either
     MathError (Matrix (Complex Double), Matrix (Complex Double))
forall a b.
Either MathError (a -> b)
-> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Int
-> Int
-> [Complex Double]
-> Either MathError (Matrix (Complex Double))
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked ([Complex Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Complex Double]
y) Int
1 [Complex Double]
y of
    Left MathError
e -> Either MathError [Complex Double]
-> IO (Either MathError [Complex Double])
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (MathError -> Either MathError [Complex Double]
forall a b. a -> Either a b
Left MathError
e)
    Right (Matrix (Complex Double)
bb, Matrix (Complex Double)
c) -> (Either MathError (Matrix (Complex Double))
 -> Either MathError [Complex Double])
-> IO (Either MathError (Matrix (Complex Double)))
-> IO (Either MathError [Complex Double])
forall a b. (a -> b) -> IO a -> IO b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ((Matrix (Complex Double) -> [Complex Double])
-> Either MathError (Matrix (Complex Double))
-> Either MathError [Complex Double]
forall a b. (a -> b) -> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap Matrix (Complex Double) -> [Complex Double]
forall a. Matrix a -> [a]
M.toList) (Backend
-> Transpose
-> Transpose
-> Complex Double
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Complex Double
-> Matrix (Complex Double)
-> IO (Either MathError (Matrix (Complex Double)))
complexGemm Backend
b Transpose
tr Transpose
NoTranspose Complex Double
alpha Matrix (Complex Double)
a Matrix (Complex Double)
bb Complex Double
beta Matrix (Complex Double)
c)

{- | Uniform scalar/transpose options, independently checked matrix triples.
Earlier calls may execute before a later invalid item; no input is mutated.
-}
batchedGemm ::
    Backend ->
    Transpose ->
    Transpose ->
    Double ->
    Double ->
    [(M.Matrix Double, M.Matrix Double, M.Matrix Double)] ->
    IO (Either MathError [M.Matrix Double])
batchedGemm :: Backend
-> Transpose
-> Transpose
-> Double
-> Double
-> [(Matrix Double, Matrix Double, Matrix Double)]
-> IO (Either MathError [Matrix Double])
batchedGemm Backend
b Transpose
ta Transpose
tb Double
alpha Double
beta [(Matrix Double, Matrix Double, Matrix Double)]
triples = ([Either MathError (Matrix Double)]
 -> Either MathError [Matrix Double])
-> IO [Either MathError (Matrix Double)]
-> IO (Either MathError [Matrix Double])
forall a b. (a -> b) -> IO a -> IO b
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap [Either MathError (Matrix Double)]
-> Either MathError [Matrix Double]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
forall (m :: * -> *) a. Monad m => [m a] -> m [a]
sequence (((Matrix Double, Matrix Double, Matrix Double)
 -> IO (Either MathError (Matrix Double)))
-> [(Matrix Double, Matrix Double, Matrix Double)]
-> IO [Either MathError (Matrix Double)]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse (\(Matrix Double
a, Matrix Double
bb, Matrix Double
c) -> Backend
-> Transpose
-> Transpose
-> Double
-> Matrix Double
-> Matrix Double
-> Double
-> Matrix Double
-> IO (Either MathError (Matrix Double))
gemm Backend
b Transpose
ta Transpose
tb Double
alpha Matrix Double
a Matrix Double
bb Double
beta Matrix Double
c) [(Matrix Double, Matrix Double, Matrix Double)]
triples)

{- | Row-major batches of A(m,k) and B(k,n). Strides count elements, may include
padding, and may be zero for broadcast. Extents are checked with Integer before
slicing; output batches are contiguous m*n matrices. No implicit transposes.
-}
stridedBatchedGemm ::
    Backend ->
    Int ->
    Int ->
    Int ->
    Int ->
    Int ->
    [Double] ->
    Int ->
    [Double] ->
    IO (Either MathError [M.Matrix Double])
stridedBatchedGemm :: Backend
-> Int
-> Int
-> Int
-> Int
-> Int
-> [Double]
-> Int
-> [Double]
-> IO (Either MathError [Matrix Double])
stridedBatchedGemm Backend
b Int
batches Int
m Int
n Int
k Int
strideA [Double]
as Int
strideB [Double]
bs = case Either MathError [(Matrix Double, Matrix Double, Matrix Double)]
prepare of
    Left MathError
e -> Either MathError [Matrix Double]
-> IO (Either MathError [Matrix Double])
forall a. a -> IO a
forall (f :: * -> *) a. Applicative f => a -> f a
pure (MathError -> Either MathError [Matrix Double]
forall a b. a -> Either a b
Left MathError
e)
    Right [(Matrix Double, Matrix Double, Matrix Double)]
triples -> Backend
-> Transpose
-> Transpose
-> Double
-> Double
-> [(Matrix Double, Matrix Double, Matrix Double)]
-> IO (Either MathError [Matrix Double])
batchedGemm Backend
b Transpose
NoTranspose Transpose
NoTranspose Double
1 Double
0 [(Matrix Double, Matrix Double, Matrix Double)]
triples
  where
    extent :: Integer -> a -> Integer
extent Integer
size a
stride = if Int
batches Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 then Integer
0 else Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int
batches Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* a -> Integer
forall a. Integral a => a -> Integer
toInteger a
stride Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
size
    prepare :: Either MathError [(Matrix Double, Matrix Double, Matrix Double)]
prepare = do
        Backend -> [Int] -> Either MathError ()
extents Backend
b [Int
batches, Int
m, Int
n, Int
k, Int
strideA, Int
strideB]
        let sa :: Integer
sa = Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
m Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
k; sb :: Integer
sb = Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
k Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n
        if (Int
strideA Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
0 Bool -> Bool -> Bool
&& Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
strideA Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
sa) Bool -> Bool -> Bool
|| (Int
strideB Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
0 Bool -> Bool -> Bool
&& Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
strideB Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
sb)
            then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"strided batch overlaps; use zero stride for broadcast")
            else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
        if Integer -> Int -> Integer
forall {a}. Integral a => Integer -> a -> Integer
extent Integer
sa Int
strideA Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
as) Bool -> Bool -> Bool
|| Integer -> Int -> Integer
forall {a}. Integral a => Integer -> a -> Integer
extent Integer
sb Int
strideB Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Int -> Integer
forall a. Integral a => a -> Integer
toInteger ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
bs)
            then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"strided batch buffer too short")
            else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
        (Int
 -> Either MathError (Matrix Double, Matrix Double, Matrix Double))
-> [Int]
-> Either MathError [(Matrix Double, Matrix Double, Matrix Double)]
forall (t :: * -> *) (f :: * -> *) a b.
(Traversable t, Applicative f) =>
(a -> f b) -> t a -> f (t b)
forall (f :: * -> *) a b.
Applicative f =>
(a -> f b) -> [a] -> f [b]
traverse
            ( \Int
i ->
                (,,)
                    (Matrix Double
 -> Matrix Double
 -> Matrix Double
 -> (Matrix Double, Matrix Double, Matrix Double))
-> Either MathError (Matrix Double)
-> Either
     MathError
     (Matrix Double
      -> Matrix Double -> (Matrix Double, Matrix Double, Matrix Double))
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int -> Int -> [Double] -> Either MathError (Matrix Double)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked Int
m Int
k (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take (Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
sa) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
strideA) [Double]
as))
                    Either
  MathError
  (Matrix Double
   -> Matrix Double -> (Matrix Double, Matrix Double, Matrix Double))
-> Either MathError (Matrix Double)
-> Either
     MathError
     (Matrix Double -> (Matrix Double, Matrix Double, Matrix Double))
forall a b.
Either MathError (a -> b)
-> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Int -> Int -> [Double] -> Either MathError (Matrix Double)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked Int
k Int
n (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take (Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
sb) (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
drop (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
strideB) [Double]
bs))
                    Either
  MathError
  (Matrix Double -> (Matrix Double, Matrix Double, Matrix Double))
-> Either MathError (Matrix Double)
-> Either MathError (Matrix Double, Matrix Double, Matrix Double)
forall a b.
Either MathError (a -> b)
-> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Int -> Int -> Either MathError (Matrix Double)
forall a. Num a => Int -> Int -> Either MathError (Matrix a)
M.zeroMatrixChecked Int
m Int
n
            )
            [Int
0 .. Int
batches Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]

{- | Given backend then real matrix, compute checked column-independent LU using the
reference pivoted LU implementation. The backend argument does not change this
routine's provider. Shape/pivot/numeric failures propagate.
-}
backendLU :: Backend -> M.Matrix Double -> Either MathError (LU.LU Double)
backendLU :: Backend -> Matrix Double -> Either MathError (LU Double)
backendLU Backend
_ = Matrix Double -> Either MathError (LU Double)
forall a. RealFloat a => Matrix a -> Either MathError (LU a)
LU.luDecompose

{- | Given backend then real matrix, compute reference column-pivoted Householder QR;
selection does not change this routine's provider. Validation/numeric errors propagate.
-}
backendQR :: Backend -> M.Matrix Double -> Either MathError (QR.QR Double)
backendQR :: Backend -> Matrix Double -> Either MathError (QR Double)
backendQR Backend
_ = Matrix Double -> Either MathError (QR Double)
forall a. RealFloat a => Matrix a -> Either MathError (QR a)
QR.qrDecompose

{- | Given backend,square real matrix,then RHS list, solve via reference LU for Pure or the
independent LAPACKE solve bridge for a native backend. Invalid dimensions, unavailable
solver, singularity or non-finite arithmetic returns Left MathError.
-}
backendSolve :: Backend -> M.Matrix Double -> [Double] -> Either MathError [Double]
backendSolve :: Backend -> Matrix Double -> [Double] -> Either MathError [Double]
backendSolve (Backend BackendInfo
_ Maybe (ForeignPtr NativeBackend)
Nothing) = Matrix Double -> [Double] -> Either MathError [Double]
forall a. RealFloat a => Matrix a -> [a] -> Either MathError [a]
LU.solve
backendSolve Backend
_ = Matrix Double -> [Double] -> Either MathError [Double]
nativeSolve

{-# NOINLINE nativeSolve #-}
nativeSolve :: M.Matrix Double -> [Double] -> Either MathError [Double]
nativeSolve :: Matrix Double -> [Double] -> Either MathError [Double]
nativeSolve Matrix Double
a [Double]
rhs = do
    let (Int
m, Int
n) = Matrix Double -> (Int, Int)
forall a. Matrix a -> (Int, Int)
M.shape Matrix Double
a
    if Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n then MathError -> Either MathError ()
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
m, Int
n)) else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    if [Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
rhs Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (Int -> Int -> MathError
RightHandSideMismatch Int
n ([Double] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Double]
rhs)) else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    if Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> CInt -> Integer
forall a. Integral a => a -> Integer
toInteger (CInt
forall a. Bounded a => a
maxBound :: CInt)
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"LAPACKE solver index exceeds Int32")
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    Matrix Double -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
I.finiteMatrix Matrix Double
a
    [Double] -> Either MathError ()
finite [Double]
rhs
    if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
        then [Double] -> Either MathError [Double]
forall a b. b -> Either a b
Right []
        else IO (Either MathError [Double]) -> Either MathError [Double]
forall a. IO a -> a
unsafePerformIO (IO (Either MathError [Double]) -> Either MathError [Double])
-> IO (Either MathError [Double]) -> Either MathError [Double]
forall a b. (a -> b) -> a -> b
$ [Double]
-> (Ptr Double -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. Storable a => [a] -> (Ptr a -> IO b) -> IO b
withArray (Matrix Double -> [Double]
forall a. Matrix a -> [a]
M.toList Matrix Double
a) ((Ptr Double -> IO (Either MathError [Double]))
 -> IO (Either MathError [Double]))
-> (Ptr Double -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. (a -> b) -> a -> b
$ \Ptr Double
pa -> [Double]
-> (Ptr Double -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. Storable a => [a] -> (Ptr a -> IO b) -> IO b
withArray [Double]
rhs ((Ptr Double -> IO (Either MathError [Double]))
 -> IO (Either MathError [Double]))
-> (Ptr Double -> IO (Either MathError [Double]))
-> IO (Either MathError [Double])
forall a b. (a -> b) -> a -> b
$ \Ptr Double
pb -> do
            code <- CInt -> Ptr Double -> Ptr Double -> IO CInt
cSolve (Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) Ptr Double
pa Ptr Double
pb
            output <- peekArray n pb
            pure $
                if code == -10000
                    then Left (DomainError "LAPACKE solve unavailable")
                    else
                        if code == -10002
                            then Left (DomainError "LAPACKE pivot allocation failed")
                            else
                                if code > 0
                                    then Left (SingularMatrix (fromIntegral code - 1))
                                    else lapackStatus "LAPACKE solve" (fromIntegral code) >> checkedList output

{- | Given backend then finite symmetric positive-definite matrix, return reference lower
Cholesky factor. Invalid shape/symmetry/pivots/numeric results fail; provider selection
is not used by this wrapper.
-}
backendCholesky :: Backend -> M.Matrix Double -> Either MathError (M.Matrix Double)
backendCholesky :: Backend -> Matrix Double -> Either MathError (Matrix Double)
backendCholesky Backend
_ = Matrix Double -> Either MathError (Matrix Double)
forall a. RealFloat a => Matrix a -> Either MathError (Matrix a)
D.cholesky

{- | Given backend,design matrix,then RHS list, solve using reference pivoted QR. Require
compatible overdetermined/full-rank input; shape/rank/numeric failures propagate.
-}
backendLeastSquares :: Backend -> M.Matrix Double -> [Double] -> Either MathError [Double]
backendLeastSquares :: Backend -> Matrix Double -> [Double] -> Either MathError [Double]
backendLeastSquares Backend
_ = Matrix Double -> [Double] -> Either MathError [Double]
forall a. RealFloat a => Matrix a -> [a] -> Either MathError [a]
QR.leastSquares

{- | Given backend then finite real matrix, use bounded pure reference SVD for Pure or
LAPACK thin SVD for native selection. Return factors or checked
shape/budget/provider/convergence errors.
-}
backendSVD :: Backend -> M.Matrix Double -> Either MathError L.SVD
backendSVD :: Backend -> Matrix Double -> Either MathError SVD
backendSVD (Backend BackendInfo
_ Maybe (ForeignPtr NativeBackend)
Nothing) = Matrix Double -> Either MathError SVD
pureSVD
backendSVD Backend
_ = Matrix Double -> Either MathError SVD
L.svd

{- | Given backend then finite symmetric square matrix, use bounded cyclic-Jacobi reference
solver for Pure or LAPACK otherwise. Return eigenvalues/vectors or
input/provider/convergence errors.
-}
backendSymmetricEigen :: Backend -> M.Matrix Double -> Either MathError L.Eigen
backendSymmetricEigen :: Backend -> Matrix Double -> Either MathError Eigen
backendSymmetricEigen (Backend BackendInfo
_ Maybe (ForeignPtr NativeBackend)
Nothing) = Matrix Double -> Either MathError Eigen
pureSymmetricEigen
backendSymmetricEigen Backend
_ = Matrix Double -> Either MathError Eigen
L.symmetricEigen

{- | Cyclic Jacobi symmetric eigensolver, at most 100 sweeps, ascending eigenvalues.
Scaled input avoids squared-norm overflow. This reference path is limited to
order 128 and reports nonconvergence rather than returning an unchecked factor.
-}
pureSymmetricEigen :: M.Matrix Double -> Either MathError L.Eigen
pureSymmetricEigen :: Matrix Double -> Either MathError Eigen
pureSymmetricEigen Matrix Double
a = do
    let m :: Int
m = (Int, Int) -> Int
forall a b. (a, b) -> a
fst (Matrix Double -> (Int, Int)
forall a. Matrix a -> (Int, Int)
M.shape Matrix Double
a)
    if Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n then MathError -> Either MathError ()
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
m, Int
n)) else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    if Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
128 then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"pure symmetric eigen: order exceeds 128") else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    Matrix Double -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
I.finiteMatrix Matrix Double
a
    if [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
or
        [ Double -> Double
forall a. Num a => a -> a
abs (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
y) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1e-12 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1 (Double -> Double -> Double
forall a. Ord a => a -> a -> a
max (Double -> Double
forall a. Num a => a -> a
abs Double
x) (Double -> Double
forall a. Num a => a -> a
abs Double
y))
        | (Double
x, Double
y) <- [Double] -> [Double] -> [(Double, Double)]
forall a b. [a] -> [b] -> [(a, b)]
zip (Matrix Double -> [Double]
forall a. Matrix a -> [a]
M.toList Matrix Double
a) (Matrix Double -> [Double]
forall a. Matrix a -> [a]
M.toList (Matrix Double -> Matrix Double
forall a. Matrix a -> Matrix a
M.transposeMatrix Matrix Double
a))
        ]
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"pure symmetric eigen requires symmetric matrix")
        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    let scale :: Double
scale = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs (Matrix Double -> [Double]
forall a. Matrix a -> [a]
M.toList Matrix Double
a))
        s :: Double
s = if Double
scale Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
1 else Double
scale
        rows :: [[Double]]
rows = ([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
s)) (Matrix Double -> [[Double]]
forall a. Matrix a -> [[a]]
I.entries Matrix Double
a)
        ident :: [[Double]]
ident = [[if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Double
1 else Double
0 | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
    (d, v) <- Int
-> [[Double]]
-> [[Double]]
-> Either MathError ([[Double]], [[Double]])
forall {a}.
RealFloat a =>
Int -> [[a]] -> [[a]] -> Either MathError ([[a]], [[a]])
sweeps Int
0 [[Double]]
rows [[Double]]
ident
    let order = (Int -> Double) -> [Int] -> [Int]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (\Int
i -> [[Double]]
d [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
        values = [[[Double]]
d [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0 | Int
i <- [Int]
order]
    finiteC values
    vectors <- M.fromListChecked n n [v !! i !! j :+ 0 | i <- [0 .. n - 1], j <- order]
    Right (L.Eigen values vectors)
  where
    n :: Int
n = (Int, Int) -> Int
forall a b. (a, b) -> b
snd (Matrix Double -> (Int, Int)
forall a. Matrix a -> (Int, Int)
M.shape Matrix Double
a)
    pairs :: [(Int, Int)]
pairs = [(Int
p, Int
q) | Int
p <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
q <- [Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
    sweeps :: Int -> [[a]] -> [[a]] -> Either MathError ([[a]], [[a]])
sweeps Int
count [[a]]
d [[a]]
v
        | ((Int, Int) -> Bool) -> [(Int, Int)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (\(Int
p, Int
q) -> a -> a
forall a. Num a => a -> a
abs ([[a]]
d [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
p [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
q) a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
2e-15) [(Int, Int)]
pairs = ([[a]], [[a]]) -> Either MathError ([[a]], [[a]])
forall a b. b -> Either a b
Right ([[a]]
d, [[a]]
v)
        | Int
count Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
100 = MathError -> Either MathError ([[a]], [[a]])
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"pure symmetric Jacobi eigen" Int
count)
        | Bool
otherwise = let ([[a]]
dd, [[a]]
vv) = (([[a]], [[a]]) -> (Int, Int) -> ([[a]], [[a]]))
-> ([[a]], [[a]]) -> [(Int, Int)] -> ([[a]], [[a]])
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl ([[a]], [[a]]) -> (Int, Int) -> ([[a]], [[a]])
forall {a}.
RealFloat a =>
([[a]], [[a]]) -> (Int, Int) -> ([[a]], [[a]])
rotate ([[a]]
d, [[a]]
v) [(Int, Int)]
pairs in Int -> [[a]] -> [[a]] -> Either MathError ([[a]], [[a]])
sweeps (Int
count Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [[a]]
dd [[a]]
vv
    rotate :: ([[a]], [[a]]) -> (Int, Int) -> ([[a]], [[a]])
rotate state :: ([[a]], [[a]])
state@([[a]]
d, [[a]]
v) (Int
p, Int
q)
        | a -> a
forall a. Num a => a -> a
abs ([[a]]
d [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
p [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
q) a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
2e-16 = ([[a]], [[a]])
state
        | Bool
otherwise =
            let ap :: a
ap = [[a]]
d [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
p [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
p
                aq :: a
aq = [[a]]
d [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
q [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
q
                b :: a
b = [[a]]
d [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
p [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
q
                tau :: a
tau = (a
aq a -> a -> a
forall a. Num a => a -> a -> a
- a
ap) a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
2 a -> a -> a
forall a. Num a => a -> a -> a
* a
b)
                t :: a
t = if a
tau a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then a
1 else a -> a
forall a. Num a => a -> a
signum a
tau a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a -> a
forall a. Num a => a -> a
abs a
tau a -> a -> a
forall a. Num a => a -> a -> a
+ [a] -> a
forall a. RealFloat a => [a] -> a
I.stableNorm [a
1, a
tau])
                c :: a
c = a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a. Floating a => a -> a
sqrt (a
1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
t a -> a -> a
forall a. Num a => a -> a -> a
* a
t)
                s :: a
s = a
t a -> a -> a
forall a. Num a => a -> a -> a
* a
c
                cols :: [[a]] -> [[a]]
cols [[a]]
x =
                    ([a] -> [a]) -> [[a]] -> [[a]]
forall a b. (a -> b) -> [a] -> [b]
map
                        ( \[a]
row ->
                            Int -> a -> [a] -> [a]
forall a. Int -> a -> [a] -> [a]
I.replace Int
q (a
s a -> a -> a
forall a. Num a => a -> a -> a
* ([a]
row [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
p) a -> a -> a
forall a. Num a => a -> a -> a
+ a
c a -> a -> a
forall a. Num a => a -> a -> a
* ([a]
row [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
q)) (Int -> a -> [a] -> [a]
forall a. Int -> a -> [a] -> [a]
I.replace Int
p (a
c a -> a -> a
forall a. Num a => a -> a -> a
* ([a]
row [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
p) a -> a -> a
forall a. Num a => a -> a -> a
- a
s a -> a -> a
forall a. Num a => a -> a -> a
* ([a]
row [a] -> Int -> a
forall a. HasCallStack => [a] -> Int -> a
!! Int
q)) [a]
row)
                        )
                        [[a]]
x
                dd :: [[a]]
dd = [[a]] -> [[a]]
forall a. [[a]] -> [[a]]
transpose ([[a]] -> [[a]]
cols ([[a]] -> [[a]]
forall a. [[a]] -> [[a]]
transpose ([[a]] -> [[a]]
cols [[a]]
d)))
                clean :: [[a]]
clean = Int -> [a] -> [[a]] -> [[a]]
forall a. Int -> a -> [a] -> [a]
I.replace Int
q (Int -> a -> [a] -> [a]
forall a. Int -> a -> [a] -> [a]
I.replace Int
p a
0 ([[a]]
dd [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
q)) (Int -> [a] -> [[a]] -> [[a]]
forall a. Int -> a -> [a] -> [a]
I.replace Int
p (Int -> a -> [a] -> [a]
forall a. Int -> a -> [a] -> [a]
I.replace Int
q a
0 ([[a]]
dd [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
p)) [[a]]
dd)
             in ([[a]]
clean, [[a]] -> [[a]]
cols [[a]]
v)

{- | Thin real SVD by one-sided Jacobi orthogonalization, at most 100 sweeps,
dimensions <=128. Unlike a normal-equation SVD it does not form A^T A.
A deterministic orthogonal completion supplies zero-singular-value columns.
-}
pureSVD :: M.Matrix Double -> Either MathError L.SVD
pureSVD :: Matrix Double -> Either MathError SVD
pureSVD Matrix Double
a = do
    Matrix Double -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
I.finiteMatrix Matrix Double
a
    let m :: Int
m = (Int, Int) -> Int
forall a b. (a, b) -> a
fst (Matrix Double -> (Int, Int)
forall a. Matrix a -> (Int, Int)
M.shape Matrix Double
a)
    if Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
m Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
128 then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"pure SVD: dimensions exceed 128") else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    if Int
m Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
n
        then do
            f <- Matrix Double -> Either MathError SVD
pureSVD (Matrix Double -> Matrix Double
forall a. Matrix a -> Matrix a
M.transposeMatrix Matrix Double
a)
            Right (L.SVD (M.transposeMatrix (L.svdVT f)) (L.singularValues f) (M.transposeMatrix (L.svdU f)))
        else
            if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
                then Matrix Double -> [Double] -> Matrix Double -> SVD
L.SVD (Matrix Double -> [Double] -> Matrix Double -> SVD)
-> Either MathError (Matrix Double)
-> Either MathError ([Double] -> Matrix Double -> SVD)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int -> Int -> [Double] -> Either MathError (Matrix Double)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked Int
m Int
0 [] Either MathError ([Double] -> Matrix Double -> SVD)
-> Either MathError [Double]
-> Either MathError (Matrix Double -> SVD)
forall a b.
Either MathError (a -> b)
-> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> [Double] -> Either MathError [Double]
forall a. a -> Either MathError a
forall (f :: * -> *) a. Applicative f => a -> f a
pure [] Either MathError (Matrix Double -> SVD)
-> Either MathError (Matrix Double) -> Either MathError SVD
forall a b.
Either MathError (a -> b)
-> Either MathError a -> Either MathError b
forall (f :: * -> *) a b. Applicative f => f (a -> b) -> f a -> f b
<*> Int -> Int -> [Double] -> Either MathError (Matrix Double)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
M.fromListChecked Int
0 Int
n []
                else do
                    let scale :: Double
scale = [Double] -> Double
forall a. Ord a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
maximum (Double
0 Double -> [Double] -> [Double]
forall a. a -> [a] -> [a]
: (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
forall a. Num a => a -> a
abs (Matrix Double -> [Double]
forall a. Matrix a -> [a]
M.toList Matrix Double
a))
                        s :: Double
s = if Double
scale Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 then Double
1 else Double
scale
                        cols :: [[Double]]
cols = [[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
transpose (([Double] -> [Double]) -> [[Double]] -> [[Double]]
forall a b. (a -> b) -> [a] -> [b]
map ((Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
s)) (Matrix Double -> [[Double]]
forall a. Matrix a -> [[a]]
I.entries Matrix Double
a))
                        eye :: [[Double]]
eye = [[if Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
j then Double
1 else Double
0 | Int
i <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]] | Int
j <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
                    if [Bool] -> Bool
forall (t :: * -> *). Foldable t => t Bool -> Bool
or
                        ((Double -> Double -> Bool) -> [Double] -> [Double] -> [Bool]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\Double
original Double
scaled -> Double
original Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0 Bool -> Bool -> Bool
&& Double
scaled Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0) (Matrix Double -> [Double]
forall a. Matrix a -> [a]
M.toList Matrix Double
a) ([[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat ([[Double]] -> [[Double]]
forall a. [[a]] -> [[a]]
transpose [[Double]]
cols)))
                        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
DomainError String
"pure SVD: input dynamic range underflows global scaling")
                        else () -> Either MathError ()
forall a b. b -> Either a b
Right ()
                    (columns, vs) <- Int
-> [[Double]]
-> [[Double]]
-> Either MathError ([[Double]], [[Double]])
forall {a}.
RealFloat a =>
Int -> [[a]] -> [[a]] -> Either MathError ([[a]], [[a]])
sweep Int
0 [[Double]]
cols [[Double]]
eye
                    let indices = (Int -> Down Double) -> [Int] -> [Int]
forall b a. Ord b => (a -> b) -> [a] -> [a]
sortOn (Double -> Down Double
forall a. a -> Down a
Down (Double -> Down Double) -> (Int -> Double) -> Int -> Down Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([Double] -> Double
forall a. RealFloat a => [a] -> a
I.stableNorm ([Double] -> Double) -> (Int -> [Double]) -> Int -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ([[Double]]
columns [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!!))) [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]
                        norms = [[Double] -> Double
forall a. RealFloat a => [a] -> a
I.stableNorm ([[Double]]
columns [[Double]] -> Int -> [Double]
forall a. HasCallStack => [a] -> Int -> a
!! Int
i) | Int
i <- [Int]
indices]
                    us <- foldM (complete m) [] (zip indices norms `zip` repeat columns)
                    let singular = (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s) [Double]
norms
                    finite singular
                    u <- M.fromListChecked m n (concat (transpose us))
                    vt <- M.fromListChecked n n (concat [vs !! i | i <- indices])
                    Right (L.SVD u singular vt)
  where
    n :: Int
n = (Int, Int) -> Int
forall a b. (a, b) -> b
snd (Matrix Double -> (Int, Int)
forall a. Matrix a -> (Int, Int)
M.shape Matrix Double
a)
    pairs :: [(Int, Int)]
pairs = [(Int
p, Int
q) | Int
p <- [Int
0 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1], Int
q <- [Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 .. Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1]]
    inner :: [a] -> [a] -> a
inner [a]
x [a]
y = [a] -> a
forall a. Num a => [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
x [a]
y)
    -- Compare normalized correlation, not squared tiny norms: distinct
    -- representable columns may have products below the subnormal range.
    converged :: [[a]] -> (Int, Int) -> Bool
converged [[a]]
cols (Int
p, Int
q) =
        let x :: [a]
x = [[a]]
cols [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
p; y :: [a]
y = [[a]]
cols [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
q; nx :: a
nx = [a] -> a
forall a. RealFloat a => [a] -> a
I.stableNorm [a]
x; ny :: a
ny = [a] -> a
forall a. RealFloat a => [a] -> a
I.stableNorm [a]
y
         in a
nx a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 Bool -> Bool -> Bool
|| a
ny a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 Bool -> Bool -> Bool
|| a -> a
forall a. Num a => a -> a
abs ([a] -> [a] -> a
forall {a}. Num a => [a] -> [a] -> a
inner ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
nx) [a]
x) ((a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
ny) [a]
y)) a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
2e-14
    sweep :: Int -> [[a]] -> [[a]] -> Either MathError ([[a]], [[a]])
sweep Int
count [[a]]
cols [[a]]
vs
        | ((Int, Int) -> Bool) -> [(Int, Int)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ([[a]] -> (Int, Int) -> Bool
forall {a}. RealFloat a => [[a]] -> (Int, Int) -> Bool
converged [[a]]
cols) [(Int, Int)]
pairs = ([[a]], [[a]]) -> Either MathError ([[a]], [[a]])
forall a b. b -> Either a b
Right ([[a]]
cols, [[a]]
vs)
        | Int
count Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
100 = MathError -> Either MathError ([[a]], [[a]])
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"pure one-sided Jacobi SVD" Int
count)
        | Bool
otherwise = let ([[a]]
cc, [[a]]
vv) = (([[a]], [[a]]) -> (Int, Int) -> ([[a]], [[a]]))
-> ([[a]], [[a]]) -> [(Int, Int)] -> ([[a]], [[a]])
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl ([[a]], [[a]]) -> (Int, Int) -> ([[a]], [[a]])
forall {a}.
RealFloat a =>
([[a]], [[a]]) -> (Int, Int) -> ([[a]], [[a]])
rotate ([[a]]
cols, [[a]]
vs) [(Int, Int)]
pairs in Int -> [[a]] -> [[a]] -> Either MathError ([[a]], [[a]])
sweep (Int
count Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) [[a]]
cc [[a]]
vv
    rotate :: ([[a]], [[a]]) -> (Int, Int) -> ([[a]], [[a]])
rotate state :: ([[a]], [[a]])
state@([[a]]
cols, [[a]]
vs) pair :: (Int, Int)
pair@(Int
p, Int
q)
        | [[a]] -> (Int, Int) -> Bool
forall {a}. RealFloat a => [[a]] -> (Int, Int) -> Bool
converged [[a]]
cols (Int, Int)
pair = ([[a]], [[a]])
state
        | Bool
otherwise =
            let x :: [a]
x = [[a]]
cols [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
p
                y :: [a]
y = [[a]]
cols [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
q
                pairScale :: a
pairScale = a -> a -> a
forall a. Ord a => a -> a -> a
max ([a] -> a
forall a. RealFloat a => [a] -> a
I.stableNorm [a]
x) ([a] -> a
forall a. RealFloat a => [a] -> a
I.stableNorm [a]
y)
                xx :: [a]
xx = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
pairScale) [a]
x
                yy :: [a]
yy = (a -> a) -> [a] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map (a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
pairScale) [a]
y
                aa :: a
aa = [a] -> [a] -> a
forall {a}. Num a => [a] -> [a] -> a
inner [a]
xx [a]
xx
                bb :: a
bb = [a] -> [a] -> a
forall {a}. Num a => [a] -> [a] -> a
inner [a]
yy [a]
yy
                ab :: a
ab = [a] -> [a] -> a
forall {a}. Num a => [a] -> [a] -> a
inner [a]
xx [a]
yy
                delta :: a
delta = (a
bb a -> a -> a
forall a. Num a => a -> a -> a
- a
aa) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
2
                -- This form avoids constructing delta/ab when a useful
                -- small rotation exists but that intermediate would overflow.
                t :: a
t = if a
delta a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 then a
1 else a
ab a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a
delta a -> a -> a
forall a. Num a => a -> a -> a
+ a -> a
forall a. Num a => a -> a
signum a
delta a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> a
forall a. RealFloat a => [a] -> a
I.stableNorm [a
delta, a
ab])
                c :: a
c = a
1 a -> a -> a
forall a. Fractional a => a -> a -> a
/ a -> a
forall a. Floating a => a -> a
sqrt (a
1 a -> a -> a
forall a. Num a => a -> a -> a
+ a
t a -> a -> a
forall a. Num a => a -> a -> a
* a
t)
                s :: a
s = a
t a -> a -> a
forall a. Num a => a -> a -> a
* a
c
                turn :: [[a]] -> [[a]]
turn [[a]]
zs =
                    Int -> [a] -> [[a]] -> [[a]]
forall a. Int -> a -> [a] -> [a]
I.replace
                        Int
q
                        ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
u a
v -> a
s a -> a -> a
forall a. Num a => a -> a -> a
* a
u a -> a -> a
forall a. Num a => a -> a -> a
+ a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
v) ([[a]]
zs [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
p) ([[a]]
zs [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
q))
                        (Int -> [a] -> [[a]] -> [[a]]
forall a. Int -> a -> [a] -> [a]
I.replace Int
p ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (\a
u a
v -> a
c a -> a -> a
forall a. Num a => a -> a -> a
* a
u a -> a -> a
forall a. Num a => a -> a -> a
- a
s a -> a -> a
forall a. Num a => a -> a -> a
* a
v) ([[a]]
zs [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
p) ([[a]]
zs [[a]] -> Int -> [a]
forall a. HasCallStack => [a] -> Int -> a
!! Int
q)) [[a]]
zs)
             in ([[a]] -> [[a]]
turn [[a]]
cols, [[a]] -> [[a]]
turn [[a]]
vs)
    complete :: p -> [[b]] -> ((Int, b), [[b]]) -> Either MathError [[b]]
complete p
m [[b]]
previous ((Int
idx, b
norm), [[b]]
cols)
        | b
norm b -> b -> Bool
forall a. Ord a => a -> a -> Bool
> b
0 = [[b]] -> Either MathError [[b]]
forall a b. b -> Either a b
Right ([[b]]
previous [[b]] -> [[b]] -> [[b]]
forall a. [a] -> [a] -> [a]
++ [(b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map (b -> b -> b
forall a. Fractional a => a -> a -> a
/ b
norm) ([[b]]
cols [[b]] -> Int -> [b]
forall a. HasCallStack => [a] -> Int -> a
!! Int
idx)])
        | Bool
otherwise =
            let project :: [b] -> [b]
project [b]
v = ([b] -> [b] -> [b]) -> [b] -> [[b]] -> [b]
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (\[b]
x [b]
u -> (b -> b -> b) -> [b] -> [b] -> [b]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (-) [b]
x ((b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map ([b] -> [b] -> b
forall {a}. Num a => [a] -> [a] -> a
inner [b]
x [b]
u b -> b -> b
forall a. Num a => a -> a -> a
*) [b]
u)) [b]
v [[b]]
previous
                candidates :: [[b]]
candidates = [[b] -> [b]
project ([b] -> [b]
project [if p
i p -> p -> Bool
forall a. Eq a => a -> a -> Bool
== p
j then b
1 else b
0 | p
i <- [p
0 .. p
m p -> p -> p
forall a. Num a => a -> a -> a
- p
1]]) | p
j <- [p
0 .. p
m p -> p -> p
forall a. Num a => a -> a -> a
- p
1]]
                good :: [[b]]
good = ([b] -> Bool) -> [[b]] -> [[b]]
forall a. (a -> Bool) -> [a] -> [a]
filter ((b -> b -> Bool
forall a. Ord a => a -> a -> Bool
> b
1e-10) (b -> Bool) -> ([b] -> b) -> [b] -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [b] -> b
forall a. RealFloat a => [a] -> a
I.stableNorm) [[b]]
candidates
             in case [[b]]
good of
                    [] -> MathError -> Either MathError [[b]]
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
"pure SVD null-space completion" ([[b]] -> Int
forall a. [a] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[b]]
previous))
                    [b]
x : [[b]]
_ -> [[b]] -> Either MathError [[b]]
forall a b. b -> Either a b
Right ([[b]]
previous [[b]] -> [[b]] -> [[b]]
forall a. [a] -> [a] -> [a]
++ [(b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map (b -> b -> b
forall a. Fractional a => a -> a -> a
/ [b] -> b
forall a. RealFloat a => [a] -> a
I.stableNorm [b]
x) [b]
x])