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

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

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


module Internal.LAPACK where

import Data.Bifunctor (first)

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

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

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

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

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

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

type F = Float
type Q = Complex Float

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


-- | 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 R -> (Matrix R, Vector R, Matrix R)
thinSVDR = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> IO CInt)
-> String -> Matrix R -> (Matrix R, Vector R, Matrix R)
forall {t} {a} {a} {a}.
(Element t, Storable a, Storable a, Storable a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a, Matrix a)
thinSVDAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgesvd String
"thinSVDR"

-- | 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 C -> (Matrix C, Vector R, Matrix C)
thinSVDC = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> IO CInt)
-> String -> Matrix C -> (Matrix C, Vector R, Matrix C)
forall {t} {a} {a} {a}.
(Element t, Storable a, Storable a, Storable a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a, Matrix a)
thinSVDAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgesvd String
"thinSVDC"

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

-- | 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 C -> (Matrix C, Vector R, Matrix C)
thinSVDCd = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> IO CInt)
-> String -> Matrix C -> (Matrix C, Vector R, Matrix C)
forall {t} {a} {a} {a}.
(Element t, Storable a, Storable a, Storable a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a, Matrix a)
thinSVDAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgesdd String
"thinSVDCdd"

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


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

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

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

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

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


-- | 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 R -> (Vector R, Matrix R)
rightSVR = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> IO CInt)
-> String -> Matrix R -> (Vector R, Matrix R)
forall {t} {a} {a} {t} {t} {t} {t} {a}.
(Element t, Storable a, Storable a, Num t, Num t, Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Vector a, Matrix a)
rightSVAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgesvd String
"rightSVR"

-- | 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 C -> (Vector R, Matrix C)
rightSVC = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> IO CInt)
-> String -> Matrix C -> (Vector R, Matrix C)
forall {t} {a} {a} {t} {t} {t} {t} {a}.
(Element t, Storable a, Storable a, Num t, Num t, Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Vector a, Matrix a)
rightSVAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgesvd String
"rightSVC"

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


-- | 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 R -> (Matrix R, Vector R)
leftSVR = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> IO CInt)
-> String -> Matrix R -> (Matrix R, Vector R)
forall {t} {a} {a} {t} {t} {t} {t} {a}.
(Element t, Storable a, Storable a, Num t, Num t, Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a)
leftSVAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgesvd String
"leftSVR"

-- | 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 C -> (Matrix C, Vector R)
leftSVC = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> IO CInt)
-> String -> Matrix C -> (Matrix C, Vector R)
forall {t} {a} {a} {t} {t} {t} {t} {a}.
(Element t, Storable a, Storable a, Num t, Num t, Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Matrix a, Vector a)
leftSVAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgesvd String
"leftSVC"

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

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

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

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


-- | 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 C -> (Vector C, Matrix C)
eigC = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> Ptr C
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> IO CInt)
-> String -> Matrix C -> (Vector C, Matrix C)
forall {t} {a} {a} {t} {t} {t} {t} {a}.
(Element t, Storable a, Storable a, Num t, Num t, Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> (Vector a, Matrix a)
eigAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgeev String
"eigC"

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

-- | 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 C -> Vector C
eigOnlyC = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> Ptr C
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> IO CInt)
-> String -> Matrix C -> Vector C
forall {t} {a} {t} {t} {t} {t} {t} {t} {t} {t} {a} {a}.
(Element t, Storable a, Num t, Num t, Num t, Num t, Num t, Num t,
 Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Vector a
eigOnlyAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgeev String
"eigOnlyC"

-- | 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 R -> (Vector C, Matrix C)
eigR Matrix R
m = (Vector C
s', Matrix C
v'')
    where (Vector C
s,Matrix R
v) = Matrix R -> (Vector C, Matrix R)
eigRaux Matrix R
m
          s' :: Vector C
s' = Vector C -> Vector C
forall {t}.
RealElement t =>
Vector (Complex t) -> Vector (Complex t)
fixeig1 Vector C
s
          v' :: [Vector R]
v' = Matrix R -> [Vector R]
forall t. Element t => Matrix t -> [Vector t]
toRows (Matrix R -> [Vector R]) -> Matrix R -> [Vector R]
forall a b. (a -> b) -> a -> b
$ Matrix R -> Matrix R
forall t. Matrix t -> Matrix t
trans Matrix R
v
          v'' :: Matrix C
v'' = [Vector C] -> Matrix C
forall t. Element t => [Vector t] -> Matrix t
fromColumns ([Vector C] -> Matrix C) -> [Vector C] -> Matrix C
forall a b. (a -> b) -> a -> b
$ [C] -> [Vector R] -> [Vector C]
forall {e} {a}.
(RealElement e, Num a, Eq a) =>
[Complex a] -> [Vector e] -> [Vector (Complex e)]
fixeig (Vector C -> [C]
forall a. Storable a => Vector a -> [a]
toList Vector C
s') [Vector R]
v'

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

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

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

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

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

-- | 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 R -> Matrix R -> (Vector C, Vector R, Matrix C)
eigG Matrix R
a Matrix R
b = (Vector C
alpha', Vector R
beta, Matrix C
v'')
  where
    (Vector C
alpha, Vector R
beta, Matrix R
v) = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> CInt
 -> Ptr C
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> IO CInt)
-> Matrix R -> Matrix R -> String -> (Vector C, Vector R, Matrix R)
forall {t} {t} {a} {a} {a} {t} {t} {t} {t} {a}.
(Element t, Element t, Storable a, Storable a, Storable a, Num t,
 Num t, Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> Matrix t -> Matrix t -> String -> (Vector a, Vector a, Matrix a)
eigGaux CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dggev Matrix R
a Matrix R
b String
"eigG"
    alpha' :: Vector C
alpha' = Vector C -> Vector C
forall {t}.
RealElement t =>
Vector (Complex t) -> Vector (Complex t)
fixeig1 Vector C
alpha
    v' :: [Vector R]
v' = Matrix R -> [Vector R]
forall t. Element t => Matrix t -> [Vector t]
toRows (Matrix R -> [Vector R]) -> Matrix R -> [Vector R]
forall a b. (a -> b) -> a -> b
$ Matrix R -> Matrix R
forall t. Matrix t -> Matrix t
trans Matrix R
v
    v'' :: Matrix C
v'' = [Vector C] -> Matrix C
forall t. Element t => [Vector t] -> Matrix t
fromColumns ([Vector C] -> Matrix C) -> [Vector C] -> Matrix C
forall a b. (a -> b) -> a -> b
$ [C] -> [Vector R] -> [Vector C]
forall {e} {a}.
(RealElement e, Ord a, Fractional a) =>
[Complex a] -> [Vector e] -> [Vector (Complex e)]
fixeigG (Vector C -> [C]
forall a. Storable a => Vector a -> [a]
toList Vector C
alpha') [Vector R]
v'

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

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

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

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

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

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

-- | 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 C -> Matrix C -> (Vector C, Vector C, Matrix C)
eigGC Matrix C
a Matrix C
b = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> Ptr C
 -> CInt
 -> Ptr C
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> IO CInt)
-> Matrix C -> Matrix C -> String -> (Vector C, Vector C, Matrix C)
forall {t} {t} {a} {a} {a} {t} {t} {t} {t} {a}.
(Element t, Element t, Storable a, Storable a, Storable a, Num t,
 Num t, Num t, Num t) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> Ptr a
 -> CInt
 -> Ptr a
 -> t
 -> t
 -> t
 -> t
 -> Ptr a
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> Matrix t -> Matrix t -> String -> (Vector a, Vector a, Matrix a)
eigGaux CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zggev Matrix C
a Matrix C
b String
"eigGC"

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

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

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

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

-- | 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 R -> (Vector R, Matrix R)
eigS Matrix R
m = (Vector R
s', Matrix R -> Matrix R
forall t. Element t => Matrix t -> Matrix t
fliprl Matrix R
v)
    where (Vector R
s,Matrix R
v) = Matrix R -> (Vector R, Matrix R)
eigS' Matrix R
m
          s' :: Vector R
s' = [R] -> Vector R
forall a. Storable a => [a] -> Vector a
fromList ([R] -> Vector R) -> (Vector R -> [R]) -> Vector R -> Vector R
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [R] -> [R]
forall a. [a] -> [a]
reverse ([R] -> [R]) -> (Vector R -> [R]) -> Vector R -> [R]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector R -> [R]
forall a. Storable a => Vector a -> [a]
toList (Vector R -> Vector R) -> Vector R -> Vector R
forall a b. (a -> b) -> a -> b
$  Vector R
s

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

-- | 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 C -> (Vector R, Matrix C)
eigH Matrix C
m = (Vector R
s', Matrix C -> Matrix C
forall t. Element t => Matrix t -> Matrix t
fliprl Matrix C
v)
  where
    (Vector R
s,Matrix C
v) = Matrix C -> (Vector R, Matrix C)
eigH' Matrix C
m
    s' :: Vector R
s' = [R] -> Vector R
forall a. Storable a => [a] -> Vector a
fromList ([R] -> Vector R) -> (Vector R -> [R]) -> Vector R -> Vector R
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [R] -> [R]
forall a. [a] -> [a]
reverse ([R] -> [R]) -> (Vector R -> [R]) -> Vector R -> [R]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector R -> [R]
forall a. Storable a => Vector a -> [a]
toList (Vector R -> Vector R) -> Vector R -> Vector R
forall a b. (a -> b) -> a -> b
$  Vector R
s

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


-- | 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 R -> Vector R
eigOnlyS = Vector R -> Vector R
vrev (Vector R -> Vector R)
-> (Matrix R -> Vector R) -> Matrix R -> Vector R
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vector R, Matrix R) -> Vector R
forall a b. (a, b) -> a
fst((Vector R, Matrix R) -> Vector R)
-> (Matrix R -> (Vector R, Matrix R)) -> Matrix R -> Vector R
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CInt -> Ptr R -> CInt -> CInt -> CInt -> CInt -> Ptr R -> IO CInt)
-> String -> Matrix R -> (Vector R, Matrix R)
forall {t} {a}.
(Element t, Storable a) =>
(CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Vector a, Matrix t)
eigSHAux (CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dsyev CInt
0) String
"eigS'"

-- | 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 C -> Vector R
eigOnlyH = Vector R -> Vector R
vrev (Vector R -> Vector R)
-> (Matrix C -> Vector R) -> Matrix C -> Vector R
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Vector R, Matrix C) -> Vector R
forall a b. (a, b) -> a
fst((Vector R, Matrix C) -> Vector R)
-> (Matrix C -> (Vector R, Matrix C)) -> Matrix C -> Vector R
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CInt -> Ptr R -> CInt -> CInt -> CInt -> CInt -> Ptr C -> IO CInt)
-> String -> Matrix C -> (Vector R, Matrix C)
forall {t} {a}.
(Element t, Storable a) =>
(CInt -> Ptr a -> CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> (Vector a, Matrix t)
eigSHAux (CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zheev CInt
0) String
"eigH'"

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

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

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

-- | 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 R -> Matrix R -> Matrix R
linearSolveR Matrix R
a Matrix R
b = (IO (Matrix R) -> IO (Matrix R))
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr R
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr R
    -> IO CInt)
-> String
-> Matrix R
-> Matrix R
-> Matrix R
forall {t} {t} {c}.
(Element t, Element t) =>
(IO (Matrix t) -> IO c)
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveSQAux IO (Matrix R) -> IO (Matrix R)
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgesv String
"linearSolveR" Matrix R
a Matrix R
b

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


-- | 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 C -> Matrix C -> Matrix C
linearSolveC Matrix C
a Matrix C
b = (IO (Matrix C) -> IO (Matrix C))
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr C
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr C
    -> IO CInt)
-> String
-> Matrix C
-> Matrix C
-> Matrix C
forall {t} {t} {c}.
(Element t, Element t) =>
(IO (Matrix t) -> IO c)
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveSQAux IO (Matrix C) -> IO (Matrix C)
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgesv String
"linearSolveC" Matrix C
a Matrix C
b

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

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


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

-- | 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 R -> Matrix R -> Matrix R
cholSolveR Matrix R
a Matrix R
b = (IO (Matrix R) -> IO (Matrix R))
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr R
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr R
    -> IO CInt)
-> String
-> Matrix R
-> Matrix R
-> Matrix R
forall {t} {t} {c}.
(Element t, Storable t) =>
(IO (Matrix t) -> IO c)
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveSQAux2 IO (Matrix R) -> IO (Matrix R)
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dpotrs String
"cholSolveR" (Matrix R -> Matrix R
forall t. Element t => Matrix t -> Matrix t
fmat Matrix R
a) Matrix R
b

-- | 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 C -> Matrix C -> Matrix C
cholSolveC Matrix C
a Matrix C
b = (IO (Matrix C) -> IO (Matrix C))
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr C
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr C
    -> IO CInt)
-> String
-> Matrix C
-> Matrix C
-> Matrix C
forall {t} {t} {c}.
(Element t, Storable t) =>
(IO (Matrix t) -> IO c)
-> (CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> CInt
    -> CInt
    -> CInt
    -> CInt
    -> Ptr t
    -> IO CInt)
-> String
-> Matrix t
-> Matrix t
-> c
linearSolveSQAux2 IO (Matrix C) -> IO (Matrix C)
forall a. a -> a
id CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zpotrs String
"cholSolveC" (Matrix C -> Matrix C
forall t. Element t => Matrix t -> Matrix t
fmat Matrix C
a) Matrix C
b

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


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

data UpLo = Lower | Upper

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

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

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

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

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

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

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

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

-- | 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 R -> Matrix R -> Matrix R
linearSolveLSR Matrix R
a Matrix R
b = (Int, Int) -> (Int, Int) -> Matrix R -> Matrix R
forall a.
Element a =>
(Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
subMatrix (Int
0,Int
0) (Matrix R -> Int
forall t. Matrix t -> Int
cols Matrix R
a, Matrix R -> Int
forall t. Matrix t -> Int
cols Matrix R
b) (Matrix R -> Matrix R) -> Matrix R -> Matrix R
forall a b. (a -> b) -> a -> b
$
                     (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> IO CInt)
-> String -> Matrix R -> Matrix R -> Matrix R
forall {t} {a}.
(Element t, Element a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Matrix a -> Matrix a
linearSolveAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgels String
"linearSolverLSR" Matrix R
a Matrix R
b

-- | 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 C -> Matrix C -> Matrix C
linearSolveLSC Matrix C
a Matrix C
b = (Int, Int) -> (Int, Int) -> Matrix C -> Matrix C
forall a.
Element a =>
(Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
subMatrix (Int
0,Int
0) (Matrix C -> Int
forall t. Matrix t -> Int
cols Matrix C
a, Matrix C -> Int
forall t. Matrix t -> Int
cols Matrix C
b) (Matrix C -> Matrix C) -> Matrix C -> Matrix C
forall a b. (a -> b) -> a -> b
$
                     (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> IO CInt)
-> String -> Matrix C -> Matrix C -> Matrix C
forall {t} {a}.
(Element t, Element a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Matrix a -> Matrix a
linearSolveAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgels String
"linearSolveLSC" Matrix C
a Matrix C
b

-- | 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 R -> Matrix R -> Matrix R -> Matrix R
linearSolveSVDR (Just R
rcond) Matrix R
a Matrix R
b = (Int, Int) -> (Int, Int) -> Matrix R -> Matrix R
forall a.
Element a =>
(Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
subMatrix (Int
0,Int
0) (Matrix R -> Int
forall t. Matrix t -> Int
cols Matrix R
a, Matrix R -> Int
forall t. Matrix t -> Int
cols Matrix R
b) (Matrix R -> Matrix R) -> Matrix R -> Matrix R
forall a b. (a -> b) -> a -> b
$
                                   (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> IO CInt)
-> String -> Matrix R -> Matrix R -> Matrix R
forall {t} {a}.
(Element t, Element a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Matrix a -> Matrix a
linearSolveAux (R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgelss R
rcond) String
"linearSolveSVDR" Matrix R
a Matrix R
b
linearSolveSVDR Maybe R
Nothing Matrix R
a Matrix R
b = Maybe R -> Matrix R -> Matrix R -> Matrix R
linearSolveSVDR (R -> Maybe R
forall a. a -> Maybe a
Just (-R
1)) Matrix R
a Matrix R
b

-- | 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 R -> Matrix C -> Matrix C -> Matrix C
linearSolveSVDC (Just R
rcond) Matrix C
a Matrix C
b = (Int, Int) -> (Int, Int) -> Matrix C -> Matrix C
forall a.
Element a =>
(Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
subMatrix (Int
0,Int
0) (Matrix C -> Int
forall t. Matrix t -> Int
cols Matrix C
a, Matrix C -> Int
forall t. Matrix t -> Int
cols Matrix C
b) (Matrix C -> Matrix C) -> Matrix C -> Matrix C
forall a b. (a -> b) -> a -> b
$
                                   (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> IO CInt)
-> String -> Matrix C -> Matrix C -> Matrix C
forall {t} {a}.
(Element t, Element a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr a
 -> IO CInt)
-> String -> Matrix t -> Matrix a -> Matrix a
linearSolveAux (R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgelss R
rcond) String
"linearSolveSVDC" Matrix C
a Matrix C
b
linearSolveSVDC Maybe R
Nothing Matrix C
a Matrix C
b = Maybe R -> Matrix C -> Matrix C -> Matrix C
linearSolveSVDC (R -> Maybe R
forall a. a -> Maybe a
Just (-R
1)) Matrix C
a Matrix C
b

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

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

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

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

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

-- | 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 C -> Maybe (Matrix C)
mbCholH = IO (Maybe (Matrix C)) -> Maybe (Matrix C)
forall a. IO a -> a
unsafePerformIO (IO (Maybe (Matrix C)) -> Maybe (Matrix C))
-> (Matrix C -> IO (Maybe (Matrix C)))
-> Matrix C
-> Maybe (Matrix C)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IO (Matrix C) -> IO (Maybe (Matrix C))
forall x. IO x -> IO (Maybe x)
mbCatch (IO (Matrix C) -> IO (Maybe (Matrix C)))
-> (Matrix C -> IO (Matrix C)) -> Matrix C -> IO (Maybe (Matrix C))
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CInt -> CInt -> CInt -> CInt -> Ptr C -> IO CInt)
-> String -> Matrix C -> IO (Matrix C)
forall {t}.
Element t =>
(CInt -> CInt -> CInt -> CInt -> Ptr t -> IO CInt)
-> String -> Matrix t -> IO (Matrix t)
cholAux CInt -> CInt -> CInt -> CInt -> Ptr C -> IO CInt
zpotrf String
"cholH"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

-- | 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 R -> [Int] -> Matrix R -> Matrix R
lusR Matrix R
a [Int]
piv Matrix R
b = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> IO CInt)
-> String -> Matrix R -> [Int] -> Matrix R -> Matrix R
forall {t} {t} {a}.
(Element t, Storable t, Integral a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> IO CInt)
-> String -> Matrix t -> [a] -> Matrix t -> Matrix t
lusAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dgetrs String
"lusR" (Matrix R -> Matrix R
forall t. Element t => Matrix t -> Matrix t
fmat Matrix R
a) [Int]
piv Matrix R
b

-- | 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 C -> [Int] -> Matrix C -> Matrix C
lusC Matrix C
a [Int]
piv Matrix C
b = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> IO CInt)
-> String -> Matrix C -> [Int] -> Matrix C -> Matrix C
forall {t} {t} {a}.
(Element t, Storable t, Integral a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> IO CInt)
-> String -> Matrix t -> [a] -> Matrix t -> Matrix t
lusAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zgetrs String
"lusC" (Matrix C -> Matrix C
forall t. Element t => Matrix t -> Matrix t
fmat Matrix C
a) [Int]
piv Matrix C
b

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

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

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

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

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

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

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

-- | 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 R -> [Int] -> Matrix R -> Matrix R
ldlsR Matrix R
a [Int]
piv Matrix R
b = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr R
 -> IO CInt)
-> String -> Matrix R -> [Int] -> Matrix R -> Matrix R
forall {t} {t} {a}.
(Element t, Storable t, Integral a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> IO CInt)
-> String -> Matrix t -> [a] -> Matrix t -> Matrix t
lusAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr R
-> IO CInt
dsytrs String
"ldlsR" (Matrix R -> Matrix R
forall t. Element t => Matrix t -> Matrix t
fmat Matrix R
a) [Int]
piv Matrix R
b

-- | 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 C -> [Int] -> Matrix C -> Matrix C
ldlsC Matrix C
a [Int]
piv Matrix C
b = (CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr C
 -> IO CInt)
-> String -> Matrix C -> [Int] -> Matrix C -> Matrix C
forall {t} {t} {a}.
(Element t, Storable t, Integral a) =>
(CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> CInt
 -> Ptr R
 -> CInt
 -> CInt
 -> CInt
 -> CInt
 -> Ptr t
 -> IO CInt)
-> String -> Matrix t -> [a] -> Matrix t -> Matrix t
lusAux CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> CInt
-> Ptr R
-> CInt
-> CInt
-> CInt
-> CInt
-> Ptr C
-> IO CInt
zsytrs String
"ldlsC" (Matrix C -> Matrix C
forall t. Element t => Matrix t -> Matrix t
fmat Matrix C
a) [Int]
piv Matrix C
b