{-# LINE 1 "src/Data/Number/Flint/Support/D/Mat/FFI.hsc" #-}
{-|
module      :  Data.Number.Flint.Support.D.Mat.FFI
copyright   :  (c) 2022 Hartmut Monien
license     :  GNU GPL, version 2 or above (see LICENSE)
maintainer  :  hmonien@uni-bonn.de
-}
module Data.Number.Flint.Support.D.Mat.FFI (
  -- * Double precision matrices
    DMat (..)
  , CDMat (..)
  -- * Smart constructors
  , newDMat
  , withDMat
  , withNewDMat
  -- * Memory management
  , d_mat_init
  , d_mat_clear
  -- * Basic assignment and manipulation
  , d_mat_set
  , d_mat_swap
  , d_mat_swap_entrywise
  , d_mat_entry
  , d_mat_get_entry
  , d_mat_entry_ptr
  , d_mat_zero
  , d_mat_one
  -- * Random matrix generation
  , d_mat_randtest
  -- * Input and output
  , d_mat_get_str
  , d_mat_fprint
  , d_mat_print       
  -- * Comparison
  , d_mat_equal
  , d_mat_approx_equal
  , d_mat_is_zero
  , d_mat_is_approx_zero
  , d_mat_is_empty
  , d_mat_is_square
  -- * Transpose
  , d_mat_transpose
  -- * Matrix multiplication
  , d_mat_mul_classical
  -- * Gram-Schmidt Orthogonalisation and QR Decomposition
  , d_mat_gso
  , d_mat_qr
) where 

-- double precision matrices ---------------------------------------------------

import Foreign.C.Types
import Foreign.C.String
import Foreign.Ptr
import Foreign.ForeignPtr
import Foreign.Storable
import Foreign.Marshal.Array

import Data.Number.Flint.Flint



-- d_mat_t ---------------------------------------------------------------------

data DMat = DMat {-# UNPACK #-} !(ForeignPtr CDMat)
data CDMat = CDMat (Ptr CDouble) CLong CLong (Ptr (Ptr CDouble))

instance Storable CDMat where
  {-# INLINE sizeOf #-}
  sizeOf :: CDMat -> Int
sizeOf CDMat
_ = (Int
32)
{-# LINE 70 "src/Data/Number/Flint/Support/D/Mat/FFI.hsc" #-}
  {-# INLINE alignment #-}
  alignment :: CDMat -> Int
alignment CDMat
_ = Int
8
{-# LINE 72 "src/Data/Number/Flint/Support/D/Mat/FFI.hsc" #-}
  peek ptr = CDMat
    <$> (\hsc_ptr -> peekByteOff hsc_ptr 0) ptr
{-# LINE 74 "src/Data/Number/Flint/Support/D/Mat/FFI.hsc" #-}
    <*> (\hsc_ptr -> peekByteOff hsc_ptr 8) ptr
{-# LINE 75 "src/Data/Number/Flint/Support/D/Mat/FFI.hsc" #-}
    <*> (\hsc_ptr -> peekByteOff hsc_ptr 16) ptr
{-# LINE 76 "src/Data/Number/Flint/Support/D/Mat/FFI.hsc" #-}
    <*> (\hsc_ptr -> peekByteOff hsc_ptr 24) ptr
{-# LINE 77 "src/Data/Number/Flint/Support/D/Mat/FFI.hsc" #-}
  poke = error "CDMat.poke: Not defined."

-- | /newDMat/ /rows/ /cols/
--
-- Construct new `DMat` with /rows/ rows and /cols/ columns.
newDMat :: CLong -> CLong -> IO DMat
newDMat CLong
rows CLong
cols = do
  ForeignPtr CDMat
x <- IO (ForeignPtr CDMat)
forall a. Storable a => IO (ForeignPtr a)
mallocForeignPtr
  ForeignPtr CDMat -> (Ptr CDMat -> IO ()) -> IO ()
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr CDMat
x ((Ptr CDMat -> IO ()) -> IO ()) -> (Ptr CDMat -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \Ptr CDMat
x -> Ptr CDMat -> CLong -> CLong -> IO ()
d_mat_init Ptr CDMat
x CLong
rows CLong
cols
  FinalizerPtr CDMat -> ForeignPtr CDMat -> IO ()
forall a. FinalizerPtr a -> ForeignPtr a -> IO ()
addForeignPtrFinalizer FinalizerPtr CDMat
p_d_mat_clear ForeignPtr CDMat
x
  DMat -> IO DMat
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return (DMat -> IO DMat) -> DMat -> IO DMat
forall a b. (a -> b) -> a -> b
$ ForeignPtr CDMat -> DMat
DMat ForeignPtr CDMat
x

-- | /withDMat/ /mat/ /f/
--
-- Apply /f/ to a /mat/.
{-# INLINE withDMat #-}
withDMat :: DMat -> (Ptr CDMat -> IO a) -> IO (DMat, a)
withDMat (DMat ForeignPtr CDMat
x) Ptr CDMat -> IO a
f = do
  ForeignPtr CDMat -> (Ptr CDMat -> IO (DMat, a)) -> IO (DMat, a)
forall a b. ForeignPtr a -> (Ptr a -> IO b) -> IO b
withForeignPtr ForeignPtr CDMat
x ((Ptr CDMat -> IO (DMat, a)) -> IO (DMat, a))
-> (Ptr CDMat -> IO (DMat, a)) -> IO (DMat, a)
forall a b. (a -> b) -> a -> b
$ \Ptr CDMat
px -> Ptr CDMat -> IO a
f Ptr CDMat
px IO a -> (a -> IO (DMat, a)) -> IO (DMat, a)
forall a b. IO a -> (a -> IO b) -> IO b
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= (DMat, a) -> IO (DMat, a)
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return ((DMat, a) -> IO (DMat, a))
-> (a -> (DMat, a)) -> a -> IO (DMat, a)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (ForeignPtr CDMat -> DMat
DMat ForeignPtr CDMat
x,)

-- | /withNewDMat/ /rows/ /cols/ /f/
--
-- Apply /f/ to a new `DMat` with /rows/ rows and /cols/ columns.
{-# INLINE withNewDMat #-}
withNewDMat :: CLong -> CLong -> (Ptr CDMat -> IO a) -> IO (DMat, a)
withNewDMat CLong
rows CLong
cols Ptr CDMat -> IO a
f = do
  DMat
x <- CLong -> CLong -> IO DMat
newDMat CLong
rows CLong
cols
  DMat -> (Ptr CDMat -> IO a) -> IO (DMat, a)
forall {a}. DMat -> (Ptr CDMat -> IO a) -> IO (DMat, a)
withDMat DMat
x Ptr CDMat -> IO a
f
  
-- Memory management -----------------------------------------------------------

-- | /d_mat_init/ /mat/ /rows/ /cols/ 
-- 
-- Initialises a matrix with the given number of rows and columns for use.
foreign import ccall "d_mat.h d_mat_init"
  d_mat_init :: Ptr CDMat -> CLong -> CLong -> IO ()

-- | /d_mat_clear/ /mat/ 
-- 
-- Clears the given matrix.
foreign import ccall "d_mat.h d_mat_clear"
  d_mat_clear :: Ptr CDMat -> IO ()

foreign import ccall "d_mat.h &d_mat_clear"
  p_d_mat_clear :: FunPtr (Ptr CDMat -> IO ())

-- Basic assignment and manipulation -------------------------------------------

-- | /d_mat_set/ /mat1/ /mat2/ 
-- 
-- Sets @mat1@ to a copy of @mat2@. The dimensions of @mat1@ and @mat2@
-- must be the same.
foreign import ccall "d_mat.h d_mat_set"
  d_mat_set :: Ptr CDMat -> Ptr CDMat -> IO ()

-- | /d_mat_swap/ /mat1/ /mat2/ 
-- 
-- Swaps two matrices. The dimensions of @mat1@ and @mat2@ are allowed to
-- be different.
foreign import ccall "d_mat.h d_mat_swap"
  d_mat_swap :: Ptr CDMat -> Ptr CDMat -> IO ()

-- | /d_mat_swap_entrywise/ /mat1/ /mat2/ 
-- 
-- Swaps two matrices by swapping the individual entries rather than
-- swapping the contents of the structs.
foreign import ccall "d_mat.h d_mat_swap_entrywise"
  d_mat_swap_entrywise :: Ptr CDMat -> Ptr CDMat -> IO ()

-- | /d_mat_entry/ /mat/ /i/ /j/ 
-- 
-- Returns the entry of @mat@ at row \(i\) and column \(j\). Both \(i\) and
-- \(j\) must not exceed the dimensions of the matrix. This function is
-- implemented as a macro.
d_mat_entry :: Ptr CDMat -> CLong -> CLong -> IO CDouble
d_mat_entry :: Ptr CDMat -> CLong -> CLong -> IO CDouble
d_mat_entry Ptr CDMat
mat CLong
i CLong
j = do
  CDMat Ptr CDouble
_ CLong
r CLong
c Ptr (Ptr CDouble)
rows <- Ptr CDMat -> IO CDMat
forall a. Storable a => Ptr a -> IO a
peek Ptr CDMat
mat
  Ptr CDouble
row_i <- Ptr (Ptr CDouble) -> IO (Ptr CDouble)
forall a. Storable a => Ptr a -> IO a
peek (Ptr (Ptr CDouble)
rows Ptr (Ptr CDouble) -> Int -> Ptr (Ptr CDouble)
forall a. Storable a => Ptr a -> Int -> Ptr a
`advancePtr` (CLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral CLong
i))
  CDouble
result <- Ptr CDouble -> IO CDouble
forall a. Storable a => Ptr a -> IO a
peek (Ptr CDouble
row_i Ptr CDouble -> Int -> Ptr CDouble
forall a. Storable a => Ptr a -> Int -> Ptr a
`advancePtr` (CLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral CLong
j))
  CDouble -> IO CDouble
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return CDouble
result
  
-- | /d_mat_get_entry/ /mat/ /i/ /j/ 
-- 
-- Returns the entry of @mat@ at row \(i\) and column \(j\). Both \(i\) and
-- \(j\) must not exceed the dimensions of the matrix.
foreign import ccall "d_mat.h d_mat_get_entry"
  d_mat_get_entry :: Ptr CDMat -> CLong -> CLong -> IO CDouble

-- | /d_mat_entry_ptr/ /mat/ /i/ /j/ 
-- 
-- Returns a pointer to the entry of @mat@ at row \(i\) and column \(j\).
-- Both \(i\) and \(j\) must not exceed the dimensions of the matrix.
foreign import ccall "d_mat.h d_mat_entry_ptr"
  d_mat_entry_ptr :: Ptr CDMat -> CLong -> CLong -> IO (Ptr CDouble)

-- | /d_mat_zero/ /mat/ 
-- 
-- Sets all entries of @mat@ to 0.
foreign import ccall "d_mat.h d_mat_zero"
  d_mat_zero :: Ptr CDMat -> IO ()

-- | /d_mat_one/ /mat/ 
-- 
-- Sets @mat@ to the unit matrix, having ones on the main diagonal and
-- zeroes elsewhere. If @mat@ is nonsquare, it is set to the truncation of
-- a unit matrix.
foreign import ccall "d_mat.h d_mat_one"
  d_mat_one :: Ptr CDMat -> IO ()

-- Random matrix generation ----------------------------------------------------

-- | /d_mat_randtest/ /mat/ /state/ /minexp/ /maxexp/ 
-- 
-- Sets the entries of @mat@ to random signed numbers with exponents
-- between @minexp@ and @maxexp@ or zero.
foreign import ccall "d_mat.h d_mat_randtest"
  d_mat_randtest :: Ptr CDMat -> Ptr CFRandState -> CLong -> CLong -> IO ()

-- Input and output ------------------------------------------------------------

-- | /d_mat_get_str/ /mat/ 
-- 
-- Returns a string representation of the given matrix.
foreign import ccall "d_mat.h d_mat_get_str"
  d_mat_get_str :: Ptr CDMat -> IO CString

-- | /d_mat_fprint/ /file/ /mat/ 
-- 
-- Prints the given matrix to the stream @stdout@.
foreign import ccall "d_mat.h d_mat_fprint"
  d_mat_fprint :: Ptr CFile -> Ptr CDMat -> IO ()

-- | /d_mat_print/ /mat/ 
-- 
-- Prints the given matrix to the stream @stdout@.
d_mat_print :: Ptr CDMat -> IO ()
d_mat_print :: Ptr CDMat -> IO ()
d_mat_print Ptr CDMat
mat = do
  (Ptr CDMat -> IO CString) -> Ptr CDMat -> IO CInt
forall a. (Ptr a -> IO CString) -> Ptr a -> IO CInt
printCStr Ptr CDMat -> IO CString
d_mat_get_str Ptr CDMat
mat
  () -> IO ()
forall a. a -> IO a
forall (m :: * -> *) a. Monad m => a -> m a
return ()

-- Comparison ------------------------------------------------------------------

-- | /d_mat_equal/ /mat1/ /mat2/ 
-- 
-- Returns a non-zero value if @mat1@ and @mat2@ have the same dimensions
-- and entries, and zero otherwise.
foreign import ccall "d_mat.h d_mat_equal"
  d_mat_equal :: Ptr CDMat -> Ptr CDMat -> IO CInt

-- | /d_mat_approx_equal/ /mat1/ /mat2/ /eps/ 
-- 
-- Returns a non-zero value if @mat1@ and @mat2@ have the same dimensions
-- and entries within @eps@ of each other, and zero otherwise.
foreign import ccall "d_mat.h d_mat_approx_equal"
  d_mat_approx_equal :: Ptr CDMat -> Ptr CDMat -> CDouble -> IO CInt

-- | /d_mat_is_zero/ /mat/ 
-- 
-- Returns a non-zero value if all entries @mat@ are zero, and otherwise
-- returns zero.
foreign import ccall "d_mat.h d_mat_is_zero"
  d_mat_is_zero :: Ptr CDMat -> IO CInt

-- | /d_mat_is_approx_zero/ /mat/ /eps/ 
-- 
-- Returns a non-zero value if all entries @mat@ are zero to within @eps@
-- and otherwise returns zero.
foreign import ccall "d_mat.h d_mat_is_approx_zero"
  d_mat_is_approx_zero :: Ptr CDMat -> CDouble -> IO CInt

-- | /d_mat_is_empty/ /mat/ 
-- 
-- Returns a non-zero value if the number of rows or the number of columns
-- in @mat@ is zero, and otherwise returns zero.
foreign import ccall "d_mat.h d_mat_is_empty"
  d_mat_is_empty :: Ptr CDMat -> IO CInt

-- | /d_mat_is_square/ /mat/ 
-- 
-- Returns a non-zero value if the number of rows is equal to the number of
-- columns in @mat@, and otherwise returns zero.
foreign import ccall "d_mat.h d_mat_is_square"
  d_mat_is_square :: Ptr CDMat -> IO CInt

-- Transpose -------------------------------------------------------------------

-- | /d_mat_transpose/ /B/ /A/ 
-- 
-- Sets \(B\) to \(A^T\), the transpose of \(A\). Dimensions must be
-- compatible. \(A\) and \(B\) are allowed to be the same object if \(A\)
-- is a square matrix.
foreign import ccall "d_mat.h d_mat_transpose"
  d_mat_transpose :: Ptr CDMat -> Ptr CDMat -> IO ()

-- Matrix multiplication -------------------------------------------------------

-- | /d_mat_mul_classical/ /C/ /A/ /B/ 
-- 
-- Sets @C@ to the matrix product \(C = A B\). The matrices must have
-- compatible dimensions for matrix multiplication (an exception is raised
-- otherwise). Aliasing is allowed.
foreign import ccall "d_mat.h d_mat_mul_classical"
  d_mat_mul_classical :: Ptr CDMat -> Ptr CDMat -> Ptr CDMat -> IO ()

-- Gram-Schmidt Orthogonalisation and QR Decomposition -------------------------

-- | /d_mat_gso/ /B/ /A/ 
-- 
-- Takes a subset of \(R^m\) \(S = {a_1, a_2, \ldots, a_n}\) (as the
-- columns of a \(m x n\) matrix @A@) and generates an orthonormal set
-- \(S' = {b_1, b_2, \ldots, b_n}\) (as the columns of the \(m x n\) matrix
-- @B@) that spans the same subspace of \(R^m\) as \(S\).
-- 
-- This uses an algorithm of Schwarz-Rutishauser. See pp. 9 of
-- <https://people.inf.ethz.ch/gander/papers/qrneu.pdf>
foreign import ccall "d_mat.h d_mat_gso"
  d_mat_gso :: Ptr CDMat -> Ptr CDMat -> IO ()

-- | /d_mat_qr/ /Q/ /R/ /A/ 
-- 
-- Computes the \(QR\) decomposition of a matrix @A@ using the Gram-Schmidt
-- process. (Sets @Q@ and @R@ such that \(A = QR\) where @R@ is an upper
-- triangular matrix and @Q@ is an orthogonal matrix.)
-- 
-- This uses an algorithm of Schwarz-Rutishauser. See pp. 9 of
-- <https://people.inf.ethz.ch/gander/papers/qrneu.pdf>
foreign import ccall "d_mat.h d_mat_qr"
  d_mat_qr :: Ptr CDMat -> Ptr CDMat -> Ptr CDMat -> IO ()