module Numeric.LAPACK.Linear.Private where

import qualified Numeric.LAPACK.Matrix.Layout.Private as Layout
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import qualified Numeric.LAPACK.Private as Private
import Numeric.LAPACK.Matrix.Layout.Private (Order(ColumnMajor))
import Numeric.LAPACK.Matrix.Private (Full)
import Numeric.LAPACK.Scalar (zero)
import Numeric.LAPACK.Private (copyToColumnMajor, peekCInt, argMsg)

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.Marshal.Alloc (alloca)
import Foreign.C.Types (CInt)
import Foreign.ForeignPtr (withForeignPtr)
import Foreign.Ptr (Ptr)

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

import Text.Printf (printf)


solver ::
   (Extent.Measure meas, Extent.C vert, Extent.C horiz,
    Shape.C height, Shape.C width, Eq height,
    Class.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 :: 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
name height
sh Int -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr CInt -> ContT () IO ()
f (Array (Layout.Full Order
order Extent meas vert horiz height width
extent) ForeignPtr a
b) =
   Full meas vert horiz height width
-> (Ptr a -> IO ()) -> Array (Full meas vert horiz height width) a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order
-> Extent meas vert horiz height width
-> Full meas vert horiz height width
forall meas vert horiz height width.
Order
-> Extent meas vert horiz height width
-> Full meas vert horiz height width
Layout.Full Order
ColumnMajor Extent meas vert horiz height width
extent) ((Ptr a -> IO ()) -> Array (Full meas vert horiz height width) a)
-> (Ptr a -> IO ()) -> Array (Full meas vert horiz height width) a
forall a b. (a -> b) -> a -> b
$
      \Ptr a
xPtr -> do
   let (height
height,width
width) = Extent meas vert horiz height width -> (height, width)
forall meas vert horiz height width.
(Measure meas, C vert, C horiz) =>
Extent meas vert horiz height width -> (height, width)
Extent.dimensions Extent meas vert horiz height width
extent
   String -> Bool -> IO ()
Call.assert (String
name String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
": height shapes mismatch") (height
sh height -> height -> Bool
forall a. Eq a => a -> a -> Bool
== height
height)
   let n :: Int
n = height -> Int
forall sh. C sh => sh -> Int
Shape.size height
height
   let nrhs :: Int
nrhs = width -> Int
forall sh. C sh => sh -> Int
Shape.size width
width
   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
nPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
n
      Ptr CInt
nrhsPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.cint Int
nrhs
      Ptr a
bPtr <- ((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
b
      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
$ Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
forall a.
Floating a =>
Order -> Int -> Int -> Ptr a -> Ptr a -> IO ()
copyToColumnMajor Order
order Int
n Int
nrhs Ptr a
bPtr Ptr a
xPtr
      Ptr CInt
ldxPtr <- Int -> FortranIO () (Ptr CInt)
forall r. Int -> FortranIO r (Ptr CInt)
Call.leadingDim Int
n
      Int -> Ptr CInt -> Ptr CInt -> Ptr a -> Ptr CInt -> ContT () IO ()
f Int
n Ptr CInt
nPtr Ptr CInt
nrhsPtr Ptr a
xPtr Ptr CInt
ldxPtr


withDeterminantInfo ::
   (Class.Floating a) =>
   String -> (Ptr CInt -> IO ()) -> IO a -> IO a
withDeterminantInfo :: forall a.
Floating a =>
String -> (Ptr CInt -> IO ()) -> IO a -> IO a
withDeterminantInfo String
name Ptr CInt -> IO ()
computation IO a
evaluation = (Ptr CInt -> IO a) -> IO a
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr CInt -> IO a) -> IO a) -> (Ptr CInt -> IO a) -> IO a
forall a b. (a -> b) -> a -> b
$ \Ptr CInt
infoPtr -> do
   Ptr CInt -> IO ()
computation Ptr CInt
infoPtr
   Int
info <- Ptr CInt -> IO Int
peekCInt Ptr CInt
infoPtr
   case Int -> Int -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Int
info (Int
0::Int) of
      Ordering
LT -> String -> IO a
forall a. HasCallStack => String -> a
error (String -> IO a) -> String -> IO a
forall a b. (a -> b) -> a -> b
$ String -> String -> Int -> String
forall r. PrintfType r => String -> r
printf String
argMsg String
name (-Int
info)
      Ordering
GT -> a -> IO a
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return a
forall a. Floating a => a
zero
      Ordering
EQ -> IO a
evaluation


withInfo :: String -> (Ptr CInt -> IO ()) -> IO ()
withInfo :: String -> (Ptr CInt -> IO ()) -> IO ()
withInfo = String -> String -> (Ptr CInt -> IO ()) -> IO ()
Private.withInfo String
diagonalMsg

diagonalMsg :: String
diagonalMsg :: String
diagonalMsg = String
"%d-th diagonal value is zero"