module Numeric.LAPACK.Linear.Private where

import qualified Numeric.LAPACK.Matrix.Shape.Private as MatrixShape
import qualified Numeric.LAPACK.Matrix.Extent.Private as Extent
import qualified Numeric.LAPACK.Private as Private
import Numeric.LAPACK.Matrix.Shape.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.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 vert horiz height width a ->
   Full vert horiz height width a
solver :: 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 ()
f (Array (MatrixShape.Full Order
order Extent vert horiz height width
extent) ForeignPtr a
b) =
   Full vert horiz height width
-> (Ptr a -> IO ()) -> Full vert horiz height width a
forall sh a.
(C sh, Storable a) =>
sh -> (Ptr a -> IO ()) -> Array sh a
Array.unsafeCreate (Order
-> Extent vert horiz height width -> Full vert horiz height width
forall vert horiz height width.
Order
-> Extent vert horiz height width -> Full vert horiz height width
MatrixShape.Full Order
ColumnMajor Extent vert horiz height width
extent) ((Ptr a -> IO ()) -> Full vert horiz height width a)
-> (Ptr a -> IO ()) -> Full vert horiz height width a
forall a b. (a -> b) -> a -> b
$
      \Ptr a
xPtr -> do
   let (height
height,width
width) = Extent vert horiz height width -> (height, width)
forall vert horiz height width.
(C vert, C horiz) =>
Extent vert horiz height width -> (height, width)
Extent.dimensions Extent 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 (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 :: 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 (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"