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

import qualified Numeric.LAPACK.Matrix.Banded.Linear as BandedLin
import qualified Numeric.LAPACK.Matrix.Banded.Basic as Banded
import qualified Numeric.LAPACK.Matrix.Symmetric.Private as Symmetric
import qualified Numeric.LAPACK.Matrix.Triangular.Private as Tri
import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import qualified Numeric.LAPACK.Vector as Vector
import Numeric.LAPACK.Linear.Private (solver, withInfo)
import Numeric.LAPACK.Matrix.Triangular.Basic
         (Triangular, Symmetric, PowerDiag, takeDiagonal, strictNonUnitDiagonal)
import Numeric.LAPACK.Matrix.Shape.Private
         (transposeFromOrder, uploFromOrder, uploOrder, charFromTriDiag)
import Numeric.LAPACK.Matrix.Modifier (Conjugation(NonConjugated))
import Numeric.LAPACK.Matrix.Private (Full)
import Numeric.LAPACK.Private (copyBlock, copyToTemp)

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 Data.Array.Comfort.Shape (triangleSize)

import System.IO.Unsafe (unsafePerformIO)

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

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


solve ::
   (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag,
    Extent.C vert, Extent.C horiz,
    Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>
   Triangular lo diag up sh a ->
   Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a
solve :: Triangular lo diag up sh a
-> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a
solve =
   MultiplyRight
  diag
  sh
  a
  (Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
  lo
  up
-> Triangular lo diag up sh a
-> Full vert horiz sh nrhs a
-> Full vert horiz sh nrhs a
forall diag sh a b lo up.
MultiplyRight diag sh a b lo up -> Triangular lo diag up sh a -> b
Tri.getMultiplyRight (MultiplyRight
   diag
   sh
   a
   (Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
   lo
   up
 -> Triangular lo diag up sh a
 -> Full vert horiz sh nrhs a
 -> Full vert horiz sh nrhs a)
-> MultiplyRight
     diag
     sh
     a
     (Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
     lo
     up
-> Triangular lo diag up sh a
-> Full vert horiz sh nrhs a
-> Full vert horiz sh nrhs a
forall a b. (a -> b) -> a -> b
$
   MultiplyRight
  diag
  sh
  a
  (Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
  Empty
  Empty
-> MultiplyRight
     diag
     sh
     a
     (Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
     Empty
     Filled
-> MultiplyRight
     diag
     sh
     a
     (Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
     Filled
     Empty
-> MultiplyRight
     diag
     sh
     a
     (Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
     Filled
     Filled
-> MultiplyRight
     diag
     sh
     a
     (Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
     lo
     up
forall lo up (f :: * -> * -> *).
(Content lo, Content up) =>
f Empty Empty
-> f Empty Filled -> f Filled Empty -> f Filled Filled -> f lo up
MatrixShape.switchDiagUpLoSym
      ((Triangular Empty diag Empty sh a
 -> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
-> MultiplyRight
     diag
     sh
     a
     (Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
     Empty
     Empty
forall diag sh a b lo up.
(Triangular lo diag up sh a -> b)
-> MultiplyRight diag sh a b lo up
Tri.MultiplyRight ((Triangular Empty diag Empty sh a
  -> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
 -> MultiplyRight
      diag
      sh
      a
      (Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
      Empty
      Empty)
-> (Triangular Empty diag Empty sh a
    -> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
-> MultiplyRight
     diag
     sh
     a
     (Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
     Empty
     Empty
forall a b. (a -> b) -> a -> b
$ Square U0 U0 sh a
-> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a
forall sub super vert horiz sh nrhs a.
(Natural sub, Natural super, C vert, C horiz, C sh, Eq sh, C nrhs,
 Floating a) =>
Square sub super sh a
-> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a
BandedLin.solve (Square U0 U0 sh a
 -> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
-> (Triangular Empty diag Empty sh a -> Square U0 U0 sh a)
-> Triangular Empty diag Empty sh a
-> Full vert horiz sh nrhs a
-> Full vert horiz sh nrhs a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Triangular Empty diag Empty sh a -> Square U0 U0 sh a
forall sh a diag.
(C sh, Floating a) =>
FlexDiagonal diag sh a -> Diagonal sh a
Banded.fromDiagonal)
      ((Triangular Empty diag Filled sh a
 -> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
-> MultiplyRight
     diag
     sh
     a
     (Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
     Empty
     Filled
forall diag sh a b lo up.
(Triangular lo diag up sh a -> b)
-> MultiplyRight diag sh a b lo up
Tri.MultiplyRight Triangular Empty diag Filled sh a
-> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a
forall lo up diag vert horiz sh nrhs a.
(UpLo lo up, TriDiag diag, C vert, C horiz, C sh, Eq sh, C nrhs,
 Floating a) =>
Triangular lo diag up sh a
-> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a
solveTriangular)
      ((Triangular Filled diag Empty sh a
 -> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
-> MultiplyRight
     diag
     sh
     a
     (Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
     Filled
     Empty
forall diag sh a b lo up.
(Triangular lo diag up sh a -> b)
-> MultiplyRight diag sh a b lo up
Tri.MultiplyRight Triangular Filled diag Empty sh a
-> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a
forall lo up diag vert horiz sh nrhs a.
(UpLo lo up, TriDiag diag, C vert, C horiz, C sh, Eq sh, C nrhs,
 Floating a) =>
Triangular lo diag up sh a
-> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a
solveTriangular)
      ((Triangular Filled diag Filled sh a
 -> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
-> MultiplyRight
     diag
     sh
     a
     (Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
     Filled
     Filled
forall diag sh a b lo up.
(Triangular lo diag up sh a -> b)
-> MultiplyRight diag sh a b lo up
Tri.MultiplyRight ((Triangular Filled diag Filled sh a
  -> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
 -> MultiplyRight
      diag
      sh
      a
      (Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
      Filled
      Filled)
-> (Triangular Filled diag Filled sh a
    -> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
-> MultiplyRight
     diag
     sh
     a
     (Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
     Filled
     Filled
forall a b. (a -> b) -> a -> b
$ Symmetric sh a
-> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a
forall vert horiz sh nrhs a.
(C vert, C horiz, C sh, Eq sh, C nrhs, Floating a) =>
Symmetric sh a
-> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a
solveSymmetric (Symmetric sh a
 -> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
-> (Triangular Filled diag Filled sh a -> Symmetric sh a)
-> Triangular Filled diag Filled sh a
-> Full vert horiz sh nrhs a
-> Full vert horiz sh nrhs a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Triangular Filled diag Filled sh a -> Symmetric sh a
forall diag lo up sh a.
TriDiag diag =>
Triangular lo diag up sh a -> Triangular lo NonUnit up sh a
strictNonUnitDiagonal)

solveTriangular ::
   (MatrixShape.UpLo lo up, MatrixShape.TriDiag diag,
    Extent.C vert, Extent.C horiz,
    Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>
   Triangular lo diag up sh a ->
   Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a
solveTriangular :: Triangular lo diag up sh a
-> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a
solveTriangular (Array (MatrixShape.Triangular diag
diag (lo, up)
uplo Order
orderA sh
shA) ForeignPtr a
a) =
   String
-> sh
-> (Int
    -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr CInt -> ContT () IO ())
-> Full vert horiz sh nrhs a
-> Full vert horiz sh nrhs 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
"Triangular.solve" sh
shA ((Int
  -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr CInt -> ContT () IO ())
 -> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a)
-> (Int
    -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr CInt -> ContT () IO ())
-> Full vert horiz sh nrhs a
-> Full 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 -> Char) -> Order -> Char
forall a b. (a -> b) -> a -> b
$ (lo, up) -> Order -> Order
forall lo up.
(Content lo, Content up) =>
(lo, up) -> Order -> Order
uploOrder (lo, up)
uplo Order
orderA
      Ptr CChar
transPtr <- 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
transposeFromOrder Order
orderA
      Ptr CChar
diagPtr <- 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
$ diag -> Char
forall diag. TriDiag diag => diag -> Char
charFromTriDiag diag
diag
      Ptr a
apPtr <- 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
      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 -> (Ptr CInt -> IO ()) -> IO ()
withInfo String
"tptrs" ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
            Ptr CChar
-> Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.tptrs Ptr CChar
uploPtr Ptr CChar
transPtr Ptr CChar
diagPtr
               Ptr CInt
nPtr Ptr CInt
nrhsPtr Ptr a
apPtr Ptr a
xPtr Ptr CInt
ldxPtr

solveSymmetric ::
   (Extent.C vert, Extent.C horiz,
    Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>
   Symmetric sh a ->
   Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a
solveSymmetric :: Symmetric sh a
-> Full vert horiz sh nrhs a -> Full vert horiz sh nrhs a
solveSymmetric (Array (MatrixShape.Triangular NonUnit
_diag (Filled, Filled)
_uplo 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
"Symmetric.solve" Conjugation
NonConjugated Order
orderA sh
shA ForeignPtr a
a


inverse ::
   (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag,
    Shape.C sh, Class.Floating a) =>
   Triangular lo diag up sh a ->
   Triangular lo (PowerDiag lo up diag) up sh a
inverse :: Triangular lo diag up sh a
-> Triangular lo (PowerDiag lo up diag) up sh a
inverse =
   Power diag sh a lo up
-> Triangular lo diag up sh a
-> Triangular lo (PowerDiag lo up diag) up sh a
forall diag sh a lo up.
Power diag sh a lo up
-> Triangular lo diag up sh a
-> Triangular lo (PowerDiag lo up diag) up sh a
Tri.getPower (Power diag sh a lo up
 -> Triangular lo diag up sh a
 -> Triangular lo (PowerDiag lo up diag) up sh a)
-> Power diag sh a lo up
-> Triangular lo diag up sh a
-> Triangular lo (PowerDiag lo up diag) up sh a
forall a b. (a -> b) -> a -> b
$
   Power diag sh a Empty Empty
-> Power diag sh a Empty Filled
-> Power diag sh a Filled Empty
-> Power diag sh a Filled Filled
-> Power diag sh a lo up
forall lo up (f :: * -> * -> *).
(Content lo, Content up) =>
f Empty Empty
-> f Empty Filled -> f Filled Empty -> f Filled Filled -> f lo up
MatrixShape.switchDiagUpLoSym
      ((Triangular Empty diag Empty sh a
 -> Triangular Empty (PowerDiag Empty Empty diag) Empty sh a)
-> Power diag sh a Empty Empty
forall diag sh a lo up.
(Triangular lo diag up sh a
 -> Triangular lo (PowerDiag lo up diag) up sh a)
-> Power diag sh a lo up
Tri.Power Triangular Empty diag Empty sh a
-> Triangular Empty (PowerDiag Empty Empty diag) Empty sh a
forall diag sh a.
(TriDiag diag, C sh, Floating a) =>
FlexDiagonal diag sh a -> FlexDiagonal diag sh a
inverseDiagonal)
      ((Triangular Empty diag Filled sh a
 -> Triangular Empty (PowerDiag Empty Filled diag) Filled sh a)
-> Power diag sh a Empty Filled
forall diag sh a lo up.
(Triangular lo diag up sh a
 -> Triangular lo (PowerDiag lo up diag) up sh a)
-> Power diag sh a lo up
Tri.Power Triangular Empty diag Filled sh a
-> Triangular Empty (PowerDiag Empty Filled diag) Filled sh a
forall lo up diag sh a.
(UpLo lo up, TriDiag diag, C sh, Floating a) =>
Triangular lo diag up sh a -> Triangular lo diag up sh a
inverseTriangular)
      ((Triangular Filled diag Empty sh a
 -> Triangular Filled (PowerDiag Filled Empty diag) Empty sh a)
-> Power diag sh a Filled Empty
forall diag sh a lo up.
(Triangular lo diag up sh a
 -> Triangular lo (PowerDiag lo up diag) up sh a)
-> Power diag sh a lo up
Tri.Power Triangular Filled diag Empty sh a
-> Triangular Filled (PowerDiag Filled Empty diag) Empty sh a
forall lo up diag sh a.
(UpLo lo up, TriDiag diag, C sh, Floating a) =>
Triangular lo diag up sh a -> Triangular lo diag up sh a
inverseTriangular)
      ((Triangular Filled diag Filled sh a
 -> Triangular Filled (PowerDiag Filled Filled diag) Filled sh a)
-> Power diag sh a Filled Filled
forall diag sh a lo up.
(Triangular lo diag up sh a
 -> Triangular lo (PowerDiag lo up diag) up sh a)
-> Power diag sh a lo up
Tri.Power ((Triangular Filled diag Filled sh a
  -> Triangular Filled (PowerDiag Filled Filled diag) Filled sh a)
 -> Power diag sh a Filled Filled)
-> (Triangular Filled diag Filled sh a
    -> Triangular Filled (PowerDiag Filled Filled diag) Filled sh a)
-> Power diag sh a Filled Filled
forall a b. (a -> b) -> a -> b
$ Symmetric sh a -> Symmetric sh a
forall sh a. (C sh, Floating a) => Symmetric sh a -> Symmetric sh a
inverseSymmetric (Symmetric sh a -> Symmetric sh a)
-> (Triangular Filled diag Filled sh a -> Symmetric sh a)
-> Triangular Filled diag Filled sh a
-> Symmetric sh a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Triangular Filled diag Filled sh a -> Symmetric sh a
forall diag lo up sh a.
TriDiag diag =>
Triangular lo diag up sh a -> Triangular lo NonUnit up sh a
strictNonUnitDiagonal)

inverseDiagonal ::
   (MatrixShape.TriDiag diag, Shape.C sh, Class.Floating a) =>
   Tri.FlexDiagonal diag sh a -> Tri.FlexDiagonal diag sh a
inverseDiagonal :: FlexDiagonal diag sh a -> FlexDiagonal diag sh a
inverseDiagonal FlexDiagonal diag sh a
a =
   diag
-> FlexDiagonal diag sh a
-> FlexDiagonal diag sh a
-> FlexDiagonal diag sh a
forall diag a. TriDiag diag => diag -> a -> a -> a
MatrixShape.caseTriDiag
      (Triangular Empty diag Empty sh -> diag
forall lo diag up size. Triangular lo diag up size -> diag
MatrixShape.triangularDiag (Triangular Empty diag Empty sh -> diag)
-> Triangular Empty diag Empty sh -> diag
forall a b. (a -> b) -> a -> b
$ FlexDiagonal diag sh a -> Triangular Empty diag Empty sh
forall sh a. Array sh a -> sh
Array.shape FlexDiagonal diag sh a
a)
      FlexDiagonal diag sh a
a (FlexDiagonal diag sh a -> FlexDiagonal diag sh a
forall sh a. (C sh, Floating a) => Vector sh a -> Vector sh a
Vector.recip FlexDiagonal diag sh a
a)

inverseTriangular ::
   (MatrixShape.UpLo lo up, MatrixShape.TriDiag diag,
    Shape.C sh, Class.Floating a) =>
   Triangular lo diag up sh a -> Triangular lo diag up sh a
inverseTriangular :: Triangular lo diag up sh a -> Triangular lo diag up sh a
inverseTriangular (Array shape :: Triangular lo diag up sh
shape@(MatrixShape.Triangular diag
diag (lo, up)
uplo Order
order sh
sh) ForeignPtr a
a) =
      Triangular lo diag up sh
-> (Int -> Ptr a -> IO ()) -> Triangular lo diag up sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize Triangular lo diag up sh
shape ((Int -> Ptr a -> IO ()) -> Triangular lo diag up sh a)
-> (Int -> Ptr a -> IO ()) -> Triangular lo diag up sh a
forall a b. (a -> b) -> a -> b
$ \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 -> Char) -> Order -> Char
forall a b. (a -> b) -> a -> b
$ (lo, up) -> Order -> Order
forall lo up.
(Content lo, Content up) =>
(lo, up) -> Order -> Order
uploOrder (lo, up)
uplo Order
order
      Ptr CChar
diagPtr <- 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
$ diag -> Char
forall diag. TriDiag diag => diag -> Char
charFromTriDiag diag
diag
      Ptr CInt
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint (Int -> FortranIO () (Ptr CInt)) -> Int -> FortranIO () (Ptr CInt)
forall a b. (a -> b) -> a -> b
$ sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh
      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 (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
         String -> (Ptr CInt -> IO ()) -> IO ()
withInfo String
"tptri" ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ Ptr CChar -> Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
forall a.
Floating a =>
Ptr CChar -> Ptr CChar -> Ptr CInt -> Ptr a -> Ptr CInt -> IO ()
LapackGen.tptri Ptr CChar
uploPtr Ptr CChar
diagPtr Ptr CInt
nPtr Ptr a
bPtr

inverseSymmetric ::
   (Shape.C sh, Class.Floating a) => Symmetric sh a -> Symmetric sh a
inverseSymmetric :: Symmetric sh a -> Symmetric sh a
inverseSymmetric (Array shape :: FlexSymmetric NonUnit sh
shape@(MatrixShape.Triangular NonUnit
_diag (Filled, Filled)
_uplo Order
order sh
sh) ForeignPtr a
a) =
   FlexSymmetric NonUnit sh
-> (Int -> Ptr a -> IO ()) -> Symmetric sh a
forall sh a.
(C sh, Storable a) =>
sh -> (Int -> Ptr a -> IO ()) -> Array sh a
Array.unsafeCreateWithSize FlexSymmetric NonUnit sh
shape ((Int -> Ptr a -> IO ()) -> Symmetric sh a)
-> (Int -> Ptr a -> IO ()) -> Symmetric 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
NonConjugated Order
order (sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh) ForeignPtr a
a


determinant ::
   (MatrixShape.Content lo, MatrixShape.Content up, MatrixShape.TriDiag diag,
    Shape.C sh, Class.Floating a) =>
   Triangular lo diag up sh a -> a
determinant :: Triangular lo diag up sh a -> a
determinant =
   MultiplyRight diag sh a a lo up -> Triangular lo diag up sh a -> a
forall diag sh a b lo up.
MultiplyRight diag sh a b lo up -> Triangular lo diag up sh a -> b
Tri.getMultiplyRight (MultiplyRight diag sh a a lo up
 -> Triangular lo diag up sh a -> a)
-> MultiplyRight diag sh a a lo up
-> Triangular lo diag up sh a
-> a
forall a b. (a -> b) -> a -> b
$
   MultiplyRight diag sh a a Empty Empty
-> MultiplyRight diag sh a a Empty Filled
-> MultiplyRight diag sh a a Filled Empty
-> MultiplyRight diag sh a a Filled Filled
-> MultiplyRight diag sh a a lo up
forall lo up (f :: * -> * -> *).
(Content lo, Content up) =>
f Empty Empty
-> f Empty Filled -> f Filled Empty -> f Filled Filled -> f lo up
MatrixShape.switchDiagUpLoSym
      ((Triangular Empty diag Empty sh a -> a)
-> MultiplyRight diag sh a a Empty Empty
forall diag sh a b lo up.
(Triangular lo diag up sh a -> b)
-> MultiplyRight diag sh a b lo up
Tri.MultiplyRight Triangular Empty diag Empty sh a -> a
forall lo up sh a diag.
(DiagUpLo lo up, C sh, Floating a) =>
Triangular lo diag up sh a -> a
determinantTriangular)
      ((Triangular Empty diag Filled sh a -> a)
-> MultiplyRight diag sh a a Empty Filled
forall diag sh a b lo up.
(Triangular lo diag up sh a -> b)
-> MultiplyRight diag sh a b lo up
Tri.MultiplyRight Triangular Empty diag Filled sh a -> a
forall lo up sh a diag.
(DiagUpLo lo up, C sh, Floating a) =>
Triangular lo diag up sh a -> a
determinantTriangular)
      ((Triangular Filled diag Empty sh a -> a)
-> MultiplyRight diag sh a a Filled Empty
forall diag sh a b lo up.
(Triangular lo diag up sh a -> b)
-> MultiplyRight diag sh a b lo up
Tri.MultiplyRight Triangular Filled diag Empty sh a -> a
forall lo up sh a diag.
(DiagUpLo lo up, C sh, Floating a) =>
Triangular lo diag up sh a -> a
determinantTriangular)
      ((Triangular Filled diag Filled sh a -> a)
-> MultiplyRight diag sh a a Filled Filled
forall diag sh a b lo up.
(Triangular lo diag up sh a -> b)
-> MultiplyRight diag sh a b lo up
Tri.MultiplyRight ((Triangular Filled diag Filled sh a -> a)
 -> MultiplyRight diag sh a a Filled Filled)
-> (Triangular Filled diag Filled sh a -> a)
-> MultiplyRight diag sh a a Filled Filled
forall a b. (a -> b) -> a -> b
$ Symmetric sh a -> a
forall sh a. (C sh, Floating a) => Symmetric sh a -> a
determinantSymmetric (Symmetric sh a -> a)
-> (Triangular Filled diag Filled sh a -> Symmetric sh a)
-> Triangular Filled diag Filled sh a
-> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Triangular Filled diag Filled sh a -> Symmetric sh a
forall diag lo up sh a.
TriDiag diag =>
Triangular lo diag up sh a -> Triangular lo NonUnit up sh a
strictNonUnitDiagonal)

determinantTriangular ::
   (MatrixShape.DiagUpLo lo up, Shape.C sh, Class.Floating a) =>
   Triangular lo diag up sh a -> a
determinantTriangular :: Triangular lo diag up sh a -> a
determinantTriangular = [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product ([a] -> a)
-> (Triangular lo diag up sh a -> [a])
-> Triangular lo diag up sh a
-> a
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Array sh a -> [a]
forall sh a. (C sh, Storable a) => Array sh a -> [a]
Array.toList (Array sh a -> [a])
-> (Triangular lo diag up sh a -> Array sh a)
-> Triangular lo diag up sh a
-> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Triangular lo diag up sh a -> Array sh a
forall lo up sh a diag.
(Content lo, Content up, C sh, Floating a) =>
Triangular lo diag up sh a -> Vector sh a
takeDiagonal

determinantSymmetric ::
   (Shape.C sh, Class.Floating a) => Symmetric sh a -> a
determinantSymmetric :: Symmetric sh a -> a
determinantSymmetric (Array (MatrixShape.Triangular NonUnit
_diag (Filled, Filled)
_uplo Order
order sh
sh) ForeignPtr a
a) =
   IO a -> a
forall a. IO a -> a
unsafePerformIO (IO a -> a) -> IO a -> a
forall a b. (a -> b) -> a -> b
$
      Conjugation
-> ((Ptr a, Maybe (Ptr a, Ptr a)) -> IO a)
-> Order
-> Int
-> ForeignPtr a
-> IO a
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
NonConjugated
         (Ptr a, Maybe (Ptr a, Ptr a)) -> IO a
forall a. Floating a => (Ptr a, Maybe (Ptr a, Ptr a)) -> IO a
peekBlockDeterminant Order
order (sh -> Int
forall sh. C sh => sh -> Int
Shape.size sh
sh) ForeignPtr a
a

peekBlockDeterminant ::
   (Class.Floating a) => (Ptr a, Maybe (Ptr a, Ptr a)) -> IO a
peekBlockDeterminant :: (Ptr a, Maybe (Ptr a, Ptr a)) -> IO a
peekBlockDeterminant (Ptr a
a0Ptr,Maybe (Ptr a, Ptr a)
ext) = do
   a
a0 <- Ptr a -> IO a
forall a. Storable a => Ptr a -> IO a
peek Ptr a
a0Ptr
   case Maybe (Ptr a, Ptr a)
ext of
      Maybe (Ptr a, Ptr a)
Nothing -> a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return a
a0
      Just (Ptr a
a1Ptr,Ptr a
bPtr) -> do
         a
a1 <- Ptr a -> IO a
forall a. Storable a => Ptr a -> IO a
peek Ptr a
a1Ptr
         a
b <- Ptr a -> IO a
forall a. Storable a => Ptr a -> IO a
peek Ptr a
bPtr
         a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (a
a0a -> a -> a
forall a. Num a => a -> a -> a
*a
a1 a -> a -> a
forall a. Num a => a -> a -> a
- a
ba -> a -> a
forall a. Num a => a -> a -> a
*a
b)