{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE ViewPatterns #-}
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
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
multiplyR :: Matrix Double -> Matrix Double -> Matrix Double
multiplyR :: Matrix R -> Matrix R -> Matrix R
multiplyR Matrix R
a Matrix R
b = {-# SCC "multiplyR" #-} (CInt
-> CInt
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String -> Matrix R -> Matrix R -> Matrix R
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 R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgemmc String
"dgemmc" Matrix R
a Matrix R
b
multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
multiplyC :: Matrix C -> Matrix C -> Matrix C
multiplyC Matrix C
a Matrix C
b = (CInt
-> CInt
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String -> Matrix C -> Matrix C -> Matrix C
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 C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgemmc String
"zgemmc" Matrix C
a Matrix C
b
multiplyF :: Matrix Float -> Matrix Float -> Matrix Float
multiplyF :: Matrix F -> Matrix F -> Matrix F
multiplyF Matrix F
a Matrix F
b = (CInt
-> CInt
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr F
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr F
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr F
-> IO CInt)
-> String -> Matrix F -> Matrix F -> Matrix F
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 F
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr F
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr F
-> IO CInt
sgemmc String
"sgemmc" Matrix F
a Matrix F
b
multiplyQ :: Matrix (Complex Float) -> Matrix (Complex Float) -> Matrix (Complex Float)
multiplyQ :: Matrix Q -> Matrix Q -> Matrix Q
multiplyQ Matrix Q
a Matrix Q
b = (CInt
-> CInt
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Q
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Q
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Q
-> IO CInt)
-> String -> Matrix Q -> Matrix Q -> Matrix Q
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 Q
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Q
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr Q
-> IO CInt
cgemmc String
"cgemmc" Matrix Q
a Matrix Q
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 ::> IO CInt)) -> IO CInt)
-> Trans (Matrix CInt) (CInt ::> (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
-> ((Z ::> (Z ::> IO CInt)) -> IO CInt)
-> Trans (Matrix Z) (Z ::> (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
svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
svdR :: Matrix R -> (Matrix R, Vector R, Matrix R)
svdR = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String -> Matrix R -> (Matrix R, Vector R, Matrix R)
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 R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgesvd String
"svdR"
svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
svdRd :: Matrix R -> (Matrix R, Vector R, Matrix R)
svdRd = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String -> Matrix R -> (Matrix R, Vector R, Matrix R)
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 R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgesdd String
"svdRdd"
svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
svdC :: Matrix C -> (Matrix C, Vector R, Matrix C)
svdC = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String -> Matrix C -> (Matrix C, Vector R, Matrix C)
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 C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgesvd String
"svdC"
svdCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
svdCd :: Matrix C -> (Matrix C, Vector R, Matrix C)
svdCd = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String -> Matrix C -> (Matrix C, Vector R, Matrix C)
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 C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> 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
thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
thinSVDR :: Matrix R -> (Matrix R, Vector R, Matrix R)
thinSVDR = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String -> Matrix R -> (Matrix R, Vector R, Matrix R)
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 R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgesvd String
"thinSVDR"
thinSVDC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
thinSVDC :: Matrix C -> (Matrix C, Vector R, Matrix C)
thinSVDC = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String -> Matrix C -> (Matrix C, Vector R, Matrix C)
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 C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgesvd String
"thinSVDC"
thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
thinSVDRd :: Matrix R -> (Matrix R, Vector R, Matrix R)
thinSVDRd = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String -> Matrix R -> (Matrix R, Vector R, Matrix R)
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 R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgesdd String
"thinSVDRdd"
thinSVDCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
thinSVDCd :: Matrix C -> (Matrix C, Vector R, Matrix C)
thinSVDCd = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String -> Matrix C -> (Matrix C, Vector R, Matrix C)
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 C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> 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
svR :: Matrix Double -> Vector Double
svR :: Matrix R -> Vector R
svR = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String -> Matrix R -> Vector R
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 R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgesvd String
"svR"
svC :: Matrix (Complex Double) -> Vector Double
svC :: Matrix C -> Vector R
svC = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String -> Matrix C -> Vector R
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 C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgesvd String
"svC"
svRd :: Matrix Double -> Vector Double
svRd :: Matrix R -> Vector R
svRd = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String -> Matrix R -> Vector R
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 R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgesdd String
"svRd"
svCd :: Matrix (Complex Double) -> Vector Double
svCd :: Matrix C -> Vector R
svCd = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String -> Matrix C -> Vector R
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 C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> 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
rightSVR :: Matrix Double -> (Vector Double, Matrix Double)
rightSVR :: Matrix R -> (Vector R, Matrix R)
rightSVR = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String -> Matrix R -> (Vector R, Matrix R)
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 R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgesvd String
"rightSVR"
rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
rightSVC :: Matrix C -> (Vector R, Matrix C)
rightSVC = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String -> Matrix C -> (Vector R, Matrix C)
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 C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> 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
leftSVR :: Matrix Double -> (Matrix Double, Vector Double)
leftSVR :: Matrix R -> (Matrix R, Vector R)
leftSVR = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String -> Matrix R -> (Matrix R, Vector R)
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 R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgesvd String
"leftSVR"
leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double)
leftSVC :: Matrix C -> (Matrix C, Vector R)
leftSVC = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String -> Matrix C -> (Matrix C, Vector R)
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 C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> 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
eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double))
eigC :: Matrix C -> (Vector C, Matrix C)
eigC = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String -> Matrix C -> (Vector C, Matrix C)
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 C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> 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
eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double)
eigOnlyC :: Matrix C -> Vector C
eigOnlyC = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String -> Matrix C -> Vector C
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 C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgeev String
"eigOnlyC"
eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double))
eigR :: Matrix R -> (Vector C, Matrix C)
eigR Matrix R
m = (Vector C
s', Matrix C
v'')
where (Vector C
s,Matrix R
v) = Matrix R -> (Vector C, Matrix R)
eigRaux Matrix R
m
s' :: Vector C
s' = Vector C -> Vector C
forall {t}.
RealElement t =>
Vector (Complex t) -> Vector (Complex t)
fixeig1 Vector C
s
v' :: [Vector R]
v' = Matrix R -> [Vector R]
forall t. Element t => Matrix t -> [Vector t]
toRows (Matrix R -> [Vector R]) -> Matrix R -> [Vector R]
forall a b. (a -> b) -> a -> b
$ Matrix R -> Matrix R
forall t. Matrix t -> Matrix t
trans Matrix R
v
v'' :: Matrix C
v'' = [Vector C] -> Matrix C
forall t. Element t => [Vector t] -> Matrix t
fromColumns ([Vector C] -> Matrix C) -> [Vector C] -> Matrix C
forall a b. (a -> b) -> a -> b
$ [C] -> [Vector R] -> [Vector C]
forall {e} {a}.
(RealElement e, Num a, Eq a) =>
[Complex a] -> [Vector e] -> [Vector (Complex e)]
fixeig (Vector C -> [C]
forall a. Storable a => Vector a -> [a]
toList Vector C
s') [Vector R]
v'
eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)
eigRaux :: Matrix R -> (Vector C, Matrix R)
eigRaux Matrix R
m = IO (Vector C, Matrix R) -> (Vector C, Matrix R)
forall a. IO a -> a
unsafePerformIO (IO (Vector C, Matrix R) -> (Vector C, Matrix R))
-> IO (Vector C, Matrix R) -> (Vector C, Matrix R)
forall a b. (a -> b) -> a -> b
$ do
Matrix R
a <- MatrixOrder -> Matrix R -> IO (Matrix R)
forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix R
m
Vector C
l <- Int -> IO (Vector C)
forall a. Storable a => Int -> IO (Vector a)
createVector Int
r
Matrix R
v <- MatrixOrder -> Int -> Int -> IO (Matrix R)
forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
r Int
r
(Matrix R
a Matrix R
-> ((C :> (CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt))
-> IO CInt)
-> Trans
(Matrix R)
(C :> (CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt))
-> IO CInt
forall {c} {b} {r}.
TransArray c =>
c -> (b -> IO r) -> Trans c b -> IO r
# Vector C
l Vector C
-> Matrix R
-> Trans (Vector C) (Trans (Matrix R) (IO CInt))
-> IO CInt
forall {c} {c} {r}.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix R
v) Trans
(Matrix R)
(C :> (CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt))
CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> C :> (CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt)
g IO CInt -> String -> IO ()
#| String
"eigR"
(Vector C, Matrix R) -> IO (Vector C, Matrix R)
forall (m :: * -> *) a. Monad m => a -> m a
return (Vector C
l,Matrix R
v)
where
r :: Int
r = Matrix R -> Int
forall t. Matrix t -> Int
rows Matrix R
m
g :: CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> C :> (CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt)
g CInt
ra CInt
ca CInt
xra CInt
xca Ptr R
pa = R
::> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> C :> (CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt))
dgeev CInt
ra CInt
ca CInt
xra CInt
xca Ptr R
pa CInt
0 CInt
0 CInt
0 CInt
0 Ptr R
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"
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"
eigOnlyR :: Matrix Double -> Vector (Complex Double)
eigOnlyR :: Matrix R -> Vector C
eigOnlyR = Vector C -> Vector C
forall {t}.
RealElement t =>
Vector (Complex t) -> Vector (Complex t)
fixeig1 (Vector C -> Vector C)
-> (Matrix R -> Vector C) -> Matrix R -> Vector C
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (R
::> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> C :> (CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt)))
-> String -> Matrix R -> Vector C
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 R
::> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> C :> (CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt))
dgeev String
"eigOnlyR"
eigG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double, Matrix (Complex Double))
eigG :: Matrix R -> Matrix R -> (Vector C, Vector R, Matrix C)
eigG Matrix R
a Matrix R
b = (Vector C
alpha', Vector R
beta, Matrix C
v'')
where
(Vector C
alpha, Vector R
beta, Matrix R
v) = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> Matrix R -> Matrix R -> String -> (Vector C, Vector R, Matrix R)
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 R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dggev Matrix R
a Matrix R
b String
"eigG"
alpha' :: Vector C
alpha' = Vector C -> Vector C
forall {t}.
RealElement t =>
Vector (Complex t) -> Vector (Complex t)
fixeig1 Vector C
alpha
v' :: [Vector R]
v' = Matrix R -> [Vector R]
forall t. Element t => Matrix t -> [Vector t]
toRows (Matrix R -> [Vector R]) -> Matrix R -> [Vector R]
forall a b. (a -> b) -> a -> b
$ Matrix R -> Matrix R
forall t. Matrix t -> Matrix t
trans Matrix R
v
v'' :: Matrix C
v'' = [Vector C] -> Matrix C
forall t. Element t => [Vector t] -> Matrix t
fromColumns ([Vector C] -> Matrix C) -> [Vector C] -> Matrix C
forall a b. (a -> b) -> a -> b
$ [C] -> [Vector R] -> [Vector C]
forall {e} {a}.
(RealElement e, Ord a, Fractional a) =>
[Complex a] -> [Vector e] -> [Vector (Complex e)]
fixeigG (Vector C -> [C]
forall a. Storable a => Vector a -> [a]
toList Vector C
alpha') [Vector R]
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
eigGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double), Matrix (Complex Double))
eigGC :: Matrix C -> Matrix C -> (Vector C, Vector C, Matrix C)
eigGC Matrix C
a Matrix C
b = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> Matrix C -> Matrix C -> String -> (Vector C, Vector C, Matrix C)
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 C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zggev Matrix C
a Matrix C
b String
"eigGC"
eigOnlyG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double)
eigOnlyG :: Matrix R -> Matrix R -> (Vector C, Vector R)
eigOnlyG Matrix R
a Matrix R
b = (Vector C -> Vector C)
-> (Vector C, Vector R) -> (Vector C, Vector R)
forall (p :: * -> * -> *) a b c.
Bifunctor p =>
(a -> b) -> p a c -> p b c
first Vector C -> Vector C
forall {t}.
RealElement t =>
Vector (Complex t) -> Vector (Complex t)
fixeig1 ((Vector C, Vector R) -> (Vector C, Vector R))
-> (Vector C, Vector R) -> (Vector C, Vector R)
forall a b. (a -> b) -> a -> b
$ (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> Matrix R -> Matrix R -> String -> (Vector C, Vector R)
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 R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dggev Matrix R
a Matrix R
b String
"eigOnlyG"
eigOnlyGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double))
eigOnlyGC :: Matrix C -> Matrix C -> (Vector C, Vector C)
eigOnlyGC Matrix C
a Matrix C
b = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> Matrix C -> Matrix C -> String -> (Vector C, Vector C)
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 C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zggev Matrix C
a Matrix C
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
eigS :: Matrix Double -> (Vector Double, Matrix Double)
eigS :: Matrix R -> (Vector R, Matrix R)
eigS Matrix R
m = (Vector R
s', Matrix R -> Matrix R
forall t. Element t => Matrix t -> Matrix t
fliprl Matrix R
v)
where (Vector R
s,Matrix R
v) = Matrix R -> (Vector R, Matrix R)
eigS' Matrix R
m
s' :: Vector R
s' = [R] -> Vector R
forall a. Storable a => [a] -> Vector a
fromList ([R] -> Vector R) -> (Vector R -> [R]) -> Vector R -> Vector R
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [R] -> [R]
forall a. [a] -> [a]
reverse ([R] -> [R]) -> (Vector R -> [R]) -> Vector R -> [R]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector R -> [R]
forall a. Storable a => Vector a -> [a]
toList (Vector R -> Vector R) -> Vector R -> Vector R
forall a b. (a -> b) -> a -> b
$ Vector R
s
eigS' :: Matrix Double -> (Vector Double, Matrix Double)
eigS' :: Matrix R -> (Vector R, Matrix R)
eigS' = (CInt -> Ptr R -> CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt)
-> String -> Matrix R -> (Vector R, Matrix R)
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 R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dsyev CInt
1) String
"eigS'"
eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
eigH :: Matrix C -> (Vector R, Matrix C)
eigH Matrix C
m = (Vector R
s', Matrix C -> Matrix C
forall t. Element t => Matrix t -> Matrix t
fliprl Matrix C
v)
where
(Vector R
s,Matrix C
v) = Matrix C -> (Vector R, Matrix C)
eigH' Matrix C
m
s' :: Vector R
s' = [R] -> Vector R
forall a. Storable a => [a] -> Vector a
fromList ([R] -> Vector R) -> (Vector R -> [R]) -> Vector R -> Vector R
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [R] -> [R]
forall a. [a] -> [a]
reverse ([R] -> [R]) -> (Vector R -> [R]) -> Vector R -> [R]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector R -> [R]
forall a. Storable a => Vector a -> [a]
toList (Vector R -> Vector R) -> Vector R -> Vector R
forall a b. (a -> b) -> a -> b
$ Vector R
s
eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
eigH' :: Matrix C -> (Vector R, Matrix C)
eigH' = (CInt -> Ptr R -> CInt -> CInt -> CInt -> CInt -> Ptr C -> IO CInt)
-> String -> Matrix C -> (Vector R, Matrix C)
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 R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zheev CInt
1) String
"eigH'"
eigOnlyS :: Matrix Double -> Vector Double
eigOnlyS :: Matrix R -> Vector R
eigOnlyS = Vector R -> Vector R
vrev (Vector R -> Vector R)
-> (Matrix R -> Vector R) -> Matrix R -> Vector R
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vector R, Matrix R) -> Vector R
forall a b. (a, b) -> a
fst((Vector R, Matrix R) -> Vector R)
-> (Matrix R -> (Vector R, Matrix R)) -> Matrix R -> Vector R
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CInt -> Ptr R -> CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt)
-> String -> Matrix R -> (Vector R, Matrix R)
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 R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dsyev CInt
0) String
"eigS'"
eigOnlyH :: Matrix (Complex Double) -> Vector Double
eigOnlyH :: Matrix C -> Vector R
eigOnlyH = Vector R -> Vector R
vrev (Vector R -> Vector R)
-> (Matrix C -> Vector R) -> Matrix C -> Vector R
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vector R, Matrix C) -> Vector R
forall a b. (a, b) -> a
fst((Vector R, Matrix C) -> Vector R)
-> (Matrix C -> (Vector R, Matrix C)) -> Matrix C -> Vector R
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CInt -> Ptr R -> CInt -> CInt -> CInt -> CInt -> Ptr C -> IO CInt)
-> String -> Matrix C -> (Vector R, Matrix C)
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 R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zheev CInt
0) String
"eigH'"
vrev :: Vector R -> Vector R
vrev = Matrix R -> Vector R
forall t. Element t => Matrix t -> Vector t
flatten (Matrix R -> Vector R)
-> (Vector R -> Matrix R) -> Vector R -> Vector R
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Matrix R -> Matrix R
forall t. Element t => Matrix t -> Matrix t
flipud (Matrix R -> Matrix R)
-> (Vector R -> Matrix R) -> Vector R -> Matrix R
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Vector R -> Matrix R
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 c)
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveSQAux IO (Matrix t) -> IO c
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 c -> c
forall a. IO a -> a
unsafePerformIO (IO c -> c) -> (IO (Matrix t) -> IO c) -> IO (Matrix t) -> c
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IO (Matrix t) -> IO c
g (IO (Matrix t) -> c) -> IO (Matrix t) -> c
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 -> c
forall a. HasCallStack => String -> a
error (String -> c) -> String -> c
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
linearSolveR :: Matrix Double -> Matrix Double -> Matrix Double
linearSolveR :: Matrix R -> Matrix R -> Matrix R
linearSolveR Matrix R
a Matrix R
b = (IO (Matrix R) -> IO (Matrix R))
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String
-> Matrix R
-> Matrix R
-> Matrix R
forall {t} {t} {c}.
(Element t, Element t) =>
(IO (Matrix t) -> IO c)
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveSQAux IO (Matrix R) -> IO (Matrix R)
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgesv String
"linearSolveR" Matrix R
a Matrix R
b
mbLinearSolveR :: Matrix Double -> Matrix Double -> Maybe (Matrix Double)
mbLinearSolveR :: Matrix R -> Matrix R -> Maybe (Matrix R)
mbLinearSolveR Matrix R
a Matrix R
b = (IO (Matrix R) -> IO (Maybe (Matrix R)))
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String
-> Matrix R
-> Matrix R
-> Maybe (Matrix R)
forall {t} {t} {c}.
(Element t, Element t) =>
(IO (Matrix t) -> IO c)
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveSQAux IO (Matrix R) -> IO (Maybe (Matrix R))
forall x. IO x -> IO (Maybe x)
mbCatch CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgesv String
"linearSolveR" Matrix R
a Matrix R
b
linearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
linearSolveC :: Matrix C -> Matrix C -> Matrix C
linearSolveC Matrix C
a Matrix C
b = (IO (Matrix C) -> IO (Matrix C))
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String
-> Matrix C
-> Matrix C
-> Matrix C
forall {t} {t} {c}.
(Element t, Element t) =>
(IO (Matrix t) -> IO c)
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveSQAux IO (Matrix C) -> IO (Matrix C)
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgesv String
"linearSolveC" Matrix C
a Matrix C
b
mbLinearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
mbLinearSolveC :: Matrix C -> Matrix C -> Maybe (Matrix C)
mbLinearSolveC Matrix C
a Matrix C
b = (IO (Matrix C) -> IO (Maybe (Matrix C)))
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String
-> Matrix C
-> Matrix C
-> Maybe (Matrix C)
forall {t} {t} {c}.
(Element t, Element t) =>
(IO (Matrix t) -> IO c)
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveSQAux IO (Matrix C) -> IO (Maybe (Matrix C))
forall x. IO x -> IO (Maybe x)
mbCatch CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgesv String
"linearSolveC" Matrix C
a Matrix C
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 c)
-> Trans
(Matrix t) (CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveSQAux2 IO (Matrix t) -> IO c
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 c -> c
forall a. IO a -> a
unsafePerformIO (IO c -> c) -> (IO (Matrix t) -> IO c) -> IO (Matrix t) -> c
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IO (Matrix t) -> IO c
g (IO (Matrix t) -> c) -> IO (Matrix t) -> c
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 -> c
forall a. HasCallStack => String -> a
error (String -> c) -> String -> c
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
cholSolveR :: Matrix Double -> Matrix Double -> Matrix Double
cholSolveR :: Matrix R -> Matrix R -> Matrix R
cholSolveR Matrix R
a Matrix R
b = (IO (Matrix R) -> IO (Matrix R))
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String
-> Matrix R
-> Matrix R
-> Matrix R
forall {t} {t} {c}.
(Element t, Storable t) =>
(IO (Matrix t) -> IO c)
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveSQAux2 IO (Matrix R) -> IO (Matrix R)
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dpotrs String
"cholSolveR" (Matrix R -> Matrix R
forall t. Element t => Matrix t -> Matrix t
fmat Matrix R
a) Matrix R
b
cholSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
cholSolveC :: Matrix C -> Matrix C -> Matrix C
cholSolveC Matrix C
a Matrix C
b = (IO (Matrix C) -> IO (Matrix C))
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String
-> Matrix C
-> Matrix C
-> Matrix C
forall {t} {t} {c}.
(Element t, Storable t) =>
(IO (Matrix t) -> IO c)
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveSQAux2 IO (Matrix C) -> IO (Matrix C)
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zpotrs String
"cholSolveC" (Matrix C -> Matrix C
forall t. Element t => Matrix t -> Matrix t
fmat Matrix C
a) Matrix C
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 c)
-> Trans
(Matrix t) (CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveTRAux2 IO (Matrix t) -> IO c
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 c -> c
forall a. IO a -> a
unsafePerformIO (IO c -> c) -> (IO (Matrix t) -> IO c) -> IO (Matrix t) -> c
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IO (Matrix t) -> IO c
g (IO (Matrix t) -> c) -> IO (Matrix t) -> c
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 -> c
forall a. HasCallStack => String -> a
error (String -> c) -> String -> c
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
triSolveR :: UpLo -> Matrix Double -> Matrix Double -> Matrix Double
triSolveR :: UpLo -> Matrix R -> Matrix R -> Matrix R
triSolveR UpLo
Lower Matrix R
a Matrix R
b = (IO (Matrix R) -> IO (Matrix R))
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String
-> Matrix R
-> Matrix R
-> Matrix R
forall {t} {t} {c}.
(Element t, Storable t) =>
(IO (Matrix t) -> IO c)
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveTRAux2 IO (Matrix R) -> IO (Matrix R)
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dtrtrs_l String
"triSolveR" (Matrix R -> Matrix R
forall t. Element t => Matrix t -> Matrix t
fmat Matrix R
a) Matrix R
b
triSolveR UpLo
Upper Matrix R
a Matrix R
b = (IO (Matrix R) -> IO (Matrix R))
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String
-> Matrix R
-> Matrix R
-> Matrix R
forall {t} {t} {c}.
(Element t, Storable t) =>
(IO (Matrix t) -> IO c)
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveTRAux2 IO (Matrix R) -> IO (Matrix R)
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dtrtrs_u String
"triSolveR" (Matrix R -> Matrix R
forall t. Element t => Matrix t -> Matrix t
fmat Matrix R
a) Matrix R
b
triSolveC :: UpLo -> Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
triSolveC :: UpLo -> Matrix C -> Matrix C -> Matrix C
triSolveC UpLo
Lower Matrix C
a Matrix C
b = (IO (Matrix C) -> IO (Matrix C))
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String
-> Matrix C
-> Matrix C
-> Matrix C
forall {t} {t} {c}.
(Element t, Storable t) =>
(IO (Matrix t) -> IO c)
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveTRAux2 IO (Matrix C) -> IO (Matrix C)
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
ztrtrs_l String
"triSolveC" (Matrix C -> Matrix C
forall t. Element t => Matrix t -> Matrix t
fmat Matrix C
a) Matrix C
b
triSolveC UpLo
Upper Matrix C
a Matrix C
b = (IO (Matrix C) -> IO (Matrix C))
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String
-> Matrix C
-> Matrix C
-> Matrix C
forall {t} {t} {c}.
(Element t, Storable t) =>
(IO (Matrix t) -> IO c)
-> (CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveTRAux2 IO (Matrix C) -> IO (Matrix C)
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
ztrtrs_u String
"triSolveC" (Matrix C -> Matrix C
forall t. Element t => Matrix t -> Matrix t
fmat Matrix C
a) Matrix C
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 c)
-> (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
-> c
linearSolveGTAux2 IO (Matrix t) -> IO c
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 c -> c
forall a. IO a -> a
unsafePerformIO (IO c -> c) -> (IO (Matrix t) -> IO c) -> IO (Matrix t) -> c
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IO (Matrix t) -> IO c
g (IO (Matrix t) -> c) -> IO (Matrix t) -> c
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 -> c
forall a. HasCallStack => String -> a
error (String -> c) -> String -> c
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
triDiagSolveR :: Vector R -> Vector R -> Vector R -> Matrix R -> Matrix R
triDiagSolveR Vector R
dl Vector R
d Vector R
du Matrix R
b = (IO (Matrix R) -> IO (Matrix R))
-> (CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String
-> Vector R
-> Vector R
-> Vector R
-> Matrix R
-> Matrix R
forall {t} {t} {t} {t} {c}.
(Element t, Element t, Element t, Storable t) =>
(IO (Matrix t) -> IO c)
-> (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
-> c
linearSolveGTAux2 IO (Matrix R) -> IO (Matrix R)
forall a. a -> a
id CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgttrs String
"triDiagSolveR" Vector R
dl Vector R
d Vector R
du Matrix R
b
triDiagSolveC :: Vector C -> Vector C -> Vector C -> Matrix C -> Matrix C
triDiagSolveC Vector C
dl Vector C
d Vector C
du Matrix C
b = (IO (Matrix C) -> IO (Matrix C))
-> (CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String
-> Vector C
-> Vector C
-> Vector C
-> Matrix C
-> Matrix C
forall {t} {t} {t} {t} {c}.
(Element t, Element t, Element t, Storable t) =>
(IO (Matrix t) -> IO c)
-> (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
-> c
linearSolveGTAux2 IO (Matrix C) -> IO (Matrix C)
forall a. a -> a
id CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgttrs String
"triDiagSolveC" Vector C
dl Vector C
d Vector C
du Matrix C
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 a
-> IO CInt)
-> String -> Matrix t -> Matrix a -> Matrix a
linearSolveAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
f String
st Matrix t
a Matrix a
b
| Int
m Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Matrix a -> Int
forall t. Matrix t -> Int
rows Matrix a
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
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 a
r <- MatrixOrder -> Int -> Int -> IO (Matrix a)
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 a -> Matrix a -> IO ()
forall a. Element a => Int -> Int -> Matrix a -> Matrix a -> IO ()
setRect Int
0 Int
0 Matrix a
b Matrix a
r
(Matrix t
a' 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
r) Trans (Matrix t) (Trans (Matrix a) (IO CInt))
CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
f IO CInt -> String -> IO ()
#| String
st
Matrix a -> IO (Matrix a)
forall (m :: * -> *) a. Monad m => a -> m a
return Matrix a
r
| Bool
otherwise = String -> Matrix a
forall a. HasCallStack => String -> a
error (String -> Matrix a) -> String -> Matrix a
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 a -> Int
forall t. Matrix t -> Int
cols Matrix a
b
linearSolveLSR :: Matrix Double -> Matrix Double -> Matrix Double
linearSolveLSR :: Matrix R -> Matrix R -> Matrix R
linearSolveLSR Matrix R
a Matrix R
b = (Int, Int) -> (Int, Int) -> Matrix R -> Matrix R
forall a.
Element a =>
(Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
subMatrix (Int
0,Int
0) (Matrix R -> Int
forall t. Matrix t -> Int
cols Matrix R
a, Matrix R -> Int
forall t. Matrix t -> Int
cols Matrix R
b) (Matrix R -> Matrix R) -> Matrix R -> Matrix R
forall a b. (a -> b) -> a -> b
$
(CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String -> Matrix R -> Matrix R -> Matrix R
forall {t} {a}.
(Element t, Element a) =>
(CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt)
-> String -> Matrix t -> Matrix a -> Matrix a
linearSolveAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgels String
"linearSolverLSR" Matrix R
a Matrix R
b
linearSolveLSC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
linearSolveLSC :: Matrix C -> Matrix C -> Matrix C
linearSolveLSC Matrix C
a Matrix C
b = (Int, Int) -> (Int, Int) -> Matrix C -> Matrix C
forall a.
Element a =>
(Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
subMatrix (Int
0,Int
0) (Matrix C -> Int
forall t. Matrix t -> Int
cols Matrix C
a, Matrix C -> Int
forall t. Matrix t -> Int
cols Matrix C
b) (Matrix C -> Matrix C) -> Matrix C -> Matrix C
forall a b. (a -> b) -> a -> b
$
(CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String -> Matrix C -> Matrix C -> Matrix C
forall {t} {a}.
(Element t, Element a) =>
(CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt)
-> String -> Matrix t -> Matrix a -> Matrix a
linearSolveAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgels String
"linearSolveLSC" Matrix C
a Matrix C
b
linearSolveSVDR :: Maybe Double
-> Matrix Double
-> Matrix Double
-> Matrix Double
linearSolveSVDR :: Maybe R -> Matrix R -> Matrix R -> Matrix R
linearSolveSVDR (Just R
rcond) Matrix R
a Matrix R
b = (Int, Int) -> (Int, Int) -> Matrix R -> Matrix R
forall a.
Element a =>
(Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
subMatrix (Int
0,Int
0) (Matrix R -> Int
forall t. Matrix t -> Int
cols Matrix R
a, Matrix R -> Int
forall t. Matrix t -> Int
cols Matrix R
b) (Matrix R -> Matrix R) -> Matrix R -> Matrix R
forall a b. (a -> b) -> a -> b
$
(CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String -> Matrix R -> Matrix R -> Matrix R
forall {t} {a}.
(Element t, Element a) =>
(CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt)
-> String -> Matrix t -> Matrix a -> Matrix a
linearSolveAux (R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgelss R
rcond) String
"linearSolveSVDR" Matrix R
a Matrix R
b
linearSolveSVDR Maybe R
Nothing Matrix R
a Matrix R
b = Maybe R -> Matrix R -> Matrix R -> Matrix R
linearSolveSVDR (R -> Maybe R
forall a. a -> Maybe a
Just (-R
1)) Matrix R
a Matrix R
b
linearSolveSVDC :: Maybe Double
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
linearSolveSVDC :: Maybe R -> Matrix C -> Matrix C -> Matrix C
linearSolveSVDC (Just R
rcond) Matrix C
a Matrix C
b = (Int, Int) -> (Int, Int) -> Matrix C -> Matrix C
forall a.
Element a =>
(Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
subMatrix (Int
0,Int
0) (Matrix C -> Int
forall t. Matrix t -> Int
cols Matrix C
a, Matrix C -> Int
forall t. Matrix t -> Int
cols Matrix C
b) (Matrix C -> Matrix C) -> Matrix C -> Matrix C
forall a b. (a -> b) -> a -> b
$
(CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String -> Matrix C -> Matrix C -> Matrix C
forall {t} {a}.
(Element t, Element a) =>
(CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt)
-> String -> Matrix t -> Matrix a -> Matrix a
linearSolveAux (R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgelss R
rcond) String
"linearSolveSVDC" Matrix C
a Matrix C
b
linearSolveSVDC Maybe R
Nothing Matrix C
a Matrix C
b = Maybe R -> Matrix C -> Matrix C -> Matrix C
linearSolveSVDC (R -> Maybe R
forall a. a -> Maybe a
Just (-R
1)) Matrix C
a Matrix C
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
cholH :: Matrix (Complex Double) -> Matrix (Complex Double)
cholH :: Matrix C -> Matrix C
cholH = IO (Matrix C) -> Matrix C
forall a. IO a -> a
unsafePerformIO (IO (Matrix C) -> Matrix C)
-> (Matrix C -> IO (Matrix C)) -> Matrix C -> Matrix C
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CInt -> CInt -> CInt -> CInt -> Ptr C -> IO CInt)
-> String -> Matrix C -> IO (Matrix C)
forall {t}.
Element t =>
(CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> IO (Matrix t)
cholAux CInt -> CInt -> CInt -> CInt -> Ptr C -> IO CInt
zpotrf String
"cholH"
cholS :: Matrix Double -> Matrix Double
cholS :: Matrix R -> Matrix R
cholS = IO (Matrix R) -> Matrix R
forall a. IO a -> a
unsafePerformIO (IO (Matrix R) -> Matrix R)
-> (Matrix R -> IO (Matrix R)) -> Matrix R -> Matrix R
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt)
-> String -> Matrix R -> IO (Matrix R)
forall {t}.
Element t =>
(CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> IO (Matrix t)
cholAux CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt
dpotrf String
"cholS"
mbCholH :: Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
mbCholH :: Matrix C -> Maybe (Matrix C)
mbCholH = IO (Maybe (Matrix C)) -> Maybe (Matrix C)
forall a. IO a -> a
unsafePerformIO (IO (Maybe (Matrix C)) -> Maybe (Matrix C))
-> (Matrix C -> IO (Maybe (Matrix C)))
-> Matrix C
-> Maybe (Matrix C)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IO (Matrix C) -> IO (Maybe (Matrix C))
forall x. IO x -> IO (Maybe x)
mbCatch (IO (Matrix C) -> IO (Maybe (Matrix C)))
-> (Matrix C -> IO (Matrix C)) -> Matrix C -> IO (Maybe (Matrix C))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CInt -> CInt -> CInt -> CInt -> Ptr C -> IO CInt)
-> String -> Matrix C -> IO (Matrix C)
forall {t}.
Element t =>
(CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> IO (Matrix t)
cholAux CInt -> CInt -> CInt -> CInt -> Ptr C -> IO CInt
zpotrf String
"cholH"
mbCholS :: Matrix Double -> Maybe (Matrix Double)
mbCholS :: Matrix R -> Maybe (Matrix R)
mbCholS = IO (Maybe (Matrix R)) -> Maybe (Matrix R)
forall a. IO a -> a
unsafePerformIO (IO (Maybe (Matrix R)) -> Maybe (Matrix R))
-> (Matrix R -> IO (Maybe (Matrix R)))
-> Matrix R
-> Maybe (Matrix R)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IO (Matrix R) -> IO (Maybe (Matrix R))
forall x. IO x -> IO (Maybe x)
mbCatch (IO (Matrix R) -> IO (Maybe (Matrix R)))
-> (Matrix R -> IO (Matrix R)) -> Matrix R -> IO (Maybe (Matrix R))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt)
-> String -> Matrix R -> IO (Matrix R)
forall {t}.
Element t =>
(CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> IO (Matrix t)
cholAux CInt -> CInt -> CInt -> CInt -> Ptr R -> 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
qrR :: Matrix Double -> (Matrix Double, Vector Double)
qrR :: Matrix R -> (Matrix R, Vector R)
qrR = (CInt -> Ptr R -> CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt)
-> String -> Matrix R -> (Matrix R, Vector R)
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 R -> CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt
dgeqr2 String
"qrR"
qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
qrC :: Matrix C -> (Matrix C, Vector C)
qrC = (CInt -> Ptr C -> CInt -> CInt -> CInt -> CInt -> Ptr C -> IO CInt)
-> String -> Matrix C -> (Matrix C, Vector C)
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 C -> CInt -> CInt -> CInt -> CInt -> Ptr C -> 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
qrgrR :: Int -> (Matrix Double, Vector Double) -> Matrix Double
qrgrR :: Int -> (Matrix R, Vector R) -> Matrix R
qrgrR = (CInt -> Ptr R -> CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt)
-> String -> Int -> (Matrix R, Vector R) -> Matrix R
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 R -> CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt
dorgqr String
"qrgrR"
qrgrC :: Int -> (Matrix (Complex Double), Vector (Complex Double)) -> Matrix (Complex Double)
qrgrC :: Int -> (Matrix C, Vector C) -> Matrix C
qrgrC = (CInt -> Ptr C -> CInt -> CInt -> CInt -> CInt -> Ptr C -> IO CInt)
-> String -> Int -> (Matrix C, Vector C) -> Matrix C
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 C -> CInt -> CInt -> CInt -> CInt -> Ptr C -> 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
hessR :: Matrix Double -> (Matrix Double, Vector Double)
hessR :: Matrix R -> (Matrix R, Vector R)
hessR = (CInt -> Ptr R -> CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt)
-> String -> Matrix R -> (Matrix R, Vector R)
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 R -> CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt
dgehrd String
"hessR"
hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
hessC :: Matrix C -> (Matrix C, Vector C)
hessC = (CInt -> Ptr C -> CInt -> CInt -> CInt -> CInt -> Ptr C -> IO CInt)
-> String -> Matrix C -> (Matrix C, Vector C)
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 C -> CInt -> CInt -> CInt -> CInt -> Ptr C -> 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
schurR :: Matrix Double -> (Matrix Double, Matrix Double)
schurR :: Matrix R -> (Matrix R, Matrix R)
schurR = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String -> Matrix R -> (Matrix R, Matrix R)
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 R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgees String
"schurR"
schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double))
schurC :: Matrix C -> (Matrix C, Matrix C)
schurC = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String -> Matrix C -> (Matrix C, Matrix C)
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 C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> 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
luR :: Matrix Double -> (Matrix Double, [Int])
luR :: Matrix R -> (Matrix R, [Int])
luR = (CInt -> Ptr R -> CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt)
-> String -> Matrix R -> (Matrix R, [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 R -> CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt
dgetrf String
"luR"
luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
luC :: Matrix C -> (Matrix C, [Int])
luC = (CInt -> Ptr R -> CInt -> CInt -> CInt -> CInt -> Ptr C -> IO CInt)
-> String -> Matrix C -> (Matrix C, [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 R -> CInt -> CInt -> CInt -> CInt -> Ptr C -> 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
lusR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double
lusR :: Matrix R -> [Int] -> Matrix R -> Matrix R
lusR Matrix R
a [Int]
piv Matrix R
b = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String -> Matrix R -> [Int] -> Matrix R -> Matrix R
forall {t} {t} {a}.
(Element t, Storable t, Integral a) =>
(CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt)
-> String -> Matrix t -> [a] -> Matrix t -> Matrix t
lusAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgetrs String
"lusR" (Matrix R -> Matrix R
forall t. Element t => Matrix t -> Matrix t
fmat Matrix R
a) [Int]
piv Matrix R
b
lusC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
lusC :: Matrix C -> [Int] -> Matrix C -> Matrix C
lusC Matrix C
a [Int]
piv Matrix C
b = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String -> Matrix C -> [Int] -> Matrix C -> Matrix C
forall {t} {t} {a}.
(Element t, Storable t, Integral a) =>
(CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt)
-> String -> Matrix t -> [a] -> Matrix t -> Matrix t
lusAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgetrs String
"lusC" (Matrix C -> Matrix C
forall t. Element t => Matrix t -> Matrix t
fmat Matrix C
a) [Int]
piv Matrix C
b
lusAux :: Trans
(Matrix t)
(CInt -> Ptr R -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> [a] -> Matrix t -> Matrix t
lusAux Trans
(Matrix t)
(CInt -> Ptr R -> 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 R -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> IO CInt)
-> Trans
(Matrix t)
(CInt -> Ptr R -> 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 R
piv' Vector R
-> Matrix t
-> Trans (Vector R) (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 R -> 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 R
piv' = [R] -> Vector R
forall a. Storable a => [a] -> Vector a
fromList ((a -> R) -> [a] -> [R]
forall a b. (a -> b) -> [a] -> [b]
map (a -> R
forall a b. (Integral a, Num b) => a -> b
fromIntegral(a -> R) -> (a -> a) -> a -> R
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
ldlR :: Matrix Double -> (Matrix Double, [Int])
ldlR :: Matrix R -> (Matrix R, [Int])
ldlR = (CInt -> Ptr R -> CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt)
-> String -> Matrix R -> (Matrix R, [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 R -> CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt
dsytrf String
"ldlR"
ldlC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
ldlC :: Matrix C -> (Matrix C, [Int])
ldlC = (CInt -> Ptr R -> CInt -> CInt -> CInt -> CInt -> Ptr C -> IO CInt)
-> String -> Matrix C -> (Matrix C, [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 R -> CInt -> CInt -> CInt -> CInt -> Ptr C -> 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
ldlsR :: Matrix Double -> [Int] -> Matrix Double -> Matrix Double
ldlsR :: Matrix R -> [Int] -> Matrix R -> Matrix R
ldlsR Matrix R
a [Int]
piv Matrix R
b = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt)
-> String -> Matrix R -> [Int] -> Matrix R -> Matrix R
forall {t} {t} {a}.
(Element t, Storable t, Integral a) =>
(CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt)
-> String -> Matrix t -> [a] -> Matrix t -> Matrix t
lusAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dsytrs String
"ldlsR" (Matrix R -> Matrix R
forall t. Element t => Matrix t -> Matrix t
fmat Matrix R
a) [Int]
piv Matrix R
b
ldlsC :: Matrix (Complex Double) -> [Int] -> Matrix (Complex Double) -> Matrix (Complex Double)
ldlsC :: Matrix C -> [Int] -> Matrix C -> Matrix C
ldlsC Matrix C
a [Int]
piv Matrix C
b = (CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt)
-> String -> Matrix C -> [Int] -> Matrix C -> Matrix C
forall {t} {t} {a}.
(Element t, Storable t, Integral a) =>
(CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> IO CInt)
-> String -> Matrix t -> [a] -> Matrix t -> Matrix t
lusAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zsytrs String
"ldlsC" (Matrix C -> Matrix C
forall t. Element t => Matrix t -> Matrix t
fmat Matrix C
a) [Int]
piv Matrix C
b