{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE GADTs #-}
module Numeric.LAPACK.Matrix.HermitianPositiveDefinite.Linear (
   solve,
   solveDecomposed,
   inverse,
   decompose,
   determinant,
   ) where

import qualified Numeric.LAPACK.Matrix.Symmetric.Unified as Symmetric
import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import qualified Numeric.LAPACK.Vector as Vector
import Numeric.LAPACK.Matrix.Hermitian.Private (Determinant(..))
import Numeric.LAPACK.Matrix.Mosaic.Private
         (withPackingLinear, label, applyFuncPair, triArg, copyTriangleToTemp)
import Numeric.LAPACK.Matrix.Mosaic.Basic (takeDiagonal)
import Numeric.LAPACK.Matrix.Layout.Private (uploFromOrder)
import Numeric.BLAS.Matrix.Modifier (Conjugation(Conjugated))
import Numeric.LAPACK.Matrix.Private (Full)
import Numeric.LAPACK.Linear.Private (solver)
import Numeric.LAPACK.Scalar (RealOf, realPart, zero)
import Numeric.LAPACK.Private
         (copySubTrapezoid, copyBlock, fill, rankMsg, definiteMsg)

import qualified Numeric.LAPACK.FFI.Generic as LapackGen
import qualified Numeric.Netlib.Utility as Call
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 Foreign.ForeignPtr (withForeignPtr)

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


type Hermitian pack sh = Array (Layout.HermitianP pack sh)
type Upper pack sh = Array (Layout.UpperTriangularP pack sh)

solve ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>
   Hermitian pack sh a ->
   Full meas vert horiz sh nrhs a -> Full meas vert horiz sh nrhs a
solve :: forall meas vert horiz sh nrhs a pack.
(Measure meas, C vert, C horiz, C sh, Eq sh, C nrhs, Floating a) =>
Hermitian pack sh a
-> Full meas vert horiz sh nrhs a -> Full meas vert horiz sh nrhs a
solve (Array shape :: HermitianP pack sh
shape@(Layout.Mosaic PackingSingleton pack
pack MirrorSingleton ConjugateMirror
_mirror UpLoSingleton Upper
_upper Order
orderA sh
shA) ForeignPtr a
a) =
   String
-> sh
-> (Int
    -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr CInt -> ContT () IO ())
-> Full meas vert horiz sh nrhs a
-> Full meas vert horiz sh nrhs a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width, Eq height,
 Floating a) =>
String
-> height
-> (Int
    -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr CInt -> ContT () IO ())
-> Full meas vert horiz height width a
-> Full meas vert horiz height width a
solver String
"Hermitian.solve" sh
shA ((Int
  -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr CInt -> ContT () IO ())
 -> Full meas vert horiz sh nrhs a
 -> Full meas vert horiz sh nrhs a)
-> (Int
    -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr CInt -> ContT () IO ())
-> Full meas vert horiz sh nrhs a
-> Full meas vert horiz sh nrhs a
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr CInt
nPtr Ptr CInt
nrhsPtr Ptr a
xPtr Ptr CInt
ldxPtr -> do
      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
orderA
      Ptr a
aPtr <- Conjugation -> Order -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r.
Floating a =>
Conjugation -> Order -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyTriangleToTemp Conjugation
Conjugated Order
orderA (HermitianP pack sh -> Int
forall sh. C sh => sh -> Int
Shape.size HermitianP pack sh
shape) ForeignPtr a
a
      String
-> PackingSingleton pack
-> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
-> ContT () IO ()
forall func pack r.
(func ~ (Ptr CInt -> IO ())) =>
String
-> PackingSingleton pack
-> Labelled2 r String func func
-> ContT r IO ()
withPackingLinear String
definiteMsg PackingSingleton pack
pack (Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
 -> ContT () IO ())
-> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
-> ContT () IO ()
forall a b. (a -> b) -> a -> b
$
         Labelled
  (FuncCont
     (Ptr CChar
      -> Ptr CInt
      -> Ptr CInt
      -> TriArg a
      -> Ptr a
      -> Ptr CInt
      -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
  String
  (FuncPacked
     (Ptr CChar
      -> Ptr CInt
      -> Ptr CInt
      -> TriArg a
      -> Ptr a
      -> Ptr CInt
      -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
-> Labelled
     (FuncCont
        (Ptr CChar
         -> Ptr CInt
         -> Ptr CInt
         -> TriArg a
         -> Ptr a
         -> Ptr CInt
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
     String
     (FuncUnpacked
        (Ptr CChar
         -> Ptr CInt
         -> Ptr CInt
         -> TriArg a
         -> Ptr a
         -> Ptr CInt
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> TriArg a
-> Ptr a
-> Ptr CInt
-> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
forall (m :: * -> *) f.
(m ~ Labelled (FuncCont f) (FuncLabel f), FunctionPair f) =>
m (FuncPacked f) -> m (FuncUnpacked f) -> f
applyFuncPair
            (String
-> (Ptr CChar
    -> Ptr CInt
    -> Ptr CInt
    -> Ptr a
    -> Ptr a
    -> Ptr CInt
    -> Ptr CInt
    -> IO ())
-> Labelled
     (FuncCont
        (Ptr CChar
         -> Ptr CInt
         -> Ptr CInt
         -> TriArg a
         -> Ptr a
         -> Ptr CInt
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
     String
     (Ptr CChar
      -> Ptr CInt
      -> Ptr CInt
      -> Ptr a
      -> Ptr a
      -> Ptr CInt
      -> Ptr CInt
      -> IO ())
forall label a r. label -> a -> Labelled r label a
label String
"ppsv" Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.ppsv)
            (String
-> (Ptr CChar
    -> Ptr CInt
    -> Ptr CInt
    -> Ptr a
    -> Ptr CInt
    -> Ptr a
    -> Ptr CInt
    -> Ptr CInt
    -> IO ())
-> Labelled
     (FuncCont
        (Ptr CChar
         -> Ptr CInt
         -> Ptr CInt
         -> TriArg a
         -> Ptr a
         -> Ptr CInt
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
     String
     (Ptr CChar
      -> Ptr CInt
      -> Ptr CInt
      -> Ptr a
      -> Ptr CInt
      -> Ptr a
      -> Ptr CInt
      -> Ptr CInt
      -> IO ())
forall label a r. label -> a -> Labelled r label a
label String
"posv" Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.posv)
            Ptr CChar
uploPtr Ptr CInt
nPtr Ptr CInt
nrhsPtr (Ptr a -> Int -> TriArg a
forall a. Ptr a -> Int -> TriArg a
triArg Ptr a
aPtr Int
n) Ptr a
xPtr Ptr CInt
ldxPtr

solveDecomposed ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>
   Upper pack sh a ->
   Full meas vert horiz sh nrhs a -> Full meas vert horiz sh nrhs a
solveDecomposed :: forall meas vert horiz sh nrhs a pack.
(Measure meas, C vert, C horiz, C sh, Eq sh, C nrhs, Floating a) =>
Upper pack sh a
-> Full meas vert horiz sh nrhs a -> Full meas vert horiz sh nrhs a
solveDecomposed
   (Array
      shape :: UpperTriangularP pack sh
shape@(Layout.Mosaic PackingSingleton pack
pack MirrorSingleton NoMirror
Layout.NoMirror UpLoSingleton Upper
_upper Order
orderA sh
shA)
      ForeignPtr a
a) =
   String
-> sh
-> (Int
    -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr CInt -> ContT () IO ())
-> Full meas vert horiz sh nrhs a
-> Full meas vert horiz sh nrhs a
forall meas vert horiz height width a.
(Measure meas, C vert, C horiz, C height, C width, Eq height,
 Floating a) =>
String
-> height
-> (Int
    -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr CInt -> ContT () IO ())
-> Full meas vert horiz height width a
-> Full meas vert horiz height width a
solver String
"Hermitian.solveDecomposed" sh
shA ((Int
  -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr CInt -> ContT () IO ())
 -> Full meas vert horiz sh nrhs a
 -> Full meas vert horiz sh nrhs a)
-> (Int
    -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr CInt -> ContT () IO ())
-> Full meas vert horiz sh nrhs a
-> Full meas vert horiz sh nrhs a
forall a b. (a -> b) -> a -> b
$ \Int
n Ptr CInt
nPtr Ptr CInt
nrhsPtr Ptr a
xPtr Ptr CInt
ldxPtr -> do
      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
orderA
      Ptr a
aPtr <- Conjugation -> Order -> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r.
Floating a =>
Conjugation -> Order -> Int -> ForeignPtr a -> ContT r IO (Ptr a)
copyTriangleToTemp Conjugation
Conjugated Order
orderA (UpperTriangularP pack sh -> Int
forall sh. C sh => sh -> Int
Shape.size UpperTriangularP pack sh
shape) ForeignPtr a
a
      String
-> PackingSingleton pack
-> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
-> ContT () IO ()
forall func pack r.
(func ~ (Ptr CInt -> IO ())) =>
String
-> PackingSingleton pack
-> Labelled2 r String func func
-> ContT r IO ()
withPackingLinear String
rankMsg PackingSingleton pack
pack (Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
 -> ContT () IO ())
-> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
-> ContT () IO ()
forall a b. (a -> b) -> a -> b
$
         Labelled
  (FuncCont
     (Ptr CChar
      -> Ptr CInt
      -> Ptr CInt
      -> TriArg a
      -> Ptr a
      -> Ptr CInt
      -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
  String
  (FuncPacked
     (Ptr CChar
      -> Ptr CInt
      -> Ptr CInt
      -> TriArg a
      -> Ptr a
      -> Ptr CInt
      -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
-> Labelled
     (FuncCont
        (Ptr CChar
         -> Ptr CInt
         -> Ptr CInt
         -> TriArg a
         -> Ptr a
         -> Ptr CInt
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
     String
     (FuncUnpacked
        (Ptr CChar
         -> Ptr CInt
         -> Ptr CInt
         -> TriArg a
         -> Ptr a
         -> Ptr CInt
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> TriArg a
-> Ptr a
-> Ptr CInt
-> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
forall (m :: * -> *) f.
(m ~ Labelled (FuncCont f) (FuncLabel f), FunctionPair f) =>
m (FuncPacked f) -> m (FuncUnpacked f) -> f
applyFuncPair
            (String
-> (Ptr CChar
    -> Ptr CInt
    -> Ptr CInt
    -> Ptr a
    -> Ptr a
    -> Ptr CInt
    -> Ptr CInt
    -> IO ())
-> Labelled
     (FuncCont
        (Ptr CChar
         -> Ptr CInt
         -> Ptr CInt
         -> TriArg a
         -> Ptr a
         -> Ptr CInt
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
     String
     (Ptr CChar
      -> Ptr CInt
      -> Ptr CInt
      -> Ptr a
      -> Ptr a
      -> Ptr CInt
      -> Ptr CInt
      -> IO ())
forall label a r. label -> a -> Labelled r label a
label String
"pptrs" Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.pptrs)
            (String
-> (Ptr CChar
    -> Ptr CInt
    -> Ptr CInt
    -> Ptr a
    -> Ptr CInt
    -> Ptr a
    -> Ptr CInt
    -> Ptr CInt
    -> IO ())
-> Labelled
     (FuncCont
        (Ptr CChar
         -> Ptr CInt
         -> Ptr CInt
         -> TriArg a
         -> Ptr a
         -> Ptr CInt
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
     String
     (Ptr CChar
      -> Ptr CInt
      -> Ptr CInt
      -> Ptr a
      -> Ptr CInt
      -> Ptr a
      -> Ptr CInt
      -> Ptr CInt
      -> IO ())
forall label a r. label -> a -> Labelled r label a
label String
"potrs" Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.potrs)
            Ptr CChar
uploPtr Ptr CInt
nPtr Ptr CInt
nrhsPtr (Ptr a -> Int -> TriArg a
forall a. Ptr a -> Int -> TriArg a
triArg Ptr a
aPtr Int
n) Ptr a
xPtr Ptr CInt
ldxPtr


inverse ::
   (Shape.C sh, Class.Floating a) => Hermitian pack sh a -> Hermitian pack sh a
inverse :: forall sh a pack.
(C sh, Floating a) =>
Hermitian pack sh a -> Hermitian pack sh a
inverse
   (Array shape :: HermitianP pack sh
shape@(Layout.Mosaic PackingSingleton pack
pack MirrorSingleton ConjugateMirror
_mirror UpLoSingleton Upper
_upper Order
order sh
sh) ForeignPtr a
a) =
      HermitianP pack sh
-> (Int -> Ptr a -> IO ()) -> Array (HermitianP pack sh) a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize HermitianP pack sh
shape ((Int -> Ptr a -> IO ()) -> Array (HermitianP pack sh) a)
-> (Int -> Ptr a -> IO ()) -> Array (HermitianP pack sh) a
forall a b. (a -> b) -> a -> b
$ \Int
triSize Ptr a
bPtr -> do

   let n :: Int
n = sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh
   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
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 CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr a
aPtr <- ((Ptr a -> IO ()) -> IO ()) -> ContT () IO (Ptr a)
forall {k} (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT (((Ptr a -> IO ()) -> IO ()) -> ContT () IO (Ptr a))
-> ((Ptr a -> IO ()) -> IO ()) -> ContT () IO (Ptr a)
forall a b. (a -> b) -> a -> b
$ ForeignPtr a -> (Ptr a -> IO ()) -> IO ()
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
a
      IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$ Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Int -> Ptr a -> Ptr a -> IO ()
copyBlock Int
triSize Ptr a
aPtr Ptr a
bPtr
      String
-> PackingSingleton pack
-> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
-> ContT () IO ()
forall func pack r.
(func ~ (Ptr CInt -> IO ())) =>
String
-> PackingSingleton pack
-> Labelled2 r String func func
-> ContT r IO ()
withPackingLinear String
definiteMsg PackingSingleton pack
pack (Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
 -> ContT () IO ())
-> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
-> ContT () IO ()
forall a b. (a -> b) -> a -> b
$
         Labelled
  (FuncCont
     (Ptr CChar
      -> Ptr CInt
      -> TriArg a
      -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
  String
  (FuncPacked
     (Ptr CChar
      -> Ptr CInt
      -> TriArg a
      -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
-> Labelled
     (FuncCont
        (Ptr CChar
         -> Ptr CInt
         -> TriArg a
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
     String
     (FuncUnpacked
        (Ptr CChar
         -> Ptr CInt
         -> TriArg a
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
-> Ptr CChar
-> Ptr CInt
-> TriArg a
-> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
forall (m :: * -> *) f.
(m ~ Labelled (FuncCont f) (FuncLabel f), FunctionPair f) =>
m (FuncPacked f) -> m (FuncUnpacked f) -> f
applyFuncPair
            (String
-> (Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ())
-> Labelled
     (FuncCont
        (Ptr CChar
         -> Ptr CInt
         -> TriArg a
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
     String
     (Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ())
forall label a r. label -> a -> Labelled r label a
label String
"pptrf" Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
LapackGen.pptrf)
            (String
-> (Ptr CChar
    -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ())
-> Labelled
     (FuncCont
        (Ptr CChar
         -> Ptr CInt
         -> TriArg a
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
     String
     (Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ())
forall label a r. label -> a -> Labelled r label a
label String
"potrf" Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ()
LapackGen.potrf)
            Ptr CChar
uploPtr Ptr CInt
nPtr (Ptr a -> Int -> TriArg a
forall a. Ptr a -> Int -> TriArg a
triArg Ptr a
bPtr Int
n)
      String
-> PackingSingleton pack
-> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
-> ContT () IO ()
forall func pack r.
(func ~ (Ptr CInt -> IO ())) =>
String
-> PackingSingleton pack
-> Labelled2 r String func func
-> ContT r IO ()
withPackingLinear String
rankMsg PackingSingleton pack
pack (Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
 -> ContT () IO ())
-> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
-> ContT () IO ()
forall a b. (a -> b) -> a -> b
$
         Labelled
  (FuncCont
     (Ptr CChar
      -> Ptr CInt
      -> TriArg a
      -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
  String
  (FuncPacked
     (Ptr CChar
      -> Ptr CInt
      -> TriArg a
      -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
-> Labelled
     (FuncCont
        (Ptr CChar
         -> Ptr CInt
         -> TriArg a
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
     String
     (FuncUnpacked
        (Ptr CChar
         -> Ptr CInt
         -> TriArg a
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
-> Ptr CChar
-> Ptr CInt
-> TriArg a
-> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
forall (m :: * -> *) f.
(m ~ Labelled (FuncCont f) (FuncLabel f), FunctionPair f) =>
m (FuncPacked f) -> m (FuncUnpacked f) -> f
applyFuncPair
            (String
-> (Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ())
-> Labelled
     (FuncCont
        (Ptr CChar
         -> Ptr CInt
         -> TriArg a
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
     String
     (Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ())
forall label a r. label -> a -> Labelled r label a
label String
"pptri" Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
LapackGen.pptri)
            (String
-> (Ptr CChar
    -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ())
-> Labelled
     (FuncCont
        (Ptr CChar
         -> Ptr CInt
         -> TriArg a
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
     String
     (Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ())
forall label a r. label -> a -> Labelled r label a
label String
"potri" Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ()
LapackGen.potri)
            Ptr CChar
uploPtr Ptr CInt
nPtr (Ptr a -> Int -> TriArg a
forall a. Ptr a -> Int -> TriArg a
triArg Ptr a
bPtr Int
n)
   PackingSingleton pack
-> Conjugation -> Order -> Int -> Ptr a -> IO ()
forall a pack.
Floating a =>
PackingSingleton pack
-> Conjugation -> Order -> Int -> Ptr a -> IO ()
Symmetric.complement PackingSingleton pack
pack Conjugation
Conjugated Order
order Int
n Ptr a
bPtr

decompose ::
   (Shape.C sh, Class.Floating a) => Hermitian pack sh a -> Upper pack sh a
decompose :: forall sh a pack.
(C sh, Floating a) =>
Hermitian pack sh a -> Upper pack sh a
decompose
   (Array (Layout.Mosaic PackingSingleton pack
pack MirrorSingleton ConjugateMirror
_mirror UpLoSingleton Upper
upper Order
order sh
sh) ForeignPtr a
a) =
      UpperTriangularP pack sh
-> (Int -> Ptr a -> IO ()) -> Array (UpperTriangularP pack sh) a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize
         (PackingSingleton pack
-> MirrorSingleton NoMirror
-> UpLoSingleton Upper
-> Order
-> sh
-> UpperTriangularP pack sh
forall pack mirror uplo size.
PackingSingleton pack
-> MirrorSingleton mirror
-> UpLoSingleton uplo
-> Order
-> size
-> Mosaic pack mirror uplo size
Layout.Mosaic PackingSingleton pack
pack MirrorSingleton NoMirror
Layout.NoMirror UpLoSingleton Upper
upper Order
order sh
sh) ((Int -> Ptr a -> IO ()) -> Array (UpperTriangularP pack sh) a)
-> (Int -> Ptr a -> IO ()) -> Array (UpperTriangularP pack sh) a
forall a b. (a -> b) -> a -> b
$
      \Int
triSize Ptr a
bPtr -> do
   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
      let uplo :: Char
uplo = Order -> Char
uploFromOrder Order
order
      Ptr CChar
uploPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
uplo
      let n :: Int
n = sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr a
aPtr <- ((Ptr a -> IO ()) -> IO ()) -> ContT () IO (Ptr a)
forall {k} (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT (((Ptr a -> IO ()) -> IO ()) -> ContT () IO (Ptr a))
-> ((Ptr a -> IO ()) -> IO ()) -> ContT () IO (Ptr a)
forall a b. (a -> b) -> a -> b
$ ForeignPtr a -> (Ptr a -> IO ()) -> IO ()
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
a
      let packed :: Bool
packed =
            case PackingSingleton pack
pack of PackingSingleton pack
Layout.Packed -> Bool
True; PackingSingleton pack
Layout.Unpacked -> Bool
False
      IO () -> ContT () IO ()
forall a. IO a -> ContT () IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO () -> ContT () IO ()) -> IO () -> ContT () IO ()
forall a b. (a -> b) -> a -> b
$
         if Bool
packed
            then Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Int -> Ptr a -> Ptr a -> IO ()
copyBlock Int
triSize Ptr a
aPtr Ptr a
bPtr
            else do
               a -> Int -> Ptr a -> IO ()
forall a. Floating a => a -> Int -> Ptr a -> IO ()
fill a
forall a. Floating a => a
zero (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n) Ptr a
bPtr
               Char -> Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
Char -> Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
copySubTrapezoid Char
uplo Int
n Int
n Int
n Ptr a
aPtr Int
n Ptr a
bPtr
      String
-> PackingSingleton pack
-> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
-> ContT () IO ()
forall func pack r.
(func ~ (Ptr CInt -> IO ())) =>
String
-> PackingSingleton pack
-> Labelled2 r String func func
-> ContT r IO ()
withPackingLinear String
definiteMsg PackingSingleton pack
pack (Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
 -> ContT () IO ())
-> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
-> ContT () IO ()
forall a b. (a -> b) -> a -> b
$
         Labelled
  (FuncCont
     (Ptr CChar
      -> Ptr CInt
      -> TriArg a
      -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
  String
  (FuncPacked
     (Ptr CChar
      -> Ptr CInt
      -> TriArg a
      -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
-> Labelled
     (FuncCont
        (Ptr CChar
         -> Ptr CInt
         -> TriArg a
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
     String
     (FuncUnpacked
        (Ptr CChar
         -> Ptr CInt
         -> TriArg a
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
-> Ptr CChar
-> Ptr CInt
-> TriArg a
-> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())
forall (m :: * -> *) f.
(m ~ Labelled (FuncCont f) (FuncLabel f), FunctionPair f) =>
m (FuncPacked f) -> m (FuncUnpacked f) -> f
applyFuncPair
            (String
-> (Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ())
-> Labelled
     (FuncCont
        (Ptr CChar
         -> Ptr CInt
         -> TriArg a
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
     String
     (Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ())
forall label a r. label -> a -> Labelled r label a
label String
"pptrf" Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
LapackGen.pptrf) (String
-> (Ptr CChar
    -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ())
-> Labelled
     (FuncCont
        (Ptr CChar
         -> Ptr CInt
         -> TriArg a
         -> Labelled2 () String (Ptr CInt -> IO ()) (Ptr CInt -> IO ())))
     String
     (Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ())
forall label a r. label -> a -> Labelled r label a
label String
"potrf" Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ()
LapackGen.potrf)
            Ptr CChar
uploPtr Ptr CInt
nPtr (Ptr a -> Int -> TriArg a
forall a. Ptr a -> Int -> TriArg a
triArg Ptr a
bPtr Int
n)


determinant :: (Shape.C sh, Class.Floating a) => Hermitian pack sh a -> RealOf a
determinant :: forall sh a pack.
(C sh, Floating a) =>
Hermitian pack sh a -> RealOf a
determinant =
   Determinant (Array (Mosaic pack ConjugateMirror Upper sh)) a
-> Array (Mosaic pack ConjugateMirror Upper sh) a -> RealOf a
forall (f :: * -> *) a. Determinant f a -> f a -> RealOf a
getDeterminant (Determinant (Array (Mosaic pack ConjugateMirror Upper sh)) a
 -> Array (Mosaic pack ConjugateMirror Upper sh) a -> RealOf a)
-> Determinant (Array (Mosaic pack ConjugateMirror Upper sh)) a
-> Array (Mosaic pack ConjugateMirror Upper sh) a
-> RealOf a
forall a b. (a -> b) -> a -> b
$
   Determinant (Array (Mosaic pack ConjugateMirror Upper sh)) Float
-> Determinant
     (Array (Mosaic pack ConjugateMirror Upper sh)) Double
-> Determinant
     (Array (Mosaic pack ConjugateMirror Upper sh)) (Complex Float)
-> Determinant
     (Array (Mosaic pack ConjugateMirror Upper sh)) (Complex Double)
-> Determinant (Array (Mosaic pack ConjugateMirror Upper sh)) a
forall a (f :: * -> *).
Floating a =>
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
forall (f :: * -> *).
f Float
-> f Double -> f (Complex Float) -> f (Complex Double) -> f a
Class.switchFloating
      ((Array (Mosaic pack ConjugateMirror Upper sh) Float
 -> RealOf Float)
-> Determinant (Array (Mosaic pack ConjugateMirror Upper sh)) Float
forall (f :: * -> *) a. (f a -> RealOf a) -> Determinant f a
Determinant Array (Mosaic pack ConjugateMirror Upper sh) Float -> Float
Array (Mosaic pack ConjugateMirror Upper sh) Float -> RealOf Float
forall sh a ar pack.
(C sh, Floating a, RealOf a ~ ar, Real ar) =>
Hermitian pack sh a -> ar
determinantAux) ((Array (Mosaic pack ConjugateMirror Upper sh) Double
 -> RealOf Double)
-> Determinant
     (Array (Mosaic pack ConjugateMirror Upper sh)) Double
forall (f :: * -> *) a. (f a -> RealOf a) -> Determinant f a
Determinant Array (Mosaic pack ConjugateMirror Upper sh) Double -> Double
Array (Mosaic pack ConjugateMirror Upper sh) Double
-> RealOf Double
forall sh a ar pack.
(C sh, Floating a, RealOf a ~ ar, Real ar) =>
Hermitian pack sh a -> ar
determinantAux)
      ((Array (Mosaic pack ConjugateMirror Upper sh) (Complex Float)
 -> RealOf (Complex Float))
-> Determinant
     (Array (Mosaic pack ConjugateMirror Upper sh)) (Complex Float)
forall (f :: * -> *) a. (f a -> RealOf a) -> Determinant f a
Determinant Array (Mosaic pack ConjugateMirror Upper sh) (Complex Float)
-> Float
Array (Mosaic pack ConjugateMirror Upper sh) (Complex Float)
-> RealOf (Complex Float)
forall sh a ar pack.
(C sh, Floating a, RealOf a ~ ar, Real ar) =>
Hermitian pack sh a -> ar
determinantAux) ((Array (Mosaic pack ConjugateMirror Upper sh) (Complex Double)
 -> RealOf (Complex Double))
-> Determinant
     (Array (Mosaic pack ConjugateMirror Upper sh)) (Complex Double)
forall (f :: * -> *) a. (f a -> RealOf a) -> Determinant f a
Determinant Array (Mosaic pack ConjugateMirror Upper sh) (Complex Double)
-> Double
Array (Mosaic pack ConjugateMirror Upper sh) (Complex Double)
-> RealOf (Complex Double)
forall sh a ar pack.
(C sh, Floating a, RealOf a ~ ar, Real ar) =>
Hermitian pack sh a -> ar
determinantAux)

determinantAux ::
   (Shape.C sh, Class.Floating a, RealOf a ~ ar, Class.Real ar) =>
   Hermitian pack sh a -> ar
determinantAux :: forall sh a ar pack.
(C sh, Floating a, RealOf a ~ ar, Real ar) =>
Hermitian pack sh a -> ar
determinantAux =
   (ar -> Int -> ar
forall a b. (Num a, Integral b) => a -> b -> a
^(Int
2::Int)) (ar -> ar)
-> (Hermitian pack sh a -> ar) -> Hermitian pack sh a -> ar
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Vector sh ar -> ar
forall sh a. (C sh, Floating a) => Vector sh a -> a
Vector.product (Vector sh ar -> ar)
-> (Hermitian pack sh a -> Vector sh ar)
-> Hermitian pack sh a
-> ar
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> ar) -> Array sh a -> Vector sh ar
forall sh a b.
(C sh, Storable a, Storable b) =>
(a -> b) -> Array sh a -> Array sh b
Array.map a -> ar
a -> RealOf a
forall a. Floating a => a -> RealOf a
realPart (Array sh a -> Vector sh ar)
-> (Hermitian pack sh a -> Array sh a)
-> Hermitian pack sh a
-> Vector sh ar
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Mosaic pack NoMirror Upper sh a -> Array sh a
forall uplo sh a pack mirror.
(UpLo uplo, C sh, Floating a) =>
Mosaic pack mirror uplo sh a -> Vector sh a
takeDiagonal (Mosaic pack NoMirror Upper sh a -> Array sh a)
-> (Hermitian pack sh a -> Mosaic pack NoMirror Upper sh a)
-> Hermitian pack sh a
-> Array sh a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Hermitian pack sh a -> Mosaic pack NoMirror Upper sh a
forall sh a pack.
(C sh, Floating a) =>
Hermitian pack sh a -> Upper pack sh a
decompose