{-# LANGUAGE ForeignFunctionInterface #-}

-- SPDX-License-Identifier: MPL-2.0

{- | Double-precision optional, dynamically loaded LAPACKE spectral routines. All foreign buffers are private
copies; no Haskell storage is exposed or mutated. Inputs must be finite and
dimensions must fit 32-bit indices; the adapter widens them for ILP64.
The default system liblapacke SONAME has its conventional LP64 ABI. A custom
HMATH_LAPACKE_LIBRARY requires an explicit HMATH_LAPACKE_INTEGER_BITS=32/64
declaration before first use. HMATH_BACKEND=pure disables native routines.
Missing native routines fail explicitly; pure alternatives live in HMath.Backend. Pure wrappers are safe
because the called routines have no observable effects outside those buffers.
LAPACKE interfaces: https://www.netlib.org/lapack/lapacke.html
-}
module HMath.Matrix.LAPACK (
    SVD (..),
    svd,
    fullSVD,
    truncatedSVD,
    randomizedSVD,
    Eigen (..),
    eigen,
    symmetricEigen,
    GeneralizedEigen (..),
    generalizedEigen,
    Schur (..),
    schur,
    QZ (..),
    qz,
    Hessenberg (..),
    hessenberg,
    Balance (..),
    balance,
    ComplexSVD (..),
    complexSVD,
    complexEigen,
) where

import Data.Complex (Complex ((:+)), imagPart, realPart)
import Foreign
import Foreign.C.Types
import HMath.Error
import HMath.Matrix
import HMath.Matrix.Internal
import HMath.Matrix.QR (qrDecompose, qrQ)
import System.IO.Unsafe (unsafePerformIO)

foreign import ccall unsafe "hmath_lapacke_zgesdd"
    cComplexSVD ::
        CInt ->
        CChar ->
        CInt ->
        CInt ->
        Ptr (Complex CDouble) ->
        CInt ->
        Ptr CDouble ->
        Ptr (Complex CDouble) ->
        CInt ->
        Ptr (Complex CDouble) ->
        CInt ->
        IO CInt
foreign import ccall unsafe "hmath_lapacke_zgeev"
    cComplexEigen ::
        CInt ->
        CChar ->
        CChar ->
        CInt ->
        Ptr (Complex CDouble) ->
        CInt ->
        Ptr (Complex CDouble) ->
        Ptr (Complex CDouble) ->
        CInt ->
        Ptr (Complex CDouble) ->
        CInt ->
        IO CInt
foreign import ccall unsafe "hmath_lapacke_dsyev"
    cSymmetric :: CInt -> CChar -> CChar -> CInt -> Ptr CDouble -> CInt -> Ptr CDouble -> IO CInt
foreign import ccall unsafe "hmath_lapacke_dgesdd"
    cSVD ::
        CInt ->
        CChar ->
        CInt ->
        CInt ->
        Ptr CDouble ->
        CInt ->
        Ptr CDouble ->
        Ptr CDouble ->
        CInt ->
        Ptr CDouble ->
        CInt ->
        IO CInt
foreign import ccall unsafe "hmath_lapacke_dgeev"
    cEigen ::
        CInt ->
        CChar ->
        CChar ->
        CInt ->
        Ptr CDouble ->
        CInt ->
        Ptr CDouble ->
        Ptr CDouble ->
        Ptr CDouble ->
        CInt ->
        Ptr CDouble ->
        CInt ->
        IO CInt
foreign import ccall unsafe "hmath_lapacke_dggev"
    cGeneralized ::
        CInt ->
        CChar ->
        CChar ->
        CInt ->
        Ptr CDouble ->
        CInt ->
        Ptr CDouble ->
        CInt ->
        Ptr CDouble ->
        Ptr CDouble ->
        Ptr CDouble ->
        Ptr CDouble ->
        CInt ->
        Ptr CDouble ->
        CInt ->
        IO CInt
foreign import ccall unsafe "hmath_lapacke_dgees"
    cSchur ::
        CInt ->
        CChar ->
        CChar ->
        FunPtr a ->
        CInt ->
        Ptr CDouble ->
        CInt ->
        Ptr CInt ->
        Ptr CDouble ->
        Ptr CDouble ->
        Ptr CDouble ->
        CInt ->
        IO CInt
foreign import ccall unsafe "hmath_lapacke_dgges"
    cQZ ::
        CInt ->
        CChar ->
        CChar ->
        CChar ->
        FunPtr a ->
        CInt ->
        Ptr CDouble ->
        CInt ->
        Ptr CDouble ->
        CInt ->
        Ptr CInt ->
        Ptr CDouble ->
        Ptr CDouble ->
        Ptr CDouble ->
        Ptr CDouble ->
        CInt ->
        Ptr CDouble ->
        CInt ->
        IO CInt
foreign import ccall unsafe "hmath_lapacke_dgehrd"
    cHessenberg :: CInt -> CInt -> CInt -> CInt -> Ptr CDouble -> CInt -> Ptr CDouble -> IO CInt
foreign import ccall unsafe "hmath_lapacke_dorghr"
    cHessenbergQ :: CInt -> CInt -> CInt -> CInt -> Ptr CDouble -> CInt -> Ptr CDouble -> IO CInt
foreign import ccall unsafe "hmath_lapacke_dgebal"
    cBalance ::
        CInt -> CChar -> CInt -> Ptr CDouble -> CInt -> Ptr CInt -> Ptr CInt -> Ptr CDouble -> IO CInt

{- | Real singular-value decomposition with U,descending nonnegative singular values,and VT.
Thin/full dimensions follow the selected factorization function.
-}
data SVD = SVD
    { SVD -> Matrix Double
svdU :: Matrix Double
    -- ^ Left singular vectors as columns of the real thin/full SVD factor.
    , SVD -> [Double]
singularValues :: [Double]
    -- ^ Nonnegative real singular values in descending order.
    , SVD -> Matrix Double
svdVT :: Matrix Double
    -- ^ Transposed right singular-vector factor in A=U*diag(s)*VT.
    }
    deriving (SVD -> SVD -> Bool
(SVD -> SVD -> Bool) -> (SVD -> SVD -> Bool) -> Eq SVD
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: SVD -> SVD -> Bool
== :: SVD -> SVD -> Bool
$c/= :: SVD -> SVD -> Bool
/= :: SVD -> SVD -> Bool
Eq, Int -> SVD -> ShowS
[SVD] -> ShowS
SVD -> String
(Int -> SVD -> ShowS)
-> (SVD -> String) -> ([SVD] -> ShowS) -> Show SVD
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> SVD -> ShowS
showsPrec :: Int -> SVD -> ShowS
$cshow :: SVD -> String
show :: SVD -> String
$cshowList :: [SVD] -> ShowS
showList :: [SVD] -> ShowS
Show)

{- | Complex eigenvalues with corresponding right eigenvector columns, including conjugate
pairs from real input.
-}
data Eigen = Eigen
    { Eigen -> [Complex Double]
eigenvalues :: [Complex Double]
    -- ^ Complex eigenvalues matching eigenvector column order.
    , Eigen -> Matrix (Complex Double)
eigenvectors :: Matrix (Complex Double)
    -- ^ Right eigenvectors as complex matrix columns.
    }
    deriving (Eigen -> Eigen -> Bool
(Eigen -> Eigen -> Bool) -> (Eigen -> Eigen -> Bool) -> Eq Eigen
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Eigen -> Eigen -> Bool
== :: Eigen -> Eigen -> Bool
$c/= :: Eigen -> Eigen -> Bool
/= :: Eigen -> Eigen -> Bool
Eq, Int -> Eigen -> ShowS
[Eigen] -> ShowS
Eigen -> String
(Int -> Eigen -> ShowS)
-> (Eigen -> String) -> ([Eigen] -> ShowS) -> Show Eigen
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Eigen -> ShowS
showsPrec :: Int -> Eigen -> ShowS
$cshow :: Eigen -> String
show :: Eigen -> String
$cshowList :: [Eigen] -> ShowS
showList :: [Eigen] -> ShowS
Show)

{- | Projective eigenvalues alpha/beta preserve infinite eigenvalues and
indeterminate singular pencils without dividing by a zero beta.
-}
data GeneralizedEigen = GeneralizedEigen
    { GeneralizedEigen -> [Complex Double]
eigenAlpha :: [Complex Double]
    -- ^ Generalized-eigenvalue numerators alpha; retain with beta to represent infinite eigenvalues.
    , GeneralizedEigen -> [Double]
eigenBeta :: [Double]
    -- ^ Generalized-eigenvalue denominators beta; zero beta can represent an infinite eigenvalue.
    , GeneralizedEigen -> Matrix (Complex Double)
generalizedEigenvectors :: Matrix (Complex Double)
    -- ^ Right eigenvector columns satisfying beta*A*v=alpha*B*v.
    }
    deriving (GeneralizedEigen -> GeneralizedEigen -> Bool
(GeneralizedEigen -> GeneralizedEigen -> Bool)
-> (GeneralizedEigen -> GeneralizedEigen -> Bool)
-> Eq GeneralizedEigen
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: GeneralizedEigen -> GeneralizedEigen -> Bool
== :: GeneralizedEigen -> GeneralizedEigen -> Bool
$c/= :: GeneralizedEigen -> GeneralizedEigen -> Bool
/= :: GeneralizedEigen -> GeneralizedEigen -> Bool
Eq, Int -> GeneralizedEigen -> ShowS
[GeneralizedEigen] -> ShowS
GeneralizedEigen -> String
(Int -> GeneralizedEigen -> ShowS)
-> (GeneralizedEigen -> String)
-> ([GeneralizedEigen] -> ShowS)
-> Show GeneralizedEigen
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> GeneralizedEigen -> ShowS
showsPrec :: Int -> GeneralizedEigen -> ShowS
$cshow :: GeneralizedEigen -> String
show :: GeneralizedEigen -> String
$cshowList :: [GeneralizedEigen] -> ShowS
showList :: [GeneralizedEigen] -> ShowS
Show)

-- | Real Schur factors A=Q*T*transpose(Q), with quasi-triangular T and complex eigenvalues.
data Schur = Schur
    { Schur -> Matrix Double
schurQ :: Matrix Double
    -- ^ Orthogonal basis in A=Q*T*transpose(Q).
    , Schur -> Matrix Double
schurT :: Matrix Double
    -- ^ Real quasi-upper-triangular Schur form, including 2-by-2 blocks for conjugate pairs.
    , Schur -> [Complex Double]
schurEigenvalues :: [Complex Double]
    -- ^ Complex eigenvalues decoded from the real Schur form.
    }
    deriving (Schur -> Schur -> Bool
(Schur -> Schur -> Bool) -> (Schur -> Schur -> Bool) -> Eq Schur
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Schur -> Schur -> Bool
== :: Schur -> Schur -> Bool
$c/= :: Schur -> Schur -> Bool
/= :: Schur -> Schur -> Bool
Eq, Int -> Schur -> ShowS
[Schur] -> ShowS
Schur -> String
(Int -> Schur -> ShowS)
-> (Schur -> String) -> ([Schur] -> ShowS) -> Show Schur
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Schur -> ShowS
showsPrec :: Int -> Schur -> ShowS
$cshow :: Schur -> String
show :: Schur -> String
$cshowList :: [Schur] -> ShowS
showList :: [Schur] -> ShowS
Show)

{- | Generalized real Schur factors for a matrix pair, retaining alpha/beta eigenvalue
representation to allow infinite eigenvalues.
-}
data QZ = QZ
    { QZ -> Matrix Double
qzQ :: Matrix Double
    -- ^ Left orthogonal QZ transform, with A=Q*S*transpose(Z).
    , QZ -> Matrix Double
qzZ :: Matrix Double
    -- ^ Right orthogonal QZ transform for both matrix factors.
    , QZ -> Matrix Double
qzS :: Matrix Double
    -- ^ First generalized real Schur factor S.
    , QZ -> Matrix Double
qzT :: Matrix Double
    -- ^ Second generalized real Schur factor T.
    , QZ -> [Complex Double]
qzAlpha :: [Complex Double]
    -- ^ Generalized Schur eigenvalue numerators alpha.
    , QZ -> [Double]
qzBeta :: [Double]
    -- ^ Generalized Schur eigenvalue denominators beta; do not divide blindly when zero.
    }
    deriving (QZ -> QZ -> Bool
(QZ -> QZ -> Bool) -> (QZ -> QZ -> Bool) -> Eq QZ
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: QZ -> QZ -> Bool
== :: QZ -> QZ -> Bool
$c/= :: QZ -> QZ -> Bool
/= :: QZ -> QZ -> Bool
Eq, Int -> QZ -> ShowS
[QZ] -> ShowS
QZ -> String
(Int -> QZ -> ShowS)
-> (QZ -> String) -> ([QZ] -> ShowS) -> Show QZ
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> QZ -> ShowS
showsPrec :: Int -> QZ -> ShowS
$cshow :: QZ -> String
show :: QZ -> String
$cshowList :: [QZ] -> ShowS
showList :: [QZ] -> ShowS
Show)

-- | Orthogonal reduction A=Q*H*transpose(Q), with H upper Hessenberg.
data Hessenberg = Hessenberg
    { Hessenberg -> Matrix Double
hessenbergQ :: Matrix Double
    -- ^ Orthogonal transform in A=Q*H*transpose(Q).
    , Hessenberg -> Matrix Double
hessenbergH :: Matrix Double
    -- ^ Upper Hessenberg form with zeros below the first subdiagonal.
    }
    deriving (Hessenberg -> Hessenberg -> Bool
(Hessenberg -> Hessenberg -> Bool)
-> (Hessenberg -> Hessenberg -> Bool) -> Eq Hessenberg
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Hessenberg -> Hessenberg -> Bool
== :: Hessenberg -> Hessenberg -> Bool
$c/= :: Hessenberg -> Hessenberg -> Bool
/= :: Hessenberg -> Hessenberg -> Bool
Eq, Int -> Hessenberg -> ShowS
[Hessenberg] -> ShowS
Hessenberg -> String
(Int -> Hessenberg -> ShowS)
-> (Hessenberg -> String)
-> ([Hessenberg] -> ShowS)
-> Show Hessenberg
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Hessenberg -> ShowS
showsPrec :: Int -> Hessenberg -> ShowS
$cshow :: Hessenberg -> String
show :: Hessenberg -> String
$cshowList :: [Hessenberg] -> ShowS
showList :: [Hessenberg] -> ShowS
Show)

-- | D^-1 A D = balancedMatrix. Scale-only balancing preserves index order.
data Balance = Balance
    { Balance -> Matrix Double
balancedMatrix :: Matrix Double
    -- ^ Matrix after LAPACK diagonal balancing.
    , Balance -> [Double]
balanceScales :: [Double]
    -- ^ Diagonal basis scaling factors used by the balancing transformation.
    }
    deriving (Balance -> Balance -> Bool
(Balance -> Balance -> Bool)
-> (Balance -> Balance -> Bool) -> Eq Balance
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: Balance -> Balance -> Bool
== :: Balance -> Balance -> Bool
$c/= :: Balance -> Balance -> Bool
/= :: Balance -> Balance -> Bool
Eq, Int -> Balance -> ShowS
[Balance] -> ShowS
Balance -> String
(Int -> Balance -> ShowS)
-> (Balance -> String) -> ([Balance] -> ShowS) -> Show Balance
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> Balance -> ShowS
showsPrec :: Int -> Balance -> ShowS
$cshow :: Balance -> String
show :: Balance -> String
$cshowList :: [Balance] -> ShowS
showList :: [Balance] -> ShowS
Show)

char :: Char -> CChar
char :: Char -> CChar
char = Int -> CChar
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CChar) -> (Char -> Int) -> Char -> CChar
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Char -> Int
forall a. Enum a => a -> Int
fromEnum
ni :: Int -> CInt
ni :: Int -> CInt
ni = Int -> CInt
forall a b. (Integral a, Num b) => a -> b
fromIntegral
withValues :: [Double] -> (Ptr CDouble -> IO a) -> IO a
withValues :: forall a. [Double] -> (Ptr CDouble -> IO a) -> IO a
withValues [Double]
xs = [CDouble] -> (Ptr CDouble -> IO a) -> IO a
forall a b. Storable a => [a] -> (Ptr a -> IO b) -> IO b
withArray ((Double -> CDouble) -> [Double] -> [CDouble]
forall a b. (a -> b) -> [a] -> [b]
map Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac [Double]
xs)
peekValues :: Int -> Ptr CDouble -> IO [Double]
peekValues :: Int -> Ptr CDouble -> IO [Double]
peekValues Int
n Ptr CDouble
p = (CDouble -> Double) -> [CDouble] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac ([CDouble] -> [Double]) -> IO [CDouble] -> IO [Double]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int -> Ptr CDouble -> IO [CDouble]
forall a. Storable a => Int -> Ptr a -> IO [a]
peekArray Int
n Ptr CDouble
p
buffer :: Int -> (Ptr CDouble -> IO a) -> IO a
buffer :: forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
n Ptr CDouble -> IO a
action = Int -> (Ptr CDouble -> IO a) -> IO a
forall a b. Storable a => Int -> (Ptr a -> IO b) -> IO b
allocaArray (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
1 Int
n) ((Ptr CDouble -> IO a) -> IO a) -> (Ptr CDouble -> IO a) -> IO a
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
p -> Ptr CDouble -> [CDouble] -> IO ()
forall a. Storable a => Ptr a -> [a] -> IO ()
pokeArray Ptr CDouble
p (Int -> CDouble -> [CDouble]
forall a. Int -> a -> [a]
replicate (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
1 Int
n) CDouble
0) IO () -> IO a -> IO a
forall a b. IO a -> IO b -> IO b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Ptr CDouble -> IO a
action Ptr CDouble
p
status :: String -> CInt -> Either MathError ()
status :: String -> CInt -> Either MathError ()
status String
name CInt
info
    | CInt
info CInt -> CInt -> Bool
forall a. Eq a => a -> a -> Bool
== CInt
0 = () -> Either MathError ()
forall a b. b -> Either a b
Right ()
    | CInt
info CInt -> CInt -> Bool
forall a. Eq a => a -> a -> Bool
== -CInt
10000 =
        MathError -> Either MathError ()
forall a b. a -> Either a b
Left
            ( String -> MathError
DomainError
                ( String
name
                    String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
": native LAPACKE unavailable; use HMath.Backend pure routines or configure an LP64/ILP64 provider"
                )
            )
    | CInt
info CInt -> CInt -> Bool
forall a. Eq a => a -> a -> Bool
== -CInt
10001 = MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument (String
name String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
": LAPACK status exceeds supported result index"))
    | CInt
info CInt -> CInt -> Bool
forall a. Ord a => a -> a -> Bool
< CInt
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
": LAPACKE rejected 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 (CInt -> Integer
forall a. Integral a => a -> Integer
toInteger CInt
info))))
    | Bool
otherwise = MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> Int -> MathError
NonConvergence String
name (CInt -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral CInt
info))
validate :: Matrix Double -> Either MathError ()
allocationSizes :: [Integer] -> Either MathError ()
allocationSizes :: [Integer] -> Either MathError ()
allocationSizes [Integer]
sizes
    | (Integer -> Bool) -> [Integer] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
any
        ( \Integer
count -> Integer
count Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 Bool -> Bool -> Bool
|| Integer
count Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger (CDouble -> Int
forall a. Storable a => a -> Int
sizeOf (CDouble
forall a. HasCallStack => a
undefined :: CDouble)) 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)
        )
        [Integer]
sizes =
        MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"LAPACKE: byte allocation overflow")
    | Bool
otherwise = () -> Either MathError ()
forall a b. b -> Either a b
Right ()
validate :: Matrix Double -> Either MathError ()
validate Matrix Double
a = do
    let (Int
m, Int
n) = Matrix Double -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix Double
a
    [Integer] -> Either MathError ()
allocationSizes [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
n]
    if Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
m 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) 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)
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"LAPACKE: dimensions exceed 32-bit integer ABI")
        else Matrix Double -> Either MathError ()
forall a. RealFloat a => Matrix a -> Either MathError ()
finiteMatrix Matrix Double
a
square :: Matrix Double -> Either MathError Int
square :: Matrix Double -> Either MathError Int
square Matrix Double
a
    | Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n = MathError -> Either MathError Int
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
m, Int
n))
    | Bool
otherwise = Matrix Double -> Either MathError ()
validate Matrix Double
a Either MathError () -> Either MathError Int -> Either MathError Int
forall a b.
Either MathError a -> Either MathError b -> Either MathError b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Int -> Either MathError Int
forall a b. b -> Either a b
Right Int
n
  where
    (Int
m, Int
n) = Matrix Double -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix Double
a

{-# NOINLINE svd #-}

{- | Compute thin real SVD A=U*diag(s)*VT via LAPACK, with descending nonnegative singular
values and k=min(m,n). Return m-by-k U and k-by-n VT. Invalid input, unavailable
backend, allocation bounds or LAPACK failure returns Left MathError.
-}
svd :: Matrix Double -> Either MathError SVD
svd :: Matrix Double -> Either MathError SVD
svd Matrix Double
a = do
    Matrix Double -> Either MathError ()
validate Matrix Double
a
    let (Int
m, Int
n) = Matrix Double -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix Double
a; k :: Int
k = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
m Int
n
    [Integer] -> Either MathError ()
allocationSizes [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, 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
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
        then Matrix Double -> [Double] -> Matrix Double -> SVD
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)
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)
fromListChecked Int
0 Int
n []
        else IO (Either MathError SVD) -> Either MathError SVD
forall a. IO a -> a
unsafePerformIO (IO (Either MathError SVD) -> Either MathError SVD)
-> IO (Either MathError SVD) -> Either MathError SVD
forall a b. (a -> b) -> a -> b
$ [Double]
-> (Ptr CDouble -> IO (Either MathError SVD))
-> IO (Either MathError SVD)
forall a. [Double] -> (Ptr CDouble -> IO a) -> IO a
withValues (Matrix Double -> [Double]
forall a. Matrix a -> [a]
toList Matrix Double
a) ((Ptr CDouble -> IO (Either MathError SVD))
 -> IO (Either MathError SVD))
-> (Ptr CDouble -> IO (Either MathError SVD))
-> IO (Either MathError SVD)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pa -> Int
-> (Ptr CDouble -> IO (Either MathError SVD))
-> IO (Either MathError SVD)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
k ((Ptr CDouble -> IO (Either MathError SVD))
 -> IO (Either MathError SVD))
-> (Ptr CDouble -> IO (Either MathError SVD))
-> IO (Either MathError SVD)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
ps -> Int
-> (Ptr CDouble -> IO (Either MathError SVD))
-> IO (Either MathError SVD)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k) ((Ptr CDouble -> IO (Either MathError SVD))
 -> IO (Either MathError SVD))
-> (Ptr CDouble -> IO (Either MathError SVD))
-> IO (Either MathError SVD)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pu -> Int
-> (Ptr CDouble -> IO (Either MathError SVD))
-> IO (Either MathError SVD)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n) ((Ptr CDouble -> IO (Either MathError SVD))
 -> IO (Either MathError SVD))
-> (Ptr CDouble -> IO (Either MathError SVD))
-> IO (Either MathError SVD)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pv -> do
            info <- CInt
-> CChar
-> CInt
-> CInt
-> Ptr CDouble
-> CInt
-> Ptr CDouble
-> Ptr CDouble
-> CInt
-> Ptr CDouble
-> CInt
-> IO CInt
cSVD CInt
101 (Char -> CChar
char Char
'S') (Int -> CInt
ni Int
m) (Int -> CInt
ni Int
n) Ptr CDouble
pa (Int -> CInt
ni Int
n) Ptr CDouble
ps Ptr CDouble
pu (Int -> CInt
ni Int
k) Ptr CDouble
pv (Int -> CInt
ni Int
n)
            s <- peekValues k ps
            u <- peekValues (m * k) pu
            v <- peekValues (k * n) pv
            pure $ do
                status "SVD" info
                finiteValues (s ++ u ++ v)
                SVD <$> fromListChecked m k u <*> pure s <*> fromListChecked k n v

{-# NOINLINE fullSVD #-}

-- | Full square U and V^T, including null-space completion.
fullSVD :: Matrix Double -> Either MathError SVD
fullSVD :: Matrix Double -> Either MathError SVD
fullSVD Matrix Double
a = do
    Matrix Double -> Either MathError ()
validate Matrix Double
a
    let (Int
m, Int
n) = Matrix Double -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix Double
a; k :: Int
k = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
m Int
n
    [Integer] -> Either MathError ()
allocationSizes [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
m, Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n]
    if 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
m 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)
        Bool -> Bool -> Bool
|| Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n 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)
        then MathError -> Either MathError SVD
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"fullSVD: full factor size overflow")
        else
            if Int
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
                then Matrix Double -> [Double] -> Matrix Double -> SVD
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 -> Either MathError (Matrix Double)
forall a. Num a => Int -> Either MathError (Matrix a)
identityChecked Int
m 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 -> Either MathError (Matrix Double)
forall a. Num a => Int -> Either MathError (Matrix a)
identityChecked Int
n
                else IO (Either MathError SVD) -> Either MathError SVD
forall a. IO a -> a
unsafePerformIO (IO (Either MathError SVD) -> Either MathError SVD)
-> IO (Either MathError SVD) -> Either MathError SVD
forall a b. (a -> b) -> a -> b
$ [Double]
-> (Ptr CDouble -> IO (Either MathError SVD))
-> IO (Either MathError SVD)
forall a. [Double] -> (Ptr CDouble -> IO a) -> IO a
withValues (Matrix Double -> [Double]
forall a. Matrix a -> [a]
toList Matrix Double
a) ((Ptr CDouble -> IO (Either MathError SVD))
 -> IO (Either MathError SVD))
-> (Ptr CDouble -> IO (Either MathError SVD))
-> IO (Either MathError SVD)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pa -> Int
-> (Ptr CDouble -> IO (Either MathError SVD))
-> IO (Either MathError SVD)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
k ((Ptr CDouble -> IO (Either MathError SVD))
 -> IO (Either MathError SVD))
-> (Ptr CDouble -> IO (Either MathError SVD))
-> IO (Either MathError SVD)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
ps -> Int
-> (Ptr CDouble -> IO (Either MathError SVD))
-> IO (Either MathError SVD)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
m) ((Ptr CDouble -> IO (Either MathError SVD))
 -> IO (Either MathError SVD))
-> (Ptr CDouble -> IO (Either MathError SVD))
-> IO (Either MathError SVD)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pu -> Int
-> (Ptr CDouble -> IO (Either MathError SVD))
-> IO (Either MathError SVD)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n) ((Ptr CDouble -> IO (Either MathError SVD))
 -> IO (Either MathError SVD))
-> (Ptr CDouble -> IO (Either MathError SVD))
-> IO (Either MathError SVD)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pv -> do
                    info <- CInt
-> CChar
-> CInt
-> CInt
-> Ptr CDouble
-> CInt
-> Ptr CDouble
-> Ptr CDouble
-> CInt
-> Ptr CDouble
-> CInt
-> IO CInt
cSVD CInt
101 (Char -> CChar
char Char
'A') (Int -> CInt
ni Int
m) (Int -> CInt
ni Int
n) Ptr CDouble
pa (Int -> CInt
ni Int
n) Ptr CDouble
ps Ptr CDouble
pu (Int -> CInt
ni Int
m) Ptr CDouble
pv (Int -> CInt
ni Int
n)
                    s <- peekValues k ps
                    u <- peekValues (m * m) pu
                    v <- peekValues (n * n) pv
                    pure $ do
                        status "full SVD" info
                        finiteValues (s ++ u ++ v)
                        SVD <$> fromListChecked m m u <*> pure s <*> fromListChecked n n v

{- | Given retained rank k then real matrix, compute thin SVD and retain its first k
components. Reject k outside the available rank dimensions; decomposition failures
propagate.
-}
truncatedSVD :: Int -> Matrix Double -> Either MathError SVD
truncatedSVD :: Int -> Matrix Double -> Either MathError SVD
truncatedSVD Int
k Matrix Double
a
    | Int
k 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 -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
m Int
n = MathError -> Either MathError SVD
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"truncatedSVD: invalid rank")
    | Bool
otherwise = do
        f <- Matrix Double -> Either MathError SVD
svd Matrix Double
a
        u <- sliceMatrixChecked (0, 0) (m, k) (svdU f)
        v <- sliceMatrixChecked (0, 0) (k, n) (svdVT f)
        Right (SVD u (take k (singularValues f)) v)
  where
    (Int
m, Int
n) = Matrix Double -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix Double
a

{- | Randomized range finder with an explicit seed, oversampling, and power
iterations, followed by a small deterministic SVD. The Gaussian test matrix
uses Box-Muller on a local LCG, so results are reproducible, not cryptographic.
-}
randomizedSVD :: Int -> Int -> Int -> Integer -> Matrix Double -> Either MathError SVD
randomizedSVD :: Int
-> Int -> Int -> Integer -> Matrix Double -> Either MathError SVD
randomizedSVD Int
k Int
oversample Int
power Integer
seed Matrix Double
a
    | Int
k 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 -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
m Int
n Bool -> Bool -> Bool
|| Int
oversample Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 Bool -> Bool -> Bool
|| Int
power Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 =
        MathError -> Either MathError SVD
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"randomizedSVD: rank, oversampling, or power budget invalid")
    | Bool
otherwise = do
        Matrix Double -> Either MathError ()
validate Matrix Double
a
        omega <- Int -> Int -> [Double] -> Either MathError (Matrix Double)
forall a. Int -> Int -> [a] -> Either MathError (Matrix a)
fromListChecked Int
n Int
l (Int -> [Double] -> [Double]
forall a. Int -> [a] -> [a]
take (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
l) (Integer -> [Double]
forall {a}. Floating a => Integer -> [a]
gaussian Integer
seed))
        y <- multiplyChecked a omega
        yFinal <- iteratePower power y
        qf <- qrDecompose yFinal
        q <- sliceMatrixChecked (0, 0) (m, l) (qrQ qf)
        b <- multiplyChecked (transposeMatrix q) a
        f <- truncatedSVD k b
        u <- multiplyChecked q (svdU f)
        Right f{svdU = u}
  where
    (Int
m, Int
n) = Matrix Double -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix Double
a
    l :: Int
l = Integer -> Int
forall a. Num a => Integer -> a
fromInteger (Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
min (Int -> Integer
forall a. Integral a => a -> Integer
toInteger (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
m Int
n)) (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
oversample))
    iteratePower :: t -> Matrix Double -> Either MathError (Matrix Double)
iteratePower t
0 Matrix Double
y = Matrix Double -> Either MathError (Matrix Double)
forall a b. b -> Either a b
Right Matrix Double
y
    iteratePower t
p Matrix Double
y = do
        qf <- Matrix Double -> Either MathError (QR Double)
forall a. RealFloat a => Matrix a -> Either MathError (QR a)
qrDecompose Matrix Double
y
        q <- sliceMatrixChecked (0, 0) (m, l) (qrQ qf)
        z <- multiplyChecked (transposeMatrix a) q
        zf <- qrDecompose z
        zq <- sliceMatrixChecked (0, 0) (n, l) (qrQ zf)
        multiplyChecked a zq >>= iteratePower (p - 1)
    gaussian :: Integer -> [a]
gaussian Integer
initial =
        [a] -> [a]
forall {a}. Floating a => [a] -> [a]
pairs
            ( (Integer -> a) -> [Integer] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map
                (\Integer
x -> (Integer -> a
forall a. Num a => Integer -> a
fromInteger Integer
x a -> a -> a
forall a. Num a => a -> a -> a
+ a
0.5) a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
4294967296)
                (Int -> [Integer] -> [Integer]
forall a. Int -> [a] -> [a]
drop Int
1 ((Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (\Integer
x -> (Integer
1664525 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
x Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1013904223) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
4294967296) (Integer
initial Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
4294967296)))
            )
    pairs :: [a] -> [a]
pairs (a
u : a
v : [a]
rest) = a -> a
forall a. Floating a => a -> a
sqrt (-a
2 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
log a
u) a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
cos (a
2 a -> a -> a
forall a. Num a => a -> a -> a
* a
forall a. Floating a => a
pi a -> a -> a
forall a. Num a => a -> a -> a
* a
v) a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a -> a
forall a. Floating a => a -> a
sqrt (-a
2 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
log a
u) a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
sin (a
2 a -> a -> a
forall a. Num a => a -> a -> a
* a
forall a. Floating a => a
pi a -> a -> a
forall a. Num a => a -> a -> a
* a
v) a -> [a] -> [a]
forall a. a -> [a] -> [a]
: [a] -> [a]
pairs [a]
rest
    pairs [a]
_ = []

complexVectors :: Int -> [Double] -> [Double] -> Either MathError (Matrix (Complex Double))
complexVectors :: Int
-> [Double]
-> [Double]
-> Either MathError (Matrix (Complex Double))
complexVectors Int
n [Double]
imag [Double]
values = Int
-> Int
-> (Int -> Int -> Complex Double)
-> Either MathError (Matrix (Complex Double))
forall a.
Int -> Int -> (Int -> Int -> a) -> Either MathError (Matrix a)
build Int
n Int
n ((Int -> Int -> Complex Double)
 -> Either MathError (Matrix (Complex Double)))
-> (Int -> Int -> Complex Double)
-> Either MathError (Matrix (Complex Double))
forall a b. (a -> b) -> a -> b
$ \Int
i Int
j ->
    if [Double]
imag [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0
        then [Double]
values [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double
0
        else
            if [Double]
imag [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! Int
j Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0
                then [Double]
values [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ [Double]
values [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
                else [Double]
values [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Double -> Double
forall a. Num a => a -> a
negate ([Double]
values [Double] -> Int -> Double
forall a. HasCallStack => [a] -> Int -> a
!! (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
j))
{-# NOINLINE eigen #-}

{- | Compute complex eigenvalues and right eigenvectors of a finite real square matrix via
LAPACK. Complex-conjugate pairs are decoded into complex columns. Invalid
shape/backend/numeric failures return Left MathError.
-}
eigen :: Matrix Double -> Either MathError Eigen
eigen :: Matrix Double -> Either MathError Eigen
eigen Matrix Double
a = do
    n <- Matrix Double -> Either MathError Int
square Matrix Double
a
    if n == 0
        then Eigen [] <$> fromListChecked 0 0 []
        else unsafePerformIO $
            withValues (toList a) $ \Ptr CDouble
pa -> Int
-> (Ptr CDouble -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
n ((Ptr CDouble -> IO (Either MathError Eigen))
 -> IO (Either MathError Eigen))
-> (Ptr CDouble -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pr -> Int
-> (Ptr CDouble -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
n ((Ptr CDouble -> IO (Either MathError Eigen))
 -> IO (Either MathError Eigen))
-> (Ptr CDouble -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pii -> Int
-> (Ptr CDouble -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
1 ((Ptr CDouble -> IO (Either MathError Eigen))
 -> IO (Either MathError Eigen))
-> (Ptr CDouble -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pl -> Int
-> (Ptr CDouble -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n) ((Ptr CDouble -> IO (Either MathError Eigen))
 -> IO (Either MathError Eigen))
-> (Ptr CDouble -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pv -> do
                info <- CInt
-> CChar
-> CChar
-> CInt
-> Ptr CDouble
-> CInt
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> CInt
-> Ptr CDouble
-> CInt
-> IO CInt
cEigen CInt
101 (Char -> CChar
char Char
'N') (Char -> CChar
char Char
'V') (Int -> CInt
ni Int
n) Ptr CDouble
pa (Int -> CInt
ni Int
n) Ptr CDouble
pr Ptr CDouble
pii Ptr CDouble
pl CInt
1 Ptr CDouble
pv (Int -> CInt
ni Int
n)
                re <- peekValues n pr
                im <- peekValues n pii
                vs <- peekValues (n * n) pv
                pure $ do
                    status "eigen" info
                    finiteValues (re ++ im ++ vs)
                    Eigen (zipWith (:+) re im) <$> complexVectors n im vs

{-# NOINLINE symmetricEigen #-}

-- | Real symmetric eigensolver with orthonormal columns, sorted ascending.
symmetricEigen :: Matrix Double -> Either MathError Eigen
symmetricEigen :: Matrix Double -> Either MathError Eigen
symmetricEigen Matrix Double
a = do
    n <- Matrix Double -> Either MathError Int
square Matrix Double
a
    if any
        (\(Double
x, Double
y) -> 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)))
        (zip (toList a) (toList (transposeMatrix a)))
        then Left (DomainError "symmetricEigen: symmetric matrix required")
        else
            if n == 0
                then Eigen [] <$> fromListChecked 0 0 []
                else unsafePerformIO $
                    withValues (toList a) $ \Ptr CDouble
pa -> Int
-> (Ptr CDouble -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
n ((Ptr CDouble -> IO (Either MathError Eigen))
 -> IO (Either MathError Eigen))
-> (Ptr CDouble -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pw -> do
                        info <- CInt
-> CChar
-> CChar
-> CInt
-> Ptr CDouble
-> CInt
-> Ptr CDouble
-> IO CInt
cSymmetric CInt
101 (Char -> CChar
char Char
'V') (Char -> CChar
char Char
'U') (Int -> CInt
ni Int
n) Ptr CDouble
pa (Int -> CInt
ni Int
n) Ptr CDouble
pw
                        ws <- peekValues n pw
                        vs <- peekValues (n * n) pa
                        pure $ do
                            status "symmetricEigen" info
                            finiteValues (ws ++ vs)
                            Eigen (map (:+ 0) ws) <$> fromListChecked n n (map (:+ 0) vs)

{-# NOINLINE generalizedEigen #-}

{- | Given equally sized square matrices A,B, compute generalized alpha/beta eigenvalue
representation and right eigenvectors. Retain alpha/beta to represent infinite
eigenvalues safely; invalid input or LAPACK failure propagates.
-}
generalizedEigen :: Matrix Double -> Matrix Double -> Either MathError GeneralizedEigen
generalizedEigen :: Matrix Double -> Matrix Double -> Either MathError GeneralizedEigen
generalizedEigen Matrix Double
a Matrix Double
b = do
    n <- Matrix Double -> Either MathError Int
square Matrix Double
a
    _ <- square b
    if shape a /= shape b
        then Left (ShapeMismatch (shape a) (shape b))
        else
            if n == 0
                then GeneralizedEigen [] [] <$> fromListChecked 0 0 []
                else unsafePerformIO $
                    withValues (toList a) $ \Ptr CDouble
pa -> [Double]
-> (Ptr CDouble -> IO (Either MathError GeneralizedEigen))
-> IO (Either MathError GeneralizedEigen)
forall a. [Double] -> (Ptr CDouble -> IO a) -> IO a
withValues (Matrix Double -> [Double]
forall a. Matrix a -> [a]
toList Matrix Double
b) ((Ptr CDouble -> IO (Either MathError GeneralizedEigen))
 -> IO (Either MathError GeneralizedEigen))
-> (Ptr CDouble -> IO (Either MathError GeneralizedEigen))
-> IO (Either MathError GeneralizedEigen)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pb -> Int
-> (Ptr CDouble -> IO (Either MathError GeneralizedEigen))
-> IO (Either MathError GeneralizedEigen)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
n ((Ptr CDouble -> IO (Either MathError GeneralizedEigen))
 -> IO (Either MathError GeneralizedEigen))
-> (Ptr CDouble -> IO (Either MathError GeneralizedEigen))
-> IO (Either MathError GeneralizedEigen)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pr -> Int
-> (Ptr CDouble -> IO (Either MathError GeneralizedEigen))
-> IO (Either MathError GeneralizedEigen)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
n ((Ptr CDouble -> IO (Either MathError GeneralizedEigen))
 -> IO (Either MathError GeneralizedEigen))
-> (Ptr CDouble -> IO (Either MathError GeneralizedEigen))
-> IO (Either MathError GeneralizedEigen)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pii -> Int
-> (Ptr CDouble -> IO (Either MathError GeneralizedEigen))
-> IO (Either MathError GeneralizedEigen)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
n ((Ptr CDouble -> IO (Either MathError GeneralizedEigen))
 -> IO (Either MathError GeneralizedEigen))
-> (Ptr CDouble -> IO (Either MathError GeneralizedEigen))
-> IO (Either MathError GeneralizedEigen)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pbe -> Int
-> (Ptr CDouble -> IO (Either MathError GeneralizedEigen))
-> IO (Either MathError GeneralizedEigen)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
1 ((Ptr CDouble -> IO (Either MathError GeneralizedEigen))
 -> IO (Either MathError GeneralizedEigen))
-> (Ptr CDouble -> IO (Either MathError GeneralizedEigen))
-> IO (Either MathError GeneralizedEigen)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pl -> Int
-> (Ptr CDouble -> IO (Either MathError GeneralizedEigen))
-> IO (Either MathError GeneralizedEigen)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n) ((Ptr CDouble -> IO (Either MathError GeneralizedEigen))
 -> IO (Either MathError GeneralizedEigen))
-> (Ptr CDouble -> IO (Either MathError GeneralizedEigen))
-> IO (Either MathError GeneralizedEigen)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pv -> do
                        info <- CInt
-> CChar
-> CChar
-> CInt
-> Ptr CDouble
-> CInt
-> Ptr CDouble
-> CInt
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> CInt
-> Ptr CDouble
-> CInt
-> IO CInt
cGeneralized CInt
101 (Char -> CChar
char Char
'N') (Char -> CChar
char Char
'V') (Int -> CInt
ni Int
n) Ptr CDouble
pa (Int -> CInt
ni Int
n) Ptr CDouble
pb (Int -> CInt
ni Int
n) Ptr CDouble
pr Ptr CDouble
pii Ptr CDouble
pbe Ptr CDouble
pl CInt
1 Ptr CDouble
pv (Int -> CInt
ni Int
n)
                        re <- peekValues n pr
                        im <- peekValues n pii
                        be <- peekValues n pbe
                        vs <- peekValues (n * n) pv
                        pure $ do
                            status "generalizedEigen" info
                            finiteValues (re ++ im ++ be ++ vs)
                            GeneralizedEigen (zipWith (:+) re im) be <$> complexVectors n im vs

{-# NOINLINE schur #-}

{- | Compute real Schur form A=Q*T*transpose(Q), returning orthogonal Q,
quasi-upper-triangular T and complex eigenvalues. Require finite square input and
available LAPACK; failures return Left MathError.
-}
schur :: Matrix Double -> Either MathError Schur
schur :: Matrix Double -> Either MathError Schur
schur Matrix Double
a = do
    n <- Matrix Double -> Either MathError Int
square Matrix Double
a
    if n == 0
        then Schur <$> identityChecked 0 <*> pure a <*> pure []
        else unsafePerformIO $
            withValues (toList a) $ \Ptr CDouble
pa -> (Ptr CInt -> IO (Either MathError Schur))
-> IO (Either MathError Schur)
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr CInt -> IO (Either MathError Schur))
 -> IO (Either MathError Schur))
-> (Ptr CInt -> IO (Either MathError Schur))
-> IO (Either MathError Schur)
forall a b. (a -> b) -> a -> b
$ \Ptr CInt
sd -> Int
-> (Ptr CDouble -> IO (Either MathError Schur))
-> IO (Either MathError Schur)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
n ((Ptr CDouble -> IO (Either MathError Schur))
 -> IO (Either MathError Schur))
-> (Ptr CDouble -> IO (Either MathError Schur))
-> IO (Either MathError Schur)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pr -> Int
-> (Ptr CDouble -> IO (Either MathError Schur))
-> IO (Either MathError Schur)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
n ((Ptr CDouble -> IO (Either MathError Schur))
 -> IO (Either MathError Schur))
-> (Ptr CDouble -> IO (Either MathError Schur))
-> IO (Either MathError Schur)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pii -> Int
-> (Ptr CDouble -> IO (Either MathError Schur))
-> IO (Either MathError Schur)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n) ((Ptr CDouble -> IO (Either MathError Schur))
 -> IO (Either MathError Schur))
-> (Ptr CDouble -> IO (Either MathError Schur))
-> IO (Either MathError Schur)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pq -> do
                info <- CInt
-> CChar
-> CChar
-> FunPtr (ZonkAny 0)
-> CInt
-> Ptr CDouble
-> CInt
-> Ptr CInt
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> CInt
-> IO CInt
forall a.
CInt
-> CChar
-> CChar
-> FunPtr a
-> CInt
-> Ptr CDouble
-> CInt
-> Ptr CInt
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> CInt
-> IO CInt
cSchur CInt
101 (Char -> CChar
char Char
'V') (Char -> CChar
char Char
'N') FunPtr (ZonkAny 0)
forall a. FunPtr a
nullFunPtr (Int -> CInt
ni Int
n) Ptr CDouble
pa (Int -> CInt
ni Int
n) Ptr CInt
sd Ptr CDouble
pr Ptr CDouble
pii Ptr CDouble
pq (Int -> CInt
ni Int
n)
                re <- peekValues n pr
                im <- peekValues n pii
                qs <- peekValues (n * n) pq
                ts <- peekValues (n * n) pa
                pure $ do
                    status "Schur" info
                    finiteValues (re ++ im ++ qs ++ ts)
                    Schur <$> fromListChecked n n qs <*> fromListChecked n n ts <*> pure (zipWith (:+) re im)

{-# NOINLINE qz #-}

{- | Compute generalized real Schur factors for a finite square matrix pair A,B, retaining
Q,Z,S,T and alpha/beta eigenvalues. Shape/backend/LAPACK errors return Left MathError.
-}
qz :: Matrix Double -> Matrix Double -> Either MathError QZ
qz :: Matrix Double -> Matrix Double -> Either MathError QZ
qz Matrix Double
a Matrix Double
b = do
    n <- Matrix Double -> Either MathError Int
square Matrix Double
a
    _ <- square b
    if shape a /= shape b
        then Left (ShapeMismatch (shape a) (shape b))
        else
            if n == 0
                then QZ <$> identityChecked 0 <*> identityChecked 0 <*> pure a <*> pure b <*> pure [] <*> pure []
                else unsafePerformIO $
                    withValues (toList a) $ \Ptr CDouble
pa -> [Double]
-> (Ptr CDouble -> IO (Either MathError QZ))
-> IO (Either MathError QZ)
forall a. [Double] -> (Ptr CDouble -> IO a) -> IO a
withValues (Matrix Double -> [Double]
forall a. Matrix a -> [a]
toList Matrix Double
b) ((Ptr CDouble -> IO (Either MathError QZ))
 -> IO (Either MathError QZ))
-> (Ptr CDouble -> IO (Either MathError QZ))
-> IO (Either MathError QZ)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pb -> (Ptr CInt -> IO (Either MathError QZ)) -> IO (Either MathError QZ)
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr CInt -> IO (Either MathError QZ))
 -> IO (Either MathError QZ))
-> (Ptr CInt -> IO (Either MathError QZ))
-> IO (Either MathError QZ)
forall a b. (a -> b) -> a -> b
$ \Ptr CInt
sd -> Int
-> (Ptr CDouble -> IO (Either MathError QZ))
-> IO (Either MathError QZ)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
n ((Ptr CDouble -> IO (Either MathError QZ))
 -> IO (Either MathError QZ))
-> (Ptr CDouble -> IO (Either MathError QZ))
-> IO (Either MathError QZ)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pr -> Int
-> (Ptr CDouble -> IO (Either MathError QZ))
-> IO (Either MathError QZ)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
n ((Ptr CDouble -> IO (Either MathError QZ))
 -> IO (Either MathError QZ))
-> (Ptr CDouble -> IO (Either MathError QZ))
-> IO (Either MathError QZ)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pii -> Int
-> (Ptr CDouble -> IO (Either MathError QZ))
-> IO (Either MathError QZ)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
n ((Ptr CDouble -> IO (Either MathError QZ))
 -> IO (Either MathError QZ))
-> (Ptr CDouble -> IO (Either MathError QZ))
-> IO (Either MathError QZ)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pbe -> Int
-> (Ptr CDouble -> IO (Either MathError QZ))
-> IO (Either MathError QZ)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n) ((Ptr CDouble -> IO (Either MathError QZ))
 -> IO (Either MathError QZ))
-> (Ptr CDouble -> IO (Either MathError QZ))
-> IO (Either MathError QZ)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pq -> Int
-> (Ptr CDouble -> IO (Either MathError QZ))
-> IO (Either MathError QZ)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n) ((Ptr CDouble -> IO (Either MathError QZ))
 -> IO (Either MathError QZ))
-> (Ptr CDouble -> IO (Either MathError QZ))
-> IO (Either MathError QZ)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pz -> do
                        info <-
                            CInt
-> CChar
-> CChar
-> CChar
-> FunPtr (ZonkAny 1)
-> CInt
-> Ptr CDouble
-> CInt
-> Ptr CDouble
-> CInt
-> Ptr CInt
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> CInt
-> Ptr CDouble
-> CInt
-> IO CInt
forall a.
CInt
-> CChar
-> CChar
-> CChar
-> FunPtr a
-> CInt
-> Ptr CDouble
-> CInt
-> Ptr CDouble
-> CInt
-> Ptr CInt
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> Ptr CDouble
-> CInt
-> Ptr CDouble
-> CInt
-> IO CInt
cQZ
                                CInt
101
                                (Char -> CChar
char Char
'V')
                                (Char -> CChar
char Char
'V')
                                (Char -> CChar
char Char
'N')
                                FunPtr (ZonkAny 1)
forall a. FunPtr a
nullFunPtr
                                (Int -> CInt
ni Int
n)
                                Ptr CDouble
pa
                                (Int -> CInt
ni Int
n)
                                Ptr CDouble
pb
                                (Int -> CInt
ni Int
n)
                                Ptr CInt
sd
                                Ptr CDouble
pr
                                Ptr CDouble
pii
                                Ptr CDouble
pbe
                                Ptr CDouble
pq
                                (Int -> CInt
ni Int
n)
                                Ptr CDouble
pz
                                (Int -> CInt
ni Int
n)
                        re <- peekValues n pr
                        im <- peekValues n pii
                        be <- peekValues n pbe
                        qs <- peekValues (n * n) pq
                        zs <- peekValues (n * n) pz
                        ss <- peekValues (n * n) pa
                        ts <- peekValues (n * n) pb
                        pure $ do
                            status "QZ" info
                            finiteValues (re ++ im ++ be ++ qs ++ zs ++ ss ++ ts)
                            QZ
                                <$> fromListChecked n n qs
                                <*> fromListChecked n n zs
                                <*> fromListChecked n n ss
                                <*> fromListChecked n n ts
                                <*> pure (zipWith (:+) re im)
                                <*> pure be

{-# NOINLINE hessenberg #-}

{- | Reduce finite real square A to upper Hessenberg H with orthogonal Q, satisfying
A=Q*H*transpose(Q). Invalid shape or backend/decomposition failures propagate.
-}
hessenberg :: Matrix Double -> Either MathError Hessenberg
hessenberg :: Matrix Double -> Either MathError Hessenberg
hessenberg Matrix Double
a = do
    n <- Matrix Double -> Either MathError Int
square Matrix Double
a
    if n == 0
        then Hessenberg <$> identityChecked 0 <*> pure a
        else unsafePerformIO $
            withValues (toList a) $ \Ptr CDouble
pa -> Int
-> (Ptr CDouble -> IO (Either MathError Hessenberg))
-> IO (Either MathError Hessenberg)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
1 (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)) ((Ptr CDouble -> IO (Either MathError Hessenberg))
 -> IO (Either MathError Hessenberg))
-> (Ptr CDouble -> IO (Either MathError Hessenberg))
-> IO (Either MathError Hessenberg)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
pt -> do
                info <- CInt
-> CInt
-> CInt
-> CInt
-> Ptr CDouble
-> CInt
-> Ptr CDouble
-> IO CInt
cHessenberg CInt
101 (Int -> CInt
ni Int
n) CInt
1 (Int -> CInt
ni Int
n) Ptr CDouble
pa (Int -> CInt
ni Int
n) Ptr CDouble
pt
                hs <- peekValues (n * n) pa
                infoQ <- cHessenbergQ 101 (ni n) 1 (ni n) pa (ni n) pt
                qs <- peekValues (n * n) pa
                pure $ do
                    status "Hessenberg" info
                    status "Hessenberg Q" infoQ
                    finiteValues (hs ++ qs)
                    h <-
                        fromListChecked
                            n
                            n
                            [if i > j + 1 then 0 else hs !! (i * n + j) | i <- [0 .. n - 1], j <- [0 .. n - 1]]
                    Hessenberg <$> fromListChecked n n qs <*> pure h

{-# NOINLINE balance #-}

{- | Balance a real square matrix by diagonal scaling, returning the balanced matrix and
scales needed to interpret its basis. Invalid finite-input/shape checks or LAPACK
failures return Left MathError.
-}
balance :: Matrix Double -> Either MathError Balance
balance :: Matrix Double -> Either MathError Balance
balance Matrix Double
a = do
    n <- Matrix Double -> Either MathError Int
square Matrix Double
a
    if n == 0
        then Right (Balance a [])
        else unsafePerformIO $
            withValues (toList a) $ \Ptr CDouble
pa -> (Ptr CInt -> IO (Either MathError Balance))
-> IO (Either MathError Balance)
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr CInt -> IO (Either MathError Balance))
 -> IO (Either MathError Balance))
-> (Ptr CInt -> IO (Either MathError Balance))
-> IO (Either MathError Balance)
forall a b. (a -> b) -> a -> b
$ \Ptr CInt
lo -> (Ptr CInt -> IO (Either MathError Balance))
-> IO (Either MathError Balance)
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr CInt -> IO (Either MathError Balance))
 -> IO (Either MathError Balance))
-> (Ptr CInt -> IO (Either MathError Balance))
-> IO (Either MathError Balance)
forall a b. (a -> b) -> a -> b
$ \Ptr CInt
hi -> Int
-> (Ptr CDouble -> IO (Either MathError Balance))
-> IO (Either MathError Balance)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
n ((Ptr CDouble -> IO (Either MathError Balance))
 -> IO (Either MathError Balance))
-> (Ptr CDouble -> IO (Either MathError Balance))
-> IO (Either MathError Balance)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
ps -> do
                info <- CInt
-> CChar
-> CInt
-> Ptr CDouble
-> CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr CDouble
-> IO CInt
cBalance CInt
101 (Char -> CChar
char Char
'S') (Int -> CInt
ni Int
n) Ptr CDouble
pa (Int -> CInt
ni Int
n) Ptr CInt
lo Ptr CInt
hi Ptr CDouble
ps
                xs <- peekValues (n * n) pa
                scales <- peekValues n ps
                pure $ do
                    status "balance" info
                    finiteValues (xs ++ scales)
                    Balance <$> fromListChecked n n xs <*> pure scales

-- | Complex128 thin SVD: A = U diag(s) V*, where complexVT already stores V*.
data ComplexSVD = ComplexSVD
    { ComplexSVD -> Matrix (Complex Double)
complexU :: Matrix (Complex Double)
    -- ^ Complex left singular vectors as columns.
    , ComplexSVD -> [Double]
complexSingularValues :: [Double]
    -- ^ Descending nonnegative real singular values of the complex matrix.
    , ComplexSVD -> Matrix (Complex Double)
complexVT :: Matrix (Complex Double)
    -- ^ Conjugate-transposed right singular vectors in complex SVD.
    }
    deriving (ComplexSVD -> ComplexSVD -> Bool
(ComplexSVD -> ComplexSVD -> Bool)
-> (ComplexSVD -> ComplexSVD -> Bool) -> Eq ComplexSVD
forall a. (a -> a -> Bool) -> (a -> a -> Bool) -> Eq a
$c== :: ComplexSVD -> ComplexSVD -> Bool
== :: ComplexSVD -> ComplexSVD -> Bool
$c/= :: ComplexSVD -> ComplexSVD -> Bool
/= :: ComplexSVD -> ComplexSVD -> Bool
Eq, Int -> ComplexSVD -> ShowS
[ComplexSVD] -> ShowS
ComplexSVD -> String
(Int -> ComplexSVD -> ShowS)
-> (ComplexSVD -> String)
-> ([ComplexSVD] -> ShowS)
-> Show ComplexSVD
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
$cshowsPrec :: Int -> ComplexSVD -> ShowS
showsPrec :: Int -> ComplexSVD -> ShowS
$cshow :: ComplexSVD -> String
show :: ComplexSVD -> String
$cshowList :: [ComplexSVD] -> ShowS
showList :: [ComplexSVD] -> ShowS
Show)

complexParts :: [Complex Double] -> [Double]
complexParts :: [Complex Double] -> [Double]
complexParts = (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])
withComplex :: [Complex Double] -> (Ptr (Complex CDouble) -> IO a) -> IO a
withComplex :: forall a.
[Complex Double] -> (Ptr (Complex CDouble) -> IO a) -> IO a
withComplex [Complex Double]
xs = [Complex CDouble] -> (Ptr (Complex CDouble) -> IO a) -> IO a
forall a b. Storable a => [a] -> (Ptr a -> IO b) -> IO b
withArray [Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Complex Double -> Double
forall a. Complex a -> a
realPart Complex Double
z) CDouble -> CDouble -> Complex CDouble
forall a. a -> a -> Complex a
:+ Double -> CDouble
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Complex Double -> Double
forall a. Complex a -> a
imagPart Complex Double
z) | Complex Double
z <- [Complex Double]
xs]
complexBuffer :: Int -> (Ptr (Complex CDouble) -> IO a) -> IO a
complexBuffer :: forall a. Int -> (Ptr (Complex CDouble) -> IO a) -> IO a
complexBuffer Int
n Ptr (Complex CDouble) -> IO a
action = Int -> (Ptr (Complex CDouble) -> IO a) -> IO a
forall a b. Storable a => Int -> (Ptr a -> IO b) -> IO b
allocaArray (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
1 Int
n) ((Ptr (Complex CDouble) -> IO a) -> IO a)
-> (Ptr (Complex CDouble) -> IO a) -> IO a
forall a b. (a -> b) -> a -> b
$ \Ptr (Complex CDouble)
p -> Ptr (Complex CDouble) -> [Complex CDouble] -> IO ()
forall a. Storable a => Ptr a -> [a] -> IO ()
pokeArray Ptr (Complex CDouble)
p (Int -> Complex CDouble -> [Complex CDouble]
forall a. Int -> a -> [a]
replicate (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
1 Int
n) Complex CDouble
0) IO () -> IO a -> IO a
forall a b. IO a -> IO b -> IO b
forall (m :: * -> *) a b. Monad m => m a -> m b -> m b
>> Ptr (Complex CDouble) -> IO a
action Ptr (Complex CDouble)
p
peekComplex :: Int -> Ptr (Complex CDouble) -> IO [Complex Double]
peekComplex :: Int -> Ptr (Complex CDouble) -> IO [Complex Double]
peekComplex Int
n Ptr (Complex CDouble)
p = (Complex CDouble -> Complex Double)
-> [Complex CDouble] -> [Complex Double]
forall a b. (a -> b) -> [a] -> [b]
map (\Complex CDouble
z -> CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Complex CDouble -> CDouble
forall a. Complex a -> a
realPart Complex CDouble
z) Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ CDouble -> Double
forall a b. (Real a, Fractional b) => a -> b
realToFrac (Complex CDouble -> CDouble
forall a. Complex a -> a
imagPart Complex CDouble
z)) ([Complex CDouble] -> [Complex Double])
-> IO [Complex CDouble] -> IO [Complex Double]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Int -> Ptr (Complex CDouble) -> IO [Complex CDouble]
forall a. Storable a => Int -> Ptr a -> IO [a]
peekArray Int
n Ptr (Complex CDouble)
p
validateComplex :: Matrix (Complex Double) -> Either MathError ()
validateComplex :: Matrix (Complex Double) -> Either MathError ()
validateComplex Matrix (Complex Double)
a = do
    let (Int
m, Int
n) = Matrix (Complex Double) -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix (Complex Double)
a
    [Integer] -> Either MathError ()
allocationSizes [Integer
2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* 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
n]
    if Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
m 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) 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)
        then MathError -> Either MathError ()
forall a b. a -> Either a b
Left (String -> MathError
InvalidArgument String
"complex LAPACKE: dimensions exceed 32-bit integer ABI")
        else
            [Double] -> Either MathError ()
forall a. RealFloat a => [a] -> Either MathError ()
finiteValues ([Complex Double] -> [Double]
complexParts (Matrix (Complex Double) -> [Complex Double]
forall a. Matrix a -> [a]
toList Matrix (Complex Double)
a))
{-# NOINLINE complexSVD #-}

{- | Compute thin complex SVD with conjugate-transposed right factor VT and real nonnegative
singular values. Require finite complex entries and available LAPACK. Invalid storage
sizes/backend/decomposition fails.
-}
complexSVD :: Matrix (Complex Double) -> Either MathError ComplexSVD
complexSVD :: Matrix (Complex Double) -> Either MathError ComplexSVD
complexSVD Matrix (Complex Double)
a = do
    Matrix (Complex Double) -> Either MathError ()
validateComplex Matrix (Complex Double)
a
    let (Int
m, Int
n) = Matrix (Complex Double) -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix (Complex Double)
a; k :: Int
k = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
m Int
n
    [Integer] -> Either MathError ()
allocationSizes [Integer
2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* 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, Integer
2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* 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
k Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
        then Matrix (Complex Double)
-> [Double] -> Matrix (Complex Double) -> ComplexSVD
ComplexSVD (Matrix (Complex Double)
 -> [Double] -> Matrix (Complex Double) -> ComplexSVD)
-> Either MathError (Matrix (Complex Double))
-> Either
     MathError ([Double] -> Matrix (Complex Double) -> ComplexSVD)
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)
fromListChecked Int
m Int
0 [] Either
  MathError ([Double] -> Matrix (Complex Double) -> ComplexSVD)
-> Either MathError [Double]
-> Either MathError (Matrix (Complex Double) -> ComplexSVD)
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 (Complex Double) -> ComplexSVD)
-> Either MathError (Matrix (Complex Double))
-> Either MathError ComplexSVD
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)
fromListChecked Int
0 Int
n []
        else IO (Either MathError ComplexSVD) -> Either MathError ComplexSVD
forall a. IO a -> a
unsafePerformIO (IO (Either MathError ComplexSVD) -> Either MathError ComplexSVD)
-> IO (Either MathError ComplexSVD) -> Either MathError ComplexSVD
forall a b. (a -> b) -> a -> b
$
            [Complex Double]
-> (Ptr (Complex CDouble) -> IO (Either MathError ComplexSVD))
-> IO (Either MathError ComplexSVD)
forall a.
[Complex Double] -> (Ptr (Complex CDouble) -> IO a) -> IO a
withComplex (Matrix (Complex Double) -> [Complex Double]
forall a. Matrix a -> [a]
toList Matrix (Complex Double)
a) ((Ptr (Complex CDouble) -> IO (Either MathError ComplexSVD))
 -> IO (Either MathError ComplexSVD))
-> (Ptr (Complex CDouble) -> IO (Either MathError ComplexSVD))
-> IO (Either MathError ComplexSVD)
forall a b. (a -> b) -> a -> b
$ \Ptr (Complex CDouble)
pa -> Int
-> (Ptr CDouble -> IO (Either MathError ComplexSVD))
-> IO (Either MathError ComplexSVD)
forall a. Int -> (Ptr CDouble -> IO a) -> IO a
buffer Int
k ((Ptr CDouble -> IO (Either MathError ComplexSVD))
 -> IO (Either MathError ComplexSVD))
-> (Ptr CDouble -> IO (Either MathError ComplexSVD))
-> IO (Either MathError ComplexSVD)
forall a b. (a -> b) -> a -> b
$ \Ptr CDouble
ps -> Int
-> (Ptr (Complex CDouble) -> IO (Either MathError ComplexSVD))
-> IO (Either MathError ComplexSVD)
forall a. Int -> (Ptr (Complex CDouble) -> IO a) -> IO a
complexBuffer (Int
m Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
k) ((Ptr (Complex CDouble) -> IO (Either MathError ComplexSVD))
 -> IO (Either MathError ComplexSVD))
-> (Ptr (Complex CDouble) -> IO (Either MathError ComplexSVD))
-> IO (Either MathError ComplexSVD)
forall a b. (a -> b) -> a -> b
$ \Ptr (Complex CDouble)
pu -> Int
-> (Ptr (Complex CDouble) -> IO (Either MathError ComplexSVD))
-> IO (Either MathError ComplexSVD)
forall a. Int -> (Ptr (Complex CDouble) -> IO a) -> IO a
complexBuffer (Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n) ((Ptr (Complex CDouble) -> IO (Either MathError ComplexSVD))
 -> IO (Either MathError ComplexSVD))
-> (Ptr (Complex CDouble) -> IO (Either MathError ComplexSVD))
-> IO (Either MathError ComplexSVD)
forall a b. (a -> b) -> a -> b
$ \Ptr (Complex CDouble)
pv -> do
                info <- CInt
-> CChar
-> CInt
-> CInt
-> Ptr (Complex CDouble)
-> CInt
-> Ptr CDouble
-> Ptr (Complex CDouble)
-> CInt
-> Ptr (Complex CDouble)
-> CInt
-> IO CInt
cComplexSVD CInt
101 (Char -> CChar
char Char
'S') (Int -> CInt
ni Int
m) (Int -> CInt
ni Int
n) Ptr (Complex CDouble)
pa (Int -> CInt
ni Int
n) Ptr CDouble
ps Ptr (Complex CDouble)
pu (Int -> CInt
ni Int
k) Ptr (Complex CDouble)
pv (Int -> CInt
ni Int
n)
                ss <- peekValues k ps
                us <- peekComplex (m * k) pu
                vs <- peekComplex (k * n) pv
                pure $ do
                    status "complex SVD" info
                    finiteValues (ss ++ complexParts us ++ complexParts vs)
                    ComplexSVD <$> fromListChecked m k us <*> pure ss <*> fromListChecked k n vs

{-# NOINLINE complexEigen #-}

{- | Compute complex eigenvalues and right eigenvectors of a finite complex square matrix
via LAPACK. Invalid shape, backend availability, bounds or decomposition errors return
Left MathError.
-}
complexEigen :: Matrix (Complex Double) -> Either MathError Eigen
complexEigen :: Matrix (Complex Double) -> Either MathError Eigen
complexEigen Matrix (Complex Double)
a = do
    Matrix (Complex Double) -> Either MathError ()
validateComplex Matrix (Complex Double)
a
    let (Int
m, Int
n) = Matrix (Complex Double) -> (Int, Int)
forall a. Matrix a -> (Int, Int)
shape Matrix (Complex Double)
a
    if Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
n
        then MathError -> Either MathError Eigen
forall a b. a -> Either a b
Left ((Int, Int) -> MathError
NonSquareMatrix (Int
m, Int
n))
        else
            if Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0
                then [Complex Double] -> Matrix (Complex Double) -> Eigen
Eigen [] (Matrix (Complex Double) -> Eigen)
-> Either MathError (Matrix (Complex Double))
-> Either MathError Eigen
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)
fromListChecked Int
0 Int
0 []
                else IO (Either MathError Eigen) -> Either MathError Eigen
forall a. IO a -> a
unsafePerformIO (IO (Either MathError Eigen) -> Either MathError Eigen)
-> IO (Either MathError Eigen) -> Either MathError Eigen
forall a b. (a -> b) -> a -> b
$
                    [Complex Double]
-> (Ptr (Complex CDouble) -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a.
[Complex Double] -> (Ptr (Complex CDouble) -> IO a) -> IO a
withComplex (Matrix (Complex Double) -> [Complex Double]
forall a. Matrix a -> [a]
toList Matrix (Complex Double)
a) ((Ptr (Complex CDouble) -> IO (Either MathError Eigen))
 -> IO (Either MathError Eigen))
-> (Ptr (Complex CDouble) -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a b. (a -> b) -> a -> b
$ \Ptr (Complex CDouble)
pa -> Int
-> (Ptr (Complex CDouble) -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a. Int -> (Ptr (Complex CDouble) -> IO a) -> IO a
complexBuffer Int
n ((Ptr (Complex CDouble) -> IO (Either MathError Eigen))
 -> IO (Either MathError Eigen))
-> (Ptr (Complex CDouble) -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a b. (a -> b) -> a -> b
$ \Ptr (Complex CDouble)
pw -> Int
-> (Ptr (Complex CDouble) -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a. Int -> (Ptr (Complex CDouble) -> IO a) -> IO a
complexBuffer Int
1 ((Ptr (Complex CDouble) -> IO (Either MathError Eigen))
 -> IO (Either MathError Eigen))
-> (Ptr (Complex CDouble) -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a b. (a -> b) -> a -> b
$ \Ptr (Complex CDouble)
pl -> Int
-> (Ptr (Complex CDouble) -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a. Int -> (Ptr (Complex CDouble) -> IO a) -> IO a
complexBuffer (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
n) ((Ptr (Complex CDouble) -> IO (Either MathError Eigen))
 -> IO (Either MathError Eigen))
-> (Ptr (Complex CDouble) -> IO (Either MathError Eigen))
-> IO (Either MathError Eigen)
forall a b. (a -> b) -> a -> b
$ \Ptr (Complex CDouble)
pv -> do
                        info <- CInt
-> CChar
-> CChar
-> CInt
-> Ptr (Complex CDouble)
-> CInt
-> Ptr (Complex CDouble)
-> Ptr (Complex CDouble)
-> CInt
-> Ptr (Complex CDouble)
-> CInt
-> IO CInt
cComplexEigen CInt
101 (Char -> CChar
char Char
'N') (Char -> CChar
char Char
'V') (Int -> CInt
ni Int
n) Ptr (Complex CDouble)
pa (Int -> CInt
ni Int
n) Ptr (Complex CDouble)
pw Ptr (Complex CDouble)
pl CInt
1 Ptr (Complex CDouble)
pv (Int -> CInt
ni Int
n)
                        ws <- peekComplex n pw
                        vs <- peekComplex (n * n) pv
                        pure $ do
                            status "complex eigen" info
                            finiteValues (complexParts ws ++ complexParts vs)
                            Eigen ws <$> fromListChecked n n vs