module Numeric.LAPACK.Matrix.Symmetric.Private where

import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import Numeric.LAPACK.Matrix.Triangular.Private
         (diagonalPointerPairs, columnMajorPointers, rowMajorPointers,
          forPointers, pack, unpackToTemp, copyTriangleToTemp)
import Numeric.LAPACK.Matrix.Shape.Private
         (Order(RowMajor,ColumnMajor), uploFromOrder)
import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated, Conjugated))
import Numeric.LAPACK.Matrix.Private (Full)
import Numeric.LAPACK.Linear.Private (solver, withDeterminantInfo, withInfo)
import Numeric.LAPACK.Scalar (zero, one)
import Numeric.LAPACK.Private (copyBlock, copyToTemp, copyCondConjugate)

import qualified Numeric.LAPACK.FFI.Generic as LapackGen
import qualified Numeric.BLAS.FFI.Generic as BlasGen
import qualified Numeric.Netlib.Utility as Call
import qualified Numeric.Netlib.Class as Class

import qualified Data.Array.Comfort.Shape as Shape
import Data.Array.Comfort.Shape (triangleSize)

import Foreign.Marshal.Array (advancePtr)
import Foreign.C.Types (CInt)
import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)
import Foreign.Ptr (Ptr)
import Foreign.Storable (Storable, peek)

import qualified System.IO.Lazy as LazyIO

import Control.Monad.Trans.Cont (ContT(ContT), evalContT)
import Control.Monad.IO.Class (liftIO)
import Control.Applicative ((<$>))


unpack :: Class.Floating a =>
   Conjugation -> Order -> Int -> Ptr a -> Ptr a -> IO ()
unpack :: Conjugation -> Order -> Int -> Ptr a -> Ptr a -> IO ()
unpack Conjugation
conj Order
order Int
n Ptr a
packedPtr Ptr a
fullPtr = 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
incxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
1
   Ptr CInt
incyPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint 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
$ case Order
order of
      Order
RowMajor ->
         [(Int, (Ptr a, Ptr a))]
-> (Ptr CInt -> (Ptr a, Ptr a) -> IO ()) -> IO ()
forall a. [(Int, a)] -> (Ptr CInt -> a -> IO ()) -> IO ()
forPointers (Int -> Ptr a -> Ptr a -> [(Int, (Ptr a, Ptr a))]
forall a.
Storable a =>
Int -> Ptr a -> Ptr a -> [(Int, (Ptr a, Ptr a))]
rowMajorPointers Int
n Ptr a
fullPtr Ptr a
packedPtr) ((Ptr CInt -> (Ptr a, Ptr a) -> IO ()) -> IO ())
-> (Ptr CInt -> (Ptr a, Ptr a) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
               \Ptr CInt
nPtr (Ptr a
dstPtr,Ptr a
srcPtr) -> do
            Conjugation
-> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Conjugation
-> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
copyCondConjugate Conjugation
conj Ptr CInt
nPtr Ptr a
srcPtr Ptr CInt
incxPtr Ptr a
dstPtr Ptr CInt
incyPtr
            Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasGen.copy Ptr CInt
nPtr Ptr a
srcPtr Ptr CInt
incxPtr Ptr a
dstPtr Ptr CInt
incxPtr
      Order
ColumnMajor ->
         [(Int, ((Ptr a, Ptr a), Ptr a))]
-> (Ptr CInt -> ((Ptr a, Ptr a), Ptr a) -> IO ()) -> IO ()
forall a. [(Int, a)] -> (Ptr CInt -> a -> IO ()) -> IO ()
forPointers (Int -> Ptr a -> Ptr a -> [(Int, ((Ptr a, Ptr a), Ptr a))]
forall a.
Storable a =>
Int -> Ptr a -> Ptr a -> [(Int, ((Ptr a, Ptr a), Ptr a))]
columnMajorPointers Int
n Ptr a
fullPtr Ptr a
packedPtr) ((Ptr CInt -> ((Ptr a, Ptr a), Ptr a) -> IO ()) -> IO ())
-> (Ptr CInt -> ((Ptr a, Ptr a), Ptr a) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
               \Ptr CInt
nPtr ((Ptr a
dstRowPtr,Ptr a
dstColumnPtr),Ptr a
srcPtr) -> do
            Conjugation
-> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Conjugation
-> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
copyCondConjugate Conjugation
conj Ptr CInt
nPtr Ptr a
srcPtr Ptr CInt
incxPtr Ptr a
dstRowPtr Ptr CInt
incyPtr
            Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
BlasGen.copy Ptr CInt
nPtr Ptr a
srcPtr Ptr CInt
incxPtr Ptr a
dstColumnPtr Ptr CInt
incxPtr


square ::
   (Class.Floating a) =>
   Conjugation -> Order -> Int -> ForeignPtr a -> Ptr a -> IO ()
square :: Conjugation -> Order -> Int -> ForeignPtr a -> Ptr a -> IO ()
square Conjugation
conj Order
order Int
n ForeignPtr a
a Ptr a
bpPtr =
   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
sidePtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'L'
      Ptr CChar
uploPtr <- Char -> FortranIO () (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char Char
'U'
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
ldPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Ptr a
aPtr <- (Int -> Ptr a -> Ptr a -> IO ())
-> Int -> ForeignPtr a -> ContT () IO (Ptr a)
forall a r.
Storable a =>
(Int -> Ptr a -> Ptr a -> IO ())
-> Int -> ForeignPtr a -> ContT r IO (Ptr a)
unpackToTemp (Conjugation -> Order -> Int -> Ptr a -> Ptr a -> IO ()
forall a.
Floating a =>
Conjugation -> Order -> Int -> Ptr a -> Ptr a -> IO ()
unpack Conjugation
conj Order
order) Int
n ForeignPtr a
a
      Ptr a
bPtr <- Int -> ContT () IO (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n)
      Ptr a
alphaPtr <- a -> ContT () IO (Ptr a)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
forall a. Floating a => a
one
      Ptr a
betaPtr <- a -> ContT () IO (Ptr a)
forall a r. Floating a => a -> FortranIO r (Ptr a)
Call.number a
forall a. Floating a => a
zero
      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
$ do
         (if Conjugation
conjConjugation -> Conjugation -> Bool
forall a. Eq a => a -> a -> Bool
==Conjugation
Conjugated then Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
BlasGen.hemm else Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> IO ()
BlasGen.symm)
            Ptr CChar
sidePtr Ptr CChar
uploPtr
            Ptr CInt
nPtr Ptr CInt
nPtr Ptr a
alphaPtr Ptr a
aPtr Ptr CInt
ldPtr
            Ptr a
aPtr Ptr CInt
ldPtr Ptr a
betaPtr Ptr a
bPtr Ptr CInt
ldPtr
         Order -> Int -> Ptr a -> Ptr a -> IO ()
forall a. Floating a => Order -> Int -> Ptr a -> Ptr a -> IO ()
pack Order
order Int
n Ptr a
bPtr Ptr a
bpPtr


solve ::
   (Extent.C vert, Extent.C horiz,
    Shape.C width, Shape.C height, Eq height, Class.Floating a) =>
   String -> Conjugation -> Order -> height -> ForeignPtr a ->
   Full vert horiz height width a ->
   Full vert horiz height width a
solve :: String
-> Conjugation
-> Order
-> height
-> ForeignPtr a
-> Full vert horiz height width a
-> Full vert horiz height width a
solve String
name Conjugation
conj Order
order height
sh ForeignPtr a
a =
   String
-> height
-> (Int
    -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr CInt -> ContT () IO ())
-> Full vert horiz height width a
-> Full vert horiz height width a
forall vert horiz height width a.
(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 vert horiz height width a
-> Full vert horiz height width a
solver String
name height
sh ((Int
  -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr CInt -> ContT () IO ())
 -> Full vert horiz height width a
 -> Full vert horiz height width a)
-> (Int
    -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr CInt -> ContT () IO ())
-> Full vert horiz height width a
-> Full vert horiz height width 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
order
      Ptr a
apPtr <- 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
conj Order
order (Int -> Int
triangleSize Int
n) ForeignPtr a
a
      Ptr CInt
ipivPtr <- Int -> FortranIO () (Ptr CInt)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray 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
$
         let (String
lapackName,Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
slv) =
               case Conjugation
conj of
                  Conjugation
Conjugated -> (String
"hpsv", 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.hpsv)
                  Conjugation
NonConjugated -> (String
"spsv", 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.spsv)
         in String -> (Ptr CInt -> IO ()) -> IO ()
withInfo String
lapackName ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
               Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
slv Ptr CChar
uploPtr Ptr CInt
nPtr Ptr CInt
nrhsPtr Ptr a
apPtr Ptr CInt
ipivPtr Ptr a
xPtr Ptr CInt
ldxPtr


inverse ::
   Class.Floating a =>
   Conjugation -> Order -> Int -> ForeignPtr a -> Int -> Ptr a -> IO ()
inverse :: Conjugation
-> Order -> Int -> ForeignPtr a -> Int -> Ptr a -> IO ()
inverse Conjugation
conj Order
order Int
n ForeignPtr a
a Int
triSize Ptr a
bPtr = 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
   Ptr CInt
ipivPtr <- Int -> FortranIO () (Ptr CInt)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
n
   Ptr a
workPtr <- Int -> ContT () IO (Ptr a)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray 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
$ do
      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
      case Conjugation
conj of
         Conjugation
Conjugated -> do
            String -> (Ptr CInt -> IO ()) -> IO ()
withInfo String
"hptrf" ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ 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.hptrf Ptr CChar
uploPtr Ptr CInt
nPtr Ptr a
bPtr Ptr CInt
ipivPtr
            String -> (Ptr CInt -> IO ()) -> IO ()
withInfo String
"hptri" ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ Ptr CChar
-> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
LapackGen.hptri Ptr CChar
uploPtr Ptr CInt
nPtr Ptr a
bPtr Ptr CInt
ipivPtr Ptr a
workPtr
         Conjugation
NonConjugated -> do
            String -> (Ptr CInt -> IO ()) -> IO ()
withInfo String
"sptrf" ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ 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.sptrf Ptr CChar
uploPtr Ptr CInt
nPtr Ptr a
bPtr Ptr CInt
ipivPtr
            String -> (Ptr CInt -> IO ()) -> IO ()
withInfo String
"sptri" ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ Ptr CChar
-> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
LapackGen.sptri Ptr CChar
uploPtr Ptr CInt
nPtr Ptr a
bPtr Ptr CInt
ipivPtr Ptr a
workPtr


blockDiagonalPointers ::
   (Storable a) =>
   Order -> [(Ptr CInt, Ptr a)] -> LazyIO.T [(Ptr a, Maybe (Ptr a, Ptr a))]
blockDiagonalPointers :: Order -> [(Ptr CInt, Ptr a)] -> T [(Ptr a, Maybe (Ptr a, Ptr a))]
blockDiagonalPointers Order
order =
   let go :: [(Ptr a, Ptr a)] -> T [(Ptr a, Maybe (Ptr a, Ptr a))]
go ((Ptr a
ipiv0Ptr,Ptr a
a0Ptr):[(Ptr a, Ptr a)]
ptrs0) = do
         a
ipiv <- IO a -> T a
forall a. IO a -> T a
LazyIO.interleave (IO a -> T a) -> IO a -> T a
forall a b. (a -> b) -> a -> b
$ Ptr a -> IO a
forall a. Storable a => Ptr a -> IO a
peek Ptr a
ipiv0Ptr
         (Maybe (Ptr a, Ptr a)
ext,[(Ptr a, Maybe (Ptr a, Ptr a))]
ptrTuples) <-
            if a
ipiv a -> a -> Bool
forall a. Ord a => a -> a -> Bool
>= a
0
               then (,) Maybe (Ptr a, Ptr a)
forall a. Maybe a
Nothing ([(Ptr a, Maybe (Ptr a, Ptr a))]
 -> (Maybe (Ptr a, Ptr a), [(Ptr a, Maybe (Ptr a, Ptr a))]))
-> T [(Ptr a, Maybe (Ptr a, Ptr a))]
-> T (Maybe (Ptr a, Ptr a), [(Ptr a, Maybe (Ptr a, Ptr a))])
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [(Ptr a, Ptr a)] -> T [(Ptr a, Maybe (Ptr a, Ptr a))]
go [(Ptr a, Ptr a)]
ptrs0
               else
                  case [(Ptr a, Ptr a)]
ptrs0 of
                     [] -> String -> T (Maybe (Ptr a, Ptr a), [(Ptr a, Maybe (Ptr a, Ptr a))])
forall a. HasCallStack => String -> a
error String
"Symmetric.determinant: incomplete 2x2 block"
                     (Ptr a
_ipiv1Ptr,Ptr a
a1Ptr):[(Ptr a, Ptr a)]
ptrs1 ->
                        let bPtr :: Ptr a
bPtr =
                              case Order
order of
                                 Order
ColumnMajor -> Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
a1Ptr (-Int
1)
                                 Order
RowMajor -> Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
a0Ptr Int
1
                        in (,) ((Ptr a, Ptr a) -> Maybe (Ptr a, Ptr a)
forall a. a -> Maybe a
Just (Ptr a
a1Ptr,Ptr a
bPtr)) ([(Ptr a, Maybe (Ptr a, Ptr a))]
 -> (Maybe (Ptr a, Ptr a), [(Ptr a, Maybe (Ptr a, Ptr a))]))
-> T [(Ptr a, Maybe (Ptr a, Ptr a))]
-> T (Maybe (Ptr a, Ptr a), [(Ptr a, Maybe (Ptr a, Ptr a))])
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [(Ptr a, Ptr a)] -> T [(Ptr a, Maybe (Ptr a, Ptr a))]
go [(Ptr a, Ptr a)]
ptrs1
         [(Ptr a, Maybe (Ptr a, Ptr a))]
-> T [(Ptr a, Maybe (Ptr a, Ptr a))]
forall (m :: * -> *) a. Monad m => a -> m a
return ([(Ptr a, Maybe (Ptr a, Ptr a))]
 -> T [(Ptr a, Maybe (Ptr a, Ptr a))])
-> [(Ptr a, Maybe (Ptr a, Ptr a))]
-> T [(Ptr a, Maybe (Ptr a, Ptr a))]
forall a b. (a -> b) -> a -> b
$ (Ptr a
a0Ptr,Maybe (Ptr a, Ptr a)
ext) (Ptr a, Maybe (Ptr a, Ptr a))
-> [(Ptr a, Maybe (Ptr a, Ptr a))]
-> [(Ptr a, Maybe (Ptr a, Ptr a))]
forall a. a -> [a] -> [a]
: [(Ptr a, Maybe (Ptr a, Ptr a))]
ptrTuples
       go [] = [(Ptr a, Maybe (Ptr a, Ptr a))]
-> T [(Ptr a, Maybe (Ptr a, Ptr a))]
forall (m :: * -> *) a. Monad m => a -> m a
return []
   in [(Ptr CInt, Ptr a)] -> T [(Ptr a, Maybe (Ptr a, Ptr a))]
forall a a.
(Storable a, Ord a, Num a, Storable a) =>
[(Ptr a, Ptr a)] -> T [(Ptr a, Maybe (Ptr a, Ptr a))]
go

determinant ::
   (Class.Floating a, Class.Floating ar) =>
   Conjugation -> ((Ptr a, Maybe (Ptr a, Ptr a)) -> IO ar) ->
   Order -> Int -> ForeignPtr a -> IO ar
determinant :: Conjugation
-> ((Ptr a, Maybe (Ptr a, Ptr a)) -> IO ar)
-> Order
-> Int
-> ForeignPtr a
-> IO ar
determinant Conjugation
conj (Ptr a, Maybe (Ptr a, Ptr a)) -> IO ar
peekBlockDeterminant Order
order Int
n ForeignPtr a
a = ContT ar IO ar -> IO ar
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT ar IO ar -> IO ar) -> ContT ar IO ar -> IO ar
forall a b. (a -> b) -> a -> b
$ do
   Ptr CChar
uploPtr <- Char -> FortranIO ar (Ptr CChar)
forall r. Char -> FortranIO r (Ptr CChar)
Call.char (Char -> FortranIO ar (Ptr CChar))
-> Char -> FortranIO ar (Ptr CChar)
forall a b. (a -> b) -> a -> b
$ Order -> Char
uploFromOrder Order
order
   Ptr CInt
nPtr <- Int -> FortranIO ar (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
   Ptr a
aPtr <- Int -> ForeignPtr a -> ContT ar 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
   Ptr CInt
ipivPtr <- Int -> FortranIO ar (Ptr CInt)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
n
   let (String
name,Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ()
trf) =
         case Conjugation
conj of
            Conjugation
Conjugated -> (String
"hptrf", 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.hptrf)
            Conjugation
NonConjugated -> (String
"sptrf", 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.sptrf)
   IO ar -> ContT ar IO ar
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO ar -> ContT ar IO ar) -> IO ar -> ContT ar IO ar
forall a b. (a -> b) -> a -> b
$ String -> (Ptr CInt -> IO ()) -> IO ar -> IO ar
forall a.
Floating a =>
String -> (Ptr CInt -> IO ()) -> IO a -> IO a
withDeterminantInfo String
name
      (Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> Ptr CInt -> IO ()
trf Ptr CChar
uploPtr Ptr CInt
nPtr Ptr a
aPtr Ptr CInt
ipivPtr)
      (((ar -> IO ar
forall (m :: * -> *) a. Monad m => a -> m a
return (ar -> IO ar) -> ar -> IO ar
forall a b. (a -> b) -> a -> b
$!) (ar -> IO ar) -> IO ar -> IO ar
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<<) (IO ar -> IO ar) -> IO ar -> IO ar
forall a b. (a -> b) -> a -> b
$
       T ar -> IO ar
forall a. T a -> IO a
LazyIO.run
         (([ar] -> ar) -> T [ar] -> T ar
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap [ar] -> ar
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product (T [ar] -> T ar) -> T [ar] -> T ar
forall a b. (a -> b) -> a -> b
$
          ((Ptr a, Maybe (Ptr a, Ptr a)) -> T ar)
-> [(Ptr a, Maybe (Ptr a, Ptr a))] -> T [ar]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
(a -> m b) -> t a -> m (t b)
mapM (IO ar -> T ar
forall a. IO a -> T a
LazyIO.interleave (IO ar -> T ar)
-> ((Ptr a, Maybe (Ptr a, Ptr a)) -> IO ar)
-> (Ptr a, Maybe (Ptr a, Ptr a))
-> T ar
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Ptr a, Maybe (Ptr a, Ptr a)) -> IO ar
peekBlockDeterminant) ([(Ptr a, Maybe (Ptr a, Ptr a))] -> T [ar])
-> T [(Ptr a, Maybe (Ptr a, Ptr a))] -> T [ar]
forall (m :: * -> *) a b. Monad m => (a -> m b) -> m a -> m b
=<<
          Order -> [(Ptr CInt, Ptr a)] -> T [(Ptr a, Maybe (Ptr a, Ptr a))]
forall a.
Storable a =>
Order -> [(Ptr CInt, Ptr a)] -> T [(Ptr a, Maybe (Ptr a, Ptr a))]
blockDiagonalPointers Order
order
            (Order -> Int -> Ptr CInt -> Ptr a -> [(Ptr CInt, Ptr a)]
forall a b.
(Storable a, Storable b) =>
Order -> Int -> Ptr a -> Ptr b -> [(Ptr a, Ptr b)]
diagonalPointerPairs Order
order Int
n Ptr CInt
ipivPtr Ptr a
aPtr)))