Linear.hs

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

import qualified Numeric.LAPACK.Matrix.Symmetric.Unified as Symmetric
import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
import Numeric.LAPACK.Matrix.Hermitian.Basic (HermitianP)
import Numeric.LAPACK.Matrix.Hermitian.Private (Determinant(..))
import Numeric.LAPACK.Scalar (RealOf, absoluteSquared)
import Numeric.LAPACK.Private (realPtr)

import qualified Numeric.Netlib.Class as Class

import qualified Data.Array.Comfort.Shape as Shape

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


determinant ::
   (Layout.Packing pack, Shape.C sh, Class.Floating a) =>
   HermitianP pack sh a -> RealOf a
determinant =
   getDeterminant $
   Class.switchFloating
      (Determinant $ Symmetric.determinant peekBlockDeterminant)
      (Determinant $ Symmetric.determinant peekBlockDeterminant)
      (Determinant $ Symmetric.determinant peekBlockDeterminant)
      (Determinant $ Symmetric.determinant peekBlockDeterminant)

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