{-
   Copyright 2016, Dominic Orchard, Andrew Rice, Mistral Contrastin, Matthew Danish

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
-}

{-
  Units of measure extension to Fortran: backend
-}

{-# LANGUAGE TupleSections #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ForeignFunctionInterface #-}

module Camfort.Specification.Units.InferenceBackendFlint where

import           Control.Monad
import           Data.List (partition)
-- import           Debug.Trace (trace, traceM, traceShowM)
import           Foreign
import           Foreign.C.Types
import           Numeric.LinearAlgebra (atIndex, rows, cols)
import qualified Numeric.LinearAlgebra as H
import           System.IO.Unsafe (unsafePerformIO)

foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_init" fmpz_mat_init :: Ptr FMPZMat -> CLong -> CLong -> IO ()
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_set" fmpz_mat_set :: Ptr FMPZMat -> Ptr FMPZMat -> IO ()
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_entry" fmpz_mat_entry :: Ptr FMPZMat -> CLong -> CLong -> IO (Ptr CLong)
foreign import ccall unsafe "flint/fmpz.h fmpz_set_si" fmpz_set_si :: Ptr CLong -> CLong -> IO ()
foreign import ccall unsafe "flint/fmpz.h fmpz_get_si" fmpz_get_si :: Ptr CLong -> IO CLong
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_clear" fmpz_mat_clear :: Ptr FMPZMat -> IO ()
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_print_pretty" fmpz_mat_print_pretty :: Ptr FMPZMat -> IO ()
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_mul" fmpz_mat_mul :: Ptr FMPZMat -> Ptr FMPZMat -> Ptr FMPZMat -> IO ()

-- fmpz_mat_window_init(fmpz_mat_t window, const fmpz_mat_t mat, slong r1, slong c1, slong r2, slong c2)
--
-- Initializes the matrix window to be an r2 - r1 by c2 - c1 submatrix
-- of mat whose (0,0) entry is the (r1, c1) entry of mat. The memory
-- for the elements of window is shared with mat.
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_window_init" fmpz_mat_window_init :: Ptr FMPZMat -> Ptr FMPZMat -> CLong -> CLong -> CLong -> CLong -> IO ()

-- Frees the window (leaving underlying matrix alone).
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_window_clear" fmpz_mat_window_clear :: Ptr FMPZMat -> IO ()

-- r <- fmp_mat_rref B den A
--
-- Uses fraction-free Gauss-Jordan elimination to set (B, den) to the
-- reduced row echelon form of A and returns the rank of A. Aliasing
-- of A and B is allowed. r is rank of A.
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_rref" fmpz_mat_rref :: Ptr FMPZMat -> Ptr CLong -> Ptr FMPZMat -> IO CLong

-- r <- fmp_mat_inv B den A
--
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_inv" fmpz_mat_inv :: Ptr FMPZMat -> Ptr CLong -> Ptr FMPZMat -> IO CLong

-- fmpz_mat_hnf H A
--
-- H is the Hermite Normal Form of A.
foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_hnf" fmpz_mat_hnf :: Ptr FMPZMat -> Ptr FMPZMat -> IO ()

foreign import ccall unsafe "flint/fmpz_mat.h fmpz_mat_rank" fmpz_mat_rank :: Ptr FMPZMat -> IO CLong

data FMPZMat

instance Storable FMPZMat where
  sizeOf :: FMPZMat -> Int
sizeOf FMPZMat
_ = Int
4 Int -> Int -> Int
forall a. Num a => a -> a -> a
* CLong -> Int
forall a. Storable a => a -> Int
sizeOf (CLong
forall a. HasCallStack => a
undefined :: CLong)
  alignment :: FMPZMat -> Int
alignment FMPZMat
_ = CLong -> Int
forall a. Storable a => a -> Int
alignment (CLong
forall a. HasCallStack => a
undefined :: CLong)
  peek :: Ptr FMPZMat -> IO FMPZMat
peek Ptr FMPZMat
_ = IO FMPZMat
forall a. HasCallStack => a
undefined
  poke :: Ptr FMPZMat -> FMPZMat -> IO ()
poke Ptr FMPZMat
_ = FMPZMat -> IO ()
forall a. HasCallStack => a
undefined

-- testFlint :: IO ()
-- testFlint = do
--   traceM "***********************8 testFlint 8***********************"
--   alloca $ \ a -> do
--     alloca $ \ b -> do
--       alloca $ \ den -> do
--         let n = 10
--         fmpz_mat_init a n n
--         fmpz_mat_init b n n
--         forM_ [0..n-1] $ \ i -> do
--           forM_ [0..n-1] $ \ j -> do
--             e <- fmpz_mat_entry a i j
--             fmpz_set_si e (2 * i + j)

--         -- fmpz_mat_mul b a a
--         r <- fmpz_mat_rref b den a
--         fmpz_mat_print_pretty a
--         fmpz_mat_print_pretty b
--         d <- peek den
--         traceM $ "r = " ++ show r ++ " den = " ++ show d
--         fmpz_mat_hnf b a
--         fmpz_mat_print_pretty b
--         fmpz_mat_clear a
--         fmpz_mat_clear b

rref :: H.Matrix Double -> (H.Matrix Double, Int, Int)
rref :: Matrix Double -> (Matrix Double, Int, Int)
rref Matrix Double
m = IO (Matrix Double, Int, Int) -> (Matrix Double, Int, Int)
forall a. IO a -> a
unsafePerformIO (IO (Matrix Double, Int, Int) -> (Matrix Double, Int, Int))
-> IO (Matrix Double, Int, Int) -> (Matrix Double, Int, Int)
forall a b. (a -> b) -> a -> b
$ do
  (Ptr FMPZMat -> IO (Matrix Double, Int, Int))
-> IO (Matrix Double, Int, Int)
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO (Matrix Double, Int, Int))
 -> IO (Matrix Double, Int, Int))
-> (Ptr FMPZMat -> IO (Matrix Double, Int, Int))
-> IO (Matrix Double, Int, Int)
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
outputM -> do
    (Ptr FMPZMat -> IO (Matrix Double, Int, Int))
-> IO (Matrix Double, Int, Int)
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO (Matrix Double, Int, Int))
 -> IO (Matrix Double, Int, Int))
-> (Ptr FMPZMat -> IO (Matrix Double, Int, Int))
-> IO (Matrix Double, Int, Int)
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
inputM -> do
      (Ptr CLong -> IO (Matrix Double, Int, Int))
-> IO (Matrix Double, Int, Int)
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr CLong -> IO (Matrix Double, Int, Int))
 -> IO (Matrix Double, Int, Int))
-> (Ptr CLong -> IO (Matrix Double, Int, Int))
-> IO (Matrix Double, Int, Int)
forall a b. (a -> b) -> a -> b
$ \ Ptr CLong
den -> do
        let numRows :: CLong
numRows = Int -> CLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CLong) -> Int -> CLong
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
m
        let numCols :: CLong
numCols = Int -> CLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CLong) -> Int -> CLong
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
m
        Ptr FMPZMat -> CLong -> CLong -> IO ()
fmpz_mat_init Ptr FMPZMat
outputM CLong
numRows CLong
numCols
        Ptr FMPZMat -> CLong -> CLong -> IO ()
fmpz_mat_init Ptr FMPZMat
inputM CLong
numRows CLong
numCols
        [CLong] -> (CLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [CLong
0..CLong
numRowsCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO ()) -> IO ()) -> (CLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ CLong
i -> do
          [CLong] -> (CLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [CLong
0..CLong
numColsCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO ()) -> IO ()) -> (CLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ CLong
j -> do
            Ptr CLong
e <- Ptr FMPZMat -> CLong -> CLong -> IO (Ptr CLong)
fmpz_mat_entry Ptr FMPZMat
inputM CLong
i CLong
j
            Ptr CLong -> CLong -> IO ()
fmpz_set_si Ptr CLong
e (Double -> CLong
forall a b. (RealFrac a, Integral b) => a -> b
floor (Matrix Double
m Matrix Double -> IndexOf Matrix -> Double
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
`atIndex` (CLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral CLong
i, CLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral CLong
j)))
        CLong
r <- Ptr FMPZMat -> Ptr CLong -> Ptr FMPZMat -> IO CLong
fmpz_mat_rref Ptr FMPZMat
outputM Ptr CLong
den Ptr FMPZMat
inputM
        CLong
d <- Ptr CLong -> IO CLong
forall a. Storable a => Ptr a -> IO a
peek Ptr CLong
den

        -- DEBUG:
        -- fmpz_mat_print_pretty outputM
        -- traceM $ "r = " ++ show r ++ " den = " ++ show d
        --

        [[Double]]
lists <- [CLong] -> (CLong -> IO [Double]) -> IO [[Double]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [CLong
0..CLong
numRowsCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO [Double]) -> IO [[Double]])
-> (CLong -> IO [Double]) -> IO [[Double]]
forall a b. (a -> b) -> a -> b
$ \ CLong
i -> do
          [CLong] -> (CLong -> IO Double) -> IO [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [CLong
0..CLong
numColsCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO Double) -> IO [Double])
-> (CLong -> IO Double) -> IO [Double]
forall a b. (a -> b) -> a -> b
$ \ CLong
j -> do
            Ptr CLong
e <- Ptr FMPZMat -> CLong -> CLong -> IO (Ptr CLong)
fmpz_mat_entry Ptr FMPZMat
outputM CLong
i CLong
j
            CLong -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (CLong -> Double) -> IO CLong -> IO Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
`fmap` Ptr CLong -> IO CLong
fmpz_get_si Ptr CLong
e
        let m' :: Matrix Double
m' = [[Double]] -> Matrix Double
forall t. Element t => [[t]] -> Matrix t
H.fromLists [[Double]]
lists
        Ptr FMPZMat -> IO ()
fmpz_mat_clear Ptr FMPZMat
inputM
        Ptr FMPZMat -> IO ()
fmpz_mat_clear Ptr FMPZMat
outputM
        (Matrix Double, Int, Int) -> IO (Matrix Double, Int, Int)
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix Double
m', CLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral CLong
d, CLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral CLong
r)

hnf :: H.Matrix Double -> H.Matrix Double
hnf :: Matrix Double -> Matrix Double
hnf Matrix Double
m = IO (Matrix Double) -> Matrix Double
forall a. IO a -> a
unsafePerformIO (IO (Matrix Double) -> Matrix Double)
-> IO (Matrix Double) -> Matrix Double
forall a b. (a -> b) -> a -> b
$ do
  (Ptr FMPZMat -> IO (Matrix Double)) -> IO (Matrix Double)
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO (Matrix Double)) -> IO (Matrix Double))
-> (Ptr FMPZMat -> IO (Matrix Double)) -> IO (Matrix Double)
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
outputM -> do
    (Ptr FMPZMat -> IO (Matrix Double)) -> IO (Matrix Double)
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO (Matrix Double)) -> IO (Matrix Double))
-> (Ptr FMPZMat -> IO (Matrix Double)) -> IO (Matrix Double)
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
inputM -> do
      let numRows :: CLong
numRows = Int -> CLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CLong) -> Int -> CLong
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
m
      let numCols :: CLong
numCols = Int -> CLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CLong) -> Int -> CLong
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
m
      Ptr FMPZMat -> CLong -> CLong -> IO ()
fmpz_mat_init Ptr FMPZMat
outputM CLong
numRows CLong
numCols
      Ptr FMPZMat -> CLong -> CLong -> IO ()
fmpz_mat_init Ptr FMPZMat
inputM CLong
numRows CLong
numCols
      [CLong] -> (CLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [CLong
0..CLong
numRowsCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO ()) -> IO ()) -> (CLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ CLong
i -> do
        [CLong] -> (CLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [CLong
0..CLong
numColsCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO ()) -> IO ()) -> (CLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ CLong
j -> do
          Ptr CLong
e <- Ptr FMPZMat -> CLong -> CLong -> IO (Ptr CLong)
fmpz_mat_entry Ptr FMPZMat
inputM CLong
i CLong
j
          Ptr CLong -> CLong -> IO ()
fmpz_set_si Ptr CLong
e (Double -> CLong
forall a b. (RealFrac a, Integral b) => a -> b
floor (Matrix Double
m Matrix Double -> IndexOf Matrix -> Double
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
`atIndex` (CLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral CLong
i, CLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral CLong
j)))
      Ptr FMPZMat -> Ptr FMPZMat -> IO ()
fmpz_mat_hnf Ptr FMPZMat
outputM Ptr FMPZMat
inputM
      CLong
r <- Ptr FMPZMat -> IO CLong
fmpz_mat_rank Ptr FMPZMat
outputM

      -- DEBUG:
      -- fmpz_mat_print_pretty outputM
      -- traceM $ "rank = " ++ show r
      --

      [[Double]]
lists <- [CLong] -> (CLong -> IO [Double]) -> IO [[Double]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [CLong
0..CLong -> CLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral CLong
rCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO [Double]) -> IO [[Double]])
-> (CLong -> IO [Double]) -> IO [[Double]]
forall a b. (a -> b) -> a -> b
$ \ CLong
i -> do
        [CLong] -> (CLong -> IO Double) -> IO [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [CLong
0..CLong
numColsCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO Double) -> IO [Double])
-> (CLong -> IO Double) -> IO [Double]
forall a b. (a -> b) -> a -> b
$ \ CLong
j -> do
          Ptr CLong
e <- Ptr FMPZMat -> CLong -> CLong -> IO (Ptr CLong)
fmpz_mat_entry Ptr FMPZMat
outputM CLong
i CLong
j
          CLong -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (CLong -> Double) -> IO CLong -> IO Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
`fmap` Ptr CLong -> IO CLong
fmpz_get_si Ptr CLong
e
      let m' :: Matrix Double
m' = [[Double]] -> Matrix Double
forall t. Element t => [[t]] -> Matrix t
H.fromLists [[Double]]
lists
      Ptr FMPZMat -> IO ()
fmpz_mat_clear Ptr FMPZMat
inputM
      Ptr FMPZMat -> IO ()
fmpz_mat_clear Ptr FMPZMat
outputM
      Matrix Double -> IO (Matrix Double)
forall (m :: * -> *) a. Monad m => a -> m a
return Matrix Double
m'

inv :: H.Matrix Double -> Maybe (H.Matrix Double, Int)
inv :: Matrix Double -> Maybe (Matrix Double, Int)
inv Matrix Double
m = IO (Maybe (Matrix Double, Int)) -> Maybe (Matrix Double, Int)
forall a. IO a -> a
unsafePerformIO (IO (Maybe (Matrix Double, Int)) -> Maybe (Matrix Double, Int))
-> IO (Maybe (Matrix Double, Int)) -> Maybe (Matrix Double, Int)
forall a b. (a -> b) -> a -> b
$ do
  (Ptr FMPZMat -> IO (Maybe (Matrix Double, Int)))
-> IO (Maybe (Matrix Double, Int))
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO (Maybe (Matrix Double, Int)))
 -> IO (Maybe (Matrix Double, Int)))
-> (Ptr FMPZMat -> IO (Maybe (Matrix Double, Int)))
-> IO (Maybe (Matrix Double, Int))
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
outputM -> do
    (Ptr FMPZMat -> IO (Maybe (Matrix Double, Int)))
-> IO (Maybe (Matrix Double, Int))
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO (Maybe (Matrix Double, Int)))
 -> IO (Maybe (Matrix Double, Int)))
-> (Ptr FMPZMat -> IO (Maybe (Matrix Double, Int)))
-> IO (Maybe (Matrix Double, Int))
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
inputM -> do
      (Ptr CLong -> IO (Maybe (Matrix Double, Int)))
-> IO (Maybe (Matrix Double, Int))
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr CLong -> IO (Maybe (Matrix Double, Int)))
 -> IO (Maybe (Matrix Double, Int)))
-> (Ptr CLong -> IO (Maybe (Matrix Double, Int)))
-> IO (Maybe (Matrix Double, Int))
forall a b. (a -> b) -> a -> b
$ \ Ptr CLong
den -> do
        let numRows :: CLong
numRows = Int -> CLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CLong) -> Int -> CLong
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
m
        let numCols :: CLong
numCols = Int -> CLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CLong) -> Int -> CLong
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
m
        Ptr FMPZMat -> CLong -> CLong -> IO ()
fmpz_mat_init Ptr FMPZMat
outputM CLong
numRows CLong
numCols
        Ptr FMPZMat -> CLong -> CLong -> IO ()
fmpz_mat_init Ptr FMPZMat
inputM CLong
numRows CLong
numCols
        [CLong] -> (CLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [CLong
0..CLong
numRowsCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO ()) -> IO ()) -> (CLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ CLong
i -> do
          [CLong] -> (CLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [CLong
0..CLong
numColsCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO ()) -> IO ()) -> (CLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ CLong
j -> do
            Ptr CLong
e <- Ptr FMPZMat -> CLong -> CLong -> IO (Ptr CLong)
fmpz_mat_entry Ptr FMPZMat
inputM CLong
i CLong
j
            Ptr CLong -> CLong -> IO ()
fmpz_set_si Ptr CLong
e (Double -> CLong
forall a b. (RealFrac a, Integral b) => a -> b
floor (Matrix Double
m Matrix Double -> IndexOf Matrix -> Double
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
`atIndex` (CLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral CLong
i, CLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral CLong
j)))
        CLong
r <- Ptr FMPZMat -> Ptr CLong -> Ptr FMPZMat -> IO CLong
fmpz_mat_inv Ptr FMPZMat
outputM Ptr CLong
den Ptr FMPZMat
inputM
        CLong
d <- Ptr CLong -> IO CLong
forall a. Storable a => Ptr a -> IO a
peek Ptr CLong
den

        -- DEBUG:
        Ptr FMPZMat -> IO ()
fmpz_mat_print_pretty Ptr FMPZMat
outputM
        -- traceM $ "r = " ++ show r ++ " den = " ++ show d
        --

        [[Double]]
lists <- [CLong] -> (CLong -> IO [Double]) -> IO [[Double]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [CLong
0..CLong
numRowsCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO [Double]) -> IO [[Double]])
-> (CLong -> IO [Double]) -> IO [[Double]]
forall a b. (a -> b) -> a -> b
$ \ CLong
i -> do
          [CLong] -> (CLong -> IO Double) -> IO [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [CLong
0..CLong
numColsCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO Double) -> IO [Double])
-> (CLong -> IO Double) -> IO [Double]
forall a b. (a -> b) -> a -> b
$ \ CLong
j -> do
            Ptr CLong
e <- Ptr FMPZMat -> CLong -> CLong -> IO (Ptr CLong)
fmpz_mat_entry Ptr FMPZMat
outputM CLong
i CLong
j
            CLong -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (CLong -> Double) -> IO CLong -> IO Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
`fmap` Ptr CLong -> IO CLong
fmpz_get_si Ptr CLong
e
        let m' :: Matrix Double
m' = [[Double]] -> Matrix Double
forall t. Element t => [[t]] -> Matrix t
H.fromLists [[Double]]
lists
        Ptr FMPZMat -> IO ()
fmpz_mat_clear Ptr FMPZMat
inputM
        Ptr FMPZMat -> IO ()
fmpz_mat_clear Ptr FMPZMat
outputM
        if CLong
r CLong -> CLong -> Bool
forall a. Eq a => a -> a -> Bool
== CLong
1 then
          Maybe (Matrix Double, Int) -> IO (Maybe (Matrix Double, Int))
forall (m :: * -> *) a. Monad m => a -> m a
return (Maybe (Matrix Double, Int) -> IO (Maybe (Matrix Double, Int)))
-> Maybe (Matrix Double, Int) -> IO (Maybe (Matrix Double, Int))
forall a b. (a -> b) -> a -> b
$ (Matrix Double, Int) -> Maybe (Matrix Double, Int)
forall a. a -> Maybe a
Just (Matrix Double
m', CLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral CLong
d)
        else
          Maybe (Matrix Double, Int) -> IO (Maybe (Matrix Double, Int))
forall (m :: * -> *) a. Monad m => a -> m a
return Maybe (Matrix Double, Int)
forall a. Maybe a
Nothing

withMatrix :: H.Matrix Double -> ((CLong, CLong, Ptr FMPZMat) -> IO b) -> IO b
withMatrix :: Matrix Double -> ((CLong, CLong, Ptr FMPZMat) -> IO b) -> IO b
withMatrix Matrix Double
m (CLong, CLong, Ptr FMPZMat) -> IO b
f = do
  (Ptr FMPZMat -> IO b) -> IO b
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO b) -> IO b) -> (Ptr FMPZMat -> IO b) -> IO b
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
outputM -> do
    let numRows :: CLong
numRows = Int -> CLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CLong) -> Int -> CLong
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
m
    let numCols :: CLong
numCols = Int -> CLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> CLong) -> Int -> CLong
forall a b. (a -> b) -> a -> b
$ Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
m
    Ptr FMPZMat -> CLong -> CLong -> IO ()
fmpz_mat_init Ptr FMPZMat
outputM CLong
numRows CLong
numCols
    [CLong] -> (CLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [CLong
0 .. CLong
numRows CLong -> CLong -> CLong
forall a. Num a => a -> a -> a
- CLong
1] ((CLong -> IO ()) -> IO ()) -> (CLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ CLong
i ->
      [CLong] -> (CLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [CLong
0 .. CLong
numCols CLong -> CLong -> CLong
forall a. Num a => a -> a -> a
- CLong
1] ((CLong -> IO ()) -> IO ()) -> (CLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ CLong
j -> do
        Ptr CLong
e <- Ptr FMPZMat -> CLong -> CLong -> IO (Ptr CLong)
fmpz_mat_entry Ptr FMPZMat
outputM CLong
i CLong
j
        Ptr CLong -> CLong -> IO ()
fmpz_set_si Ptr CLong
e (Double -> CLong
forall a b. (RealFrac a, Integral b) => a -> b
floor (Matrix Double
m Matrix Double -> IndexOf Matrix -> Double
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
`atIndex` (CLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral CLong
i, CLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral CLong
j)) :: CLong)
    b
x <- (CLong, CLong, Ptr FMPZMat) -> IO b
f (CLong
numRows, CLong
numCols, Ptr FMPZMat
outputM)
    Ptr FMPZMat -> IO ()
fmpz_mat_clear Ptr FMPZMat
outputM
    b -> IO b
forall (m :: * -> *) a. Monad m => a -> m a
return b
x

withBlankMatrix :: CLong -> CLong -> (Ptr FMPZMat -> IO b) -> IO b
withBlankMatrix :: CLong -> CLong -> (Ptr FMPZMat -> IO b) -> IO b
withBlankMatrix CLong
numRows CLong
numCols Ptr FMPZMat -> IO b
f = do
  (Ptr FMPZMat -> IO b) -> IO b
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO b) -> IO b) -> (Ptr FMPZMat -> IO b) -> IO b
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
outputM -> do
    Ptr FMPZMat -> CLong -> CLong -> IO ()
fmpz_mat_init Ptr FMPZMat
outputM (CLong -> CLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral CLong
numRows) (CLong -> CLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral CLong
numCols)
    b
x <- Ptr FMPZMat -> IO b
f Ptr FMPZMat
outputM
    Ptr FMPZMat -> IO ()
fmpz_mat_clear Ptr FMPZMat
outputM
    b -> IO b
forall (m :: * -> *) a. Monad m => a -> m a
return b
x

withWindow :: Ptr FMPZMat -> CLong -> CLong -> CLong -> CLong -> (Ptr FMPZMat -> IO b) -> IO b
withWindow :: Ptr FMPZMat
-> CLong
-> CLong
-> CLong
-> CLong
-> (Ptr FMPZMat -> IO b)
-> IO b
withWindow Ptr FMPZMat
underM CLong
r1 CLong
c1 CLong
r2 CLong
c2 Ptr FMPZMat -> IO b
f = do
  (Ptr FMPZMat -> IO b) -> IO b
forall a b. Storable a => (Ptr a -> IO b) -> IO b
alloca ((Ptr FMPZMat -> IO b) -> IO b) -> (Ptr FMPZMat -> IO b) -> IO b
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
window -> do
    Ptr FMPZMat
-> Ptr FMPZMat -> CLong -> CLong -> CLong -> CLong -> IO ()
fmpz_mat_window_init Ptr FMPZMat
window Ptr FMPZMat
underM CLong
r1 CLong
c1 CLong
r2 CLong
c2
    b
x <- Ptr FMPZMat -> IO b
f Ptr FMPZMat
window
    Ptr FMPZMat -> IO ()
fmpz_mat_window_clear Ptr FMPZMat
window
    b -> IO b
forall (m :: * -> *) a. Monad m => a -> m a
return b
x

flintToHMatrix :: CLong -> CLong -> Ptr FMPZMat -> IO (H.Matrix Double)
flintToHMatrix :: CLong -> CLong -> Ptr FMPZMat -> IO (Matrix Double)
flintToHMatrix CLong
numRows CLong
numCols Ptr FMPZMat
flintM = do
  [[Double]]
lists <- [CLong] -> (CLong -> IO [Double]) -> IO [[Double]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [CLong
0..CLong
numRowsCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO [Double]) -> IO [[Double]])
-> (CLong -> IO [Double]) -> IO [[Double]]
forall a b. (a -> b) -> a -> b
$ \ CLong
i -> do
    [CLong] -> (CLong -> IO Double) -> IO [Double]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [CLong
0..CLong
numColsCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO Double) -> IO [Double])
-> (CLong -> IO Double) -> IO [Double]
forall a b. (a -> b) -> a -> b
$ \ CLong
j -> do
      Ptr CLong
e <- Ptr FMPZMat -> CLong -> CLong -> IO (Ptr CLong)
fmpz_mat_entry Ptr FMPZMat
flintM CLong
i CLong
j
      CLong -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (CLong -> Double) -> IO CLong -> IO Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
`fmap` Ptr CLong -> IO CLong
fmpz_get_si Ptr CLong
e
  Matrix Double -> IO (Matrix Double)
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix Double -> IO (Matrix Double))
-> Matrix Double -> IO (Matrix Double)
forall a b. (a -> b) -> a -> b
$ [[Double]] -> Matrix Double
forall t. Element t => [[t]] -> Matrix t
H.fromLists [[Double]]
lists

pokeM :: Ptr FMPZMat -> CLong -> CLong -> CLong -> IO ()
pokeM :: Ptr FMPZMat -> CLong -> CLong -> CLong -> IO ()
pokeM Ptr FMPZMat
flintM CLong
i CLong
j CLong
v = do
  Ptr CLong
e <- Ptr FMPZMat -> CLong -> CLong -> IO (Ptr CLong)
fmpz_mat_entry Ptr FMPZMat
flintM CLong
i CLong
j
  Ptr CLong -> CLong -> IO ()
fmpz_set_si Ptr CLong
e CLong
v

peekM :: Ptr FMPZMat -> CLong -> CLong -> IO CLong
peekM :: Ptr FMPZMat -> CLong -> CLong -> IO CLong
peekM Ptr FMPZMat
flintM CLong
i CLong
j = do
  Ptr CLong
e <- Ptr FMPZMat -> CLong -> CLong -> IO (Ptr CLong)
fmpz_mat_entry Ptr FMPZMat
flintM CLong
i CLong
j
  Ptr CLong -> IO CLong
fmpz_get_si Ptr CLong
e

copyMatrix :: Ptr FMPZMat -> Ptr FMPZMat -> CLong -> CLong -> CLong -> CLong -> IO ()
copyMatrix :: Ptr FMPZMat
-> Ptr FMPZMat -> CLong -> CLong -> CLong -> CLong -> IO ()
copyMatrix Ptr FMPZMat
m1 Ptr FMPZMat
m2 CLong
r1 CLong
c1 CLong
r2 CLong
c2 =
  [CLong] -> (CLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [CLong
r1 .. CLong
r2CLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO ()) -> IO ()) -> (CLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ CLong
i ->
    [CLong] -> (CLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [CLong
c1 .. CLong
c2CLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO ()) -> IO ()) -> (CLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ CLong
j ->
      Ptr FMPZMat -> CLong -> CLong -> IO CLong
peekM Ptr FMPZMat
m2 CLong
i CLong
j IO CLong -> (CLong -> IO ()) -> IO ()
forall (m :: * -> *) a b. Monad m => m a -> (a -> m b) -> m b
>>= Ptr FMPZMat -> CLong -> CLong -> CLong -> IO ()
pokeM Ptr FMPZMat
m1 CLong
i CLong
j


-- copyMatrix' m1 m2 r1 c1 r2 c2 =
--   withWindow m2 r1 c1 r2 c2 $ \ w2 ->
--     withWindow m1 r1 c1 r2 c2 $ \ w1 -> do
--       fmpz_mat_set w1 w2

--------------------------------------------------
-- 'normalising' Hermite Normal Form, for lack of a better name
--
-- The problem with HNF is that it will happily return a matrix where
-- the leading-coefficient diagonal contains integers > 1.
--
-- In some cases the leading-coefficient does not divide some of the
-- remaining numbers in the row.
--
-- This corresponds to the case where one of the solutions would be
-- fractional if we were dealing with rational matrices.
--
-- But we don't want fractional solutions. We need to bump the matrix
-- so that this situation does not arise.
--
-- This tends to happen due to implicit polymorphism.
--
-- We do this by adding another column that copies the column with the
-- problematic leading-coefficient.
--
-- We then set up a new row expressing the constraint that the old
-- column is equal to the new column, 1:1.
--
-- This prevents a 'fractional' solution.
--
-- Along the way we look for leading-coefficients > 1 that do divide
-- their entire row. Then we simply apply elementary row scaling to
-- the row so that the leading-coefficients is 1.
--
-- The result is a matrix where the leading-coefficients of all the
-- rows that matter are 1. It's not quite a 'reduced' form though
-- because it can be bigger than the original matrix.
--
-- As a result we also return a list of columns that were cloned this
-- way.
--
-- The resulting matrix can have non-zeroes above the
-- leading-coefficients in the same column, so I apply some elementary
-- row operations to fix that and bring things into RREF.
--
-- Running time is computation of Hermite Normal Form twice, plus
-- construction of a slightly larger matrix (possibly), plus scanning
-- through the matrix for leading-coefficients, and then again to
-- divide them out sometimes.

normHNF :: H.Matrix Double -> (H.Matrix Double, [Int])
normHNF :: Matrix Double -> (Matrix Double, [Int])
normHNF Matrix Double
m
  | (Matrix Double -> Int
forall t. Matrix t -> Int
rows Matrix Double
m, Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
m) (Int, Int) -> (Int, Int) -> Bool
forall a. Eq a => a -> a -> Bool
== (Int
1, Int
1) = (Int -> Matrix Double
forall a. (Num a, Element a) => Int -> Matrix a
H.ident (if Matrix Double -> IndexOf Matrix -> Double
forall (c :: * -> *) e. Container c e => c e -> IndexOf c -> e
H.atIndex Matrix Double
m (Int
0, Int
0) Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
/= Double
0 then Int
1 else Int
0), [])
  | Bool
otherwise = (Matrix Double
m', [Int]
indices)
  where
    numCols :: Int
numCols = Matrix Double -> Int
forall t. Matrix t -> Int
cols Matrix Double
m
    indexLookup :: Int -> Int
indexLookup Int
j | Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
numCols = Int
j
                  | Bool
otherwise = [Int]
indices [Int] -> Int -> Int
forall a. [a] -> Int -> a
!! (Int
j Int -> Int -> Int
forall a. Integral a => a -> a -> a
`mod` Int
numCols)
    indices :: [Int]
indices = (Int -> Int) -> [Int] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Int
indexLookup ([Int] -> [Int]) -> [Int] -> [Int]
forall a b. (a -> b) -> a -> b
$ ((Matrix Double, [Int]) -> [Int])
-> [(Matrix Double, [Int])] -> [Int]
forall (t :: * -> *) a b. Foldable t => (a -> [b]) -> t a -> [b]
concatMap (Matrix Double, [Int]) -> [Int]
forall a b. (a, b) -> b
snd [(Matrix Double, [Int])]
rs1
    ([(Matrix Double, [Int])]
rs1, (Matrix Double
m',[]):[(Matrix Double, [Int])]
_) = ((Matrix Double, [Int]) -> Bool)
-> [(Matrix Double, [Int])]
-> ([(Matrix Double, [Int])], [(Matrix Double, [Int])])
forall a. (a -> Bool) -> [a] -> ([a], [a])
break ([Int] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null ([Int] -> Bool)
-> ((Matrix Double, [Int]) -> [Int])
-> (Matrix Double, [Int])
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Matrix Double, [Int]) -> [Int]
forall a b. (a, b) -> b
snd) [(Matrix Double, [Int])]
results
    results :: [(Matrix Double, [Int])]
results = [(Matrix Double, [Int])] -> [(Matrix Double, [Int])]
forall a. [a] -> [a]
tail ([(Matrix Double, [Int])] -> [(Matrix Double, [Int])])
-> [(Matrix Double, [Int])] -> [(Matrix Double, [Int])]
forall a b. (a -> b) -> a -> b
$ ((Matrix Double, [Int]) -> (Matrix Double, [Int]))
-> (Matrix Double, [Int]) -> [(Matrix Double, [Int])]
forall a. (a -> a) -> a -> [a]
iterate (Matrix Double -> (Matrix Double, [Int])
normHNF' (Matrix Double -> (Matrix Double, [Int]))
-> ((Matrix Double, [Int]) -> Matrix Double)
-> (Matrix Double, [Int])
-> (Matrix Double, [Int])
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Matrix Double, [Int]) -> Matrix Double
forall a b. (a, b) -> a
fst) (Matrix Double
m, [])

normHNF' :: H.Matrix Double -> (H.Matrix Double, [Int])
normHNF' :: Matrix Double -> (Matrix Double, [Int])
normHNF' Matrix Double
m = ([CLong] -> [Int])
-> (Matrix Double, [CLong]) -> (Matrix Double, [Int])
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ((CLong -> Int) -> [CLong] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map CLong -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral) ((Matrix Double, [CLong]) -> (Matrix Double, [Int]))
-> (IO (Matrix Double, [CLong]) -> (Matrix Double, [CLong]))
-> IO (Matrix Double, [CLong])
-> (Matrix Double, [Int])
forall b c a. (b -> c) -> (a -> b) -> a -> c
. IO (Matrix Double, [CLong]) -> (Matrix Double, [CLong])
forall a. IO a -> a
unsafePerformIO (IO (Matrix Double, [CLong]) -> (Matrix Double, [Int]))
-> IO (Matrix Double, [CLong]) -> (Matrix Double, [Int])
forall a b. (a -> b) -> a -> b
$ Matrix Double
-> ((CLong, CLong, Ptr FMPZMat) -> IO (Matrix Double, [CLong]))
-> IO (Matrix Double, [CLong])
forall b.
Matrix Double -> ((CLong, CLong, Ptr FMPZMat) -> IO b) -> IO b
withMatrix Matrix Double
m (CLong, CLong, Ptr FMPZMat) -> IO (Matrix Double, [CLong])
normhnf

normhnf :: (CLong, CLong, Ptr FMPZMat) -> IO (H.Matrix Double, [CLong])
normhnf :: (CLong, CLong, Ptr FMPZMat) -> IO (Matrix Double, [CLong])
normhnf (CLong
numRows, CLong
numCols, Ptr FMPZMat
inputM) = do
  CLong
-> CLong
-> (Ptr FMPZMat -> IO (Matrix Double, [CLong]))
-> IO (Matrix Double, [CLong])
forall b. CLong -> CLong -> (Ptr FMPZMat -> IO b) -> IO b
withBlankMatrix CLong
numRows CLong
numCols ((Ptr FMPZMat -> IO (Matrix Double, [CLong]))
 -> IO (Matrix Double, [CLong]))
-> (Ptr FMPZMat -> IO (Matrix Double, [CLong]))
-> IO (Matrix Double, [CLong])
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
outputM -> do
    Ptr FMPZMat -> Ptr FMPZMat -> IO ()
fmpz_mat_hnf Ptr FMPZMat
outputM Ptr FMPZMat
inputM
    CLong
rank <- Ptr FMPZMat -> IO CLong
fmpz_mat_rank Ptr FMPZMat
outputM
    -- scale the rows so that it is a RREF, returning the indices
    -- where leading-coefficients had to be scaled to 1.
    [(CLong, CLong)]
indices <- Ptr FMPZMat -> CLong -> CLong -> IO [(CLong, CLong)]
elemrowscale Ptr FMPZMat
outputM CLong
rank CLong
numCols
    -- HNF allows non-zeroes above the leading-coefficients that
    -- were greater than 1, so also fix that to bring into RREF.
    Ptr FMPZMat -> CLong -> CLong -> [(CLong, CLong)] -> IO ()
elemrowadds Ptr FMPZMat
outputM CLong
rank CLong
numCols [(CLong, CLong)]
indices

    -- column indices of leading co-efficients > 1
    [((CLong, CLong), [CLong])]
lcoefs <- (((CLong, CLong), [CLong]) -> Bool)
-> [((CLong, CLong), [CLong])] -> [((CLong, CLong), [CLong])]
forall a. (a -> Bool) -> [a] -> [a]
filter ((CLong -> CLong -> Bool
forall a. Ord a => a -> a -> Bool
> CLong
1) (CLong -> Bool)
-> (((CLong, CLong), [CLong]) -> CLong)
-> ((CLong, CLong), [CLong])
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [CLong] -> CLong
forall a. [a] -> a
head ([CLong] -> CLong)
-> (((CLong, CLong), [CLong]) -> [CLong])
-> ((CLong, CLong), [CLong])
-> CLong
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((CLong, CLong), [CLong]) -> [CLong]
forall a b. (a, b) -> b
snd) ([((CLong, CLong), [CLong])] -> [((CLong, CLong), [CLong])])
-> IO [((CLong, CLong), [CLong])] -> IO [((CLong, CLong), [CLong])]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [CLong]
-> (CLong -> IO ((CLong, CLong), [CLong]))
-> IO [((CLong, CLong), [CLong])]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [CLong
0 .. CLong
rankCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] (\ CLong
i -> do
      [(CLong, CLong)]
cs <- [CLong] -> [CLong] -> [(CLong, CLong)]
forall a b. [a] -> [b] -> [(a, b)]
zip [CLong
0..] ([CLong] -> [(CLong, CLong)]) -> IO [CLong] -> IO [(CLong, CLong)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
`fmap` [IO CLong] -> IO [CLong]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence [ Ptr FMPZMat -> CLong -> CLong -> IO CLong
peekM Ptr FMPZMat
outputM CLong
i CLong
j | CLong
j <- [CLong
0 .. CLong
numColsCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ]
      let ([(CLong, CLong)]
_, (CLong
j, CLong
lcoef):[(CLong, CLong)]
rest) = ((CLong, CLong) -> Bool)
-> [(CLong, CLong)] -> ([(CLong, CLong)], [(CLong, CLong)])
forall a. (a -> Bool) -> [a] -> ([a], [a])
span ((CLong -> CLong -> Bool
forall a. Eq a => a -> a -> Bool
== CLong
0) (CLong -> Bool)
-> ((CLong, CLong) -> CLong) -> (CLong, CLong) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CLong, CLong) -> CLong
forall a b. (a, b) -> b
snd) [(CLong, CLong)]
cs
      ((CLong, CLong), [CLong]) -> IO ((CLong, CLong), [CLong])
forall (m :: * -> *) a. Monad m => a -> m a
return ((CLong
i, CLong
j), CLong
lcoefCLong -> [CLong] -> [CLong]
forall a. a -> [a] -> [a]
:((CLong, CLong) -> CLong) -> [(CLong, CLong)] -> [CLong]
forall a b. (a -> b) -> [a] -> [b]
map (CLong, CLong) -> CLong
forall a b. (a, b) -> b
snd [(CLong, CLong)]
rest))

    -- split the identified rows into two categories: those that can
    -- be divided out by the leading co-efficient, and those that
    -- cannot.
    let ([((CLong, CLong), [CLong])]
multCands, [((CLong, CLong), [CLong])]
consCands) = (((CLong, CLong), [CLong]) -> Bool)
-> [((CLong, CLong), [CLong])]
-> ([((CLong, CLong), [CLong])], [((CLong, CLong), [CLong])])
forall a. (a -> Bool) -> [a] -> ([a], [a])
partition (\ ((CLong, CLong)
_, CLong
lcoef:[CLong]
rest) -> (CLong -> Bool) -> [CLong] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ((CLong -> CLong -> Bool
forall a. Eq a => a -> a -> Bool
== CLong
0) (CLong -> Bool) -> (CLong -> CLong) -> CLong -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CLong -> CLong -> CLong
forall a. Integral a => a -> a -> a
`rem` CLong
lcoef)) [CLong]
rest) [((CLong, CLong), [CLong])]
lcoefs

    -- apply elementary row scaling
    [((CLong, CLong), [CLong])]
-> (((CLong, CLong), [CLong]) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [((CLong, CLong), [CLong])]
multCands ((((CLong, CLong), [CLong]) -> IO ()) -> IO ())
-> (((CLong, CLong), [CLong]) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ ((CLong
i, CLong
j), CLong
lcoef:[CLong]
_) ->
      [CLong] -> (CLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [CLong
j..CLong
numCols CLong -> CLong -> CLong
forall a. Num a => a -> a -> a
- CLong
1] ((CLong -> IO ()) -> IO ()) -> (CLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ CLong
j' -> do
        CLong
x <- Ptr FMPZMat -> CLong -> CLong -> IO CLong
peekM Ptr FMPZMat
outputM CLong
i CLong
j'
        Ptr FMPZMat -> CLong -> CLong -> CLong -> IO ()
pokeM Ptr FMPZMat
outputM CLong
i CLong
j' (CLong -> IO ()) -> CLong -> IO ()
forall a b. (a -> b) -> a -> b
$ CLong
x CLong -> CLong -> CLong
forall a. Integral a => a -> a -> a
`div` CLong
lcoef

    -- identify columns that need additional constraints generated and
    -- their associated value d, which is the value of the non-leading
    -- co-efficient divided by the GCD of the row.
    let genColCons :: ((a, a), [b]) -> (a, b)
genColCons ((a
_, a
j), b
lcoef:[b]
rest) = (a
j, b
minNLcoef b -> b -> b
forall a. Integral a => a -> a -> a
`div` b -> b -> b
forall a. Integral a => a -> a -> a
gcd b
lcoef b
minNLcoef)
          where
            restABS :: [b]
restABS   = (b -> b) -> [b] -> [b]
forall a b. (a -> b) -> [a] -> [b]
map b -> b
forall a. Num a => a -> a
abs [b]
rest
            minNLcoef :: b
minNLcoef = [b] -> b
forall (t :: * -> *) a. (Foldable t, Ord a) => t a -> a
minimum ((b -> Bool) -> [b] -> [b]
forall a. (a -> Bool) -> [a] -> [a]
filter (b -> b -> Bool
forall a. Eq a => a -> a -> Bool
/= b
0) [b]
restABS)
        genColCons ((a, a), [b])
_ = [Char] -> (a, b)
forall a. HasCallStack => [Char] -> a
error [Char]
"normhnf: genColCons: impossible: missing leading co-efficient"

    let consCols :: [(CLong, CLong)]
consCols = (((CLong, CLong), [CLong]) -> (CLong, CLong))
-> [((CLong, CLong), [CLong])] -> [(CLong, CLong)]
forall a b. (a -> b) -> [a] -> [b]
map ((CLong, CLong), [CLong]) -> (CLong, CLong)
forall b a a. Integral b => ((a, a), [b]) -> (a, b)
genColCons [((CLong, CLong), [CLong])]
consCands

    -- generate operations that poke 1.0 and -d into columns of the
    -- given row (matrix parameter supplied later); d is the value of
    -- the non-leading co-efficient that isn't divisible by the
    -- leading-coefficient, divided by the GCD of the row.
    let ops :: [Ptr FMPZMat -> CLong -> IO ()]
ops = [ \ Ptr FMPZMat
flintM CLong
i -> do Ptr FMPZMat -> CLong -> CLong -> CLong -> IO ()
pokeM Ptr FMPZMat
flintM CLong
i CLong
j CLong
1
                                 Ptr FMPZMat -> CLong -> CLong -> CLong -> IO ()
pokeM Ptr FMPZMat
flintM CLong
i (CLong
numCols CLong -> CLong -> CLong
forall a. Num a => a -> a -> a
+ CLong
k) (-CLong
d)
              | ((CLong
j, CLong
d), CLong
k) <- [(CLong, CLong)] -> [CLong] -> [((CLong, CLong), CLong)]
forall a b. [a] -> [b] -> [(a, b)]
zip [(CLong, CLong)]
consCols [CLong
0..] ]
    let numOps :: CLong
numOps = Int -> CLong
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Ptr FMPZMat -> CLong -> IO ()] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Ptr FMPZMat -> CLong -> IO ()]
ops)
    let numRows' :: CLong
numRows' = CLong
numOps CLong -> CLong -> CLong
forall a. Num a => a -> a -> a
+ CLong
rank
    let numCols' :: CLong
numCols' = CLong
numOps CLong -> CLong -> CLong
forall a. Num a => a -> a -> a
+ CLong
numCols

    CLong
-> CLong
-> (Ptr FMPZMat -> IO (Matrix Double, [CLong]))
-> IO (Matrix Double, [CLong])
forall b. CLong -> CLong -> (Ptr FMPZMat -> IO b) -> IO b
withBlankMatrix CLong
numRows' CLong
numCols' ((Ptr FMPZMat -> IO (Matrix Double, [CLong]))
 -> IO (Matrix Double, [CLong]))
-> (Ptr FMPZMat -> IO (Matrix Double, [CLong]))
-> IO (Matrix Double, [CLong])
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
outputM' -> do
      -- create larger matrix containing outputM as submatrix
      Ptr FMPZMat
-> Ptr FMPZMat -> CLong -> CLong -> CLong -> CLong -> IO ()
copyMatrix Ptr FMPZMat
outputM' Ptr FMPZMat
outputM CLong
0 CLong
0 CLong
rank CLong
numCols

      -- apply the operations on the extra space to write the
      -- additional rows and columns providing the additional
      -- constraints
      [(CLong, Ptr FMPZMat -> CLong -> IO ())]
-> ((CLong, Ptr FMPZMat -> CLong -> IO ()) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ ([CLong]
-> [Ptr FMPZMat -> CLong -> IO ()]
-> [(CLong, Ptr FMPZMat -> CLong -> IO ())]
forall a b. [a] -> [b] -> [(a, b)]
zip [CLong
rank..] [Ptr FMPZMat -> CLong -> IO ()]
ops) (((CLong, Ptr FMPZMat -> CLong -> IO ()) -> IO ()) -> IO ())
-> ((CLong, Ptr FMPZMat -> CLong -> IO ()) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ (CLong
i, Ptr FMPZMat -> CLong -> IO ()
op) -> Ptr FMPZMat -> CLong -> IO ()
op Ptr FMPZMat
outputM' CLong
i

      -- re-run HNF
      CLong
-> CLong
-> (Ptr FMPZMat -> IO (Matrix Double, [CLong]))
-> IO (Matrix Double, [CLong])
forall b. CLong -> CLong -> (Ptr FMPZMat -> IO b) -> IO b
withBlankMatrix CLong
numRows' CLong
numCols' ((Ptr FMPZMat -> IO (Matrix Double, [CLong]))
 -> IO (Matrix Double, [CLong]))
-> (Ptr FMPZMat -> IO (Matrix Double, [CLong]))
-> IO (Matrix Double, [CLong])
forall a b. (a -> b) -> a -> b
$ \ Ptr FMPZMat
outputM'' -> do
        Ptr FMPZMat -> Ptr FMPZMat -> IO ()
fmpz_mat_hnf Ptr FMPZMat
outputM'' Ptr FMPZMat
outputM'
        CLong
rank' <- Ptr FMPZMat -> IO CLong
fmpz_mat_rank Ptr FMPZMat
outputM''
        -- scale the rows so that it is a RREF, returning the indices
        -- where leading-coefficients had to be scaled to 1.
        [(CLong, CLong)]
indices' <- Ptr FMPZMat -> CLong -> CLong -> IO [(CLong, CLong)]
elemrowscale Ptr FMPZMat
outputM'' CLong
rank' CLong
numCols'
        -- HNF allows non-zeroes above the leading-coefficients that
        -- were greater than 1, so also fix that to bring into RREF.
        Ptr FMPZMat -> CLong -> CLong -> [(CLong, CLong)] -> IO ()
elemrowadds Ptr FMPZMat
outputM'' CLong
rank' CLong
numCols' [(CLong, CLong)]
indices'
        -- convert back to HMatrix form
        Matrix Double
h <- CLong -> CLong -> Ptr FMPZMat -> IO (Matrix Double)
flintToHMatrix CLong
rank' CLong
numCols' Ptr FMPZMat
outputM''
        (Matrix Double, [CLong]) -> IO (Matrix Double, [CLong])
forall (m :: * -> *) a. Monad m => a -> m a
return (Matrix Double
h, ((CLong, CLong) -> CLong) -> [(CLong, CLong)] -> [CLong]
forall a b. (a -> b) -> [a] -> [b]
map (CLong, CLong) -> CLong
forall a b. (a, b) -> a
fst [(CLong, CLong)]
consCols)

-- find leading-coefficients that are greater than 1 and scale those
-- rows accordingly to reach RREF
--
-- precondition: matrix outputM is in HNF
elemrowscale :: Ptr FMPZMat -> CLong -> CLong -> IO [(CLong, CLong)]
elemrowscale :: Ptr FMPZMat -> CLong -> CLong -> IO [(CLong, CLong)]
elemrowscale Ptr FMPZMat
outputM CLong
rank CLong
numCols = do
  -- column indices of leading co-efficients > 1
  [((CLong, CLong), [CLong])]
lcoefs <- (((CLong, CLong), [CLong]) -> Bool)
-> [((CLong, CLong), [CLong])] -> [((CLong, CLong), [CLong])]
forall a. (a -> Bool) -> [a] -> [a]
filter ((CLong -> CLong -> Bool
forall a. Ord a => a -> a -> Bool
> CLong
1) (CLong -> Bool)
-> (((CLong, CLong), [CLong]) -> CLong)
-> ((CLong, CLong), [CLong])
-> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [CLong] -> CLong
forall a. [a] -> a
head ([CLong] -> CLong)
-> (((CLong, CLong), [CLong]) -> [CLong])
-> ((CLong, CLong), [CLong])
-> CLong
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((CLong, CLong), [CLong]) -> [CLong]
forall a b. (a, b) -> b
snd) ([((CLong, CLong), [CLong])] -> [((CLong, CLong), [CLong])])
-> IO [((CLong, CLong), [CLong])] -> IO [((CLong, CLong), [CLong])]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> [CLong]
-> (CLong -> IO ((CLong, CLong), [CLong]))
-> IO [((CLong, CLong), [CLong])]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [CLong
0 .. CLong
rankCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] (\ CLong
i -> do
    [(CLong, CLong)]
cs <- [CLong] -> [CLong] -> [(CLong, CLong)]
forall a b. [a] -> [b] -> [(a, b)]
zip [CLong
0..] ([CLong] -> [(CLong, CLong)]) -> IO [CLong] -> IO [(CLong, CLong)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
`fmap` [IO CLong] -> IO [CLong]
forall (t :: * -> *) (m :: * -> *) a.
(Traversable t, Monad m) =>
t (m a) -> m (t a)
sequence [ Ptr FMPZMat -> CLong -> CLong -> IO CLong
peekM Ptr FMPZMat
outputM CLong
i CLong
j | CLong
j <- [CLong
0 .. CLong
numColsCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ]
    let ([(CLong, CLong)]
_, (CLong
j, CLong
lcoef):[(CLong, CLong)]
rest) = ((CLong, CLong) -> Bool)
-> [(CLong, CLong)] -> ([(CLong, CLong)], [(CLong, CLong)])
forall a. (a -> Bool) -> [a] -> ([a], [a])
span ((CLong -> CLong -> Bool
forall a. Eq a => a -> a -> Bool
== CLong
0) (CLong -> Bool)
-> ((CLong, CLong) -> CLong) -> (CLong, CLong) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CLong, CLong) -> CLong
forall a b. (a, b) -> b
snd) [(CLong, CLong)]
cs
    ((CLong, CLong), [CLong]) -> IO ((CLong, CLong), [CLong])
forall (m :: * -> *) a. Monad m => a -> m a
return ((CLong
i, CLong
j), CLong
lcoefCLong -> [CLong] -> [CLong]
forall a. a -> [a] -> [a]
:((CLong, CLong) -> CLong) -> [(CLong, CLong)] -> [CLong]
forall a b. (a -> b) -> [a] -> [b]
map (CLong, CLong) -> CLong
forall a b. (a, b) -> b
snd [(CLong, CLong)]
rest))

  let multCands :: [((CLong, CLong), [CLong])]
multCands = (((CLong, CLong), [CLong]) -> Bool)
-> [((CLong, CLong), [CLong])] -> [((CLong, CLong), [CLong])]
forall a. (a -> Bool) -> [a] -> [a]
filter (\ ((CLong, CLong)
_, CLong
lcoef:[CLong]
rest) -> (CLong -> Bool) -> [CLong] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all ((CLong -> CLong -> Bool
forall a. Eq a => a -> a -> Bool
== CLong
0) (CLong -> Bool) -> (CLong -> CLong) -> CLong -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CLong -> CLong -> CLong
forall a. Integral a => a -> a -> a
`rem` CLong
lcoef)) [CLong]
rest) [((CLong, CLong), [CLong])]
lcoefs

  -- apply elementary row scaling
  [((CLong, CLong), [CLong])]
-> (((CLong, CLong), [CLong]) -> IO (CLong, CLong))
-> IO [(CLong, CLong)]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [((CLong, CLong), [CLong])]
multCands ((((CLong, CLong), [CLong]) -> IO (CLong, CLong))
 -> IO [(CLong, CLong)])
-> (((CLong, CLong), [CLong]) -> IO (CLong, CLong))
-> IO [(CLong, CLong)]
forall a b. (a -> b) -> a -> b
$ \ ((CLong
i, CLong
j), CLong
lcoef:[CLong]
_) -> do
    [CLong] -> (CLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [CLong
j..CLong
numCols CLong -> CLong -> CLong
forall a. Num a => a -> a -> a
- CLong
1] ((CLong -> IO ()) -> IO ()) -> (CLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ CLong
j' -> do
      CLong
x <- Ptr FMPZMat -> CLong -> CLong -> IO CLong
peekM Ptr FMPZMat
outputM CLong
i CLong
j'
      Ptr FMPZMat -> CLong -> CLong -> CLong -> IO ()
pokeM Ptr FMPZMat
outputM CLong
i CLong
j' (CLong -> IO ()) -> CLong -> IO ()
forall a b. (a -> b) -> a -> b
$ CLong
x CLong -> CLong -> CLong
forall a. Integral a => a -> a -> a
`div` CLong
lcoef
    (CLong, CLong) -> IO (CLong, CLong)
forall (m :: * -> *) a. Monad m => a -> m a
return (CLong
i, CLong
j)

-- use indices to guide elementary row additions to reach RREF
--
-- precondition: matrix outputM is in HNF save for work done by
-- elemrowscale, and indices is a list of coordinates where we have
-- just scaled the leading-coefficient to 1 and now we must look to
-- see if there are any non-zeroes in the column above the
-- leading-coefficient, because that is allowed by HNF.
elemrowadds :: Ptr FMPZMat -> CLong -> CLong -> [(CLong, CLong)] -> IO ()
elemrowadds :: Ptr FMPZMat -> CLong -> CLong -> [(CLong, CLong)] -> IO ()
elemrowadds Ptr FMPZMat
outputM CLong
_ CLong
numCols [(CLong, CLong)]
indices = do
  -- look for non-zero members of the columns above the
  -- leading-coefficient and wipe them out.
  [(CLong, CLong)] -> ((CLong, CLong) -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [(CLong, CLong)]
indices (((CLong, CLong) -> IO ()) -> IO ())
-> ((CLong, CLong) -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ (CLong
lcI, CLong
lcJ) -> do
    let j :: CLong
j = CLong
lcJ
    [CLong] -> (CLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [CLong
0..CLong
lcICLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO ()) -> IO ()) -> (CLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ CLong
i -> do
      -- (i, j) ranges over the elements of the column above leading
      -- co-efficient (lcI, lcJ).
      CLong
x <- Ptr FMPZMat -> CLong -> CLong -> IO CLong
peekM Ptr FMPZMat
outputM CLong
i CLong
j
      if CLong
x CLong -> CLong -> Bool
forall a. Eq a => a -> a -> Bool
== CLong
0 then
        () -> IO ()
forall (f :: * -> *) a. Applicative f => a -> f a
pure () -- nothing to do at row i
      else do
        -- (i, j) is non-zero and row i must be cancelled x times
        let sf :: CLong
sf = CLong
x -- scaling factor = non-zero magnitude
        [CLong] -> (CLong -> IO ()) -> IO ()
forall (t :: * -> *) (m :: * -> *) a b.
(Foldable t, Monad m) =>
t a -> (a -> m b) -> m ()
forM_ [CLong
lcJ..CLong
numColsCLong -> CLong -> CLong
forall a. Num a => a -> a -> a
-CLong
1] ((CLong -> IO ()) -> IO ()) -> (CLong -> IO ()) -> IO ()
forall a b. (a -> b) -> a -> b
$ \ CLong
j' -> do
          -- (i,   j') ranges over the row where we discovered a non-zero
          -- (lcI, j') ranges over the row with the leading co-efficient
          CLong
x1 <- Ptr FMPZMat -> CLong -> CLong -> IO CLong
peekM Ptr FMPZMat
outputM CLong
i CLong
j'
          CLong
x2 <- Ptr FMPZMat -> CLong -> CLong -> IO CLong
peekM Ptr FMPZMat
outputM CLong
lcI CLong
j'
          -- add lower row scaled by sf to upper row
          Ptr FMPZMat -> CLong -> CLong -> CLong -> IO ()
pokeM Ptr FMPZMat
outputM CLong
i CLong
j' (CLong
x1 CLong -> CLong -> CLong
forall a. Num a => a -> a -> a
- CLong
x2 CLong -> CLong -> CLong
forall a. Num a => a -> a -> a
* CLong
sf)

--------------------------------------------------

-- m1 :: H.Matrix Double
-- m1 = (8><6)
--  [ 1.0, 0.0, 0.0, -1.0,  0.0,  0.0
--  , 0.0, 1.0, 0.0,  0.0, -1.0,  0.0
--  , 0.0, 0.0, 1.0,  0.0,  0.0, -1.0
--  , 1.0, 0.0, 0.0, -1.0,  0.0,  0.0
--  , 0.0, 1.0, 0.0,  0.0, -1.0,  0.0
--  , 0.0, 0.0, 1.0,  0.0,  0.0, -1.0
--  , 0.0, 0.0, 0.0,  1.0, -4.0,  0.0
--  , 0.0, 0.0, 0.0,  0.0,  4.0, -3.0 ]

-- m2 :: H.Matrix Double
-- m2 = (8><6)
--  [ 1.0, 0.0, 0.0, -1.0,  0.0,  0.0
--  , 0.0, 1.0, 0.0,  0.0, -1.0,  0.0
--  , 0.0, 0.0, 1.0,  0.0,  0.0, -1.0
--  , 1.0, 0.0, 0.0, -1.0,  0.0,  0.0
--  , 0.0, 1.0, 0.0,  0.0, -1.0,  0.0
--  , 0.0, 0.0, 1.0,  0.0,  0.0, -1.0
--  , 0.0, 0.0, 0.0,  1.0, -6.0,  0.0
--  , 0.0, 0.0, 0.0,  0.0,  6.0, -4.0 ]