{-# LANGUAGE TypeFamilies #-}
module Numeric.LAPACK.Matrix.Banded.Linear (
   solve,
   solveColumnMajor,
   solveTriangular,
   determinant,
   ) where

import qualified Numeric.LAPACK.Matrix.Banded.Basic as Banded
import qualified Numeric.LAPACK.Matrix.Shape.Omni as Omni
import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import qualified Numeric.LAPACK.Permutation.Private as Perm
import qualified Numeric.LAPACK.Private as Private
import Numeric.LAPACK.Linear.Private (solver, withDeterminantInfo, withInfo)
import Numeric.LAPACK.Matrix.Layout.Private
         (Order(RowMajor,ColumnMajor), transposeFromOrder)
import Numeric.LAPACK.Matrix.Private (Full)
import Numeric.LAPACK.Private (copySubMatrix)

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

import qualified Type.Data.Num.Unary as Unary
import Type.Data.Num (integralFromProxy)

import qualified Data.Array.Comfort.Shape as Shape
import Data.Array.Comfort.Storable.Unchecked (Array(Array))

import System.IO.Unsafe (unsafePerformIO)

import Foreign.Marshal.Array (peekArray, advancePtr)
import Foreign.ForeignPtr (withForeignPtr)

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


solve ::
   (Unary.Natural sub, Unary.Natural super,
    Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>
   Banded.Square sub super sh a ->
   Full meas vert horiz sh nrhs a -> Full meas vert horiz sh nrhs a
solve :: Square sub super sh a
-> Full meas vert horiz sh nrhs a -> Full meas vert horiz sh nrhs a
solve (Array (Layout.Banded (UnaryProxy sub, UnaryProxy super)
numOff Order
order Extent Shape Small Small sh sh
extent) 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
"Banded.solve" (Extent Shape Small Small sh sh -> sh
forall shape. Square shape -> shape
Extent.squareSize Extent Shape Small Small sh sh
extent) ((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
      let (Int
kl,Int
ku) = Order -> (UnaryProxy sub, UnaryProxy super) -> (Int, Int)
forall sub super.
(Natural sub, Natural super) =>
Order -> (UnaryProxy sub, UnaryProxy super) -> (Int, Int)
Layout.numOffDiagonals Order
order (UnaryProxy sub, UnaryProxy super)
numOff
      let k :: Int
k = Int
klInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
ku
      let ldab :: Int
ldab = Int
klInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
k
      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
order
      Ptr CInt
klPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
kl
      Ptr CInt
kuPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
ku
      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 a
abPtr <- 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
ldab)
      Ptr CInt
ldabPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
ldab
      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
$ do
         Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
copySubMatrix Int
k Int
n Int
k Ptr a
aPtr Int
ldab (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
abPtr Int
kl)
         String -> (Ptr CInt -> IO ()) -> IO ()
withInfo String
"gbtrf" ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
            Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.gbtrf Ptr CInt
nPtr Ptr CInt
nPtr Ptr CInt
klPtr Ptr CInt
kuPtr Ptr a
abPtr Ptr CInt
ldabPtr Ptr CInt
ipivPtr
         String -> (Ptr CInt -> IO ()) -> IO ()
withInfo String
"gbtrs" ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
            Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.gbtrs Ptr CChar
transPtr Ptr CInt
nPtr Ptr CInt
klPtr Ptr CInt
kuPtr Ptr CInt
nrhsPtr
               Ptr a
abPtr Ptr CInt
ldabPtr Ptr CInt
ipivPtr Ptr a
xPtr Ptr CInt
ldxPtr

solveColumnMajor ::
   (Unary.Natural sub, Unary.Natural super,
    Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>
   Banded.Square sub super sh a ->
   Full meas vert horiz sh nrhs a -> Full meas vert horiz sh nrhs a
solveColumnMajor :: Square sub super sh a
-> Full meas vert horiz sh nrhs a -> Full meas vert horiz sh nrhs a
solveColumnMajor
      (Array (Layout.Banded (UnaryProxy sub
sub,UnaryProxy super
super) Order
ColumnMajor Extent Shape Small Small sh sh
extent) 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
"Banded.solve" (Extent Shape Small Small sh sh -> sh
forall shape. Square shape -> shape
Extent.squareSize Extent Shape Small Small sh sh
extent) ((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
      let kl :: Int
kl = UnaryProxy sub -> Int
forall x y. (Integer x, Num y) => Proxy x -> y
integralFromProxy UnaryProxy sub
sub
      let ku :: Int
ku = UnaryProxy super -> Int
forall x y. (Integer x, Num y) => Proxy x -> y
integralFromProxy UnaryProxy super
super
      let k :: Int
k = Int
klInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
ku
      let ldab :: Int
ldab = Int
klInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
k
      Ptr CInt
klPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
kl
      Ptr CInt
kuPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
ku
      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 a
abPtr <- 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
ldab)
      Ptr CInt
ldabPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
ldab
      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
$ do
         Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
copySubMatrix Int
k Int
n Int
k Ptr a
aPtr Int
ldab (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
abPtr Int
kl)
         String -> (Ptr CInt -> IO ()) -> IO ()
withInfo String
"gbsv" ((Ptr CInt -> IO ()) -> IO ()) -> (Ptr CInt -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$
            Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.gbsv Ptr CInt
nPtr Ptr CInt
klPtr Ptr CInt
kuPtr Ptr CInt
nrhsPtr
               Ptr a
abPtr Ptr CInt
ldabPtr Ptr CInt
ipivPtr Ptr a
xPtr Ptr CInt
ldxPtr
solveColumnMajor (Array (Layout.Banded (UnaryProxy sub, UnaryProxy super)
_ Order
RowMajor Extent Shape Small Small sh sh
_) ForeignPtr a
_) =
   String
-> Full meas vert horiz sh nrhs a -> Full meas vert horiz sh nrhs a
forall a. HasCallStack => String -> a
error String
"Linear.Banded.solveColumnMajor: RowMajor intentionally unimplemented"

solveTriangular ::
   (Omni.BandedTriangular sub super, Omni.TriDiag diag,
    Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C sh, Eq sh, Shape.C nrhs, Class.Floating a) =>
   Omni.DiagSingleton diag ->
   Banded.Square sub super sh a ->
   Full meas vert horiz sh nrhs a -> Full meas vert horiz sh nrhs a
solveTriangular :: DiagSingleton diag
-> Square sub super sh a
-> Full meas vert horiz sh nrhs a
-> Full meas vert horiz sh nrhs a
solveTriangular DiagSingleton diag
diag (Array (Layout.Banded (UnaryProxy sub, UnaryProxy super)
numOff Order
order Extent Shape Small Small sh sh
extent) 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
"Banded.solveTriangular" (Extent Shape Small Small sh sh -> sh
forall shape. Square shape -> shape
Extent.squareSize Extent Shape Small Small sh sh
extent) ((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
      let (Int
kl,Int
ku) = Order -> (UnaryProxy sub, UnaryProxy super) -> (Int, Int)
forall sub super.
(Natural sub, Natural super) =>
Order -> (UnaryProxy sub, UnaryProxy super) -> (Int, Int)
Layout.numOffDiagonals Order
order (UnaryProxy sub, UnaryProxy super)
numOff
      let k :: Int
k = Int
klInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
ku
      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
$ if Int
klInt -> Int -> Bool
forall a. Eq a => a -> a -> Bool
==Int
0 then Char
'U' else Char
'L'
      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
Layout.transposeFromOrder 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
$ DiagSingleton diag -> Char
forall diag. TriDiag diag => DiagSingleton diag -> Char
Omni.charFromTriDiag DiagSingleton diag
diag
      Ptr CInt
kPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
k
      Ptr a
abPtr <- ((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
ldabPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
      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
"tbtrs" ((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 CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CChar
-> Ptr CChar
-> Ptr CChar
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.tbtrs Ptr CChar
uploPtr Ptr CChar
transPtr Ptr CChar
diagPtr Ptr CInt
nPtr Ptr CInt
kPtr Ptr CInt
nrhsPtr
               Ptr a
abPtr Ptr CInt
ldabPtr Ptr a
xPtr Ptr CInt
ldxPtr

determinant ::
   (Unary.Natural sub, Unary.Natural super, Shape.C sh, Class.Floating a) =>
   Banded.Square sub super sh a -> a
determinant :: Square sub super sh a -> a
determinant (Array (Layout.Banded (UnaryProxy sub, UnaryProxy super)
numOff Order
order Extent Shape Small Small sh sh
extent) 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
$ do
   let n :: Int
n = sh -> Int
forall sh. C sh => sh -> Int
Shape.size (sh -> Int) -> sh -> Int
forall a b. (a -> b) -> a -> b
$ Extent Shape Small Small sh sh -> sh
forall shape. Square shape -> shape
Extent.squareSize Extent Shape Small Small sh sh
extent
   ContT a IO a -> IO a
forall (m :: * -> *) r. Monad m => ContT r m r -> m r
evalContT (ContT a IO a -> IO a) -> ContT a IO a -> IO a
forall a b. (a -> b) -> a -> b
$ do
      let (Int
kl,Int
ku) = Order -> (UnaryProxy sub, UnaryProxy super) -> (Int, Int)
forall sub super.
(Natural sub, Natural super) =>
Order -> (UnaryProxy sub, UnaryProxy super) -> (Int, Int)
Layout.numOffDiagonals Order
order (UnaryProxy sub, UnaryProxy super)
numOff
      let k :: Int
k = Int
klInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
ku
      let ldab :: Int
ldab = Int
klInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
k
      Ptr CInt
nPtr <- Int -> FortranIO a (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
klPtr <- Int -> FortranIO a (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
kl
      Ptr CInt
kuPtr <- Int -> FortranIO a (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
ku
      Ptr a
aPtr <- ((Ptr a -> IO a) -> IO a) -> ContT a IO (Ptr a)
forall k (r :: k) (m :: k -> *) a.
((a -> m r) -> m r) -> ContT r m a
ContT (((Ptr a -> IO a) -> IO a) -> ContT a IO (Ptr a))
-> ((Ptr a -> IO a) -> IO a) -> ContT a IO (Ptr a)
forall a b. (a -> b) -> a -> b
$ ForeignPtr a -> (Ptr a -> IO a) -> IO a
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr a
a
      Ptr a
abPtr <- Int -> ContT a 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
ldab)
      Ptr CInt
ldabPtr <- Int -> FortranIO a (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
ldab
      Ptr CInt
ipivPtr <- Int -> FortranIO a (Ptr CInt)
forall a r. Storable a => Int -> FortranIO r (Ptr a)
Call.allocaArray Int
n
      IO a -> ContT a IO a
forall (m :: * -> *) a. MonadIO m => IO a -> m a
liftIO (IO a -> ContT a IO a) -> IO a -> ContT a IO a
forall a b. (a -> b) -> a -> b
$ do
         Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
forall a.
Floating a =>
Int -> Int -> Int -> Ptr a -> Int -> Ptr a -> IO ()
copySubMatrix Int
k Int
n Int
k Ptr a
aPtr Int
ldab (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
abPtr Int
kl)
         String -> (Ptr CInt -> IO ()) -> IO a -> IO a
forall a.
Floating a =>
String -> (Ptr CInt -> IO ()) -> IO a -> IO a
withDeterminantInfo String
"gbtrf"
            (Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> IO ()
forall a.
Floating a =>
Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> Ptr a
-> Ptr CInt
-> Ptr CInt
-> Ptr CInt
-> IO ()
LapackGen.gbtrf Ptr CInt
nPtr Ptr CInt
nPtr Ptr CInt
klPtr Ptr CInt
kuPtr Ptr a
abPtr Ptr CInt
ldabPtr Ptr CInt
ipivPtr)
            (do
               a
det <- Int -> Ptr a -> Int -> IO a
forall a. Floating a => Int -> Ptr a -> Int -> IO a
Private.product Int
n (Ptr a -> Int -> Ptr a
forall a. Storable a => Ptr a -> Int -> Ptr a
advancePtr Ptr a
abPtr (Int
klInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
ku)) Int
ldab
               [CInt]
ipiv <- Int -> Ptr CInt -> IO [CInt]
forall a. Storable a => Int -> Ptr a -> IO [a]
peekArray Int
n Ptr CInt
ipivPtr
               a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (a -> IO a) -> a -> IO a
forall a b. (a -> b) -> a -> b
$ [CInt] -> a -> a
forall a. Floating a => [CInt] -> a -> a
Perm.condNegate [CInt]
ipiv a
det)