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

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

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


module Internal.LAPACK where

import Data.Bifunctor (first)

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

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

infixr 1 #
c
a # :: c -> (b -> IO r) -> Trans c b -> IO r
# b -> IO r
b = 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 forall {c} {b} {r}.
TransArray c =>
c -> (b -> IO r) -> Trans c b -> IO r
# c
b forall {c} {b} {r}.
TransArray c =>
c -> (b -> IO r) -> Trans c b -> 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 -> a
isT (forall t. Matrix t -> Bool
rowOrder -> Bool
False) = a
0
isT Matrix t
_                   = a
1

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

-- | Matrix product based on BLAS's /dgemm/.
multiplyR :: Matrix Double -> Matrix Double -> Matrix Double
multiplyR :: Matrix Double -> Matrix Double -> Matrix Double
multiplyR Matrix Double
a Matrix Double
b = {-# SCC "multiplyR" #-} 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 -> TMMM Double
dgemmc String
"dgemmc" Matrix Double
a Matrix Double
b

-- | Matrix product based on BLAS's /zgemm/.
multiplyC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
multiplyC :: Matrix (Complex Double)
-> Matrix (Complex Double) -> Matrix (Complex Double)
multiplyC Matrix (Complex Double)
a Matrix (Complex Double)
b = 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 -> TMMM (Complex Double)
zgemmc String
"zgemmc" Matrix (Complex Double)
a Matrix (Complex Double)
b

-- | Matrix product based on BLAS's /sgemm/.
multiplyF :: Matrix Float -> Matrix Float -> Matrix Float
multiplyF :: Matrix Float -> Matrix Float -> Matrix Float
multiplyF Matrix Float
a Matrix Float
b = 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 -> TMMM Float
sgemmc String
"sgemmc" Matrix Float
a Matrix Float
b

-- | Matrix product based on BLAS's /cgemm/.
multiplyQ :: Matrix (Complex Float) -> Matrix (Complex Float) -> Matrix (Complex Float)
multiplyQ :: Matrix (Complex Float)
-> Matrix (Complex Float) -> Matrix (Complex Float)
multiplyQ Matrix (Complex Float)
a Matrix (Complex Float)
b = 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 -> TMMM (Complex Float)
cgemmc String
"cgemmc" Matrix (Complex Float)
a Matrix (Complex Float)
b

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

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

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

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

-- | Full SVD of a real matrix using LAPACK's /dgesvd/.
svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
svdR = 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 TSVD Double
dgesvd String
"svdR"

-- | Full SVD of a real matrix using LAPACK's /dgesdd/.
svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
svdRd = 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 TSVD Double
dgesdd String
"svdRdd"

-- | Full SVD of a complex matrix using LAPACK's /zgesvd/.
svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
svdC :: Matrix (Complex Double)
-> (Matrix (Complex Double), Vector Double,
    Matrix (Complex Double))
svdC = 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 TSVD (Complex Double)
zgesvd String
"svdC"

-- | Full SVD of a complex matrix using LAPACK's /zgesdd/.
svdCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
svdCd :: Matrix (Complex Double)
-> (Matrix (Complex Double), Vector Double,
    Matrix (Complex Double))
svdCd = 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 TSVD (Complex Double)
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 = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$ do
    Matrix t
a <- forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
x
    Matrix a
u <- forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
r Int
r
    Vector a
s <- forall a. Storable a => Int -> IO (Vector a)
createVector (forall a. Ord a => a -> a -> a
min Int
r Int
c)
    Matrix a
v <- forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
c Int
c
    (Matrix t
a forall {c} {b} {r}.
TransArray c =>
c -> (b -> IO r) -> Trans c b -> IO r
# Matrix a
u forall {c} {b} {r}.
TransArray c =>
c -> (b -> IO r) -> Trans c b -> IO r
# Vector a
s forall {c} {c} {r}.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix a
v) 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
    forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix a
u,Vector a
s,Matrix a
v)
  where
    r :: Int
r = forall t. Matrix t -> Int
rows Matrix t
x
    c :: Int
c = forall t. Matrix t -> Int
cols Matrix t
x


-- | Thin SVD of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'S\'.
thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
thinSVDR = 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 TSVD Double
dgesvd String
"thinSVDR"

-- | Thin SVD of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'S\'.
thinSVDC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
thinSVDC :: Matrix (Complex Double)
-> (Matrix (Complex Double), Vector Double,
    Matrix (Complex Double))
thinSVDC = 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 TSVD (Complex Double)
zgesvd String
"thinSVDC"

-- | Thin SVD of a real matrix, using LAPACK's /dgesdd/ with jobz == \'S\'.
thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
thinSVDRd = 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 TSVD Double
dgesdd String
"thinSVDRdd"

-- | Thin SVD of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'S\'.
thinSVDCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
thinSVDCd :: Matrix (Complex Double)
-> (Matrix (Complex Double), Vector Double,
    Matrix (Complex Double))
thinSVDCd = 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 TSVD (Complex Double)
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 = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$ do
    Matrix t
a <- forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
x
    Matrix a
u <- forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
r Int
q
    Vector a
s <- forall a. Storable a => Int -> IO (Vector a)
createVector Int
q
    Matrix a
v <- forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
q Int
c
    (Matrix t
a forall {c} {b} {r}.
TransArray c =>
c -> (b -> IO r) -> Trans c b -> IO r
# Matrix a
u forall {c} {b} {r}.
TransArray c =>
c -> (b -> IO r) -> Trans c b -> IO r
# Vector a
s forall {c} {c} {r}.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix a
v) 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
    forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix a
u,Vector a
s,Matrix a
v)
  where
    r :: Int
r = forall t. Matrix t -> Int
rows Matrix t
x
    c :: Int
c = forall t. Matrix t -> Int
cols Matrix t
x
    q :: Int
q = forall a. Ord a => a -> a -> a
min Int
r Int
c


-- | Singular values of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'N\'.
svR :: Matrix Double -> Vector Double
svR :: Matrix Double -> Vector Double
svR = 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 TSVD Double
dgesvd String
"svR"

-- | Singular values of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'N\'.
svC :: Matrix (Complex Double) -> Vector Double
svC :: Matrix (Complex Double) -> Vector Double
svC = 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 TSVD (Complex Double)
zgesvd String
"svC"

-- | Singular values of a real matrix, using LAPACK's /dgesdd/ with jobz == \'N\'.
svRd :: Matrix Double -> Vector Double
svRd :: Matrix Double -> Vector Double
svRd = 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 TSVD Double
dgesdd String
"svRd"

-- | Singular values of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'N\'.
svCd :: Matrix (Complex Double) -> Vector Double
svCd :: Matrix (Complex Double) -> Vector Double
svCd = 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 TSVD (Complex Double)
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 = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$ do
    Matrix t
a <- forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
x
    Vector a
s <- forall a. Storable a => Int -> IO (Vector a)
createVector Int
q
    (Matrix t
a forall {c} {c} {r}.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Vector a
s) CInt -> CInt -> CInt -> CInt -> Ptr t -> CInt -> Ptr a -> IO CInt
g IO CInt -> String -> IO ()
#| String
st
    forall (m :: * -> *) a. Monad m => a -> m a
return Vector a
s
  where
    r :: Int
r = forall t. Matrix t -> Int
rows Matrix t
x
    c :: Int
c = forall t. Matrix t -> Int
cols Matrix t
x
    q :: Int
q = 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 forall a. Ptr a
nullPtr CInt
nb Ptr a
pb t
0 t
0 t
0 t
0 forall a. Ptr a
nullPtr


-- | Singular values and all right singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'N\' and jobvt == \'A\'.
rightSVR :: Matrix Double -> (Vector Double, Matrix Double)
rightSVR :: Matrix Double -> (Vector Double, Matrix Double)
rightSVR = 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 TSVD Double
dgesvd String
"rightSVR"

-- | Singular values and all right singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'N\' and jobvt == \'A\'.
rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
rightSVC = 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 TSVD (Complex Double)
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 = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$ do
    Matrix t
a <- forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
x
    Vector a
s <- forall a. Storable a => Int -> IO (Vector a)
createVector Int
q
    Matrix a
v <- forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
c Int
c
    (Matrix t
a forall {c} {b} {r}.
TransArray c =>
c -> (b -> IO r) -> Trans c b -> IO r
# Vector a
s forall {c} {c} {r}.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix a
v) CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
g IO CInt -> String -> IO ()
#| String
st
    forall (m :: * -> *) a. Monad m => a -> m a
return (Vector a
s,Matrix a
v)
  where
    r :: Int
r = forall t. Matrix t -> Int
rows Matrix t
x
    c :: Int
c = forall t. Matrix t -> Int
cols Matrix t
x
    q :: Int
q = 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 forall a. Ptr a
nullPtr


-- | Singular values and all left singular vectors of a real matrix, using LAPACK's /dgesvd/  with jobu == \'A\' and jobvt == \'N\'.
leftSVR :: Matrix Double -> (Matrix Double, Vector Double)
leftSVR :: Matrix Double -> (Matrix Double, Vector Double)
leftSVR = 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 TSVD Double
dgesvd String
"leftSVR"

-- | Singular values and all left singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'A\' and jobvt == \'N\'.
leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double)
leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double)
leftSVC = 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 TSVD (Complex Double)
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 = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$ do
    Matrix t
a <- forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
x
    Matrix a
u <- forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
r Int
r
    Vector a
s <- forall a. Storable a => Int -> IO (Vector a)
createVector Int
q
    (Matrix t
a forall {c} {b} {r}.
TransArray c =>
c -> (b -> IO r) -> Trans c b -> IO r
# Matrix a
u forall {c} {c} {r}.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Vector a
s) CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> CInt
-> Ptr a
-> IO CInt
g IO CInt -> String -> IO ()
#| String
st
    forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix a
u,Vector a
s)
  where
    r :: Int
r = forall t. Matrix t -> Int
rows Matrix t
x
    c :: Int
c = forall t. Matrix t -> Int
cols Matrix t
x
    q :: Int
q = 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 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 = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$ do
    Matrix t
a <- forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
m
    Vector a
l <- forall a. Storable a => Int -> IO (Vector a)
createVector Int
r
    Matrix a
v <- forall a. Storable a => MatrixOrder -> Int -> Int -> IO (Matrix a)
createMatrix MatrixOrder
ColumnMajor Int
r Int
r
    (Matrix t
a forall {c} {b} {r}.
TransArray c =>
c -> (b -> IO r) -> Trans c b -> IO r
# Vector a
l forall {c} {c} {r}.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix a
v) CInt
-> CInt
-> CInt
-> CInt
-> Ptr t
-> CInt
-> Ptr a
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr a
-> IO CInt
g IO CInt -> String -> IO ()
#| String
st
    forall (m :: * -> *) a. Monad m => a -> m a
return (Vector a
l,Matrix a
v)
  where
    r :: Int
r = 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 forall a. Ptr a
nullPtr


-- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/.
-- The eigenvectors are the columns of v. The eigenvalues are not sorted.
eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double))
eigC :: Matrix (Complex Double)
-> (Vector (Complex Double), Matrix (Complex Double))
eigC = 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 Complex Double
::> (Complex Double
     ::> (Complex Double :> (Complex Double ::> IO CInt)))
zgeev String
"eigC"

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

-- | Eigenvalues of a general complex matrix, using LAPACK's /zgeev/ with jobz == \'N\'.
-- The eigenvalues are not sorted.
eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double)
eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double)
eigOnlyC = 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 Complex Double
::> (Complex Double
     ::> (Complex Double :> (Complex Double ::> IO CInt)))
zgeev String
"eigOnlyC"

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

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

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

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

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

-- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'.
-- The eigenvalues are not sorted.
eigOnlyR :: Matrix Double -> Vector (Complex Double)
eigOnlyR :: Matrix Double -> Vector (Complex Double)
eigOnlyR = forall {e}.
RealElement e =>
Vector (Complex e) -> Vector (Complex e)
fixeig1 forall b c a. (b -> c) -> (a -> b) -> a -> 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 Double
::> (CInt
     -> CInt
     -> CInt
     -> CInt
     -> Ptr Double
     -> Complex Double :> (Double ::> IO CInt))
dgeev String
"eigOnlyR"

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

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

    (Matrix t
a forall {c} {b} {r}.
TransArray c =>
c -> (b -> IO r) -> Trans c b -> IO r
# Matrix t
b forall {c} {b} {r}.
TransArray c =>
c -> (b -> IO r) -> Trans c b -> IO r
# Vector a
alpha forall {c} {b} {r}.
TransArray c =>
c -> (b -> IO r) -> Trans c b -> IO r
# Vector a
beta forall {c} {c} {r}.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Matrix a
vr) 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

    forall (m :: * -> *) a. Monad m => a -> m a
return (Vector a
alpha, Vector a
beta, Matrix a
vr)
  where
    r :: Int
r = 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 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 = forall a. IO a -> a
unsafePerformIO forall a b. (a -> b) -> a -> b
$ do
    Matrix t
a <- forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
ma
    Matrix t
b <- forall t. Element t => MatrixOrder -> Matrix t -> IO (Matrix t)
copy MatrixOrder
ColumnMajor Matrix t
mb
    Vector a
alpha <- forall a. Storable a => Int -> IO (Vector a)
createVector Int
r
    Vector a
beta <- forall a. Storable a => Int -> IO (Vector a)
createVector Int
r

    (Matrix t
a forall {c} {b} {r}.
TransArray c =>
c -> (b -> IO r) -> Trans c b -> IO r
# Matrix t
b forall {c} {b} {r}.
TransArray c =>
c -> (b -> IO r) -> Trans c b -> IO r
# Vector a
alpha forall {c} {c} {r}.
(TransArray c, TransArray c) =>
c -> c -> Trans c (Trans c (IO r)) -> IO r
#! Vector a
beta) 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

    forall (m :: * -> *) a. Monad m => a -> m a
return (Vector a
alpha, Vector a
beta)
  where
    r :: Int
r = 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 forall a. Ptr a
nullPtr t
0 t
0 t
0 t
0 forall a. Ptr a
nullPtr

-- | Generalized eigenvalues and right eigenvectors of a pair of complex matrices, using LAPACK's /zggev/.
-- The eigenvectors are the columns of v. The eigenvalues are represented as alphas / betas and not sorted.
eigGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double), Matrix (Complex Double))
eigGC :: Matrix (Complex Double)
-> Matrix (Complex Double)
-> (Vector (Complex Double), Vector (Complex Double),
    Matrix (Complex Double))
eigGC Matrix (Complex Double)
a Matrix (Complex Double)
b = 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 Complex Double
::> (Complex Double
     ::> (Complex Double
          :> (Complex Double
              :> (Complex Double ::> (Complex Double ::> IO CInt)))))
zggev Matrix (Complex Double)
a Matrix (Complex Double)
b String
"eigGC"

eigOnlyG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double)
eigOnlyG :: Matrix Double
-> Matrix Double -> (Vector (Complex Double), Vector Double)
eigOnlyG Matrix Double
a Matrix Double
b = forall (p :: * -> * -> *) a b c.
Bifunctor p =>
(a -> b) -> p a c -> p b c
first forall {e}.
RealElement e =>
Vector (Complex e) -> Vector (Complex e)
fixeig1 forall a b. (a -> b) -> a -> b
$ 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 Double
::> (Double
     ::> (Complex Double
          :> (Double :> (Double ::> (Double ::> IO CInt)))))
dggev Matrix Double
a Matrix Double
b String
"eigOnlyG"

eigOnlyGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double))
eigOnlyGC :: Matrix (Complex Double)
-> Matrix (Complex Double)
-> (Vector (Complex Double), Vector (Complex Double))
eigOnlyGC Matrix (Complex Double)
a Matrix (Complex Double)
b = 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 Complex Double
::> (Complex Double
     ::> (Complex Double
          :> (Complex Double
              :> (Complex Double ::> (Complex Double ::> IO CInt)))))
zggev Matrix (Complex Double)
a Matrix (Complex Double)
b String
"eigOnlyGC"

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

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

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

-- | 'eigS' in ascending order
eigS' :: Matrix Double -> (Vector Double, Matrix Double)
eigS' :: Matrix Double -> (Vector Double, Matrix Double)
eigS' = 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 -> Double :> (Double ::> IO CInt)
dsyev CInt
1) String
"eigS'"

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

-- | 'eigH' in ascending order
eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
eigH' = 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 -> Double :> (Complex Double ::> IO CInt)
zheev CInt
1) String
"eigH'"


-- | Eigenvalues of a symmetric real matrix, using LAPACK's /dsyev/ with jobz == \'N\'.
-- The eigenvalues are sorted in descending order.
eigOnlyS :: Matrix Double -> Vector Double
eigOnlyS :: Matrix Double -> Vector Double
eigOnlyS = Vector Double -> Vector Double
vrev forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a, b) -> a
fstforall b c a. (b -> c) -> (a -> b) -> a -> 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 -> Double :> (Double ::> IO CInt)
dsyev CInt
0) String
"eigS'"

-- | Eigenvalues of a hermitian complex matrix, using LAPACK's /zheev/ with jobz == \'N\'.
-- The eigenvalues are sorted in descending order.
eigOnlyH :: Matrix (Complex Double) -> Vector Double
eigOnlyH :: Matrix (Complex Double) -> Vector Double
eigOnlyH = Vector Double -> Vector Double
vrev forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall a b. (a, b) -> a
fstforall b c a. (b -> c) -> (a -> b) -> a -> 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 -> Double :> (Complex Double ::> IO CInt)
zheev CInt
0) String
"eigH'"

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

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

mbLinearSolveR :: Matrix Double -> Matrix Double -> Maybe (Matrix Double)
mbLinearSolveR :: Matrix Double -> Matrix Double -> Maybe (Matrix Double)
mbLinearSolveR Matrix Double
a Matrix Double
b = 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 forall x. IO x -> IO (Maybe x)
mbCatch Double ::> (Double ::> IO CInt)
dgesv String
"linearSolveR" Matrix Double
a Matrix Double
b


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

mbLinearSolveC :: Matrix (Complex Double) -> Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
mbLinearSolveC :: Matrix (Complex Double)
-> Matrix (Complex Double) -> Maybe (Matrix (Complex Double))
mbLinearSolveC Matrix (Complex Double)
a Matrix (Complex Double)
b = 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 forall x. IO x -> IO (Maybe x)
mbCatch Complex Double ::> (Complex Double ::> IO CInt)
zgesv String
"linearSolveC" Matrix (Complex Double)
a Matrix (Complex Double)
b

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


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

-- | Solves a symmetric positive definite system of linear equations using a precomputed Cholesky factorization obtained by 'cholS'.
cholSolveR :: Matrix Double -> Matrix Double -> Matrix Double
cholSolveR :: Matrix Double -> Matrix Double -> Matrix Double
cholSolveR Matrix Double
a Matrix Double
b = 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 forall a. a -> a
id Double ::> (Double ::> IO CInt)
dpotrs String
"cholSolveR" (forall t. Element t => Matrix t -> Matrix t
fmat Matrix Double
a) Matrix Double
b

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

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


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

data UpLo = Lower | Upper

-- | Solves a triangular system of linear equations.
triSolveR :: UpLo -> Matrix Double -> Matrix Double -> Matrix Double
triSolveR :: UpLo -> Matrix Double -> Matrix Double -> Matrix Double
triSolveR UpLo
Lower Matrix Double
a Matrix Double
b = 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 forall a. a -> a
id Double ::> (Double ::> IO CInt)
dtrtrs_l String
"triSolveR" (forall t. Element t => Matrix t -> Matrix t
fmat Matrix Double
a) Matrix Double
b
triSolveR UpLo
Upper Matrix Double
a Matrix Double
b = 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 forall a. a -> a
id Double ::> (Double ::> IO CInt)
dtrtrs_u String
"triSolveR" (forall t. Element t => Matrix t -> Matrix t
fmat Matrix Double
a) Matrix Double
b

-- | Solves a triangular system of linear equations.
triSolveC :: UpLo -> Matrix (Complex Double) -> Matrix (Complex Double) -> Matrix (Complex Double)
triSolveC :: UpLo
-> Matrix (Complex Double)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
triSolveC UpLo
Lower Matrix (Complex Double)
a Matrix (Complex Double)
b = 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 forall a. a -> a
id Complex Double ::> (Complex Double ::> IO CInt)
ztrtrs_l String
"triSolveC" (forall t. Element t => Matrix t -> Matrix t
fmat Matrix (Complex Double)
a) Matrix (Complex Double)
b
triSolveC UpLo
Upper Matrix (Complex Double)
a Matrix (Complex Double)
b = 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 forall a. a -> a
id Complex Double ::> (Complex Double ::> IO CInt)
ztrtrs_u String
"triSolveC" (forall t. Element t => Matrix t -> Matrix t
fmat Matrix (Complex Double)
a) Matrix (Complex Double)
b

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

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

-- | Solves a tridiagonal system of linear equations.
triDiagSolveR :: Vector Double
-> Vector Double -> Vector Double -> Matrix Double -> Matrix Double
triDiagSolveR Vector Double
dl Vector Double
d Vector Double
du Matrix Double
b = 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 forall a. a -> a
id Double :> (Double :> (Double :> (Double ::> IO CInt)))
dgttrs String
"triDiagSolveR" Vector Double
dl Vector Double
d Vector Double
du Matrix Double
b
triDiagSolveC :: Vector (Complex Double)
-> Vector (Complex Double)
-> Vector (Complex Double)
-> Matrix (Complex Double)
-> Matrix (Complex Double)
triDiagSolveC Vector (Complex Double)
dl Vector (Complex Double)
d Vector (Complex Double)
du Matrix (Complex Double)
b = 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 forall a. a -> a
id Complex Double
:> (Complex Double
    :> (Complex Double :> (Complex Double ::> IO CInt)))
zgttrs String
"triDiagSolveC" Vector (Complex Double)
dl Vector (Complex Double)
d Vector (Complex Double)
du Matrix (Complex Double)
b

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

-- | QR factorization of a real matrix, using LAPACK's /dgeqr2/.
qrR :: Matrix Double -> (Matrix Double, Vector Double)
qrR :: Matrix Double -> (Matrix Double, Vector Double)
qrR = 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 Double :> (Double ::> IO CInt)
dgeqr2 String
"qrR"

-- | QR factorization of a complex matrix, using LAPACK's /zgeqr2/.
qrC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
qrC :: Matrix (Complex Double)
-> (Matrix (Complex Double), Vector (Complex Double))
qrC = 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 Complex Double :> (Complex Double ::> IO CInt)
zgeqr2 String
"qrC"

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

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

-- | build rotation from reflectors
qrgrR :: Int -> (Matrix Double, Vector Double) -> Matrix Double
qrgrR :: Int -> (Matrix Double, Vector Double) -> Matrix Double
qrgrR = 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 Double :> (Double ::> IO CInt)
dorgqr String
"qrgrR"
-- | build rotation from reflectors
qrgrC :: Int -> (Matrix (Complex Double), Vector (Complex Double)) -> Matrix (Complex Double)
qrgrC :: Int
-> (Matrix (Complex Double), Vector (Complex Double))
-> Matrix (Complex Double)
qrgrC = 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 Complex Double :> (Complex Double ::> IO CInt)
zungqr String
"qrgrC"

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

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

-- | Hessenberg factorization of a square real matrix, using LAPACK's /dgehrd/.
hessR :: Matrix Double -> (Matrix Double, Vector Double)
hessR :: Matrix Double -> (Matrix Double, Vector Double)
hessR = 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 Double :> (Double ::> IO CInt)
dgehrd String
"hessR"

-- | Hessenberg factorization of a square complex matrix, using LAPACK's /zgehrd/.
hessC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector (Complex Double))
hessC :: Matrix (Complex Double)
-> (Matrix (Complex Double), Vector (Complex Double))
hessC = 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 Complex Double :> (Complex Double ::> IO CInt)
zgehrd String
"hessC"

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

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

-- | Schur factorization of a square real matrix, using LAPACK's /dgees/.
schurR :: Matrix Double -> (Matrix Double, Matrix Double)
schurR :: Matrix Double -> (Matrix Double, Matrix Double)
schurR = 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 Double ::> (Double ::> IO CInt)
dgees String
"schurR"

-- | Schur factorization of a square complex matrix, using LAPACK's /zgees/.
schurC :: Matrix (Complex Double) -> (Matrix (Complex Double), Matrix (Complex Double))
schurC :: Matrix (Complex Double)
-> (Matrix (Complex Double), Matrix (Complex Double))
schurC = 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 Complex Double ::> (Complex Double ::> IO CInt)
zgees String
"schurC"

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

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

-- | LU factorization of a general real matrix, using LAPACK's /dgetrf/.
luR :: Matrix Double -> (Matrix Double, [Int])
luR :: Matrix Double -> (Matrix Double, [Int])
luR = 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 Double :> (Double ::> IO CInt)
dgetrf String
"luR"

-- | LU factorization of a general complex matrix, using LAPACK's /zgetrf/.
luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
luC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
luC = 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 Double :> (Complex Double ::> IO CInt)
zgetrf String
"luC"

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

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

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

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

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

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

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

-- | LDL factorization of a symmetric real matrix, using LAPACK's /dsytrf/.
ldlR :: Matrix Double -> (Matrix Double, [Int])
ldlR :: Matrix Double -> (Matrix Double, [Int])
ldlR = 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 Double :> (Double ::> IO CInt)
dsytrf String
"ldlR"

-- | LDL factorization of a hermitian complex matrix, using LAPACK's /zhetrf/.
ldlC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
ldlC :: Matrix (Complex Double) -> (Matrix (Complex Double), [Int])
ldlC = 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 Double :> (Complex Double ::> IO CInt)
zhetrf String
"ldlC"

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

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

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

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

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