{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE ViewPatterns #-}

{-# OPTIONS_GHC -fno-warn-missing-signatures #-}

-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.LinearAlgebra.LAPACK
-- Copyright   :  (c) Alberto Ruiz 2006-14
-- License     :  BSD3
-- Maintainer  :  Alberto Ruiz
-- Stability   :  provisional
--
-- Functional interface to selected LAPACK functions (<http://www.netlib.org/lapack>).
--
-----------------------------------------------------------------------------


module Internal.LAPACK where

import Data.Bifunctor (first)

import Internal.Devel
import Internal.Vector
import Internal.Matrix hiding ((#), (#!))
import Internal.Conversion
import Internal.Element
import Foreign.Ptr(nullPtr)
import Foreign.C.Types
import Control.Monad(when)
import System.IO.Unsafe(unsafePerformIO)

-----------------------------------------------------------------------------------

infixr 1 #
c
a # :: c -> (b -> IO r) -> Trans c b -> IO r
# b -> IO r
b = c -> (b -> IO r) -> Trans c b -> IO r
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
apply c
a b -> IO r
b
{-# INLINE (#) #-}

c
a #! :: c -> c -> Trans c (Trans c (IO r)) -> IO r
#! c
b = c
a c -> (Trans c (IO r) -> IO r) -> Trans c (Trans c (IO r)) -> IO r
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# c
b c -> (IO r -> IO r) -> Trans c (IO r) -> IO r
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# IO r -> IO r
forall a. a -> a
id
{-# INLINE (#!) #-}

-----------------------------------------------------------------------------------

type TMMM t = t ::> t ::> t ::> Ok

type F = Float
type Q = Complex Float

foreign import ccall unsafe "multiplyR" dgemmc :: CInt -> CInt -> TMMM R
foreign import ccall unsafe "multiplyC" zgemmc :: CInt -> CInt -> TMMM C
foreign import ccall unsafe "multiplyF" sgemmc :: CInt -> CInt -> TMMM F
foreign import ccall unsafe "multiplyQ" cgemmc :: CInt -> CInt -> TMMM Q
foreign import ccall unsafe "multiplyI" c_multiplyI :: I -> TMMM I
foreign import ccall unsafe "multiplyL" c_multiplyL :: Z -> TMMM Z

isT :: Matrix t -> p
isT (Matrix t -> Bool
forall t. Matrix t -> Bool
rowOrder -> Bool
False) = p
0
isT Matrix t
_                   = p
1

tt :: Matrix t -> Matrix t
tt x :: Matrix t
x@(Matrix t -> Bool
forall t. Matrix t -> Bool
rowOrder -> Bool
False) = Matrix t
x
tt Matrix t
x                     = Matrix t -> Matrix t
forall t. Matrix t -> Matrix t
trans Matrix t
x

multiplyAux :: (t
 -> t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Matrix t -> Matrix a
multiplyAux t
-> t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
f String
st Matrix t
a Matrix t
b = IO (Matrix a) -> Matrix a
forall a. IO a -> a
unsafePerformIO (IO (Matrix a) -> Matrix a) -> IO (Matrix a) -> Matrix a
forall a b. (a -> b) -> a -> b
$ do
    Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
b) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ String -> IO ()
forall a. HasCallStack => String -> a
error (String -> IO ()) -> String -> IO ()
forall a b. (a -> b) -> a -> b
$ String
"inconsistent dimensions in matrix product "String -> String -> String
forall a. [a] -> [a] -> [a]
++
                                       (Int, Int) -> String
forall a. Show a => a -> String
show (Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
a,Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
a) String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
" x " String -> String -> String
forall a. [a] -> [a] -> [a]
++ (Int, Int) -> String
forall a. Show a => a -> String
show (Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
b, Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
b)
    Matrix a
s <- MatrixOrder -> Int -> Int -> IO (Matrix a)
forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor (Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
a) (Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
b)
    ((Matrix t -> Matrix t
forall t. Matrix t -> Matrix t
tt Matrix t
a) Matrix t
-> ((CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr t
     -> CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr a
     -> IO CInt)
    -> IO CInt)
-> Trans
     (Matrix t)
     (CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr t
      -> CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr a
      -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# (Matrix t -> Matrix t
forall t. Matrix t -> Matrix t
tt Matrix t
b) Matrix t
-> Matrix a
-> Trans (Matrix t) (Trans (Matrix a) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix a
s) (t
-> t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
f (Matrix t -> t
forall p t. Num p => Matrix t -> p
isT Matrix t
a) (Matrix t -> t
forall p t. Num p => Matrix t -> p
isT Matrix t
b)) IO CInt -> String -> IO ()
#| String
st
    Matrix a -> IO (Matrix a)
forall (m :: * -> *) a. Monad m => a -> m a
return Matrix a
s

-- | Matrix product based on BLAS's /dgemm/.
multiplyR :: Matrix Double -> Matrix Double -> Matrix Double
multiplyR :: Matrix Double -> Matrix Double -> Matrix Double
multiplyR Matrix Double
a Matrix Double
b = {-# SCC "multiplyR" #-} (CInt
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String -> Matrix Double -> Matrix Double -> Matrix Double
forall a t t t t.
(Storable a, Storable t, Storable t, Num t, Num t) =>
(t
 -> t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Matrix t -> Matrix a
multiplyAux CInt
-> CInt
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgemmc String
"dgemmc" Matrix Double
a Matrix Double
b

-- | Matrix product based on BLAS's /zgemm/.
multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
multiplyC :: Matrix (Complex Double)
-> Matrix (Complex Double) -> Matrix (Complex Double)
multiplyC Matrix (Complex Double)
a Matrix (Complex Double)
b = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
forall a t t t t.
(Storable a, Storable t, Storable t, Num t, Num t) =>
(t
 -> t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Matrix t -> Matrix a
multiplyAux CInt
-> CInt
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgemmc String
"zgemmc" Matrix (Complex Double)
a Matrix (Complex Double)
b

-- | Matrix product based on BLAS's /sgemm/.
multiplyF :: Matrix Float -> Matrix Float -> Matrix Float
multiplyF :: Matrix Float -> Matrix Float -> Matrix Float
multiplyF Matrix Float
a Matrix Float
b = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Float
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Float
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Float
 -> IO CInt)
-> String -> Matrix Float -> Matrix Float -> Matrix Float
forall a t t t t.
(Storable a, Storable t, Storable t, Num t, Num t) =>
(t
 -> t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Matrix t -> Matrix a
multiplyAux CInt
-> CInt
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Float
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Float
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Float
-> IO CInt
sgemmc String
"sgemmc" Matrix Float
a Matrix Float
b

-- | Matrix product based on BLAS's /cgemm/.
multiplyQ :: Matrix (Complex Float) -> Matrix (Complex Float) -> Matrix (Complex Float)
multiplyQ :: Matrix (Complex Float)
-> Matrix (Complex Float) -> Matrix (Complex Float)
multiplyQ Matrix (Complex Float)
a Matrix (Complex Float)
b = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Float)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Float)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Float)
 -> IO CInt)
-> String
-> Matrix (Complex Float)
-> Matrix (Complex Float)
-> Matrix (Complex Float)
forall a t t t t.
(Storable a, Storable t, Storable t, Num t, Num t) =>
(t
 -> t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Matrix t -> Matrix a
multiplyAux CInt
-> CInt
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Float)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Float)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Float)
-> IO CInt
cgemmc String
"cgemmc" Matrix (Complex Float)
a Matrix (Complex Float)
b

multiplyI :: I -> Matrix CInt -> Matrix CInt -> Matrix CInt
multiplyI :: CInt -> Matrix CInt -> Matrix CInt -> Matrix CInt
multiplyI CInt
m Matrix CInt
a Matrix CInt
b = IO (Matrix CInt) -> Matrix CInt
forall a. IO a -> a
unsafePerformIO (IO (Matrix CInt) -> Matrix CInt)
-> IO (Matrix CInt) -> Matrix CInt
forall a b. (a -> b) -> a -> b
$ do
    Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Matrix CInt -> Int
forall t. Matrix t -> Int
cols Matrix CInt
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Matrix CInt -> Int
forall t. Matrix t -> Int
rows Matrix CInt
b) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ String -> IO ()
forall a. HasCallStack => String -> a
error (String -> IO ()) -> String -> IO ()
forall a b. (a -> b) -> a -> b
$
        String
"inconsistent dimensions in matrix product "String -> String -> String
forall a. [a] -> [a] -> [a]
++ Matrix CInt -> String
forall t. Matrix t -> String
shSize Matrix CInt
a String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
" x " String -> String -> String
forall a. [a] -> [a] -> [a]
++ Matrix CInt -> String
forall t. Matrix t -> String
shSize Matrix CInt
b
    Matrix CInt
s <- MatrixOrder -> Int -> Int -> IO (Matrix CInt)
forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor (Matrix CInt -> Int
forall t. Matrix t -> Int
rows Matrix CInt
a) (Matrix CInt -> Int
forall t. Matrix t -> Int
cols Matrix CInt
b)
    (Matrix CInt
a Matrix CInt
-> ((CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr CInt
     -> CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr CInt
     -> IO CInt)
    -> IO CInt)
-> Trans
     (Matrix CInt)
     (CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr CInt
      -> CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr CInt
      -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Matrix CInt
b Matrix CInt
-> Matrix CInt
-> Trans (Matrix CInt) (Trans (Matrix CInt) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix CInt
s) (CInt -> TMMM CInt
c_multiplyI CInt
m) IO CInt -> String -> IO ()
#|String
"c_multiplyI"
    Matrix CInt -> IO (Matrix CInt)
forall (m :: * -> *) a. Monad m => a -> m a
return Matrix CInt
s

multiplyL :: Z -> Matrix Z -> Matrix Z -> Matrix Z
multiplyL :: Z -> Matrix Z -> Matrix Z -> Matrix Z
multiplyL Z
m Matrix Z
a Matrix Z
b = IO (Matrix Z) -> Matrix Z
forall a. IO a -> a
unsafePerformIO (IO (Matrix Z) -> Matrix Z) -> IO (Matrix Z) -> Matrix Z
forall a b. (a -> b) -> a -> b
$ do
    Bool -> IO () -> IO ()
forall (f :: * -> *). Applicative f => Bool -> f () -> f ()
when (Matrix Z -> Int
forall t. Matrix t -> Int
cols Matrix Z
a Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Matrix Z -> Int
forall t. Matrix t -> Int
rows Matrix Z
b) (IO () -> IO ()) -> IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ String -> IO ()
forall a. HasCallStack => String -> a
error (String -> IO ()) -> String -> IO ()
forall a b. (a -> b) -> a -> b
$
        String
"inconsistent dimensions in matrix product "String -> String -> String
forall a. [a] -> [a] -> [a]
++ Matrix Z -> String
forall t. Matrix t -> String
shSize Matrix Z
a String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
" x " String -> String -> String
forall a. [a] -> [a] -> [a]
++ Matrix Z -> String
forall t. Matrix t -> String
shSize Matrix Z
b
    Matrix Z
s <- MatrixOrder -> Int -> Int -> IO (Matrix Z)
forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor (Matrix Z -> Int
forall t. Matrix t -> Int
rows Matrix Z
a) (Matrix Z -> Int
forall t. Matrix t -> Int
cols Matrix Z
b)
    (Matrix Z
a Matrix Z
-> ((CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr Z
     -> CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr Z
     -> IO CInt)
    -> IO CInt)
-> Trans
     (Matrix Z)
     (CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr Z
      -> CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr Z
      -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Matrix Z
b Matrix Z
-> Matrix Z
-> Trans (Matrix Z) (Trans (Matrix Z) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix Z
s) (Z -> TMMM Z
c_multiplyL Z
m) IO CInt -> String -> IO ()
#|String
"c_multiplyL"
    Matrix Z -> IO (Matrix Z)
forall (m :: * -> *) a. Monad m => a -> m a
return Matrix Z
s

-----------------------------------------------------------------------------

type TSVD t = t ::> t ::> R :> t ::> Ok

foreign import ccall unsafe "svd_l_R" dgesvd :: TSVD R
foreign import ccall unsafe "svd_l_C" zgesvd :: TSVD C
foreign import ccall unsafe "svd_l_Rdd" dgesdd :: TSVD R
foreign import ccall unsafe "svd_l_Cdd" zgesdd :: TSVD C

-- | Full SVD of a real matrix using LAPACK's /dgesvd/.
svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
svdR = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String
-> Matrix Double
-> (Matrix Double, Vector Double, Matrix Double)
forall t a a a.
(Element t, Storable a, Storable a, Storable a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a, Matrix a)
svdAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgesvd String
"svdR"

-- | Full SVD of a real matrix using LAPACK's /dgesdd/.
svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
svdRd = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String
-> Matrix Double
-> (Matrix Double, Vector Double, Matrix Double)
forall t a a a.
(Element t, Storable a, Storable a, Storable a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a, Matrix a)
svdAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgesdd String
"svdRdd"

-- | Full SVD of a complex matrix using LAPACK's /zgesvd/.
svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
svdC :: Matrix (Complex Double)
-> (Matrix (Complex Double), Vector Double,
    Matrix (Complex Double))
svdC = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> (Matrix (Complex Double), Vector Double,
    Matrix (Complex Double))
forall t a a a.
(Element t, Storable a, Storable a, Storable a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a, Matrix a)
svdAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgesvd String
"svdC"

-- | Full SVD of a complex matrix using LAPACK's /zgesdd/.
svdCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
svdCd :: Matrix (Complex Double)
-> (Matrix (Complex Double), Vector Double,
    Matrix (Complex Double))
svdCd = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> (Matrix (Complex Double), Vector Double,
    Matrix (Complex Double))
forall t a a a.
(Element t, Storable a, Storable a, Storable a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a, Matrix a)
svdAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgesdd String
"svdCdd"

svdAux :: (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a, Matrix a)
svdAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
f String
st Matrix t
x = IO (Matrix a, Vector a, Matrix a) -> (Matrix a, Vector a, Matrix a)
forall a. IO a -> a
unsafePerformIO (IO (Matrix a, Vector a, Matrix a)
 -> (Matrix a, Vector a, Matrix a))
-> IO (Matrix a, Vector a, Matrix a)
-> (Matrix a, Vector a, Matrix a)
forall a b. (a -> b) -> a -> b
$ do
    Matrix t
a <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
x
    Matrix a
u <- MatrixOrder -> Int -> Int -> IO (Matrix a)
forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
r Int
r
    Vector a
s <- Int -> IO (Vector a)
forall a. Storable a => Int -> IO (Vector a)
createVector (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
r Int
c)
    Matrix a
v <- MatrixOrder -> Int -> Int -> IO (Matrix a)
forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
c Int
c
    (Matrix t
a Matrix t
-> ((CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr a
     -> CInt
     -> Ptr a
     -> CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr a
     -> IO CInt)
    -> IO CInt)
-> Trans
     (Matrix t)
     (CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr a
      -> CInt
      -> Ptr a
      -> CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr a
      -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Matrix a
u Matrix a
-> ((CInt
     -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)
    -> IO CInt)
-> Trans
     (Matrix a)
     (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Vector a
s Vector a
-> Matrix a
-> Trans (Vector a) (Trans (Matrix a) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix a
v) Trans
  (Matrix t)
  (CInt
   -> CInt
   -> CInt
   -> CInt
   -> Ptr a
   -> CInt
   -> Ptr a
   -> CInt
   -> CInt
   -> CInt
   -> CInt
   -> Ptr a
   -> IO CInt)
CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
f IO CInt -> String -> IO ()
#| String
st
    (Matrix a, Vector a, Matrix a) -> IO (Matrix a, Vector a, Matrix a)
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix a
u,Vector a
s,Matrix a
v)
  where
    r :: Int
r = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
x
    c :: Int
c = Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
x


-- | Thin SVD of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'S\'.
thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
thinSVDR = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String
-> Matrix Double
-> (Matrix Double, Vector Double, Matrix Double)
forall t a a a.
(Element t, Storable a, Storable a, Storable a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a, Matrix a)
thinSVDAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgesvd String
"thinSVDR"

-- | Thin SVD of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'S\'.
thinSVDC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
thinSVDC :: Matrix (Complex Double)
-> (Matrix (Complex Double), Vector Double,
    Matrix (Complex Double))
thinSVDC = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> (Matrix (Complex Double), Vector Double,
    Matrix (Complex Double))
forall t a a a.
(Element t, Storable a, Storable a, Storable a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a, Matrix a)
thinSVDAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgesvd String
"thinSVDC"

-- | Thin SVD of a real matrix, using LAPACK's /dgesdd/ with jobz == \'S\'.
thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
thinSVDRd = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String
-> Matrix Double
-> (Matrix Double, Vector Double, Matrix Double)
forall t a a a.
(Element t, Storable a, Storable a, Storable a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a, Matrix a)
thinSVDAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgesdd String
"thinSVDRdd"

-- | Thin SVD of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'S\'.
thinSVDCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
thinSVDCd :: Matrix (Complex Double)
-> (Matrix (Complex Double), Vector Double,
    Matrix (Complex Double))
thinSVDCd = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> (Matrix (Complex Double), Vector Double,
    Matrix (Complex Double))
forall t a a a.
(Element t, Storable a, Storable a, Storable a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a, Matrix a)
thinSVDAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgesdd String
"thinSVDCdd"

thinSVDAux :: (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a, Matrix a)
thinSVDAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
f String
st Matrix t
x = IO (Matrix a, Vector a, Matrix a) -> (Matrix a, Vector a, Matrix a)
forall a. IO a -> a
unsafePerformIO (IO (Matrix a, Vector a, Matrix a)
 -> (Matrix a, Vector a, Matrix a))
-> IO (Matrix a, Vector a, Matrix a)
-> (Matrix a, Vector a, Matrix a)
forall a b. (a -> b) -> a -> b
$ do
    Matrix t
a <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
x
    Matrix a
u <- MatrixOrder -> Int -> Int -> IO (Matrix a)
forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
r Int
q
    Vector a
s <- Int -> IO (Vector a)
forall a. Storable a => Int -> IO (Vector a)
createVector Int
q
    Matrix a
v <- MatrixOrder -> Int -> Int -> IO (Matrix a)
forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
q Int
c
    (Matrix t
a Matrix t
-> ((CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr a
     -> CInt
     -> Ptr a
     -> CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr a
     -> IO CInt)
    -> IO CInt)
-> Trans
     (Matrix t)
     (CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr a
      -> CInt
      -> Ptr a
      -> CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr a
      -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Matrix a
u Matrix a
-> ((CInt
     -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)
    -> IO CInt)
-> Trans
     (Matrix a)
     (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Vector a
s Vector a
-> Matrix a
-> Trans (Vector a) (Trans (Matrix a) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix a
v) Trans
  (Matrix t)
  (CInt
   -> CInt
   -> CInt
   -> CInt
   -> Ptr a
   -> CInt
   -> Ptr a
   -> CInt
   -> CInt
   -> CInt
   -> CInt
   -> Ptr a
   -> IO CInt)
CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
f IO CInt -> String -> IO ()
#| String
st
    (Matrix a, Vector a, Matrix a) -> IO (Matrix a, Vector a, Matrix a)
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix a
u,Vector a
s,Matrix a
v)
  where
    r :: Int
r = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
x
    c :: Int
c = Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
x
    q :: Int
q = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
r Int
c


-- | Singular values of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'N\'.
svR :: Matrix Double -> Vector Double
svR :: Matrix Double -> Vector Double
svR = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String -> Matrix Double -> Vector Double
forall t a t t t t t t t t a a.
(Element t, Storable a, Num t, Num t, Num t, Num t, Num t, Num t,
 Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Vector a
svAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgesvd String
"svR"

-- | Singular values of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'N\'.
svC :: Matrix (Complex Double) -> Vector Double
svC :: Matrix (Complex Double) -> Vector Double
svC = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String -> Matrix (Complex Double) -> Vector Double
forall t a t t t t t t t t a a.
(Element t, Storable a, Num t, Num t, Num t, Num t, Num t, Num t,
 Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Vector a
svAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgesvd String
"svC"

-- | Singular values of a real matrix, using LAPACK's /dgesdd/ with jobz == \'N\'.
svRd :: Matrix Double -> Vector Double
svRd :: Matrix Double -> Vector Double
svRd = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String -> Matrix Double -> Vector Double
forall t a t t t t t t t t a a.
(Element t, Storable a, Num t, Num t, Num t, Num t, Num t, Num t,
 Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Vector a
svAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgesdd String
"svRd"

-- | Singular values of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'N\'.
svCd :: Matrix (Complex Double) -> Vector Double
svCd :: Matrix (Complex Double) -> Vector Double
svCd = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String -> Matrix (Complex Double) -> Vector Double
forall t a t t t t t t t t a a.
(Element t, Storable a, Num t, Num t, Num t, Num t, Num t, Num t,
 Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Vector a
svAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgesdd String
"svCd"

svAux :: (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Vector a
svAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> t
-> t
-> t
-> t
-> Ptr a
-> CInt
-> Ptr a
-> t
-> t
-> t
-> t
-> Ptr a
-> IO CInt
f String
st Matrix t
x = IO (Vector a) -> Vector a
forall a. IO a -> a
unsafePerformIO (IO (Vector a) -> Vector a) -> IO (Vector a) -> Vector a
forall a b. (a -> b) -> a -> b
$ do
    Matrix t
a <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
x
    Vector a
s <- Int -> IO (Vector a)
forall a. Storable a => Int -> IO (Vector a)
createVector Int
q
    (Matrix t
a Matrix t
-> Vector a
-> Trans (Matrix t) (Trans (Vector a) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Vector a
s) Trans (Matrix t) (Trans (Vector a) (IO CInt))
CInt -> CInt -> CInt -> CInt -> Ptr t -> CInt -> Ptr a -> IO CInt
g IO CInt -> String -> IO ()
#| String
st
    Vector a -> IO (Vector a)
forall (m :: * -> *) a. Monad m => a -> m a
return Vector a
s
  where
    r :: Int
r = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
x
    c :: Int
c = Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
x
    q :: Int
q = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
r Int
c
    g :: CInt -> CInt -> CInt -> CInt -> Ptr t -> CInt -> Ptr a -> IO CInt
g CInt
ra CInt
ca CInt
xra CInt
xca Ptr t
pa CInt
nb Ptr a
pb = CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> t
-> t
-> t
-> t
-> Ptr a
-> CInt
-> Ptr a
-> t
-> t
-> t
-> t
-> Ptr a
-> IO CInt
f CInt
ra CInt
ca CInt
xra CInt
xca Ptr t
pa t
0 t
0 t
0 t
0 Ptr a
forall a. Ptr a
nullPtr CInt
nb Ptr a
pb t
0 t
0 t
0 t
0 Ptr a
forall a. Ptr a
nullPtr


-- | Singular values and all right singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'N\' and jobvt == \'A\'.
rightSVR :: Matrix Double -> (Vector Double, Matrix Double)
rightSVR :: Matrix Double -> (Vector Double, Matrix Double)
rightSVR = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String -> Matrix Double -> (Vector Double, Matrix Double)
forall t a a t t t t a.
(Element t, Storable a, Storable a, Num t, Num t, Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Vector a, Matrix a)
rightSVAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgesvd String
"rightSVR"

-- | Singular values and all right singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'N\' and jobvt == \'A\'.
rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
rightSVC = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> (Vector Double, Matrix (Complex Double))
forall t a a t t t t a.
(Element t, Storable a, Storable a, Num t, Num t, Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Vector a, Matrix a)
rightSVAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgesvd String
"rightSVC"

rightSVAux :: (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Vector a, Matrix a)
rightSVAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> t
-> t
-> t
-> t
-> Ptr a
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
f String
st Matrix t
x = IO (Vector a, Matrix a) -> (Vector a, Matrix a)
forall a. IO a -> a
unsafePerformIO (IO (Vector a, Matrix a) -> (Vector a, Matrix a))
-> IO (Vector a, Matrix a) -> (Vector a, Matrix a)
forall a b. (a -> b) -> a -> b
$ do
    Matrix t
a <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
x
    Vector a
s <- Int -> IO (Vector a)
forall a. Storable a => Int -> IO (Vector a)
createVector Int
q
    Matrix a
v <- MatrixOrder -> Int -> Int -> IO (Matrix a)
forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
c Int
c
    (Matrix t
a Matrix t
-> ((CInt
     -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)
    -> IO CInt)
-> Trans
     (Matrix t)
     (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Vector a
s Vector a
-> Matrix a
-> Trans (Vector a) (Trans (Matrix a) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix a
v) Trans
  (Matrix t)
  (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)
CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
g IO CInt -> String -> IO ()
#| String
st
    (Vector a, Matrix a) -> IO (Vector a, Matrix a)
forall (m :: * -> *) a. Monad m => a -> m a
return (Vector a
s,Matrix a
v)
  where
    r :: Int
r = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
x
    c :: Int
c = Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
x
    q :: Int
q = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
r Int
c
    g :: CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
g CInt
ra CInt
ca CInt
xra CInt
xca Ptr t
pa = CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> t
-> t
-> t
-> t
-> Ptr a
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
f CInt
ra CInt
ca CInt
xra CInt
xca Ptr t
pa t
0 t
0 t
0 t
0 Ptr a
forall a. Ptr a
nullPtr


-- | Singular values and all left singular vectors of a real matrix, using LAPACK's /dgesvd/  with jobu == \'A\' and jobvt == \'N\'.
leftSVR :: Matrix Double -> (Matrix Double, Vector Double)
leftSVR :: Matrix Double -> (Matrix Double, Vector Double)
leftSVR = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String -> Matrix Double -> (Matrix Double, Vector Double)
forall t a a t t t t a.
(Element t, Storable a, Storable a, Num t, Num t, Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a)
leftSVAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgesvd String
"leftSVR"

-- | Singular values and all left singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'A\' and jobvt == \'N\'.
leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double)
leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double)
leftSVC = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> (Matrix (Complex Double), Vector Double)
forall t a a t t t t a.
(Element t, Storable a, Storable a, Num t, Num t, Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a)
leftSVAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgesvd String
"leftSVC"

leftSVAux :: (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a)
leftSVAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> t
-> t
-> t
-> t
-> Ptr a
-> IO CInt
f String
st Matrix t
x = IO (Matrix a, Vector a) -> (Matrix a, Vector a)
forall a. IO a -> a
unsafePerformIO (IO (Matrix a, Vector a) -> (Matrix a, Vector a))
-> IO (Matrix a, Vector a) -> (Matrix a, Vector a)
forall a b. (a -> b) -> a -> b
$ do
    Matrix t
a <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
x
    Matrix a
u <- MatrixOrder -> Int -> Int -> IO (Matrix a)
forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
r Int
r
    Vector a
s <- Int -> IO (Vector a)
forall a. Storable a => Int -> IO (Vector a)
createVector Int
q
    (Matrix t
a Matrix t
-> ((CInt
     -> CInt -> CInt -> CInt -> Ptr a -> CInt -> Ptr a -> IO CInt)
    -> IO CInt)
-> Trans
     (Matrix t)
     (CInt -> CInt -> CInt -> CInt -> Ptr a -> CInt -> Ptr a -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Matrix a
u Matrix a
-> Vector a
-> Trans (Matrix a) (Trans (Vector a) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Vector a
s) Trans
  (Matrix t)
  (CInt -> CInt -> CInt -> CInt -> Ptr a -> CInt -> Ptr a -> IO CInt)
CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> IO CInt
g IO CInt -> String -> IO ()
#| String
st
    (Matrix a, Vector a) -> IO (Matrix a, Vector a)
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix a
u,Vector a
s)
  where
    r :: Int
r = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
x
    c :: Int
c = Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
x
    q :: Int
q = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
r Int
c
    g :: CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> IO CInt
g CInt
ra CInt
ca CInt
xra CInt
xca Ptr t
pa CInt
ru CInt
cu CInt
xru CInt
xcu Ptr a
pu CInt
nb Ptr a
pb = CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> t
-> t
-> t
-> t
-> Ptr a
-> IO CInt
f CInt
ra CInt
ca CInt
xra CInt
xca Ptr t
pa CInt
ru CInt
cu CInt
xru CInt
xcu Ptr a
pu CInt
nb Ptr a
pb t
0 t
0 t
0 t
0 Ptr a
forall a. Ptr a
nullPtr

-----------------------------------------------------------------------------

foreign import ccall unsafe "eig_l_R" dgeev :: R ::> R ::> C :> R ::> Ok
foreign import ccall unsafe "eig_l_G" dggev :: R ::> R ::> C :> R :> R ::> R ::> Ok
foreign import ccall unsafe "eig_l_C" zgeev :: C ::> C ::> C :> C ::> Ok
foreign import ccall unsafe "eig_l_GC" zggev :: C ::> C ::> C :> C :> C ::> C ::> Ok
foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> R :> R ::> Ok
foreign import ccall unsafe "eig_l_H" zheev :: CInt -> R :> C ::> Ok

eigAux :: (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Vector a, Matrix a)
eigAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> t
-> t
-> t
-> t
-> Ptr a
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
f String
st Matrix t
m = IO (Vector a, Matrix a) -> (Vector a, Matrix a)
forall a. IO a -> a
unsafePerformIO (IO (Vector a, Matrix a) -> (Vector a, Matrix a))
-> IO (Vector a, Matrix a) -> (Vector a, Matrix a)
forall a b. (a -> b) -> a -> b
$ do
    Matrix t
a <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
m
    Vector a
l <- Int -> IO (Vector a)
forall a. Storable a => Int -> IO (Vector a)
createVector Int
r
    Matrix a
v <- MatrixOrder -> Int -> Int -> IO (Matrix a)
forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
r Int
r
    (Matrix t
a Matrix t
-> ((CInt
     -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)
    -> IO CInt)
-> Trans
     (Matrix t)
     (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Vector a
l Vector a
-> Matrix a
-> Trans (Vector a) (Trans (Matrix a) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix a
v) Trans
  (Matrix t)
  (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)
CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
g IO CInt -> String -> IO ()
#| String
st
    (Vector a, Matrix a) -> IO (Vector a, Matrix a)
forall (m :: * -> *) a. Monad m => a -> m a
return (Vector a
l,Matrix a
v)
  where
    r :: Int
r = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
m
    g :: CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
g CInt
ra CInt
ca CInt
xra CInt
xca Ptr t
pa = CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> t
-> t
-> t
-> t
-> Ptr a
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
f CInt
ra CInt
ca CInt
xra CInt
xca Ptr t
pa t
0 t
0 t
0 t
0 Ptr a
forall a. Ptr a
nullPtr


-- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/.
-- The eigenvectors are the columns of v. The eigenvalues are not sorted.
eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double))
eigC :: Matrix (Complex Double)
-> (Vector (Complex Double), Matrix (Complex Double))
eigC = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> (Vector (Complex Double), Matrix (Complex Double))
forall t a a t t t t a.
(Element t, Storable a, Storable a, Num t, Num t, Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Vector a, Matrix a)
eigAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgeev String
"eigC"

eigOnlyAux :: (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Vector a
eigOnlyAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> t
-> t
-> t
-> t
-> Ptr a
-> CInt
-> Ptr a
-> t
-> t
-> t
-> t
-> Ptr a
-> IO CInt
f String
st Matrix t
m = IO (Vector a) -> Vector a
forall a. IO a -> a
unsafePerformIO (IO (Vector a) -> Vector a) -> IO (Vector a) -> Vector a
forall a b. (a -> b) -> a -> b
$ do
    Matrix t
a <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
m
    Vector a
l <- Int -> IO (Vector a)
forall a. Storable a => Int -> IO (Vector a)
createVector Int
r
    (Matrix t
a Matrix t
-> Vector a
-> Trans (Matrix t) (Trans (Vector a) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Vector a
l) Trans (Matrix t) (Trans (Vector a) (IO CInt))
CInt -> CInt -> CInt -> CInt -> Ptr t -> CInt -> Ptr a -> IO CInt
g IO CInt -> String -> IO ()
#| String
st
    Vector a -> IO (Vector a)
forall (m :: * -> *) a. Monad m => a -> m a
return Vector a
l
  where
    r :: Int
r = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
m
    g :: CInt -> CInt -> CInt -> CInt -> Ptr t -> CInt -> Ptr a -> IO CInt
g CInt
ra CInt
ca CInt
xra CInt
xca Ptr t
pa CInt
nl Ptr a
pl = CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> t
-> t
-> t
-> t
-> Ptr a
-> CInt
-> Ptr a
-> t
-> t
-> t
-> t
-> Ptr a
-> IO CInt
f CInt
ra CInt
ca CInt
xra CInt
xca Ptr t
pa t
0 t
0 t
0 t
0 Ptr a
forall a. Ptr a
nullPtr CInt
nl Ptr a
pl t
0 t
0 t
0 t
0 Ptr a
forall a. Ptr a
nullPtr

-- | Eigenvalues of a general complex matrix, using LAPACK's /zgeev/ with jobz == \'N\'.
-- The eigenvalues are not sorted.
eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double)
eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double)
eigOnlyC = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String -> Matrix (Complex Double) -> Vector (Complex Double)
forall t a t t t t t t t t a a.
(Element t, Storable a, Num t, Num t, Num t, Num t, Num t, Num t,
 Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Vector a
eigOnlyAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgeev String
"eigOnlyC"

-- | Eigenvalues and right eigenvectors of a general real matrix, using LAPACK's /dgeev/.
-- The eigenvectors are the columns of v. The eigenvalues are not sorted.
eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double))
eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double))
eigR Matrix Double
m = (Vector (Complex Double)
s', Matrix (Complex Double)
v'')
    where (Vector (Complex Double)
s,Matrix Double
v) = Matrix Double -> (Vector (Complex Double), Matrix Double)
eigRaux Matrix Double
m
          s' :: Vector (Complex Double)
s' = Vector (Complex Double) -> Vector (Complex Double)
forall t. RealElement t => Vector (Complex t) -> Vector (Complex t)
fixeig1 Vector (Complex Double)
s
          v' :: [Vector Double]
v' = Matrix Double -> [Vector Double]
forall t. Element t => Matrix t -> [Vector t]
toRows (Matrix Double -> [Vector Double])
-> Matrix Double -> [Vector Double]
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Matrix Double
forall t. Matrix t -> Matrix t
trans Matrix Double
v
          v'' :: Matrix (Complex Double)
v'' = [Vector (Complex Double)] -> Matrix (Complex Double)
forall t. Element t => [Vector t] -> Matrix t
fromColumns ([Vector (Complex Double)] -> Matrix (Complex Double))
-> [Vector (Complex Double)] -> Matrix (Complex Double)
forall a b. (a -> b) -> a -> b
$ [Complex Double] -> [Vector Double] -> [Vector (Complex Double)]
forall e a.
(RealElement e, Num a, Eq a) =>
[Complex a] -> [Vector e] -> [Vector (Complex e)]
fixeig (Vector (Complex Double) -> [Complex Double]
forall a. Storable a => Vector a -> [a]
toList Vector (Complex Double)
s') [Vector Double]
v'

eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)
eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)
eigRaux Matrix Double
m = IO (Vector (Complex Double), Matrix Double)
-> (Vector (Complex Double), Matrix Double)
forall a. IO a -> a
unsafePerformIO (IO (Vector (Complex Double), Matrix Double)
 -> (Vector (Complex Double), Matrix Double))
-> IO (Vector (Complex Double), Matrix Double)
-> (Vector (Complex Double), Matrix Double)
forall a b. (a -> b) -> a -> b
$ do
    Matrix Double
a <- MatrixOrder -> Matrix Double -> IO (Matrix Double)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix Double
m
    Vector (Complex Double)
l <- Int -> IO (Vector (Complex Double))
forall a. Storable a => Int -> IO (Vector a)
createVector Int
r
    Matrix Double
v <- MatrixOrder -> Int -> Int -> IO (Matrix Double)
forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
r Int
r
    (Matrix Double
a Matrix Double
-> ((CInt
     -> Ptr (Complex Double)
     -> CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr Double
     -> IO CInt)
    -> IO CInt)
-> Trans
     (Matrix Double)
     (CInt
      -> Ptr (Complex Double)
      -> CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr Double
      -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Vector (Complex Double)
l Vector (Complex Double)
-> Matrix Double
-> Trans
     (Vector (Complex Double)) (Trans (Matrix Double) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix Double
v) Trans
  (Matrix Double)
  (CInt
   -> Ptr (Complex Double)
   -> CInt
   -> CInt
   -> CInt
   -> CInt
   -> Ptr Double
   -> IO CInt)
CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
g IO CInt -> String -> IO ()
#| String
"eigR"
    (Vector (Complex Double), Matrix Double)
-> IO (Vector (Complex Double), Matrix Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Vector (Complex Double)
l,Matrix Double
v)
  where
    r :: Int
r = Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
m
    g :: CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
g CInt
ra CInt
ca CInt
xra CInt
xca Ptr Double
pa = Double
::> (CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr Double
     -> CInt
     -> Ptr (Complex Double)
     -> CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr Double
     -> IO CInt)
dgeev CInt
ra CInt
ca CInt
xra CInt
xca Ptr Double
pa CInt
0 CInt
0 CInt
0 CInt
0 Ptr Double
forall a. Ptr a
nullPtr

fixeig1 :: Vector (Complex t) -> Vector (Complex t)
fixeig1 Vector (Complex t)
s = (Vector t, Vector t) -> Vector (Complex t)
forall (c :: * -> *) e.
(Complexable c, RealElement e) =>
(c e, c e) -> c (Complex e)
toComplex' (Int -> Int -> Vector t -> Vector t
forall t. Storable t => Int -> Int -> Vector t -> Vector t
subVector Int
0 Int
r (Vector (Complex t) -> Vector t
forall a.
(RealFloat a, Storable a) =>
Vector (Complex a) -> Vector a
asReal Vector (Complex t)
s), Int -> Int -> Vector t -> Vector t
forall t. Storable t => Int -> Int -> Vector t -> Vector t
subVector Int
r Int
r (Vector (Complex t) -> Vector t
forall a.
(RealFloat a, Storable a) =>
Vector (Complex a) -> Vector a
asReal Vector (Complex t)
s))
    where r :: Int
r = Vector (Complex t) -> Int
forall t. Storable t => Vector t -> Int
dim Vector (Complex t)
s

fixeig :: [Complex a] -> [Vector e] -> [Vector (Complex e)]
fixeig  []  [Vector e]
_ =  []
fixeig [Complex a
_] [Vector e
v] = [Vector e -> Vector (Complex e)
forall (c :: * -> *) e.
(Complexable c, RealElement e) =>
c e -> c (Complex e)
comp' Vector e
v]
fixeig ((a
r1:+a
i1):(a
r2:+a
i2):[Complex a]
r) (Vector e
v1:Vector e
v2:[Vector e]
vs)
    | a
r1 a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
r2 Bool -> Bool -> Bool
&& a
i1 a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== (-a
i2) = (Vector e, Vector e) -> Vector (Complex e)
forall (c :: * -> *) e.
(Complexable c, RealElement e) =>
(c e, c e) -> c (Complex e)
toComplex' (Vector e
v1,Vector e
v2) Vector (Complex e) -> [Vector (Complex e)] -> [Vector (Complex e)]
forall a. a -> [a] -> [a]
: (Vector e, Vector e) -> Vector (Complex e)
forall (c :: * -> *) e.
(Complexable c, RealElement e) =>
(c e, c e) -> c (Complex e)
toComplex' (Vector e
v1, (e -> e) -> Vector e -> Vector e
forall a b.
(Storable a, Storable b) =>
(a -> b) -> Vector a -> Vector b
mapVector e -> e
forall a. Num a => a -> a
negate Vector e
v2) Vector (Complex e) -> [Vector (Complex e)] -> [Vector (Complex e)]
forall a. a -> [a] -> [a]
: [Complex a] -> [Vector e] -> [Vector (Complex e)]
fixeig [Complex a]
r [Vector e]
vs
    | Bool
otherwise = Vector e -> Vector (Complex e)
forall (c :: * -> *) e.
(Complexable c, RealElement e) =>
c e -> c (Complex e)
comp' Vector e
v1 Vector (Complex e) -> [Vector (Complex e)] -> [Vector (Complex e)]
forall a. a -> [a] -> [a]
: [Complex a] -> [Vector e] -> [Vector (Complex e)]
fixeig ((a
r2a -> a -> Complex a
forall a. a -> a -> Complex a
:+a
i2)Complex a -> [Complex a] -> [Complex a]
forall a. a -> [a] -> [a]
:[Complex a]
r) (Vector e
v2Vector e -> [Vector e] -> [Vector e]
forall a. a -> [a] -> [a]
:[Vector e]
vs)
fixeig [Complex a]
_ [Vector e]
_ = String -> [Vector (Complex e)]
forall a. HasCallStack => String -> a
error String
"fixeig with impossible inputs"

-- For dggev alpha(i) / beta(i), alpha(i+1) / beta(i+1) form a complex conjugate pair when Im alpha(i) != 0.
-- However, this does not lead to Re alpha(i) == Re alpha(i+1), since beta(i) and beta(i+1)
-- can be different. Therefore old 'fixeig' would fail for 'eigG'.
fixeigG :: [Complex a] -> [Vector e] -> [Vector (Complex e)]
fixeigG  []  [Vector e]
_  = []
fixeigG [Complex a
_] [Vector e
v] = [Vector e -> Vector (Complex e)
forall (c :: * -> *) e.
(Complexable c, RealElement e) =>
c e -> c (Complex e)
comp' Vector e
v]
fixeigG ((a
_:+a
ai1) : Complex a
an : [Complex a]
as) (Vector e
v1:Vector e
v2:[Vector e]
vs)
    | a -> a
forall a. Num a => a -> a
abs a
ai1 a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
1e-13 = (Vector e, Vector e) -> Vector (Complex e)
forall (c :: * -> *) e.
(Complexable c, RealElement e) =>
(c e, c e) -> c (Complex e)
toComplex' (Vector e
v1, Vector e
v2) Vector (Complex e) -> [Vector (Complex e)] -> [Vector (Complex e)]
forall a. a -> [a] -> [a]
: (Vector e, Vector e) -> Vector (Complex e)
forall (c :: * -> *) e.
(Complexable c, RealElement e) =>
(c e, c e) -> c (Complex e)
toComplex' (Vector e
v1, (e -> e) -> Vector e -> Vector e
forall a b.
(Storable a, Storable b) =>
(a -> b) -> Vector a -> Vector b
mapVector e -> e
forall a. Num a => a -> a
negate Vector e
v2) Vector (Complex e) -> [Vector (Complex e)] -> [Vector (Complex e)]
forall a. a -> [a] -> [a]
: [Complex a] -> [Vector e] -> [Vector (Complex e)]
fixeigG [Complex a]
as [Vector e]
vs
    | Bool
otherwise = Vector e -> Vector (Complex e)
forall (c :: * -> *) e.
(Complexable c, RealElement e) =>
c e -> c (Complex e)
comp' Vector e
v1 Vector (Complex e) -> [Vector (Complex e)] -> [Vector (Complex e)]
forall a. a -> [a] -> [a]
: [Complex a] -> [Vector e] -> [Vector (Complex e)]
fixeigG (Complex a
anComplex a -> [Complex a] -> [Complex a]
forall a. a -> [a] -> [a]
:[Complex a]
as) (Vector e
v2Vector e -> [Vector e] -> [Vector e]
forall a. a -> [a] -> [a]
:[Vector e]
vs)
fixeigG [Complex a]
_ [Vector e]
_ = String -> [Vector (Complex e)]
forall a. HasCallStack => String -> a
error String
"fixeigG with impossible inputs"

-- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'.
-- The eigenvalues are not sorted.
eigOnlyR :: Matrix Double -> Vector (Complex Double)
eigOnlyR :: Matrix Double -> Vector (Complex Double)
eigOnlyR = Vector (Complex Double) -> Vector (Complex Double)
forall t. RealElement t => Vector (Complex t) -> Vector (Complex t)
fixeig1 (Vector (Complex Double) -> Vector (Complex Double))
-> (Matrix Double -> Vector (Complex Double))
-> Matrix Double
-> Vector (Complex Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double
 ::> (CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr Double
      -> CInt
      -> Ptr (Complex Double)
      -> CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr Double
      -> IO CInt))
-> String -> Matrix Double -> Vector (Complex Double)
forall t a t t t t t t t t a a.
(Element t, Storable a, Num t, Num t, Num t, Num t, Num t, Num t,
 Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Vector a
eigOnlyAux Double
::> (CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr Double
     -> CInt
     -> Ptr (Complex Double)
     -> CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr Double
     -> IO CInt)
dgeev String
"eigOnlyR"

-- | Generalized eigenvalues and right eigenvectors of a pair of real matrices, using LAPACK's /dggev/.
-- The eigenvectors are the columns of v. The eigenvalues are represented as alphas / betas and not sorted.
eigG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double, Matrix (Complex Double))
eigG :: Matrix Double
-> Matrix Double
-> (Vector (Complex Double), Vector Double,
    Matrix (Complex Double))
eigG Matrix Double
a Matrix Double
b = (Vector (Complex Double)
alpha', Vector Double
beta, Matrix (Complex Double)
v'')
  where
    (Vector (Complex Double)
alpha, Vector Double
beta, Matrix Double
v) = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> Matrix Double
-> Matrix Double
-> String
-> (Vector (Complex Double), Vector Double, Matrix Double)
forall t t a a a t t t t a.
(Element t, Element t, Storable a, Storable a, Storable a, Num t,
 Num t, Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> Matrix t -> Matrix t -> String -> (Vector a, Vector a, Matrix a)
eigGaux CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dggev Matrix Double
a Matrix Double
b String
"eigG"
    alpha' :: Vector (Complex Double)
alpha' = Vector (Complex Double) -> Vector (Complex Double)
forall t. RealElement t => Vector (Complex t) -> Vector (Complex t)
fixeig1 Vector (Complex Double)
alpha
    v' :: [Vector Double]
v' = Matrix Double -> [Vector Double]
forall t. Element t => Matrix t -> [Vector t]
toRows (Matrix Double -> [Vector Double])
-> Matrix Double -> [Vector Double]
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Matrix Double
forall t. Matrix t -> Matrix t
trans Matrix Double
v
    v'' :: Matrix (Complex Double)
v'' = [Vector (Complex Double)] -> Matrix (Complex Double)
forall t. Element t => [Vector t] -> Matrix t
fromColumns ([Vector (Complex Double)] -> Matrix (Complex Double))
-> [Vector (Complex Double)] -> Matrix (Complex Double)
forall a b. (a -> b) -> a -> b
$ [Complex Double] -> [Vector Double] -> [Vector (Complex Double)]
forall e a.
(RealElement e, Ord a, Fractional a) =>
[Complex a] -> [Vector e] -> [Vector (Complex e)]
fixeigG (Vector (Complex Double) -> [Complex Double]
forall a. Storable a => Vector a -> [a]
toList Vector (Complex Double)
alpha') [Vector Double]
v'

eigGaux :: (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> Matrix t -> Matrix t -> String -> (Vector a, Vector a, Matrix a)
eigGaux CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> t
-> t
-> t
-> t
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
f Matrix t
ma Matrix t
mb String
st = IO (Vector a, Vector a, Matrix a) -> (Vector a, Vector a, Matrix a)
forall a. IO a -> a
unsafePerformIO (IO (Vector a, Vector a, Matrix a)
 -> (Vector a, Vector a, Matrix a))
-> IO (Vector a, Vector a, Matrix a)
-> (Vector a, Vector a, Matrix a)
forall a b. (a -> b) -> a -> b
$ do
    Matrix t
a <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
ma
    Matrix t
b <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
mb
    Vector a
alpha <- Int -> IO (Vector a)
forall a. Storable a => Int -> IO (Vector a)
createVector Int
r
    Vector a
beta <- Int -> IO (Vector a)
forall a. Storable a => Int -> IO (Vector a)
createVector Int
r
    Matrix a
vr <- MatrixOrder -> Int -> Int -> IO (Matrix a)
forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
r Int
r

    (Matrix t
a Matrix t
-> ((CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr t
     -> CInt
     -> Ptr a
     -> CInt
     -> Ptr a
     -> CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr a
     -> IO CInt)
    -> IO CInt)
-> Trans
     (Matrix t)
     (CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr t
      -> CInt
      -> Ptr a
      -> CInt
      -> Ptr a
      -> CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr a
      -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Matrix t
b Matrix t
-> ((CInt
     -> Ptr a
     -> CInt
     -> Ptr a
     -> CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr a
     -> IO CInt)
    -> IO CInt)
-> Trans
     (Matrix t)
     (CInt
      -> Ptr a
      -> CInt
      -> Ptr a
      -> CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr a
      -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Vector a
alpha Vector a
-> ((CInt
     -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)
    -> IO CInt)
-> Trans
     (Vector a)
     (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr a -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Vector a
beta Vector a
-> Matrix a
-> Trans (Vector a) (Trans (Matrix a) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix a
vr) Trans
  (Matrix t)
  (CInt
   -> CInt
   -> CInt
   -> CInt
   -> Ptr t
   -> CInt
   -> Ptr a
   -> CInt
   -> Ptr a
   -> CInt
   -> CInt
   -> CInt
   -> CInt
   -> Ptr a
   -> IO CInt)
CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
g IO CInt -> String -> IO ()
#| String
st

    (Vector a, Vector a, Matrix a) -> IO (Vector a, Vector a, Matrix a)
forall (m :: * -> *) a. Monad m => a -> m a
return (Vector a
alpha, Vector a
beta, Matrix a
vr)
  where
    r :: Int
r = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
ma
    g :: CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
g CInt
ar CInt
ac CInt
xra CInt
xca Ptr t
pa CInt
br CInt
bc CInt
xrb CInt
xcb Ptr t
pb CInt
alphan Ptr a
palpha CInt
betan Ptr a
pbeta = CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> t
-> t
-> t
-> t
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
f CInt
ar CInt
ac CInt
xra CInt
xca Ptr t
pa CInt
br CInt
bc CInt
xrb CInt
xcb Ptr t
pb CInt
alphan Ptr a
palpha CInt
betan Ptr a
pbeta t
0 t
0 t
0 t
0 Ptr a
forall a. Ptr a
nullPtr 

eigGOnlyAux :: (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> Matrix t -> Matrix t -> String -> (Vector a, Vector a)
eigGOnlyAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> t
-> t
-> t
-> t
-> Ptr a
-> t
-> t
-> t
-> t
-> Ptr a
-> IO CInt
f Matrix t
ma Matrix t
mb String
st = IO (Vector a, Vector a) -> (Vector a, Vector a)
forall a. IO a -> a
unsafePerformIO (IO (Vector a, Vector a) -> (Vector a, Vector a))
-> IO (Vector a, Vector a) -> (Vector a, Vector a)
forall a b. (a -> b) -> a -> b
$ do
    Matrix t
a <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
ma
    Matrix t
b <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
mb
    Vector a
alpha <- Int -> IO (Vector a)
forall a. Storable a => Int -> IO (Vector a)
createVector Int
r
    Vector a
beta <- Int -> IO (Vector a)
forall a. Storable a => Int -> IO (Vector a)
createVector Int
r

    (Matrix t
a Matrix t
-> ((CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr t
     -> CInt
     -> Ptr a
     -> CInt
     -> Ptr a
     -> IO CInt)
    -> IO CInt)
-> Trans
     (Matrix t)
     (CInt
      -> CInt
      -> CInt
      -> CInt
      -> Ptr t
      -> CInt
      -> Ptr a
      -> CInt
      -> Ptr a
      -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Matrix t
b Matrix t
-> ((CInt -> Ptr a -> CInt -> Ptr a -> IO CInt) -> IO CInt)
-> Trans (Matrix t) (CInt -> Ptr a -> CInt -> Ptr a -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Vector a
alpha Vector a
-> Vector a
-> Trans (Vector a) (Trans (Vector a) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Vector a
beta) Trans
  (Matrix t)
  (CInt
   -> CInt
   -> CInt
   -> CInt
   -> Ptr t
   -> CInt
   -> Ptr a
   -> CInt
   -> Ptr a
   -> IO CInt)
CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> IO CInt
g IO CInt -> String -> IO ()
#| String
st

    (Vector a, Vector a) -> IO (Vector a, Vector a)
forall (m :: * -> *) a. Monad m => a -> m a
return (Vector a
alpha, Vector a
beta)
  where
    r :: Int
r = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
ma
    g :: CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> IO CInt
g CInt
ar CInt
ac CInt
xra CInt
xca Ptr t
pa CInt
br CInt
bc CInt
xrb CInt
xcb Ptr t
pb CInt
alphan Ptr a
palpha CInt
betan Ptr a
pbeta = CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> t
-> t
-> t
-> t
-> Ptr a
-> t
-> t
-> t
-> t
-> Ptr a
-> IO CInt
f CInt
ar CInt
ac CInt
xra CInt
xca Ptr t
pa CInt
br CInt
bc CInt
xrb CInt
xcb Ptr t
pb CInt
alphan Ptr a
palpha CInt
betan Ptr a
pbeta t
0 t
0 t
0 t
0 Ptr a
forall a. Ptr a
nullPtr t
0 t
0 t
0 t
0 Ptr a
forall a. Ptr a
nullPtr

-- | Generalized eigenvalues and right eigenvectors of a pair of complex matrices, using LAPACK's /zggev/.
-- The eigenvectors are the columns of v. The eigenvalues are represented as alphas / betas and not sorted.
eigGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double), Matrix (Complex Double))
eigGC :: Matrix (Complex Double)
-> Matrix (Complex Double)
-> (Vector (Complex Double), Vector (Complex Double),
    Matrix (Complex Double))
eigGC Matrix (Complex Double)
a Matrix (Complex Double)
b = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> String
-> (Vector (Complex Double), Vector (Complex Double),
    Matrix (Complex Double))
forall t t a a a t t t t a.
(Element t, Element t, Storable a, Storable a, Storable a, Num t,
 Num t, Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> Matrix t -> Matrix t -> String -> (Vector a, Vector a, Matrix a)
eigGaux CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zggev Matrix (Complex Double)
a Matrix (Complex Double)
b String
"eigGC"

eigOnlyG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double)
eigOnlyG :: Matrix Double
-> Matrix Double -> (Vector (Complex Double), Vector Double)
eigOnlyG Matrix Double
a Matrix Double
b = (Vector (Complex Double) -> Vector (Complex Double))
-> (Vector (Complex Double), Vector Double)
-> (Vector (Complex Double), Vector Double)
forall (p :: * -> * -> *) a b c.
Bifunctor p =>
(a -> b) -> p a c -> p b c
first Vector (Complex Double) -> Vector (Complex Double)
forall t. RealElement t => Vector (Complex t) -> Vector (Complex t)
fixeig1 ((Vector (Complex Double), Vector Double)
 -> (Vector (Complex Double), Vector Double))
-> (Vector (Complex Double), Vector Double)
-> (Vector (Complex Double), Vector Double)
forall a b. (a -> b) -> a -> b
$ (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> Matrix Double
-> Matrix Double
-> String
-> (Vector (Complex Double), Vector Double)
forall t t a a t t t t t t t t a a.
(Element t, Element t, Storable a, Storable a, Num t, Num t, Num t,
 Num t, Num t, Num t, Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> Matrix t -> Matrix t -> String -> (Vector a, Vector a)
eigGOnlyAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dggev Matrix Double
a Matrix Double
b String
"eigOnlyG"

eigOnlyGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double))
eigOnlyGC :: Matrix (Complex Double)
-> Matrix (Complex Double)
-> (Vector (Complex Double), Vector (Complex Double))
eigOnlyGC Matrix (Complex Double)
a Matrix (Complex Double)
b = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> String
-> (Vector (Complex Double), Vector (Complex Double))
forall t t a a t t t t t t t t a a.
(Element t, Element t, Storable a, Storable a, Num t, Num t, Num t,
 Num t, Num t, Num t, Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> Matrix t -> Matrix t -> String -> (Vector a, Vector a)
eigGOnlyAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zggev Matrix (Complex Double)
a Matrix (Complex Double)
b String
"eigOnlyGC"

-----------------------------------------------------------------------------

eigSHAux :: (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Vector a, Matrix t)
eigSHAux CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt
f String
st Matrix t
m = IO (Vector a, Matrix t) -> (Vector a, Matrix t)
forall a. IO a -> a
unsafePerformIO (IO (Vector a, Matrix t) -> (Vector a, Matrix t))
-> IO (Vector a, Matrix t) -> (Vector a, Matrix t)
forall a b. (a -> b) -> a -> b
$ do
    Vector a
l <- Int -> IO (Vector a)
forall a. Storable a => Int -> IO (Vector a)
createVector Int
r
    Matrix t
v <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
m
    (Vector a
l Vector a
-> Matrix t
-> Trans (Vector a) (Trans (Matrix t) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix t
v) Trans (Vector a) (Trans (Matrix t) (IO CInt))
CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt
f IO CInt -> String -> IO ()
#| String
st
    (Vector a, Matrix t) -> IO (Vector a, Matrix t)
forall (m :: * -> *) a. Monad m => a -> m a
return (Vector a
l,Matrix t
v)
  where
    r :: Int
r = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
m

-- | Eigenvalues and right eigenvectors of a symmetric real matrix, using LAPACK's /dsyev/.
-- The eigenvectors are the columns of v.
-- The eigenvalues are sorted in descending order (use 'eigS'' for ascending order).
eigS :: Matrix Double -> (Vector Double, Matrix Double)
eigS :: Matrix Double -> (Vector Double, Matrix Double)
eigS Matrix Double
m = (Vector Double
s', Matrix Double -> Matrix Double
forall t. Element t => Matrix t -> Matrix t
fliprl Matrix Double
v)
    where (Vector Double
s,Matrix Double
v) = Matrix Double -> (Vector Double, Matrix Double)
eigS' Matrix Double
m
          s' :: Vector Double
s' = [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList ([Double] -> Vector Double)
-> (Vector Double -> [Double]) -> Vector Double -> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> [Double]
forall a. [a] -> [a]
reverse ([Double] -> [Double])
-> (Vector Double -> [Double]) -> Vector Double -> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> [Double]
forall a. Storable a => Vector a -> [a]
toList (Vector Double -> Vector Double) -> Vector Double -> Vector Double
forall a b. (a -> b) -> a -> b
$  Vector Double
s

-- | 'eigS' in ascending order
eigS' :: Matrix Double -> (Vector Double, Matrix Double)
eigS' :: Matrix Double -> (Vector Double, Matrix Double)
eigS' = (CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String -> Matrix Double -> (Vector Double, Matrix Double)
forall t a.
(Element t, Storable a) =>
(CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Vector a, Matrix t)
eigSHAux (CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dsyev CInt
1) String
"eigS'"

-- | Eigenvalues and right eigenvectors of a hermitian complex matrix, using LAPACK's /zheev/.
-- The eigenvectors are the columns of v.
-- The eigenvalues are sorted in descending order (use 'eigH'' for ascending order).
eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
eigH Matrix (Complex Double)
m = (Vector Double
s', Matrix (Complex Double) -> Matrix (Complex Double)
forall t. Element t => Matrix t -> Matrix t
fliprl Matrix (Complex Double)
v)
  where
    (Vector Double
s,Matrix (Complex Double)
v) = Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
eigH' Matrix (Complex Double)
m
    s' :: Vector Double
s' = [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList ([Double] -> Vector Double)
-> (Vector Double -> [Double]) -> Vector Double -> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Double] -> [Double]
forall a. [a] -> [a]
reverse ([Double] -> [Double])
-> (Vector Double -> [Double]) -> Vector Double -> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector Double -> [Double]
forall a. Storable a => Vector a -> [a]
toList (Vector Double -> Vector Double) -> Vector Double -> Vector Double
forall a b. (a -> b) -> a -> b
$  Vector Double
s

-- | 'eigH' in ascending order
eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
eigH' = (CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> (Vector Double, Matrix (Complex Double))
forall t a.
(Element t, Storable a) =>
(CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Vector a, Matrix t)
eigSHAux (CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zheev CInt
1) String
"eigH'"


-- | Eigenvalues of a symmetric real matrix, using LAPACK's /dsyev/ with jobz == \'N\'.
-- The eigenvalues are sorted in descending order.
eigOnlyS :: Matrix Double -> Vector Double
eigOnlyS :: Matrix Double -> Vector Double
eigOnlyS = Vector Double -> Vector Double
vrev (Vector Double -> Vector Double)
-> (Matrix Double -> Vector Double)
-> Matrix Double
-> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vector Double, Matrix Double) -> Vector Double
forall a b. (a, b) -> a
fst((Vector Double, Matrix Double) -> Vector Double)
-> (Matrix Double -> (Vector Double, Matrix Double))
-> Matrix Double
-> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String -> Matrix Double -> (Vector Double, Matrix Double)
forall t a.
(Element t, Storable a) =>
(CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Vector a, Matrix t)
eigSHAux (CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dsyev CInt
0) String
"eigS'"

-- | Eigenvalues of a hermitian complex matrix, using LAPACK's /zheev/ with jobz == \'N\'.
-- The eigenvalues are sorted in descending order.
eigOnlyH :: Matrix (Complex Double) -> Vector Double
eigOnlyH :: Matrix (Complex Double) -> Vector Double
eigOnlyH = Vector Double -> Vector Double
vrev (Vector Double -> Vector Double)
-> (Matrix (Complex Double) -> Vector Double)
-> Matrix (Complex Double)
-> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vector Double, Matrix (Complex Double)) -> Vector Double
forall a b. (a, b) -> a
fst((Vector Double, Matrix (Complex Double)) -> Vector Double)
-> (Matrix (Complex Double)
    -> (Vector Double, Matrix (Complex Double)))
-> Matrix (Complex Double)
-> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> (Vector Double, Matrix (Complex Double))
forall t a.
(Element t, Storable a) =>
(CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Vector a, Matrix t)
eigSHAux (CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zheev CInt
0) String
"eigH'"

vrev :: Vector Double -> Vector Double
vrev = Matrix Double -> Vector Double
forall t. Element t => Matrix t -> Vector t
flatten (Matrix Double -> Vector Double)
-> (Vector Double -> Matrix Double)
-> Vector Double
-> Vector Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix Double -> Matrix Double
forall t. Element t => Matrix t -> Matrix t
flipud (Matrix Double -> Matrix Double)
-> (Vector Double -> Matrix Double)
-> Vector Double
-> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Vector Double -> Matrix Double
forall t. Storable t => Int -> Vector t -> Matrix t
reshape Int
1

-----------------------------------------------------------------------------
foreign import ccall unsafe "linearSolveR_l" dgesv :: R ::> R ::> Ok
foreign import ccall unsafe "linearSolveC_l" zgesv :: C ::> C ::> Ok

linearSolveSQAux :: (IO (Matrix t) -> IO p)
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> p
linearSolveSQAux IO (Matrix t) -> IO p
g CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt
f String
st Matrix t
a Matrix t
b
    | Int
n1Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
n2 Bool -> Bool -> Bool
&& Int
n1Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
r = IO p -> p
forall a. IO a -> a
unsafePerformIO (IO p -> p) -> (IO (Matrix t) -> IO p) -> IO (Matrix t) -> p
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IO (Matrix t) -> IO p
g (IO (Matrix t) -> p) -> IO (Matrix t) -> p
forall a b. (a -> b) -> a -> b
$ do
        Matrix t
a' <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
a
        Matrix t
s  <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
b
        (Matrix t
a' Matrix t
-> Matrix t
-> Trans (Matrix t) (Trans (Matrix t) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix t
s) Trans (Matrix t) (Trans (Matrix t) (IO CInt))
CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt
f IO CInt -> String -> IO ()
#| String
st
        Matrix t -> IO (Matrix t)
forall (m :: * -> *) a. Monad m => a -> m a
return Matrix t
s
    | Bool
otherwise = String -> p
forall a. HasCallStack => String -> a
error (String -> p) -> String -> p
forall a b. (a -> b) -> a -> b
$ String
st String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
" of nonsquare matrix"
  where
    n1 :: Int
n1 = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
a
    n2 :: Int
n2 = Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
a
    r :: Int
r  = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
b

-- | Solve a real linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /dgesv/. For underconstrained or overconstrained systems use 'linearSolveLSR' or 'linearSolveSVDR'. See also 'lusR'.
linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double
linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double
linearSolveR Matrix Double
a Matrix Double
b = (IO (Matrix Double) -> IO (Matrix Double))
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr Double
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr Double
    -> IO CInt)
-> String
-> Matrix Double
-> Matrix Double
-> Matrix Double
forall t t p.
(Element t, Element t) =>
(IO (Matrix t) -> IO p)
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> p
linearSolveSQAux IO (Matrix Double) -> IO (Matrix Double)
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgesv String
"linearSolveR" Matrix Double
a Matrix Double
b

mbLinearSolveR :: Matrix Double -> Matrix Double -> Maybe (Matrix Double)
mbLinearSolveR :: Matrix Double -> Matrix Double -> Maybe (Matrix Double)
mbLinearSolveR Matrix Double
a Matrix Double
b = (IO (Matrix Double) -> IO (Maybe (Matrix Double)))
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr Double
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr Double
    -> IO CInt)
-> String
-> Matrix Double
-> Matrix Double
-> Maybe (Matrix Double)
forall t t p.
(Element t, Element t) =>
(IO (Matrix t) -> IO p)
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> p
linearSolveSQAux IO (Matrix Double) -> IO (Maybe (Matrix Double))
forall x. IO x -> IO (Maybe x)
mbCatch CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgesv String
"linearSolveR" Matrix Double
a Matrix Double
b


-- | Solve a complex linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition, based on LAPACK's /zgesv/. For underconstrained or overconstrained systems use 'linearSolveLSC' or 'linearSolveSVDC'. See also 'lusC'.
linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
linearSolveC :: Matrix (Complex Double)
-> Matrix (Complex Double) -> Matrix (Complex Double)
linearSolveC Matrix (Complex Double)
a Matrix (Complex Double)
b = (IO (Matrix (Complex Double)) -> IO (Matrix (Complex Double)))
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr (Complex Double)
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr (Complex Double)
    -> IO CInt)
-> String
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
forall t t p.
(Element t, Element t) =>
(IO (Matrix t) -> IO p)
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> p
linearSolveSQAux IO (Matrix (Complex Double)) -> IO (Matrix (Complex Double))
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgesv String
"linearSolveC" Matrix (Complex Double)
a Matrix (Complex Double)
b

mbLinearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
mbLinearSolveC :: Matrix (Complex Double)
-> Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
mbLinearSolveC Matrix (Complex Double)
a Matrix (Complex Double)
b = (IO (Matrix (Complex Double))
 -> IO (Maybe (Matrix (Complex Double))))
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr (Complex Double)
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr (Complex Double)
    -> IO CInt)
-> String
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Maybe (Matrix (Complex Double))
forall t t p.
(Element t, Element t) =>
(IO (Matrix t) -> IO p)
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> p
linearSolveSQAux IO (Matrix (Complex Double))
-> IO (Maybe (Matrix (Complex Double)))
forall x. IO x -> IO (Maybe x)
mbCatch CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgesv String
"linearSolveC" Matrix (Complex Double)
a Matrix (Complex Double)
b

--------------------------------------------------------------------------------
foreign import ccall unsafe "cholSolveR_l" dpotrs  :: R ::> R ::> Ok
foreign import ccall unsafe "cholSolveC_l" zpotrs  :: C ::> C ::> Ok


linearSolveSQAux2 :: (IO (Matrix t) -> IO p)
-> Trans
     (Matrix t) (CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> p
linearSolveSQAux2 IO (Matrix t) -> IO p
g Trans (Matrix t) (CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
f String
st Matrix t
a Matrix t
b
    | Int
n1Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
n2 Bool -> Bool -> Bool
&& Int
n1Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
r = IO p -> p
forall a. IO a -> a
unsafePerformIO (IO p -> p) -> (IO (Matrix t) -> IO p) -> IO (Matrix t) -> p
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IO (Matrix t) -> IO p
g (IO (Matrix t) -> p) -> IO (Matrix t) -> p
forall a b. (a -> b) -> a -> b
$ do
        Matrix t
s <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
b
        (Matrix t
a Matrix t
-> Matrix t
-> Trans (Matrix t) (Trans (Matrix t) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix t
s) Trans (Matrix t) (Trans (Matrix t) (IO CInt))
Trans (Matrix t) (CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
f IO CInt -> String -> IO ()
#| String
st
        Matrix t -> IO (Matrix t)
forall (m :: * -> *) a. Monad m => a -> m a
return Matrix t
s
    | Bool
otherwise = String -> p
forall a. HasCallStack => String -> a
error (String -> p) -> String -> p
forall a b. (a -> b) -> a -> b
$ String
st String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
" of nonsquare matrix"
  where
    n1 :: Int
n1 = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
a
    n2 :: Int
n2 = Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
a
    r :: Int
r  = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
b

-- | Solves a symmetric positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholS'.
cholSolveR :: Matrix Double -> Matrix Double -> Matrix Double
cholSolveR :: Matrix Double -> Matrix Double -> Matrix Double
cholSolveR Matrix Double
a Matrix Double
b = (IO (Matrix Double) -> IO (Matrix Double))
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr Double
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr Double
    -> IO CInt)
-> String
-> Matrix Double
-> Matrix Double
-> Matrix Double
forall t t p.
(Element t, Storable t) =>
(IO (Matrix t) -> IO p)
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> p
linearSolveSQAux2 IO (Matrix Double) -> IO (Matrix Double)
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dpotrs String
"cholSolveR" (Matrix Double -> Matrix Double
forall t. Element t => Matrix t -> Matrix t
fmat Matrix Double
a) Matrix Double
b

-- | Solves a Hermitian positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholH'.
cholSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
cholSolveC :: Matrix (Complex Double)
-> Matrix (Complex Double) -> Matrix (Complex Double)
cholSolveC Matrix (Complex Double)
a Matrix (Complex Double)
b = (IO (Matrix (Complex Double)) -> IO (Matrix (Complex Double)))
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr (Complex Double)
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr (Complex Double)
    -> IO CInt)
-> String
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
forall t t p.
(Element t, Storable t) =>
(IO (Matrix t) -> IO p)
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> p
linearSolveSQAux2 IO (Matrix (Complex Double)) -> IO (Matrix (Complex Double))
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zpotrs String
"cholSolveC" (Matrix (Complex Double) -> Matrix (Complex Double)
forall t. Element t => Matrix t -> Matrix t
fmat Matrix (Complex Double)
a) Matrix (Complex Double)
b

--------------------------------------------------------------------------------
foreign import ccall unsafe "triSolveR_l_u" dtrtrs_u  :: R ::> R ::> Ok
foreign import ccall unsafe "triSolveC_l_u" ztrtrs_u  :: C ::> C ::> Ok
foreign import ccall unsafe "triSolveR_l_l" dtrtrs_l  :: R ::> R ::> Ok
foreign import ccall unsafe "triSolveC_l_l" ztrtrs_l  :: C ::> C ::> Ok


linearSolveTRAux2 :: (IO (Matrix t) -> IO p)
-> Trans
     (Matrix t) (CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> p
linearSolveTRAux2 IO (Matrix t) -> IO p
g Trans (Matrix t) (CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
f String
st Matrix t
a Matrix t
b
    | Int
n1Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
n2 Bool -> Bool -> Bool
&& Int
n1Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
r = IO p -> p
forall a. IO a -> a
unsafePerformIO (IO p -> p) -> (IO (Matrix t) -> IO p) -> IO (Matrix t) -> p
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IO (Matrix t) -> IO p
g (IO (Matrix t) -> p) -> IO (Matrix t) -> p
forall a b. (a -> b) -> a -> b
$ do
        Matrix t
s <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
b
        (Matrix t
a Matrix t
-> Matrix t
-> Trans (Matrix t) (Trans (Matrix t) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix t
s) Trans (Matrix t) (Trans (Matrix t) (IO CInt))
Trans (Matrix t) (CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
f IO CInt -> String -> IO ()
#| String
st
        Matrix t -> IO (Matrix t)
forall (m :: * -> *) a. Monad m => a -> m a
return Matrix t
s
    | Bool
otherwise = String -> p
forall a. HasCallStack => String -> a
error (String -> p) -> String -> p
forall a b. (a -> b) -> a -> b
$ String
st String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
" of nonsquare matrix"
  where
    n1 :: Int
n1 = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
a
    n2 :: Int
n2 = Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
a
    r :: Int
r  = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
b

data UpLo = Lower | Upper

-- | Solves a triangular system of linear equations.
triSolveR :: UpLo -> Matrix Double -> Matrix Double -> Matrix Double
triSolveR :: UpLo -> Matrix Double -> Matrix Double -> Matrix Double
triSolveR UpLo
Lower Matrix Double
a Matrix Double
b = (IO (Matrix Double) -> IO (Matrix Double))
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr Double
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr Double
    -> IO CInt)
-> String
-> Matrix Double
-> Matrix Double
-> Matrix Double
forall t t p.
(Element t, Storable t) =>
(IO (Matrix t) -> IO p)
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> p
linearSolveTRAux2 IO (Matrix Double) -> IO (Matrix Double)
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dtrtrs_l String
"triSolveR" (Matrix Double -> Matrix Double
forall t. Element t => Matrix t -> Matrix t
fmat Matrix Double
a) Matrix Double
b
triSolveR UpLo
Upper Matrix Double
a Matrix Double
b = (IO (Matrix Double) -> IO (Matrix Double))
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr Double
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr Double
    -> IO CInt)
-> String
-> Matrix Double
-> Matrix Double
-> Matrix Double
forall t t p.
(Element t, Storable t) =>
(IO (Matrix t) -> IO p)
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> p
linearSolveTRAux2 IO (Matrix Double) -> IO (Matrix Double)
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dtrtrs_u String
"triSolveR" (Matrix Double -> Matrix Double
forall t. Element t => Matrix t -> Matrix t
fmat Matrix Double
a) Matrix Double
b

-- | Solves a triangular system of linear equations.
triSolveC :: UpLo -> Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
triSolveC :: UpLo
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
triSolveC UpLo
Lower Matrix (Complex Double)
a Matrix (Complex Double)
b = (IO (Matrix (Complex Double)) -> IO (Matrix (Complex Double)))
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr (Complex Double)
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr (Complex Double)
    -> IO CInt)
-> String
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
forall t t p.
(Element t, Storable t) =>
(IO (Matrix t) -> IO p)
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> p
linearSolveTRAux2 IO (Matrix (Complex Double)) -> IO (Matrix (Complex Double))
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
ztrtrs_l String
"triSolveC" (Matrix (Complex Double) -> Matrix (Complex Double)
forall t. Element t => Matrix t -> Matrix t
fmat Matrix (Complex Double)
a) Matrix (Complex Double)
b
triSolveC UpLo
Upper Matrix (Complex Double)
a Matrix (Complex Double)
b = (IO (Matrix (Complex Double)) -> IO (Matrix (Complex Double)))
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr (Complex Double)
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr (Complex Double)
    -> IO CInt)
-> String
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
forall t t p.
(Element t, Storable t) =>
(IO (Matrix t) -> IO p)
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> p
linearSolveTRAux2 IO (Matrix (Complex Double)) -> IO (Matrix (Complex Double))
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
ztrtrs_u String
"triSolveC" (Matrix (Complex Double) -> Matrix (Complex Double)
forall t. Element t => Matrix t -> Matrix t
fmat Matrix (Complex Double)
a) Matrix (Complex Double)
b

--------------------------------------------------------------------------------
foreign import ccall unsafe "triDiagSolveR_l" dgttrs  :: R :> R :> R :> R ::> Ok
foreign import ccall unsafe "triDiagSolveC_l" zgttrs  :: C :> C :> C :> C ::> Ok

linearSolveGTAux2 :: (IO (Matrix t) -> IO p)
-> (CInt
    -> Ptr t
    -> Trans
         (Vector t)
         (CInt
          -> Ptr t -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt))
-> String
-> Vector t
-> Vector t
-> Vector t
-> Matrix t
-> p
linearSolveGTAux2 IO (Matrix t) -> IO p
g CInt
-> Ptr t
-> Trans
     (Vector t)
     (CInt -> Ptr t -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
f String
st Vector t
dl Vector t
d Vector t
du Matrix t
b
    | Int
ndl  Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
nd Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1 Bool -> Bool -> Bool
&&
      Int
ndu  Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
nd Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1 Bool -> Bool -> Bool
&&
      Int
nd   Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
r = IO p -> p
forall a. IO a -> a
unsafePerformIO (IO p -> p) -> (IO (Matrix t) -> IO p) -> IO (Matrix t) -> p
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IO (Matrix t) -> IO p
g (IO (Matrix t) -> p) -> IO (Matrix t) -> p
forall a b. (a -> b) -> a -> b
$ do
        Vector t
dl' <- [Vector t] -> Vector t
forall a. [a] -> a
head ([Vector t] -> Vector t)
-> (Matrix t -> [Vector t]) -> Matrix t -> Vector t
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix t -> [Vector t]
forall t. Element t => Matrix t -> [Vector t]
toRows (Matrix t -> Vector t) -> IO (Matrix t) -> IO (Vector t)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor ([Vector t] -> Matrix t
forall t. Element t => [Vector t] -> Matrix t
fromRows [Vector t
dl])
        Vector t
du' <- [Vector t] -> Vector t
forall a. [a] -> a
head ([Vector t] -> Vector t)
-> (Matrix t -> [Vector t]) -> Matrix t -> Vector t
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix t -> [Vector t]
forall t. Element t => Matrix t -> [Vector t]
toRows (Matrix t -> Vector t) -> IO (Matrix t) -> IO (Vector t)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor ([Vector t] -> Matrix t
forall t. Element t => [Vector t] -> Matrix t
fromRows [Vector t
du])
        Matrix t
s <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
b
        (Vector t
dl' Vector t
-> (Trans
      (Vector t)
      (CInt -> Ptr t -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
    -> IO CInt)
-> Trans
     (Vector t)
     (Trans
        (Vector t)
        (CInt
         -> Ptr t -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt))
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Vector t
d Vector t
-> ((CInt
     -> Ptr t -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
    -> IO CInt)
-> Trans
     (Vector t)
     (CInt -> Ptr t -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Vector t
du' Vector t
-> Matrix t
-> Trans (Vector t) (Trans (Matrix t) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix t
s) Trans
  (Vector t)
  (Trans
     (Vector t)
     (CInt
      -> Ptr t -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt))
CInt
-> Ptr t
-> Trans
     (Vector t)
     (CInt -> Ptr t -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
f IO CInt -> String -> IO ()
#| String
st
        Matrix t -> IO (Matrix t)
forall (m :: * -> *) a. Monad m => a -> m a
return Matrix t
s
    | Bool
otherwise = String -> p
forall a. HasCallStack => String -> a
error (String -> p) -> String -> p
forall a b. (a -> b) -> a -> b
$ String
st String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
" of nonsquare matrix"
  where
    ndl :: Int
ndl  = Vector t -> Int
forall t. Storable t => Vector t -> Int
dim Vector t
dl
    nd :: Int
nd   = Vector t -> Int
forall t. Storable t => Vector t -> Int
dim Vector t
d
    ndu :: Int
ndu  = Vector t -> Int
forall t. Storable t => Vector t -> Int
dim Vector t
du
    r :: Int
r    = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
b

-- | Solves a tridiagonal system of linear equations.
triDiagSolveR :: Vector Double
-> Vector Double -> Vector Double -> Matrix Double -> Matrix Double
triDiagSolveR Vector Double
dl Vector Double
d Vector Double
du Matrix Double
b = (IO (Matrix Double) -> IO (Matrix Double))
-> (CInt
    -> Ptr Double
    -> CInt
    -> Ptr Double
    -> CInt
    -> Ptr Double
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr Double
    -> IO CInt)
-> String
-> Vector Double
-> Vector Double
-> Vector Double
-> Matrix Double
-> Matrix Double
forall t t t t p.
(Element t, Element t, Element t, Storable t) =>
(IO (Matrix t) -> IO p)
-> (CInt
    -> Ptr t
    -> CInt
    -> Ptr t
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Vector t
-> Vector t
-> Vector t
-> Matrix t
-> p
linearSolveGTAux2 IO (Matrix Double) -> IO (Matrix Double)
forall a. a -> a
id CInt
-> Ptr Double
-> CInt
-> Ptr Double
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgttrs String
"triDiagSolveR" Vector Double
dl Vector Double
d Vector Double
du Matrix Double
b
triDiagSolveC :: Vector (Complex Double)
-> Vector (Complex Double)
-> Vector (Complex Double)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
triDiagSolveC Vector (Complex Double)
dl Vector (Complex Double)
d Vector (Complex Double)
du Matrix (Complex Double)
b = (IO (Matrix (Complex Double)) -> IO (Matrix (Complex Double)))
-> (CInt
    -> Ptr (Complex Double)
    -> CInt
    -> Ptr (Complex Double)
    -> CInt
    -> Ptr (Complex Double)
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr (Complex Double)
    -> IO CInt)
-> String
-> Vector (Complex Double)
-> Vector (Complex Double)
-> Vector (Complex Double)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
forall t t t t p.
(Element t, Element t, Element t, Storable t) =>
(IO (Matrix t) -> IO p)
-> (CInt
    -> Ptr t
    -> CInt
    -> Ptr t
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Vector t
-> Vector t
-> Vector t
-> Matrix t
-> p
linearSolveGTAux2 IO (Matrix (Complex Double)) -> IO (Matrix (Complex Double))
forall a. a -> a
id CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgttrs String
"triDiagSolveC" Vector (Complex Double)
dl Vector (Complex Double)
d Vector (Complex Double)
du Matrix (Complex Double)
b

-----------------------------------------------------------------------------------

foreign import ccall unsafe "linearSolveLSR_l"   dgels ::           R ::> R ::> Ok
foreign import ccall unsafe "linearSolveLSC_l"   zgels ::           C ::> C ::> Ok
foreign import ccall unsafe "linearSolveSVDR_l" dgelss :: Double -> R ::> R ::> Ok
foreign import ccall unsafe "linearSolveSVDC_l" zgelss :: Double -> C ::> C ::> Ok

linearSolveAux :: (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> IO CInt)
-> String -> Matrix t -> Matrix t -> Matrix t
linearSolveAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt
f String
st Matrix t
a Matrix t
b
    | Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
b = IO (Matrix t) -> Matrix t
forall a. IO a -> a
unsafePerformIO (IO (Matrix t) -> Matrix t) -> IO (Matrix t) -> Matrix t
forall a b. (a -> b) -> a -> b
$ do
        Matrix t
a' <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
a
        Matrix t
r  <- MatrixOrder -> Int -> Int -> IO (Matrix t)
forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
m Int
n) Int
nrhs
        Int -> Int -> Matrix t -> Matrix t -> IO ()
forall a. Element a => Int -> Int -> Matrix a -> Matrix a -> IO ()
setRect Int
0 Int
0 Matrix t
b Matrix t
r
        (Matrix t
a' Matrix t
-> Matrix t
-> Trans (Matrix t) (Trans (Matrix t) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix t
r) Trans (Matrix t) (Trans (Matrix t) (IO CInt))
CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt
f IO CInt -> String -> IO ()
#| String
st
        Matrix t -> IO (Matrix t)
forall (m :: * -> *) a. Monad m => a -> m a
return Matrix t
r
    | Bool
otherwise = String -> Matrix t
forall a. HasCallStack => String -> a
error (String -> Matrix t) -> String -> Matrix t
forall a b. (a -> b) -> a -> b
$ String
"different number of rows in linearSolve ("String -> String -> String
forall a. [a] -> [a] -> [a]
++String
stString -> String -> String
forall a. [a] -> [a] -> [a]
++String
")"
  where
    m :: Int
m = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
a
    n :: Int
n = Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
a
    nrhs :: Int
nrhs = Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
b

-- | Least squared error solution of an overconstrained real linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /dgels/. For rank-deficient systems use 'linearSolveSVDR'.
linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
linearSolveLSR Matrix Double
a Matrix Double
b = (Int, Int) -> (Int, Int) -> Matrix Double -> Matrix Double
forall a.
Element a =>
(Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
subMatrix (Int
0,Int
0) (Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
a, Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
b) (Matrix Double -> Matrix Double) -> Matrix Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$
                     (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String -> Matrix Double -> Matrix Double -> Matrix Double
forall t t.
(Element t, Element t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> IO CInt)
-> String -> Matrix t -> Matrix t -> Matrix t
linearSolveAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgels String
"linearSolverLSR" Matrix Double
a Matrix Double
b

-- | Least squared error solution of an overconstrained complex linear system, or the minimum norm solution of an underconstrained system, using LAPACK's /zgels/. For rank-deficient systems use 'linearSolveSVDC'.
linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
linearSolveLSC :: Matrix (Complex Double)
-> Matrix (Complex Double) -> Matrix (Complex Double)
linearSolveLSC Matrix (Complex Double)
a Matrix (Complex Double)
b = (Int, Int)
-> (Int, Int) -> Matrix (Complex Double) -> Matrix (Complex Double)
forall a.
Element a =>
(Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
subMatrix (Int
0,Int
0) (Matrix (Complex Double) -> Int
forall t. Matrix t -> Int
cols Matrix (Complex Double)
a, Matrix (Complex Double) -> Int
forall t. Matrix t -> Int
cols Matrix (Complex Double)
b) (Matrix (Complex Double) -> Matrix (Complex Double))
-> Matrix (Complex Double) -> Matrix (Complex Double)
forall a b. (a -> b) -> a -> b
$
                     (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
forall t t.
(Element t, Element t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> IO CInt)
-> String -> Matrix t -> Matrix t -> Matrix t
linearSolveAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgels String
"linearSolveLSC" Matrix (Complex Double)
a Matrix (Complex Double)
b

-- | Minimum norm solution of a general real linear least squares problem Ax=B using the SVD, based on LAPACK's /dgelss/. Admits rank-deficient systems but it is slower than 'linearSolveLSR'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used.
linearSolveSVDR :: Maybe Double   -- ^ rcond
                -> Matrix Double  -- ^ coefficient matrix
                -> Matrix Double  -- ^ right hand sides (as columns)
                -> Matrix Double  -- ^ solution vectors (as columns)
linearSolveSVDR :: Maybe Double -> Matrix Double -> Matrix Double -> Matrix Double
linearSolveSVDR (Just Double
rcond) Matrix Double
a Matrix Double
b = (Int, Int) -> (Int, Int) -> Matrix Double -> Matrix Double
forall a.
Element a =>
(Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
subMatrix (Int
0,Int
0) (Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
a, Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
b) (Matrix Double -> Matrix Double) -> Matrix Double -> Matrix Double
forall a b. (a -> b) -> a -> b
$
                                   (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String -> Matrix Double -> Matrix Double -> Matrix Double
forall t t.
(Element t, Element t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> IO CInt)
-> String -> Matrix t -> Matrix t -> Matrix t
linearSolveAux (Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgelss Double
rcond) String
"linearSolveSVDR" Matrix Double
a Matrix Double
b
linearSolveSVDR Maybe Double
Nothing Matrix Double
a Matrix Double
b = Maybe Double -> Matrix Double -> Matrix Double -> Matrix Double
linearSolveSVDR (Double -> Maybe Double
forall a. a -> Maybe a
Just (-Double
1)) Matrix Double
a Matrix Double
b

-- | Minimum norm solution of a general complex linear least squares problem Ax=B using the SVD, based on LAPACK's /zgelss/. Admits rank-deficient systems but it is slower than 'linearSolveLSC'. The effective rank of A is determined by treating as zero those singular valures which are less than rcond times the largest singular value. If rcond == Nothing machine precision is used.
linearSolveSVDC :: Maybe Double            -- ^ rcond
                -> Matrix (Complex Double) -- ^ coefficient matrix
                -> Matrix (Complex Double) -- ^ right hand sides (as columns)
                -> Matrix (Complex Double) -- ^ solution vectors (as columns)
linearSolveSVDC :: Maybe Double
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
linearSolveSVDC (Just Double
rcond) Matrix (Complex Double)
a Matrix (Complex Double)
b = (Int, Int)
-> (Int, Int) -> Matrix (Complex Double) -> Matrix (Complex Double)
forall a.
Element a =>
(Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
subMatrix (Int
0,Int
0) (Matrix (Complex Double) -> Int
forall t. Matrix t -> Int
cols Matrix (Complex Double)
a, Matrix (Complex Double) -> Int
forall t. Matrix t -> Int
cols Matrix (Complex Double)
b) (Matrix (Complex Double) -> Matrix (Complex Double))
-> Matrix (Complex Double) -> Matrix (Complex Double)
forall a b. (a -> b) -> a -> b
$
                                   (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
forall t t.
(Element t, Element t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> IO CInt)
-> String -> Matrix t -> Matrix t -> Matrix t
linearSolveAux (Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgelss Double
rcond) String
"linearSolveSVDC" Matrix (Complex Double)
a Matrix (Complex Double)
b
linearSolveSVDC Maybe Double
Nothing Matrix (Complex Double)
a Matrix (Complex Double)
b = Maybe Double
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
linearSolveSVDC (Double -> Maybe Double
forall a. a -> Maybe a
Just (-Double
1)) Matrix (Complex Double)
a Matrix (Complex Double)
b

-----------------------------------------------------------------------------------

foreign import ccall unsafe "chol_l_H" zpotrf :: C ::> Ok
foreign import ccall unsafe "chol_l_S" dpotrf :: R ::> Ok

cholAux :: (CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> IO (Matrix t)
cholAux CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt
f String
st Matrix t
a = do
    Matrix t
r <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
a
    (Matrix t
r Matrix t
-> (IO CInt -> IO CInt) -> Trans (Matrix t) (IO CInt) -> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# IO CInt -> IO CInt
forall a. a -> a
id) Trans (Matrix t) (IO CInt)
CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt
f IO CInt -> String -> IO ()
#| String
st
    Matrix t -> IO (Matrix t)
forall (m :: * -> *) a. Monad m => a -> m a
return Matrix t
r

-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/.
cholH :: Matrix (Complex Double) -> Matrix (Complex Double)
cholH :: Matrix (Complex Double) -> Matrix (Complex Double)
cholH = IO (Matrix (Complex Double)) -> Matrix (Complex Double)
forall a. IO a -> a
unsafePerformIO (IO (Matrix (Complex Double)) -> Matrix (Complex Double))
-> (Matrix (Complex Double) -> IO (Matrix (Complex Double)))
-> Matrix (Complex Double)
-> Matrix (Complex Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CInt -> CInt -> CInt -> CInt -> Ptr (Complex Double) -> IO CInt)
-> String
-> Matrix (Complex Double)
-> IO (Matrix (Complex Double))
forall t.
Element t =>
(CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> IO (Matrix t)
cholAux CInt -> CInt -> CInt -> CInt -> Ptr (Complex Double) -> IO CInt
zpotrf String
"cholH"

-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/.
cholS :: Matrix Double -> Matrix Double
cholS :: Matrix Double -> Matrix Double
cholS =  IO (Matrix Double) -> Matrix Double
forall a. IO a -> a
unsafePerformIO (IO (Matrix Double) -> Matrix Double)
-> (Matrix Double -> IO (Matrix Double))
-> Matrix Double
-> Matrix Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CInt -> CInt -> CInt -> CInt -> Ptr Double -> IO CInt)
-> String -> Matrix Double -> IO (Matrix Double)
forall t.
Element t =>
(CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> IO (Matrix t)
cholAux CInt -> CInt -> CInt -> CInt -> Ptr Double -> IO CInt
dpotrf String
"cholS"

-- | Cholesky factorization of a complex Hermitian positive definite matrix, using LAPACK's /zpotrf/ ('Maybe' version).
mbCholH :: Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
mbCholH :: Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
mbCholH = IO (Maybe (Matrix (Complex Double)))
-> Maybe (Matrix (Complex Double))
forall a. IO a -> a
unsafePerformIO (IO (Maybe (Matrix (Complex Double)))
 -> Maybe (Matrix (Complex Double)))
-> (Matrix (Complex Double)
    -> IO (Maybe (Matrix (Complex Double))))
-> Matrix (Complex Double)
-> Maybe (Matrix (Complex Double))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IO (Matrix (Complex Double))
-> IO (Maybe (Matrix (Complex Double)))
forall x. IO x -> IO (Maybe x)
mbCatch (IO (Matrix (Complex Double))
 -> IO (Maybe (Matrix (Complex Double))))
-> (Matrix (Complex Double) -> IO (Matrix (Complex Double)))
-> Matrix (Complex Double)
-> IO (Maybe (Matrix (Complex Double)))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CInt -> CInt -> CInt -> CInt -> Ptr (Complex Double) -> IO CInt)
-> String
-> Matrix (Complex Double)
-> IO (Matrix (Complex Double))
forall t.
Element t =>
(CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> IO (Matrix t)
cholAux CInt -> CInt -> CInt -> CInt -> Ptr (Complex Double) -> IO CInt
zpotrf String
"cholH"

-- | Cholesky factorization of a real symmetric positive definite matrix, using LAPACK's /dpotrf/  ('Maybe' version).
mbCholS :: Matrix Double -> Maybe (Matrix Double)
mbCholS :: Matrix Double -> Maybe (Matrix Double)
mbCholS =  IO (Maybe (Matrix Double)) -> Maybe (Matrix Double)
forall a. IO a -> a
unsafePerformIO (IO (Maybe (Matrix Double)) -> Maybe (Matrix Double))
-> (Matrix Double -> IO (Maybe (Matrix Double)))
-> Matrix Double
-> Maybe (Matrix Double)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IO (Matrix Double) -> IO (Maybe (Matrix Double))
forall x. IO x -> IO (Maybe x)
mbCatch (IO (Matrix Double) -> IO (Maybe (Matrix Double)))
-> (Matrix Double -> IO (Matrix Double))
-> Matrix Double
-> IO (Maybe (Matrix Double))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CInt -> CInt -> CInt -> CInt -> Ptr Double -> IO CInt)
-> String -> Matrix Double -> IO (Matrix Double)
forall t.
Element t =>
(CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> IO (Matrix t)
cholAux CInt -> CInt -> CInt -> CInt -> Ptr Double -> IO CInt
dpotrf String
"cholS"

-----------------------------------------------------------------------------------

type TMVM t = t ::> t :> t ::> Ok

foreign import ccall unsafe "qr_l_R" dgeqr2 :: R :> R ::> Ok
foreign import ccall unsafe "qr_l_C" zgeqr2 :: C :> C ::> Ok

-- | QR factorization of a real matrix, using LAPACK's /dgeqr2/.
qrR :: Matrix Double -> (Matrix Double, Vector Double)
qrR :: Matrix Double -> (Matrix Double, Vector Double)
qrR = (CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String -> Matrix Double -> (Matrix Double, Vector Double)
forall t a.
(Element t, Storable a) =>
(CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Matrix t, Vector a)
qrAux CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgeqr2 String
"qrR"

-- | QR factorization of a complex matrix, using LAPACK's /zgeqr2/.
qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
qrC :: Matrix (Complex Double)
-> (Matrix (Complex Double), Vector (Complex Double))
qrC = (CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> (Matrix (Complex Double), Vector (Complex Double))
forall t a.
(Element t, Storable a) =>
(CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Matrix t, Vector a)
qrAux CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgeqr2 String
"qrC"

qrAux :: (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Matrix t, Vector a)
qrAux CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt
f String
st Matrix t
a = IO (Matrix t, Vector a) -> (Matrix t, Vector a)
forall a. IO a -> a
unsafePerformIO (IO (Matrix t, Vector a) -> (Matrix t, Vector a))
-> IO (Matrix t, Vector a) -> (Matrix t, Vector a)
forall a b. (a -> b) -> a -> b
$ do
    Matrix t
r <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
a
    Vector a
tau <- Int -> IO (Vector a)
forall a. Storable a => Int -> IO (Vector a)
createVector Int
mn
    (Vector a
tau Vector a
-> Matrix t
-> Trans (Vector a) (Trans (Matrix t) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix t
r) Trans (Vector a) (Trans (Matrix t) (IO CInt))
CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt
f IO CInt -> String -> IO ()
#| String
st
    (Matrix t, Vector a) -> IO (Matrix t, Vector a)
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix t
r,Vector a
tau)
  where
    m :: Int
m = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
a
    n :: Int
n = Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
a
    mn :: Int
mn = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
m Int
n

foreign import ccall unsafe "c_dorgqr" dorgqr :: R :> R ::> Ok
foreign import ccall unsafe "c_zungqr" zungqr :: C :> C ::> Ok

-- | build rotation from reflectors
qrgrR :: Int -> (Matrix Double, Vector Double) -> Matrix Double
qrgrR :: Int -> (Matrix Double, Vector Double) -> Matrix Double
qrgrR = (CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String -> Int -> (Matrix Double, Vector Double) -> Matrix Double
forall t t.
(Element t, Element t, Num t) =>
(CInt -> Ptr t -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Int -> (Matrix t, Vector t) -> Matrix t
qrgrAux CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dorgqr String
"qrgrR"
-- | build rotation from reflectors
qrgrC :: Int -> (Matrix (Complex Double), Vector (Complex Double)) -> Matrix (Complex Double)
qrgrC :: Int
-> (Matrix (Complex Double), Vector (Complex Double))
-> Matrix (Complex Double)
qrgrC = (CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Int
-> (Matrix (Complex Double), Vector (Complex Double))
-> Matrix (Complex Double)
forall t t.
(Element t, Element t, Num t) =>
(CInt -> Ptr t -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Int -> (Matrix t, Vector t) -> Matrix t
qrgrAux CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zungqr String
"qrgrC"

qrgrAux :: (CInt -> Ptr t -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Int -> (Matrix t, Vector t) -> Matrix t
qrgrAux CInt -> Ptr t -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt
f String
st Int
n (Matrix t
a, Vector t
tau) = IO (Matrix t) -> Matrix t
forall a. IO a -> a
unsafePerformIO (IO (Matrix t) -> Matrix t) -> IO (Matrix t) -> Matrix t
forall a b. (a -> b) -> a -> b
$ do
    Matrix t
res <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor ((Int, Int) -> (Int, Int) -> Matrix t -> Matrix t
forall a.
Element a =>
(Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
subMatrix (Int
0,Int
0) (Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
a,Int
n) Matrix t
a)
    ((Int -> Int -> Vector t -> Vector t
forall t. Storable t => Int -> Int -> Vector t -> Vector t
subVector Int
0 Int
n Vector t
tau') Vector t
-> Matrix t
-> Trans (Vector t) (Trans (Matrix t) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix t
res) Trans (Vector t) (Trans (Matrix t) (IO CInt))
CInt -> Ptr t -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt
f IO CInt -> String -> IO ()
#| String
st
    Matrix t -> IO (Matrix t)
forall (m :: * -> *) a. Monad m => a -> m a
return Matrix t
res
  where
    tau' :: Vector t
tau' = [Vector t] -> Vector t
forall t. Storable t => [Vector t] -> Vector t
vjoin [Vector t
tau, t -> Int -> Vector t
forall a. Element a => a -> Int -> Vector a
constantD t
0 Int
n]

-----------------------------------------------------------------------------------
foreign import ccall unsafe "hess_l_R" dgehrd :: R :> R ::> Ok
foreign import ccall unsafe "hess_l_C" zgehrd :: C :> C ::> Ok

-- | Hessenberg factorization of a square real matrix, using LAPACK's /dgehrd/.
hessR :: Matrix Double -> (Matrix Double, Vector Double)
hessR :: Matrix Double -> (Matrix Double, Vector Double)
hessR = (CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String -> Matrix Double -> (Matrix Double, Vector Double)
forall t a.
(Element t, Storable a) =>
(CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Matrix t, Vector a)
hessAux CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgehrd String
"hessR"

-- | Hessenberg factorization of a square complex matrix, using LAPACK's /zgehrd/.
hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
hessC :: Matrix (Complex Double)
-> (Matrix (Complex Double), Vector (Complex Double))
hessC = (CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> (Matrix (Complex Double), Vector (Complex Double))
forall t a.
(Element t, Storable a) =>
(CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Matrix t, Vector a)
hessAux CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgehrd String
"hessC"

hessAux :: (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Matrix t, Vector a)
hessAux CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt
f String
st Matrix t
a = IO (Matrix t, Vector a) -> (Matrix t, Vector a)
forall a. IO a -> a
unsafePerformIO (IO (Matrix t, Vector a) -> (Matrix t, Vector a))
-> IO (Matrix t, Vector a) -> (Matrix t, Vector a)
forall a b. (a -> b) -> a -> b
$ do
    Matrix t
r <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
a
    Vector a
tau <- Int -> IO (Vector a)
forall a. Storable a => Int -> IO (Vector a)
createVector (Int
mnInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
    (Vector a
tau Vector a
-> Matrix t
-> Trans (Vector a) (Trans (Matrix t) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix t
r) Trans (Vector a) (Trans (Matrix t) (IO CInt))
CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt
f IO CInt -> String -> IO ()
#| String
st
    (Matrix t, Vector a) -> IO (Matrix t, Vector a)
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix t
r,Vector a
tau)
  where
    m :: Int
m = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
a
    n :: Int
n = Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
a
    mn :: Int
mn = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
m Int
n

-----------------------------------------------------------------------------------
foreign import ccall unsafe "schur_l_R" dgees :: R ::> R ::> Ok
foreign import ccall unsafe "schur_l_C" zgees :: C ::> C ::> Ok

-- | Schur factorization of a square real matrix, using LAPACK's /dgees/.
schurR :: Matrix Double -> (Matrix Double, Matrix Double)
schurR :: Matrix Double -> (Matrix Double, Matrix Double)
schurR = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String -> Matrix Double -> (Matrix Double, Matrix Double)
forall t a.
(Element t, Storable a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Matrix t)
schurAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgees String
"schurR"

-- | Schur factorization of a square complex matrix, using LAPACK's /zgees/.
schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double))
schurC :: Matrix (Complex Double)
-> (Matrix (Complex Double), Matrix (Complex Double))
schurC = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> (Matrix (Complex Double), Matrix (Complex Double))
forall t a.
(Element t, Storable a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Matrix t)
schurAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgees String
"schurC"

schurAux :: (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Matrix t)
schurAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt
f String
st Matrix t
a = IO (Matrix a, Matrix t) -> (Matrix a, Matrix t)
forall a. IO a -> a
unsafePerformIO (IO (Matrix a, Matrix t) -> (Matrix a, Matrix t))
-> IO (Matrix a, Matrix t) -> (Matrix a, Matrix t)
forall a b. (a -> b) -> a -> b
$ do
    Matrix a
u <- MatrixOrder -> Int -> Int -> IO (Matrix a)
forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
n Int
n
    Matrix t
s <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
a
    (Matrix a
u Matrix a
-> Matrix t
-> Trans (Matrix a) (Trans (Matrix t) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix t
s) Trans (Matrix a) (Trans (Matrix t) (IO CInt))
CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt
f IO CInt -> String -> IO ()
#| String
st
    (Matrix a, Matrix t) -> IO (Matrix a, Matrix t)
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix a
u,Matrix t
s)
  where
    n :: Int
n = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
a

-----------------------------------------------------------------------------------
foreign import ccall unsafe "lu_l_R" dgetrf :: R :> R ::> Ok
foreign import ccall unsafe "lu_l_C" zgetrf :: R :> C ::> Ok

-- | LU factorization of a general real matrix, using LAPACK's /dgetrf/.
luR :: Matrix Double -> (Matrix Double, [Int])
luR :: Matrix Double -> (Matrix Double, [Int])
luR = (CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String -> Matrix Double -> (Matrix Double, [Int])
forall t a b.
(Element t, Storable a, RealFrac a, Integral b) =>
(CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Matrix t, [b])
luAux CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgetrf String
"luR"

-- | LU factorization of a general complex matrix, using LAPACK's /zgetrf/.
luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
luC = (CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> (Matrix (Complex Double), [Int])
forall t a b.
(Element t, Storable a, RealFrac a, Integral b) =>
(CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Matrix t, [b])
luAux CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgetrf String
"luC"

luAux :: (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Matrix t, [b])
luAux CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt
f String
st Matrix t
a = IO (Matrix t, [b]) -> (Matrix t, [b])
forall a. IO a -> a
unsafePerformIO (IO (Matrix t, [b]) -> (Matrix t, [b]))
-> IO (Matrix t, [b]) -> (Matrix t, [b])
forall a b. (a -> b) -> a -> b
$ do
    Matrix t
lu <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
a
    Vector a
piv <- Int -> IO (Vector a)
forall a. Storable a => Int -> IO (Vector a)
createVector (Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
n Int
m)
    (Vector a
piv Vector a
-> Matrix t
-> Trans (Vector a) (Trans (Matrix t) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix t
lu) Trans (Vector a) (Trans (Matrix t) (IO CInt))
CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt
f IO CInt -> String -> IO ()
#| String
st
    (Matrix t, [b]) -> IO (Matrix t, [b])
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix t
lu, (a -> b) -> [a] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map (b -> b
forall a. Enum a => a -> a
pred(b -> b) -> (a -> b) -> a -> b
forall b c a. (b -> c) -> (a -> b) -> a -> c
.a -> b
forall a b. (RealFrac a, Integral b) => a -> b
round) (Vector a -> [a]
forall a. Storable a => Vector a -> [a]
toList Vector a
piv))
  where
    n :: Int
n = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
a
    m :: Int
m = Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
a

-----------------------------------------------------------------------------------

foreign import ccall unsafe "luS_l_R" dgetrs :: R ::> R :> R ::> Ok
foreign import ccall unsafe "luS_l_C" zgetrs :: C ::> R :> C ::> Ok

-- | Solve a real linear system from a precomputed LU decomposition ('luR'), using LAPACK's /dgetrs/.
lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double
lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double
lusR Matrix Double
a [Int]
piv Matrix Double
b = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String
-> Matrix Double
-> [Int]
-> Matrix Double
-> Matrix Double
forall t t a.
(Element t, Storable t, Integral a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> IO CInt)
-> String -> Matrix t -> [a] -> Matrix t -> Matrix t
lusAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dgetrs String
"lusR" (Matrix Double -> Matrix Double
forall t. Element t => Matrix t -> Matrix t
fmat Matrix Double
a) [Int]
piv Matrix Double
b

-- | Solve a complex linear system from a precomputed LU decomposition ('luC'), using LAPACK's /zgetrs/.
lusC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
lusC :: Matrix (Complex Double)
-> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
lusC Matrix (Complex Double)
a [Int]
piv Matrix (Complex Double)
b = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> [Int]
-> Matrix (Complex Double)
-> Matrix (Complex Double)
forall t t a.
(Element t, Storable t, Integral a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> IO CInt)
-> String -> Matrix t -> [a] -> Matrix t -> Matrix t
lusAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zgetrs String
"lusC" (Matrix (Complex Double) -> Matrix (Complex Double)
forall t. Element t => Matrix t -> Matrix t
fmat Matrix (Complex Double)
a) [Int]
piv Matrix (Complex Double)
b

lusAux :: Trans
  (Matrix t)
  (CInt
   -> Ptr Double -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> [a] -> Matrix t -> Matrix t
lusAux Trans
  (Matrix t)
  (CInt
   -> Ptr Double -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
f String
st Matrix t
a [a]
piv Matrix t
b
    | Int
n1Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
n2 Bool -> Bool -> Bool
&& Int
n2Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
n =IO (Matrix t) -> Matrix t
forall a. IO a -> a
unsafePerformIO (IO (Matrix t) -> Matrix t) -> IO (Matrix t) -> Matrix t
forall a b. (a -> b) -> a -> b
$ do
         Matrix t
x <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
b
         (Matrix t
a Matrix t
-> ((CInt
     -> Ptr Double -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
    -> IO CInt)
-> Trans
     (Matrix t)
     (CInt
      -> Ptr Double -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> IO CInt
forall c b r. TransArray c => c -> (b -> IO r) -> Trans c b -> IO r
# Vector Double
piv' Vector Double
-> Matrix t
-> Trans (Vector Double) (Trans (Matrix t) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix t
x) Trans
  (Matrix t)
  (CInt
   -> Ptr Double -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
f IO CInt -> String -> IO ()
#| String
st
         Matrix t -> IO (Matrix t)
forall (m :: * -> *) a. Monad m => a -> m a
return Matrix t
x
    | Bool
otherwise = String -> Matrix t
forall a. HasCallStack => String -> a
error String
st
  where
    n1 :: Int
n1 = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
a
    n2 :: Int
n2 = Matrix t -> Int
forall t. Matrix t -> Int
cols Matrix t
a
    n :: Int
n = Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
b
    piv' :: Vector Double
piv' = [Double] -> Vector Double
forall a. Storable a => [a] -> Vector a
fromList ((a -> Double) -> [a] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map (a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral(a -> Double) -> (a -> a) -> a -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
.a -> a
forall a. Enum a => a -> a
succ) [a]
piv) :: Vector Double

-----------------------------------------------------------------------------------
foreign import ccall unsafe "ldl_R" dsytrf :: R :> R ::> Ok
foreign import ccall unsafe "ldl_C" zhetrf :: R :> C ::> Ok

-- | LDL factorization of a symmetric real matrix, using LAPACK's /dsytrf/.
ldlR :: Matrix Double -> (Matrix Double, [Int])
ldlR :: Matrix Double -> (Matrix Double, [Int])
ldlR = (CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String -> Matrix Double -> (Matrix Double, [Int])
forall t a b.
(Element t, Storable a, RealFrac a, Integral b) =>
(CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Matrix t, [b])
ldlAux CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dsytrf String
"ldlR"

-- | LDL factorization of a hermitian complex matrix, using LAPACK's /zhetrf/.
ldlC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
ldlC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
ldlC = (CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> (Matrix (Complex Double), [Int])
forall t a b.
(Element t, Storable a, RealFrac a, Integral b) =>
(CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Matrix t, [b])
ldlAux CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zhetrf String
"ldlC"

ldlAux :: (CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Matrix t, [b])
ldlAux CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt
f String
st Matrix t
a = IO (Matrix t, [b]) -> (Matrix t, [b])
forall a. IO a -> a
unsafePerformIO (IO (Matrix t, [b]) -> (Matrix t, [b]))
-> IO (Matrix t, [b]) -> (Matrix t, [b])
forall a b. (a -> b) -> a -> b
$ do
    Matrix t
ldl <- MatrixOrder -> Matrix t -> IO (Matrix t)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
a
    Vector a
piv <- Int -> IO (Vector a)
forall a. Storable a => Int -> IO (Vector a)
createVector (Matrix t -> Int
forall t. Matrix t -> Int
rows Matrix t
a)
    (Vector a
piv Vector a
-> Matrix t
-> Trans (Vector a) (Trans (Matrix t) (IO CInt))
-> IO CInt
forall c c r.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix t
ldl) Trans (Vector a) (Trans (Matrix t) (IO CInt))
CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt
f IO CInt -> String -> IO ()
#| String
st
    (Matrix t, [b]) -> IO (Matrix t, [b])
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix t
ldl, (a -> b) -> [a] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map (b -> b
forall a. Enum a => a -> a
pred(b -> b) -> (a -> b) -> a -> b
forall b c a. (b -> c) -> (a -> b) -> a -> c
.a -> b
forall a b. (RealFrac a, Integral b) => a -> b
round) (Vector a -> [a]
forall a. Storable a => Vector a -> [a]
toList Vector a
piv))

-----------------------------------------------------------------------------------

foreign import ccall unsafe "ldl_S_R" dsytrs :: R ::> R :> R ::> Ok
foreign import ccall unsafe "ldl_S_C" zsytrs :: C ::> R :> C ::> Ok

-- | Solve a real linear system from a precomputed LDL decomposition ('ldlR'), using LAPACK's /dsytrs/.
ldlsR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double
ldlsR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double
ldlsR Matrix Double
a [Int]
piv Matrix Double
b = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr Double
 -> IO CInt)
-> String
-> Matrix Double
-> [Int]
-> Matrix Double
-> Matrix Double
forall t t a.
(Element t, Storable t, Integral a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> IO CInt)
-> String -> Matrix t -> [a] -> Matrix t -> Matrix t
lusAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Double
-> IO CInt
dsytrs String
"ldlsR" (Matrix Double -> Matrix Double
forall t. Element t => Matrix t -> Matrix t
fmat Matrix Double
a) [Int]
piv Matrix Double
b

-- | Solve a complex linear system from a precomputed LDL decomposition ('ldlC'), using LAPACK's /zsytrs/.
ldlsC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
ldlsC :: Matrix (Complex Double)
-> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
ldlsC Matrix (Complex Double)
a [Int]
piv Matrix (Complex Double)
b = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr (Complex Double)
 -> IO CInt)
-> String
-> Matrix (Complex Double)
-> [Int]
-> Matrix (Complex Double)
-> Matrix (Complex Double)
forall t t a.
(Element t, Storable t, Integral a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> Ptr Double
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> IO CInt)
-> String -> Matrix t -> [a] -> Matrix t -> Matrix t
lusAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> CInt
-> Ptr Double
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr (Complex Double)
-> IO CInt
zsytrs String
"ldlsC" (Matrix (Complex Double) -> Matrix (Complex Double)
forall t. Element t => Matrix t -> Matrix t
fmat Matrix (Complex Double)
a) [Int]
piv Matrix (Complex Double)
b