{-# LANGUAGE TypeFamilies #-}
module Numeric.LAPACK.Matrix.Hermitian.Eigen (
   values,
   decompose,
   ) where

import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
import Numeric.LAPACK.Matrix.Hermitian.Basic (Hermitian)
import Numeric.LAPACK.Matrix.Hermitian.Private (TakeDiagonal(..))
import Numeric.LAPACK.Matrix.Square.Basic (Square)
import Numeric.LAPACK.Matrix.Shape.Private (Order(ColumnMajor), uploFromOrder)
import Numeric.LAPACK.Matrix.Modifier (conjugatedOnRowMajor)
import Numeric.LAPACK.Vector (Vector)
import Numeric.LAPACK.Scalar (RealOf)
import Numeric.LAPACK.Private
         (copyToTemp, copyCondConjugateToTemp, withInfo, eigenMsg)

import qualified Numeric.LAPACK.FFI.Complex as LapackComplex
import qualified Numeric.LAPACK.FFI.Real as LapackReal
import qualified Numeric.Netlib.Utility as Call
import qualified Numeric.Netlib.Class as Class

import qualified Data.Array.Comfort.Storable.Unchecked.Monadic as ArrayIO
import qualified Data.Array.Comfort.Storable.Unchecked as Array
import qualified Data.Array.Comfort.Shape as Shape
import Data.Array.Comfort.Storable.Unchecked (Array(Array))
import Data.Array.Comfort.Shape (triangleSize)

import Foreign.C.Types (CInt, CChar)
import Foreign.Ptr (Ptr, nullPtr)
import Foreign.Storable (Storable)

import Control.Monad.Trans.Cont (evalContT)
import Control.Monad.IO.Class (liftIO)

import Data.Complex (Complex)


values ::
   (Shape.C sh, Class.Floating a) =>
   Hermitian sh a -> Vector sh (RealOf a)
values :: Hermitian sh a -> Vector sh (RealOf a)
values =
   TakeDiagonal (Array (Hermitian sh)) (Array sh) a
-> Hermitian sh a -> Vector sh (RealOf a)
forall (f :: * -> *) (g :: * -> *) a.
TakeDiagonal f g a -> f a -> g (RealOf a)
runTakeDiagonal (TakeDiagonal (Array (Hermitian sh)) (Array sh) a
 -> Hermitian sh a -> Vector sh (RealOf a))
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) a
-> Hermitian sh a
-> Vector sh (RealOf a)
forall a b. (a -> b) -> a -> b
$
   TakeDiagonal (Array (Hermitian sh)) (Array sh) Float
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) Double
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) (Complex Float)
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) (Complex Double)
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      ((Array (Hermitian sh) Float -> Array sh (RealOf Float))
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) Float
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (RealOf a)) -> TakeDiagonal f g a
TakeDiagonal Array (Hermitian sh) Float -> Array sh (RealOf Float)
forall sh a ar.
(C sh, Floating a, RealOf a ~ ar, Storable ar) =>
Hermitian sh a -> Vector sh ar
valuesAux) ((Array (Hermitian sh) Double -> Array sh (RealOf Double))
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) Double
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (RealOf a)) -> TakeDiagonal f g a
TakeDiagonal Array (Hermitian sh) Double -> Array sh (RealOf Double)
forall sh a ar.
(C sh, Floating a, RealOf a ~ ar, Storable ar) =>
Hermitian sh a -> Vector sh ar
valuesAux)
      ((Array (Hermitian sh) (Complex Float)
 -> Array sh (RealOf (Complex Float)))
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) (Complex Float)
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (RealOf a)) -> TakeDiagonal f g a
TakeDiagonal Array (Hermitian sh) (Complex Float)
-> Array sh (RealOf (Complex Float))
forall sh a ar.
(C sh, Floating a, RealOf a ~ ar, Storable ar) =>
Hermitian sh a -> Vector sh ar
valuesAux) ((Array (Hermitian sh) (Complex Double)
 -> Array sh (RealOf (Complex Double)))
-> TakeDiagonal (Array (Hermitian sh)) (Array sh) (Complex Double)
forall (f :: * -> *) (g :: * -> *) a.
(f a -> g (RealOf a)) -> TakeDiagonal f g a
TakeDiagonal Array (Hermitian sh) (Complex Double)
-> Array sh (RealOf (Complex Double))
forall sh a ar.
(C sh, Floating a, RealOf a ~ ar, Storable ar) =>
Hermitian sh a -> Vector sh ar
valuesAux)

valuesAux ::
   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>
   Hermitian sh a -> Vector sh ar
valuesAux :: Hermitian sh a -> Vector sh ar
valuesAux (Array (MatrixShape.Hermitian Order
order sh
size) ForeignPtr a
a) =
   sh -> (Int -> Ptr ar -> IO ()) -> Vector sh ar
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize sh
size ((Int -> Ptr ar -> IO ()) -> Vector sh ar)
-> (Int -> Ptr ar -> IO ()) -> Vector sh ar
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr ar
wPtr ->
   ContT () IO () -> IO ()
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT () IO () -> IO ()) -> ContT () IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
      Ptr CChar
jobzPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'N'
      Ptr CChar
uploPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char (Char -> FortranIO () (Ptr CChar))
-> Char -> FortranIO () (Ptr CChar)
forall a b. (a -> b) -> a -> b
$ Order -> Char
uploFromOrder Order
order
      Ptr a
aPtr <- Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r. Storable a => Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyToTemp (Int -> Int
triangleSize Int
n) ForeignPtr a
a
      let zPtr :: Ptr a
zPtr = Ptr a
forall a. Ptr a
nullPtr
      Ptr CInt
ldzPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      IO () -> ContT () IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ String -> String -> (Ptr CInt -> IO ()) -> IO ()
withInfo String
eigenMsg String
"hpev" ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
         HPEV_ (RealOf a) a
forall a. Floating a => HPEV_ (RealOf a) a
hpev Ptr CChar
jobzPtr Ptr CChar
uploPtr Int
n Ptr a
aPtr Ptr ar
Ptr (RealOf a)
wPtr Ptr a
forall a. Ptr a
zPtr Ptr CInt
ldzPtr


decompose ::
   (Shape.C sh, Class.Floating a) =>
   Hermitian sh a -> (Square sh a, Vector sh (RealOf a))
decompose :: Hermitian sh a -> (Square sh a, Vector sh (RealOf a))
decompose =
   Decompose sh a
-> Hermitian sh a -> (Square sh a, Vector sh (RealOf a))
forall sh a. Decompose sh a -> Decompose_ sh a
getDecompose (Decompose sh a
 -> Hermitian sh a -> (Square sh a, Vector sh (RealOf a)))
-> Decompose sh a
-> Hermitian sh a
-> (Square sh a, Vector sh (RealOf a))
forall a b. (a -> b) -> a -> b
$
   Decompose sh Float
-> Decompose sh Double
-> Decompose sh (Complex Float)
-> Decompose sh (Complex Double)
-> Decompose sh a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (Decompose_ sh Float -> Decompose sh Float
forall sh a. Decompose_ sh a -> Decompose sh a
Decompose Decompose_ sh Float
forall sh a ar.
(C sh, Floating a, RealOf a ~ ar, Storable ar) =>
Decompose_ sh a
decomposeAux) (Decompose_ sh Double -> Decompose sh Double
forall sh a. Decompose_ sh a -> Decompose sh a
Decompose Decompose_ sh Double
forall sh a ar.
(C sh, Floating a, RealOf a ~ ar, Storable ar) =>
Decompose_ sh a
decomposeAux)
      (Decompose_ sh (Complex Float) -> Decompose sh (Complex Float)
forall sh a. Decompose_ sh a -> Decompose sh a
Decompose Decompose_ sh (Complex Float)
forall sh a ar.
(C sh, Floating a, RealOf a ~ ar, Storable ar) =>
Decompose_ sh a
decomposeAux) (Decompose_ sh (Complex Double) -> Decompose sh (Complex Double)
forall sh a. Decompose_ sh a -> Decompose sh a
Decompose Decompose_ sh (Complex Double)
forall sh a ar.
(C sh, Floating a, RealOf a ~ ar, Storable ar) =>
Decompose_ sh a
decomposeAux)

type Decompose_ sh a = Hermitian sh a -> (Square sh a, Vector sh (RealOf a))

newtype Decompose sh a = Decompose {Decompose sh a -> Decompose_ sh a
getDecompose :: Decompose_ sh a}

decomposeAux ::
   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Storable ar) =>
   Decompose_ sh a
decomposeAux :: Decompose_ sh a
decomposeAux (Array (MatrixShape.Hermitian Order
order sh
size) ForeignPtr a
a) =
   Full Small Small sh sh
-> (Int -> Ptr a -> IO (Array sh ar))
-> (Array (Full Small Small sh sh) a, Array sh ar)
forall sh a b.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO b) -> (Array sh a, b)
Array.unsafeCreateWithSizeAndResult (Order -> sh -> Full Small Small sh sh
forall sh. Order -> sh -> Square sh
MatrixShape.square Order
ColumnMajor sh
size) ((Int -> Ptr a -> IO (Array sh ar))
 -> (Array (Full Small Small sh sh) a, Array sh ar))
-> (Int -> Ptr a -> IO (Array sh ar))
-> (Array (Full Small Small sh sh) a, Array sh ar)
forall a b. (a -> b) -> a -> b
$
      \Int
_ Ptr a
zPtr ->
   sh -> (Int -> Ptr ar -> IO ()) -> IO (Array sh ar)
forall (m :: * -> *) sh a.
(PrimMonad m, C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> m (Array sh a)
ArrayIO.unsafeCreateWithSize sh
size ((Int -> Ptr ar -> IO ()) -> IO (Array sh ar))
-> (Int -> Ptr ar -> IO ()) -> IO (Array sh ar)
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr ar
wPtr ->
   ContT () IO () -> IO ()
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT () IO () -> IO ()) -> ContT () IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
      Ptr CChar
jobzPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'V'
      Ptr CChar
uploPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char (Char -> FortranIO () (Ptr CChar))
-> Char -> FortranIO () (Ptr CChar)
forall a b. (a -> b) -> a -> b
$ Order -> Char
uploFromOrder Order
order
      Ptr a
aPtr <-
         Conjugation -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r.
Floating a =>
Conjugation -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyCondConjugateToTemp (Order -> Conjugation
conjugatedOnRowMajor Order
order) (Int -> Int
triangleSize Int
n) ForeignPtr a
a
      Ptr CInt
ldzPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      IO () -> ContT () IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ String -> String -> (Ptr CInt -> IO ()) -> IO ()
withInfo String
eigenMsg String
"hpev" ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
         HPEV_ (RealOf a) a
forall a. Floating a => HPEV_ (RealOf a) a
hpev Ptr CChar
jobzPtr Ptr CChar
uploPtr Int
n Ptr a
aPtr Ptr ar
Ptr (RealOf a)
wPtr Ptr a
zPtr Ptr CInt
ldzPtr


type HPEV_ ar a =
   Ptr CChar -> Ptr CChar -> Int -> Ptr a -> Ptr ar ->
   Ptr a -> Ptr CInt -> Ptr CInt -> IO ()

newtype HPEV a = HPEV {HPEV a -> HPEV_ (RealOf a) a
getHPEV :: HPEV_ (RealOf a) a}

hpev :: Class.Floating a => HPEV_ (RealOf a) a
hpev :: HPEV_ (RealOf a) a
hpev =
   HPEV a -> HPEV_ (RealOf a) a
forall a. HPEV a -> HPEV_ (RealOf a) a
getHPEV (HPEV a -> HPEV_ (RealOf a) a) -> HPEV a -> HPEV_ (RealOf a) a
forall a b. (a -> b) -> a -> b
$
   HPEV Float
-> HPEV Double
-> HPEV (Complex Float)
-> HPEV (Complex Double)
-> HPEV a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      (HPEV_ (RealOf Float) Float -> HPEV Float
forall a. HPEV_ (RealOf a) a -> HPEV a
HPEV HPEV_ (RealOf Float) Float
forall a. Real a => HPEV_ a a
spevReal) (HPEV_ (RealOf Double) Double -> HPEV Double
forall a. HPEV_ (RealOf a) a -> HPEV a
HPEV HPEV_ (RealOf Double) Double
forall a. Real a => HPEV_ a a
spevReal) (HPEV_ (RealOf (Complex Float)) (Complex Float)
-> HPEV (Complex Float)
forall a. HPEV_ (RealOf a) a -> HPEV a
HPEV HPEV_ (RealOf (Complex Float)) (Complex Float)
forall a. Real a => HPEV_ a (Complex a)
hpevComplex) (HPEV_ (RealOf (Complex Double)) (Complex Double)
-> HPEV (Complex Double)
forall a. HPEV_ (RealOf a) a -> HPEV a
HPEV HPEV_ (RealOf (Complex Double)) (Complex Double)
forall a. Real a => HPEV_ a (Complex a)
hpevComplex)

spevReal :: Class.Real a => HPEV_ a a
spevReal :: HPEV_ a a
spevReal Ptr CChar
jobzPtr Ptr CChar
uploPtr Int
n Ptr a
apPtr Ptr a
wPtr Ptr a
zPtr Ptr CInt
ldzPtr Ptr CInt
infoPtr =
   ContT () IO () -> IO ()
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT () IO () -> IO ()) -> ContT () IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr a
workPtr <- Int -> FortranIO () (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray (Int
3Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n)
      IO () -> ContT () IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> IO ()
forall a.
Real a =>
Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> IO ()
LapackReal.spev
            Ptr CChar
jobzPtr Ptr CChar
uploPtr Ptr CInt
nPtr Ptr a
apPtr Ptr a
wPtr Ptr a
zPtr Ptr CInt
ldzPtr Ptr a
workPtr Ptr CInt
infoPtr

hpevComplex :: Class.Real a => HPEV_ a (Complex a)
hpevComplex :: HPEV_ a (Complex a)
hpevComplex Ptr CChar
jobzPtr Ptr CChar
uploPtr Int
n Ptr (Complex a)
apPtr Ptr a
wPtr Ptr (Complex a)
zPtr Ptr CInt
ldzPtr Ptr CInt
infoPtr =
   ContT () IO () -> IO ()
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT () IO () -> IO ()) -> ContT () IO () -> IO ()
forall a b. (a -> b) -> a -> b
$ do
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr (Complex a)
workPtr <- Int -> FortranIO () (Ptr (Complex a))
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
1 (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1))
      Ptr a
rworkPtr <- Int -> FortranIO () (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
1 (Int
3Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2))
      IO () -> ContT () IO ()
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$
         Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr a
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr a
-> Ptr CInt
-> IO ()
forall a.
Real a =>
Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr a
-> Ptr (Complex a)
-> Ptr CInt
-> Ptr (Complex a)
-> Ptr a
-> Ptr CInt
-> IO ()
LapackComplex.hpev
            Ptr CChar
jobzPtr Ptr CChar
uploPtr Ptr CInt
nPtr Ptr (Complex a)
apPtr Ptr a
wPtr Ptr (Complex a)
zPtr Ptr CInt
ldzPtr Ptr (Complex a)
workPtr Ptr a
rworkPtr Ptr CInt
infoPtr