{-# LANGUAGE TypeFamilies #-}
module Numeric.LAPACK.Matrix.Hermitian.Linear (
   solve,
   inverse,
   determinant,
   ) where

import qualified Numeric.LAPACK.Matrix.Symmetric.Private as Symmetric
import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import Numeric.LAPACK.Matrix.Hermitian.Basic (Hermitian)
import Numeric.LAPACK.Matrix.Hermitian.Private (Determinant(..))
import Numeric.LAPACK.Matrix.Modifier (Conjugation(Conjugated))
import Numeric.LAPACK.Matrix.Private (Full)
import Numeric.LAPACK.Scalar (RealOf, absoluteSquared)
import Numeric.LAPACK.Private (realPtr)

import qualified Numeric.Netlib.Class as Class

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 System.IO.Unsafe (unsafePerformIO)

import Foreign.Ptr (Ptr)
import Foreign.Storable (peek)


solve ::
   (Extent.C vert, Extent.C horiz,
    Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>
   Hermitian sh a -> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a
solve :: Hermitian sh a
-> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a
solve (Array (MatrixShape.Hermitian Order
orderA sh
shA) ForeignPtr a
a) =
   String
-> Conjugation
-> Order
-> sh
-> ForeignPtr a
-> Full vert horiz sh nrhs a
-> Full vert horiz sh nrhs a
forall vert horiz width height a.
(C vert, C horiz, C width, C height, Eq height, Floating a) =>
String
-> Conjugation
-> Order
-> height
-> ForeignPtr a
-> Full vert horiz height width a
-> Full vert horiz height width a
Symmetric.solve String
"Hermitian.solve" Conjugation
Conjugated Order
orderA sh
shA ForeignPtr a
a


inverse :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> Hermitian sh a
inverse :: Hermitian sh a -> Hermitian sh a
inverse (Array shape :: Hermitian sh
shape@(MatrixShape.Hermitian Order
order sh
sh) ForeignPtr a
a) =
   Hermitian sh -> (Int -> Ptr a -> IO ()) -> Hermitian sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize Hermitian sh
shape ((Int -> Ptr a -> IO ()) -> Hermitian sh a)
-> (Int -> Ptr a -> IO ()) -> Hermitian sh a
forall a b. (a -> b) -> a -> b
$
      Conjugation
-> Order -> Int -> ForeignPtr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
Conjugation
-> Order -> Int -> ForeignPtr a -> Int -> Ptr a -> IO ()
Symmetric.inverse Conjugation
Conjugated Order
order (sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh) ForeignPtr a
a


determinant :: (Shape.C sh, Class.Floating a) => Hermitian sh a -> RealOf a
determinant :: Hermitian sh a -> RealOf a
determinant =
   Determinant (Array (Hermitian sh)) a -> Hermitian sh a -> RealOf a
forall (f :: * -> *) a. Determinant f a -> f a -> RealOf a
getDeterminant (Determinant (Array (Hermitian sh)) a
 -> Hermitian sh a -> RealOf a)
-> Determinant (Array (Hermitian sh)) a
-> Hermitian sh a
-> RealOf a
forall a b. (a -> b) -> a -> b
$
   Determinant (Array (Hermitian sh)) Float
-> Determinant (Array (Hermitian sh)) Double
-> Determinant (Array (Hermitian sh)) (Complex Float)
-> Determinant (Array (Hermitian sh)) (Complex Double)
-> Determinant (Array (Hermitian 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 -> RealOf Float)
-> Determinant (Array (Hermitian sh)) Float
forall (f :: * -> *) a. (f a -> RealOf a) -> Determinant f a
Determinant Array (Hermitian sh) Float -> RealOf Float
forall sh a ar.
(C sh, Floating a, RealOf a ~ ar, Real ar) =>
Hermitian sh a -> ar
determinantAux) ((Array (Hermitian sh) Double -> RealOf Double)
-> Determinant (Array (Hermitian sh)) Double
forall (f :: * -> *) a. (f a -> RealOf a) -> Determinant f a
Determinant Array (Hermitian sh) Double -> RealOf Double
forall sh a ar.
(C sh, Floating a, RealOf a ~ ar, Real ar) =>
Hermitian sh a -> ar
determinantAux)
      ((Array (Hermitian sh) (Complex Float) -> RealOf (Complex Float))
-> Determinant (Array (Hermitian sh)) (Complex Float)
forall (f :: * -> *) a. (f a -> RealOf a) -> Determinant f a
Determinant Array (Hermitian sh) (Complex Float) -> RealOf (Complex Float)
forall sh a ar.
(C sh, Floating a, RealOf a ~ ar, Real ar) =>
Hermitian sh a -> ar
determinantAux) ((Array (Hermitian sh) (Complex Double) -> RealOf (Complex Double))
-> Determinant (Array (Hermitian sh)) (Complex Double)
forall (f :: * -> *) a. (f a -> RealOf a) -> Determinant f a
Determinant Array (Hermitian sh) (Complex Double) -> RealOf (Complex Double)
forall sh a ar.
(C sh, Floating a, RealOf a ~ ar, Real ar) =>
Hermitian sh a -> ar
determinantAux)

determinantAux ::
   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Hermitian sh a -> ar
determinantAux :: Hermitian sh a -> ar
determinantAux (Array (MatrixShape.Hermitian Order
order sh
sh) ForeignPtr a
a) =
   IO ar -> ar
forall a. IO a -> a
unsafePerformIO (IO ar -> ar) -> IO ar -> ar
forall a b. (a -> b) -> a -> b
$
      Conjugation
-> ((Ptr a, Maybe (Ptr a, Ptr a)) -> IO ar)
-> Order
-> Int
-> ForeignPtr a
-> IO ar
forall a ar.
(Floating a, Floating ar) =>
Conjugation
-> ((Ptr a, Maybe (Ptr a, Ptr a)) -> IO ar)
-> Order
-> Int
-> ForeignPtr a
-> IO ar
Symmetric.determinant Conjugation
Conjugated
         (Ptr a, Maybe (Ptr a, Ptr a)) -> IO ar
forall a ar.
(Floating a, RealOf a ~ ar, Real ar) =>
(Ptr a, Maybe (Ptr a, Ptr a)) -> IO ar
peekBlockDeterminant Order
order (sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh) ForeignPtr a
a

peekBlockDeterminant ::
   (Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   (Ptr a, Maybe (Ptr a, Ptr a)) -> IO ar
peekBlockDeterminant :: (Ptr a, Maybe (Ptr a, Ptr a)) -> IO ar
peekBlockDeterminant (Ptr a
a0Ptr,Maybe (Ptr a, Ptr a)
ext) = do
   let peekReal :: Ptr a -> IO (RealOf a)
peekReal = Ptr (RealOf a) -> IO (RealOf a)
forall a. Storable a => Ptr a -> IO a
peek (Ptr (RealOf a) -> IO (RealOf a))
-> (Ptr a -> Ptr (RealOf a)) -> Ptr a -> IO (RealOf a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Ptr a -> Ptr (RealOf a)
forall a. Ptr a -> Ptr (RealOf a)
realPtr
   ar
a0 <- Ptr a -> IO (RealOf a)
peekReal Ptr a
a0Ptr
   case Maybe (Ptr a, Ptr a)
ext of
      Maybe (Ptr a, Ptr a)
Nothing -> ar -> IO ar
forall (m :: * -> *) a. Monad m => a -> m a
return ar
a0
      Just (Ptr a
a1Ptr,Ptr a
bPtr) -> do
         ar
a1 <- Ptr a -> IO (RealOf a)
peekReal Ptr a
a1Ptr
         a
b <- Ptr a -> IO a
forall a. Storable a => Ptr a -> IO a
peek Ptr a
bPtr
         ar -> IO ar
forall (m :: * -> *) a. Monad m => a -> m a
return (ar
a0ar -> ar -> ar
forall a. Num a => a -> a -> a
*ar
a1 ar -> ar -> ar
forall a. Num a => a -> a -> a
- a -> RealOf a
forall a. Floating a => a -> RealOf a
absoluteSquared a
b)