eigen 1.2.5 → 2.0.0

raw patch · 11 files changed

+575/−404 lines, 11 filesPVP ok

version bump matches the API change (PVP)

API changes (from Hackage documentation)

Files

• Data/Eigen/Internal.hs +114/−34
• Data/Eigen/LA.hs +13/−13
• Data/Eigen/Matrix.hs +182/−146
• Data/Eigen/Matrix/Mutable.hs +33/−17
• Data/Eigen/Parallel.hs +5/−5
• cbits/eigen-proxy.cpp +188/−149
• cbits/eigen-proxy.h +31/−34
• eigen.cabal +2/−2
• test/rank.hs +2/−2
• test/regression.hs +4/−2
• test/solve.hs +1/−0

@@ -1,26 +1,46 @@-{-# LANGUAGE MultiParamTypeClasses, ForeignFunctionInterface, EmptyDataDecls #-}+{-# LANGUAGE MultiParamTypeClasses, ForeignFunctionInterface, ScopedTypeVariables, FunctionalDependencies, FlexibleInstances #-} module Data.Eigen.Internal where  import Foreign.Ptr+import Foreign.Storable import Foreign.C.Types import Foreign.C.String import Control.Monad+import Control.Applicative import System.IO.Unsafe+import Data.Complex +class (Num a, Ord a, Cast a b, Cast b a, Storable b, Code b) => Elem a b | a -> b where++instance Elem Float CFloat where+instance Elem Double CDouble where+ class Cast a b where     cast :: a -> b -instance Cast CDouble Double where-    cast (CDouble x) = x+data CComplex a = CComplex !a !a -instance Cast Double CDouble where-    cast = CDouble+instance Storable a => Storable (CComplex a) where+    sizeOf _ = sizeOf (undefined :: a) * 2+    alignment _ = alignment (undefined :: a)+    poke p (CComplex x y) = do+        pokeElemOff (castPtr p) 0 x+        pokeElemOff (castPtr p) 1 y+    peek p = CComplex+        <$> peekElemOff (castPtr p) 0+        <*> peekElemOff (castPtr p) 1 -instance Cast CInt Int where-    cast = fromIntegral+instance Cast CInt Int where; cast = fromIntegral+instance Cast Int CInt where; cast = fromIntegral+instance Cast CFloat Float where; cast (CFloat x) = x+instance Cast Float CFloat where; cast = CFloat+instance Cast CDouble Double where; cast (CDouble x) = x+instance Cast Double CDouble where; cast = CDouble+instance Cast (CComplex CFloat) (Complex Float) where; cast (CComplex x y) = cast x :+ cast y+instance Cast (Complex Float) (CComplex CFloat) where; cast (x :+ y) = CComplex (cast x) (cast y)+instance Cast (CComplex CDouble) (Complex Double) where; cast (CComplex x y) = cast x :+ cast y+instance Cast (Complex Double) (CComplex CDouble) where; cast (x :+ y) = CComplex (cast x) (cast y) -instance Cast Int CInt where-    cast = fromIntegral  performIO :: IO a -> a performIO = unsafeDupablePerformIO@@ -31,32 +51,92 @@ call func = func >>= \c_str -> when (c_str /= nullPtr) $     peekCString c_str >>= \str -> c_freeString c_str >> fail str +foreign import ccall "eigen-proxy.h free" free :: Ptr a -> IO ()+ foreign import ccall "eigen-proxy.h eigen_setNbThreads" c_setNbThreads :: CInt -> IO () foreign import ccall "eigen-proxy.h eigen_getNbThreads" c_getNbThreads :: IO CInt -foreign import ccall "eigen-proxy.h eigen_random"      c_random      :: Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_add"         c_add         :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_sub"         c_sub         :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_mul"         c_mul         :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_diagonal"    c_diagonal    :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_transpose"   c_transpose   :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_inverse"     c_inverse     :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_adjoint"     c_adjoint     :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_conjugate"   c_conjugate   :: Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_normalize"   c_normalize   :: Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_sum"         c_sum         :: Ptr CDouble -> CInt -> CInt -> IO CDouble-foreign import ccall "eigen-proxy.h eigen_prod"        c_prod        :: Ptr CDouble -> CInt -> CInt -> IO CDouble-foreign import ccall "eigen-proxy.h eigen_mean"        c_mean        :: Ptr CDouble -> CInt -> CInt -> IO CDouble-foreign import ccall "eigen-proxy.h eigen_norm"        c_norm        :: Ptr CDouble -> CInt -> CInt -> IO CDouble-foreign import ccall "eigen-proxy.h eigen_trace"       c_trace       :: Ptr CDouble -> CInt -> CInt -> IO CDouble-foreign import ccall "eigen-proxy.h eigen_squaredNorm" c_squaredNorm :: Ptr CDouble -> CInt -> CInt -> IO CDouble-foreign import ccall "eigen-proxy.h eigen_blueNorm"    c_blueNorm    :: Ptr CDouble -> CInt -> CInt -> IO CDouble-foreign import ccall "eigen-proxy.h eigen_hypotNorm"   c_hypotNorm   :: Ptr CDouble -> CInt -> CInt -> IO CDouble-foreign import ccall "eigen-proxy.h eigen_determinant" c_determinant :: Ptr CDouble -> CInt -> CInt -> IO CDouble+foreign import ccall "eigen-proxy.h eigen_random"      c_random      :: CInt -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_add"         c_add         :: CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_sub"         c_sub         :: CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_mul"         c_mul         :: CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_diagonal"    c_diagonal    :: CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_transpose"   c_transpose   :: CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_inverse"     c_inverse     :: CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_adjoint"     c_adjoint     :: CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_conjugate"   c_conjugate   :: CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_normalize"   c_normalize   :: CInt -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_sum"         c_sum         :: CInt -> Ptr a -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_prod"        c_prod        :: CInt -> Ptr a -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_mean"        c_mean        :: CInt -> Ptr a -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_norm"        c_norm        :: CInt -> Ptr a -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_trace"       c_trace       :: CInt -> Ptr a -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_squaredNorm" c_squaredNorm :: CInt -> Ptr a -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_blueNorm"    c_blueNorm    :: CInt -> Ptr a -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_hypotNorm"   c_hypotNorm   :: CInt -> Ptr a -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_determinant" c_determinant :: CInt -> Ptr a -> Ptr a -> CInt -> CInt -> IO CString -foreign import ccall "eigen-proxy.h eigen_rank"         c_rank       :: CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_image"        c_image      :: CInt -> Ptr (Ptr CDouble) -> Ptr CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_kernel"       c_kernel     :: CInt -> Ptr (Ptr CDouble) -> Ptr CInt -> Ptr CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h free"               c_free       :: Ptr a -> IO ()-foreign import ccall "eigen-proxy.h eigen_solve"        c_solve       :: CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_relativeError" c_relativeError :: Ptr CDouble -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_rank"         c_rank       :: CInt -> CInt -> Ptr CInt -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_image"        c_image      :: CInt -> CInt -> Ptr (Ptr a) -> Ptr CInt -> Ptr CInt -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_kernel"       c_kernel     :: CInt -> CInt -> Ptr (Ptr a) -> Ptr CInt -> Ptr CInt -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_solve"        c_solve       :: CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+foreign import ccall "eigen-proxy.h eigen_relativeError" c_relativeError :: CInt -> Ptr a -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+++class Code a where; code :: a -> CInt+instance Code CFloat where; code _ = 0+instance Code CDouble where; code _ = 1+instance Code (CComplex CFloat) where; code _ = 2+instance Code (CComplex CDouble) where; code _ = 3++random      :: forall a . Code a => Ptr a -> CInt -> CInt -> IO CString+add         :: forall a . Code a => Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+sub         :: forall a . Code a => Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+mul         :: forall a . Code a => Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+diagonal    :: forall a . Code a => Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+transpose   :: forall a . Code a => Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+inverse     :: forall a . Code a => Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+adjoint     :: forall a . Code a => Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+conjugate   :: forall a . Code a => Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+normalize   :: forall a . Code a => Ptr a -> CInt -> CInt -> IO CString+sum         :: forall a . Code a => Ptr a -> Ptr a -> CInt -> CInt -> IO CString+prod        :: forall a . Code a => Ptr a -> Ptr a -> CInt -> CInt -> IO CString+mean        :: forall a . Code a => Ptr a -> Ptr a -> CInt -> CInt -> IO CString+norm        :: forall a . Code a => Ptr a -> Ptr a -> CInt -> CInt -> IO CString+trace       :: forall a . Code a => Ptr a -> Ptr a -> CInt -> CInt -> IO CString+squaredNorm :: forall a . Code a => Ptr a -> Ptr a -> CInt -> CInt -> IO CString+blueNorm    :: forall a . Code a => Ptr a -> Ptr a -> CInt -> CInt -> IO CString+hypotNorm   :: forall a . Code a => Ptr a -> Ptr a -> CInt -> CInt -> IO CString+determinant :: forall a . Code a => Ptr a -> Ptr a -> CInt -> CInt -> IO CString++rank       :: forall a . Code a => CInt -> Ptr CInt -> Ptr a -> CInt -> CInt -> IO CString+image      :: forall a . Code a => CInt -> Ptr (Ptr a) -> Ptr CInt -> Ptr CInt -> Ptr a -> CInt -> CInt -> IO CString+kernel     :: forall a . Code a => CInt -> Ptr (Ptr a) -> Ptr CInt -> Ptr CInt -> Ptr a -> CInt -> CInt -> IO CString+solve      :: forall a . Code a => CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString+relativeError :: forall a . Code a => Ptr a -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString++random      = c_random      (code (undefined :: a))+add         = c_add         (code (undefined :: a))+sub         = c_sub         (code (undefined :: a))+mul         = c_mul         (code (undefined :: a))+diagonal    = c_diagonal    (code (undefined :: a))+transpose   = c_transpose   (code (undefined :: a))+inverse     = c_inverse     (code (undefined :: a))+adjoint     = c_adjoint     (code (undefined :: a))+conjugate   = c_conjugate   (code (undefined :: a))+normalize   = c_normalize   (code (undefined :: a))+sum         = c_sum         (code (undefined :: a))+prod        = c_prod        (code (undefined :: a))+mean        = c_mean        (code (undefined :: a))+norm        = c_norm        (code (undefined :: a))+trace       = c_trace       (code (undefined :: a))+squaredNorm = c_squaredNorm (code (undefined :: a))+blueNorm    = c_blueNorm    (code (undefined :: a))+hypotNorm   = c_hypotNorm   (code (undefined :: a))+determinant = c_determinant (code (undefined :: a))++rank        = c_rank (code (undefined :: a))+image       = c_image (code (undefined :: a))+kernel      = c_kernel (code (undefined :: a))+solve       = c_solve (code (undefined :: a))+relativeError = c_relativeError (code (undefined :: a))

@@ -151,68 +151,68 @@   -- | [x = solve d a b] finds a solution @x@ of @ax = b@ equation using decomposition @d@-solve :: Decomposition -> Matrix -> Matrix -> Matrix+solve :: I.Elem a b => Decomposition -> Matrix a b -> Matrix a b -> Matrix a b solve d a b = I.performIO $ do-    x <- M.new (m_cols a) 1+    x <- M.new (cols a) 1     M.unsafeWith x $ \x_vals x_rows x_cols ->         unsafeWith a $ \a_vals a_rows a_cols ->             unsafeWith b $ \b_vals b_rows b_cols ->-                I.call $ I.c_solve (I.cast $ fromEnum d)+                I.call $ I.solve (I.cast $ fromEnum d)                     x_vals x_rows x_cols                     a_vals a_rows a_cols                     b_vals b_rows b_cols     unsafeFreeze x  -- | [e = relativeError x a b] computes @norm (ax - b) / norm b@ where @norm@ is L2 norm-relativeError :: Matrix -> Matrix -> Matrix -> Double+relativeError :: I.Elem a b => Matrix a b -> Matrix a b -> Matrix a b -> a relativeError x a b = I.performIO $     unsafeWith x $ \x_vals x_rows x_cols ->         unsafeWith a $ \a_vals a_rows a_cols ->             unsafeWith b $ \b_vals b_rows b_cols ->                 alloca $ \pe -> do-                    I.call $ I.c_relativeError pe+                    I.call $ I.relativeError pe                         x_vals x_rows x_cols                         a_vals a_rows a_cols                         b_vals b_rows b_cols                     I.cast <$> peek pe  -- | The rank of the matrix-rank :: Decomposition -> Matrix -> Int+rank :: I.Elem a b => Decomposition -> Matrix a b -> Int rank d m = I.performIO $ alloca $ \pr -> do-    I.call $ unsafeWith m $ I.c_rank (I.cast $ fromEnum d) pr+    I.call $ unsafeWith m $ I.rank (I.cast $ fromEnum d) pr     I.cast <$> peek pr  -- | Return matrix whose columns form a basis of the null-space of @A@-kernel :: Decomposition -> Matrix -> Matrix+kernel :: I.Elem a b => Decomposition -> Matrix a b -> Matrix a b kernel d m1 = I.performIO $     alloca $ \pvals ->     alloca $ \prows ->     alloca $ \pcols ->         unsafeWith m1 $ \vals1 rows1 cols1 -> do-            I.call $ I.c_kernel (I.cast $ fromEnum d)+            I.call $ I.kernel (I.cast $ fromEnum d)                 pvals prows pcols                 vals1 rows1 cols1             vals <- peek pvals             rows <- I.cast <$> peek prows             cols <- I.cast <$> peek pcols-            fp <- FC.newForeignPtr vals $ I.c_free vals+            fp <- FC.newForeignPtr vals $ I.free vals             return $ Matrix rows cols $ VS.unsafeFromForeignPtr0 fp $ rows * cols   -- | Return a matrix whose columns form a basis of the column-space of @A@-image :: Decomposition -> Matrix -> Matrix+image :: I.Elem a b => Decomposition -> Matrix a b -> Matrix a b image d m1 = I.performIO $     alloca $ \pvals ->     alloca $ \prows ->     alloca $ \pcols ->         unsafeWith m1 $ \vals1 rows1 cols1 -> do-            I.call $ I.c_image (I.cast $ fromEnum d)+            I.call $ I.image (I.cast $ fromEnum d)                 pvals prows pcols                 vals1 rows1 cols1             vals <- peek pvals             rows <- I.cast <$> peek prows             cols <- I.cast <$> peek pcols-            fp <- FC.newForeignPtr vals $ I.c_free vals+            fp <- FC.newForeignPtr vals $ I.free vals             return $ Matrix rows cols $ VS.unsafeFromForeignPtr0 fp $ rows * cols  

@@ -1,7 +1,11 @@-{-# LANGUAGE RecordWildCards, MultiParamTypeClasses, Rank2Types #-}+{-# LANGUAGE RecordWildCards, MultiParamTypeClasses, Rank2Types, ScopedTypeVariables, FunctionalDependencies, GADTs #-} module Data.Eigen.Matrix (     -- * Matrix type     Matrix(..),+    MatrixXf,+    MatrixXd,+    MatrixXcf,+    MatrixXcd,     valid,     -- * Matrix conversions     fromList,@@ -19,6 +23,7 @@     -- * Accessing matrix data     cols,     rows,+    dims,     (!),     coeff,     unsafeCoeff,@@ -46,6 +51,8 @@     fold',     ifold,     ifold',+    fold1,+    fold1',     -- * Boolean reductions     all,     any,@@ -67,23 +74,27 @@     conjugate,     normalize,     modify,-    upperTriangule,-    lowerTriangule,+    upperTriangle,+    lowerTriangle,     -- * Mutable matrices     thaw,     freeze,     unsafeThaw,     unsafeFreeze,+    -- * Raw pointers     unsafeWith, ) where  import qualified Prelude as P import Prelude hiding (null, sum, all, any, map, filter)-import Data.List (intercalate)+import Data.List (intercalate, foldl1') import Data.Tuple+import Data.Complex hiding (conjugate) import Foreign.Ptr import Foreign.C.Types import Foreign.C.String+import Foreign.Storable+import Foreign.Marshal.Alloc import Text.Printf import Control.Monad import Control.Monad.ST@@ -95,20 +106,28 @@ import qualified Data.Eigen.Matrix.Mutable as M  -- | Matrix to be used in pure computations, uses column major memory layout, features copy-free FFI with C++ <http://eigen.tuxfamily.org Eigen> library.-data Matrix = Matrix {-    m_rows :: Int,-    m_cols :: Int,-    m_vals :: VS.Vector CDouble-}; +data Matrix a b where +    Matrix :: I.Elem a b => !Int -> !Int -> !(VS.Vector b) -> Matrix a b++-- | Alias for single precision matrix+type MatrixXf = Matrix Float CFloat+-- | Alias for double precision matrix+type MatrixXd = Matrix Double CDouble+-- | Alias for single previsiom matrix of complex numbers+type MatrixXcf = Matrix (Complex Float) (I.CComplex CFloat)+-- | Alias for double prevision matrix of complex numbers+type MatrixXcd = Matrix (Complex Double) (I.CComplex CDouble)+ -- | Pretty prints the matrix-instance Show Matrix where-    show m@Matrix{..} = concat [-        "Matrix ", show m_rows, "x", show m_cols, +instance (I.Elem a b, Show a) => Show (Matrix a b) where+    show m@(Matrix rows cols _) = concat [+        "Matrix ", show rows, "x", show cols,         "\n", intercalate "\n" $ P.map (intercalate "\t" . P.map show) $ toList m, "\n"] + -- | Shortcuts for basic matrix math-instance Num Matrix where+instance I.Elem a b => Num (Matrix a b) where     (*) = mul     (+) = add     (-) = sub@@ -117,132 +136,138 @@     abs = map abs  -- | Empty 0x0 matrix-empty :: Matrix+empty :: I.Elem a b => Matrix a b empty = Matrix 0 0 VS.empty  -- | Is matrix empty?-null :: Matrix -> Bool-null Matrix{..} = m_rows == 0 && m_cols == 0+null :: I.Elem a b => Matrix a b -> Bool+null (Matrix rows cols _) = rows == 0 && cols == 0  -- | Is matrix square?-square :: Matrix -> Bool-square Matrix{..} = m_rows == m_cols+square :: I.Elem a b => Matrix a b -> Bool+square (Matrix rows cols _) = rows == cols  -- | Matrix where all coeffs are filled with given value-constant :: Int -> Int -> Double -> Matrix+constant :: I.Elem a b => Int -> Int -> a -> Matrix a b constant rows cols val = Matrix rows cols $ VS.replicate (rows * cols) (I.cast val)  -- | Matrix where all coeff are 0-zero :: Int -> Int -> Matrix+zero :: I.Elem a b => Int -> Int -> Matrix a b zero rows cols = constant rows cols 0  -- | Matrix where all coeff are 1-ones :: Int -> Int -> Matrix+ones :: I.Elem a b => Int -> Int -> Matrix a b ones rows cols = constant rows cols 1  -- | Square matrix with 1 on main diagonal and 0 elsewhere-identity :: Int -> Matrix+identity :: forall a b . I.Elem a b => Int -> Matrix a b identity size = Matrix size size $ VS.create $ do-    vm <- VSM.replicate (size * size) 0+    vm <- VSM.replicate (size * size) (I.cast (0::a))     forM_ [0..pred size] $ \n ->-        VSM.write vm (n * size + n) 1+        VSM.write vm (n * size + n) (I.cast (1::a))     return vm  -- | The random matrix of a given size-random :: Int -> Int -> IO Matrix+random :: I.Elem a b => Int -> Int -> IO (Matrix a b) random rows cols = do     m <- M.new rows cols-    I.call $ M.unsafeWith m I.c_random+    I.call $ M.unsafeWith m I.random     unsafeFreeze m  -- | Number of rows for the matrix-rows :: Matrix -> Int-rows = m_rows+rows :: I.Elem a b => Matrix a b -> Int+rows (Matrix rows _ _) = rows  -- | Number of columns for the matrix-cols :: Matrix -> Int-cols = m_cols+cols :: I.Elem a b => Matrix a b -> Int+cols (Matrix _ cols _) = cols +vals :: I.Elem a b => Matrix a b -> VS.Vector b+vals (Matrix _ _ vals) = vals+ -- | Mtrix size as (rows, cols) pair-dims :: Matrix -> (Int, Int)-dims Matrix{..} = (m_rows, m_cols)+dims :: I.Elem a b => Matrix a b -> (Int, Int)+dims (Matrix rows cols _) = (rows, cols)  -- | Matrix coefficient at specific row and col-(!) :: Matrix -> (Int, Int) -> Double+(!) :: I.Elem a b => Matrix a b -> (Int, Int) -> a (!) m (row,col) = coeff row col m  -- | Matrix coefficient at specific row and col-coeff :: Int -> Int -> Matrix -> Double-coeff row col m@Matrix{..}+coeff :: I.Elem a b => Int -> Int -> Matrix a b -> a+coeff row col m@(Matrix rows cols _)     | not (valid m) = error "matrix is not valid"-    | row < 0 || row >= m_rows = error $ printf "Matrix.coeff: row %d is out of bounds [0..%d)" row m_rows-    | col < 0 || col >= m_cols = error $ printf "Matrix.coeff: col %d is out of bounds [0..%d)" col m_cols+    | row < 0 || row >= rows = error $ printf "Matrix.coeff: row %d is out of bounds [0..%d)" row rows+    | col < 0 || col >= cols = error $ printf "Matrix.coeff: col %d is out of bounds [0..%d)" col cols     | otherwise = unsafeCoeff row col m  -- | Unsafe version of coeff function. No bounds check performed so SEGFAULT possible-unsafeCoeff :: Int -> Int -> Matrix -> Double-unsafeCoeff row col Matrix{..} = I.cast $ VS.unsafeIndex m_vals $ col * m_rows + row+unsafeCoeff :: I.Elem a b => Int -> Int -> Matrix a b -> a+unsafeCoeff row col (Matrix rows _ vals) = I.cast $ VS.unsafeIndex vals $ col * rows + row  -- | List of coefficients for the given col-col :: Int -> Matrix -> [Double]-col c m@Matrix{..} = [coeff r c m | r <- [0..pred m_rows]]+col :: I.Elem a b => Int -> Matrix a b -> [a]+col c m@(Matrix rows _ _) = [coeff r c m | r <- [0..pred rows]]  -- | List of coefficients for the given row-row :: Int -> Matrix -> [Double]-row r m@Matrix{..} = [coeff r c m | c <- [0..pred m_cols]]+row :: I.Elem a b => Int -> Matrix a b -> [a]+row r m@(Matrix _ cols _) = [coeff r c m | c <- [0..pred cols]]  -- | Extract rectangular block from matrix defined by startRow startCol blockRows blockCols-block :: Int -> Int -> Int -> Int -> Matrix -> Matrix+block :: I.Elem a b => Int -> Int -> Int -> Int -> Matrix a b -> Matrix a b block startRow startCol blockRows blockCols m =     generate blockRows blockCols $ \row col ->         coeff (startRow + row) (startCol + col) m --- | Verify matrix dimensions and memory layout-valid :: Matrix -> Bool-valid Matrix{..} = m_rows >= 0 && m_cols >= 0 && VS.length m_vals == m_rows * m_cols+-- | Verify matrix dimensions and values memory layout+valid :: I.Elem a b => Matrix a b -> Bool+valid (Matrix rows cols vals) = rows >= 0 && cols >= 0 && VS.length vals == rows * cols  -- | The maximum coefficient of the matrix-maxCoeff :: Matrix -> Double-maxCoeff Matrix{..} = I.cast $ VS.maximum m_vals+maxCoeff :: I.Elem a b => Matrix a b -> a+maxCoeff = fold1' max  -- | The minimum coefficient of the matrix-minCoeff :: Matrix -> Double-minCoeff Matrix{..} = I.cast $ VS.minimum m_vals+minCoeff :: I.Elem a b => Matrix a b -> a+minCoeff = fold1' min  -- | Top @N@ rows of matrix-topRows :: Int -> Matrix -> Matrix-topRows rows m@Matrix{..} = block 0 0 rows m_cols m+topRows :: I.Elem a b => Int -> Matrix a b -> Matrix a b+topRows n m@(Matrix _ cols _) = block 0 0 n cols m  -- | Bottom @N@ rows of matrix-bottomRows :: Int -> Matrix -> Matrix-bottomRows rows m@Matrix{..} = block (m_rows - rows) 0 rows m_cols m+bottomRows :: I.Elem a b => Int -> Matrix a b -> Matrix a b+bottomRows n m@(Matrix rows cols _) = block (rows - n) 0 n cols m  -- | Left @N@ columns of matrix-leftCols :: Int -> Matrix -> Matrix-leftCols cols m@Matrix{..} = block 0 0 m_rows cols m+leftCols :: I.Elem a b => Int -> Matrix a b -> Matrix a b+leftCols n m@(Matrix rows _ _) = block 0 0 rows n m  -- | Right @N@ columns of matrix-rightCols :: Int -> Matrix -> Matrix-rightCols cols m@Matrix{..} = block 0 (m_cols - cols) m_rows cols m+rightCols :: I.Elem a b => Int -> Matrix a b -> Matrix a b+rightCols n m@(Matrix rows cols _) = block 0 (cols - n) rows n m  -- | Construct matrix from a list of rows, column count is detected as maximum row length. Missing values are filled with 0-fromList :: [[Double]] -> Matrix+fromList :: forall a b . I.Elem a b => [[a]] -> Matrix a b fromList list = Matrix rows cols vals where     rows = length list     cols = maximum $ P.map length list     vals = VS.create $ do-        vm <- VSM.replicate (rows * cols) 0+        vm <- VSM.replicate (rows * cols) (I.cast (0::a))         forM_ (zip [0..] list) $ \(row, vals) ->             forM_ (zip [0..] vals) $ \(col, val) ->                 VSM.write vm (col * rows + row) (I.cast val)         return vm  -- | Convert matrix to a list of rows-toList :: Matrix -> [[Double]]-toList Matrix{..} = [[I.cast $ m_vals VS.! (col * m_rows + row) | col <- [0..pred m_cols]] | row <- [0..pred m_rows]]+toList :: I.Elem a b => Matrix a b -> [[a]]+toList (Matrix rows cols vals) = [[I.cast $ vals VS.! (col * rows + row) | col <- [0..pred cols]] | row <- [0..pred rows]] --- | Create matrix using generator function f :: row -> col -> val-generate :: Int -> Int -> (Int -> Int -> Double) -> Matrix+-- | [generate rows cols (λ row col -> val)]+--+-- Create matrix using generator function @f row col val@+--+generate :: I.Elem a b => Int -> Int -> (Int -> Int -> a) -> Matrix a b generate rows cols f = Matrix rows cols $ VS.create $ do     vals <- VSM.new (rows * cols)     forM_ [0..pred rows] $ \row ->@@ -251,74 +276,74 @@     return vals  -- | The sum of all coefficients of the matrix-sum :: Matrix -> Double-sum = _prop I.c_sum+sum :: I.Elem a b => Matrix a b -> a+sum = _prop I.sum  -- | The product of all coefficients of the matrix-prod :: Matrix -> Double-prod = _prop I.c_prod+prod :: I.Elem a b => Matrix a b -> a+prod = _prop I.prod  -- | The mean of all coefficients of the matrix-mean :: Matrix -> Double-mean = _prop I.c_prod+mean :: I.Elem a b => Matrix a b -> a+mean = _prop I.prod  -- | The trace of a matrix is the sum of the diagonal coefficients and can also be computed as sum (diagonal m)-trace :: Matrix -> Double-trace = _prop I.c_trace+trace :: I.Elem a b => Matrix a b -> a+trace = _prop I.trace  -- | Applied to a predicate and a matrix, all determines if all elements of the matrix satisfies the predicate-all :: (Double -> Bool) -> Matrix -> Bool-all f = VS.all (f . I.cast) . m_vals+all :: I.Elem a b => (a -> Bool) -> Matrix a b -> Bool+all f = VS.all (f . I.cast) . vals  -- | Applied to a predicate and a matrix, any determines if any element of the matrix satisfies the predicate-any :: (Double -> Bool) -> Matrix -> Bool-any f = VS.any (f . I.cast) . m_vals+any :: I.Elem a b => (a -> Bool) -> Matrix a b -> Bool+any f = VS.any (f . I.cast) . vals  -- | Returns the number of coefficients in a given matrix that evaluate to true-count :: (Double -> Bool) -> Matrix -> Int-count f = VS.foldl' (\n x -> if f (I.cast x) then succ n else n) 0 . m_vals+count :: I.Elem a b => (a -> Bool) -> Matrix a b -> Int+count f = VS.foldl' (\n x -> if f (I.cast x) then succ n else n) 0 . vals  {-| For vectors, the l2 norm, and for matrices the Frobenius norm.     In both cases, it consists in the square root of the sum of the square of all the matrix entries.     For vectors, this is also equals to the square root of the dot product of this with itself. -}-norm :: Matrix -> Double-norm = _prop I.c_norm+norm :: I.Elem a b => Matrix a b -> a+norm = _prop I.norm  -- | For vectors, the squared l2 norm, and for matrices the Frobenius norm. In both cases, it consists in the sum of the square of all the matrix entries. For vectors, this is also equals to the dot product of this with itself.-squaredNorm :: Matrix -> Double-squaredNorm = _prop I.c_squaredNorm+squaredNorm :: I.Elem a b => Matrix a b -> a+squaredNorm = _prop I.squaredNorm  -- | The l2 norm of the matrix using the Blue's algorithm. A Portable Fortran Program to Find the Euclidean Norm of a Vector, ACM TOMS, Vol 4, Issue 1, 1978.-blueNorm :: Matrix -> Double-blueNorm = _prop I.c_blueNorm+blueNorm :: I.Elem a b => Matrix a b -> a+blueNorm = _prop I.blueNorm  -- | The l2 norm of the matrix avoiding undeflow and overflow. This version use a concatenation of hypot calls, and it is very slow.-hypotNorm :: Matrix -> Double-hypotNorm = _prop I.c_hypotNorm+hypotNorm :: I.Elem a b => Matrix a b -> a+hypotNorm = _prop I.hypotNorm  -- | The determinant of the matrix-determinant :: Matrix -> Double+determinant :: I.Elem a b => Matrix a b -> a determinant m-    | square m = _prop I.c_determinant m+    | square m = _prop I.determinant m     | otherwise = error "Matrix.determinant: non-square matrix"  -- | Adding two matrices by adding the corresponding entries together. You can use @(+)@ function as well.-add :: Matrix -> Matrix -> Matrix+add :: I.Elem a b => Matrix a b -> Matrix a b -> Matrix a b add m1 m2-    | dims m1 == dims m2 = _binop const I.c_add m1 m2+    | dims m1 == dims m2 = _binop const I.add m1 m2     | otherwise = error "Matrix.add: matrices should have the same size"  -- | Subtracting two matrices by subtracting the corresponding entries together. You can use @(-)@ function as well.-sub :: Matrix -> Matrix -> Matrix+sub :: I.Elem a b => Matrix a b -> Matrix a b -> Matrix a b sub m1 m2-    | dims m1 == dims m2 = _binop const I.c_sub m1 m2+    | dims m1 == dims m2 = _binop const I.sub m1 m2     | otherwise = error "Matrix.add: matrices should have the same size"  -- | Matrix multiplication. You can use @(*)@ function as well.-mul :: Matrix -> Matrix -> Matrix+mul :: I.Elem a b => Matrix a b -> Matrix a b -> Matrix a b mul m1 m2-    | m_cols m1 == m_rows m2 = _binop (\(rows, _) (_, cols) -> (rows, cols)) I.c_mul m1 m2+    | cols m1 == rows m2 = _binop (\(rows, _) (_, cols) -> (rows, cols)) I.mul m1 m2     | otherwise = error "Matrix.mul: number of columns for lhs matrix should be the same as number of rows for rhs matrix"  {- | Apply a given function to each element of the matrix.@@ -338,8 +363,8 @@ 30.0    40.0  -}-map :: (Double -> Double) -> Matrix -> Matrix-map f m@Matrix{..} = m { m_vals = VS.map (I.cast . f . I.cast) m_vals }+map :: I.Elem a b => (a -> a) -> Matrix a b -> Matrix a b+map f (Matrix rows cols vals) = Matrix rows cols (VS.map (I.cast . f . I.cast) vals)  {- | Apply a given function to each element of the matrix. @@ -361,104 +386,115 @@  -} -imap :: (Int -> Int -> Double -> Double) -> Matrix -> Matrix-imap f m@Matrix{..} = m { m_vals = VS.imap (\n -> let (c, r) = divMod n m_rows in I.cast . f r c . I.cast) m_vals }+imap :: I.Elem a b => (Int -> Int -> a -> a) -> Matrix a b -> Matrix a b+imap f (Matrix rows cols vals) = Matrix rows cols (VS.imap (\n -> let (c, r) = divMod n rows in I.cast . f r c . I.cast) vals)  -- | Upper trinagle of the matrix-upperTriangule :: Matrix -> Matrix-upperTriangule = imap (\row col val -> if row <= col then val else 0)+upperTriangle :: I.Elem a b => Matrix a b -> Matrix a b+upperTriangle = imap (\row col val -> if row <= col then val else 0)  -- | Lower trinagle of the matrix-lowerTriangule :: Matrix -> Matrix-lowerTriangule = imap (\row col val -> if row >= col then val else 0)+lowerTriangle :: I.Elem a b => Matrix a b -> Matrix a b+lowerTriangle = imap (\row col val -> if row >= col then val else 0)  -- | Filter elements in the matrix. Filtered elements will be replaced by 0-filter :: (Double -> Bool) -> Matrix -> Matrix+filter :: I.Elem a b => (a -> Bool) -> Matrix a b -> Matrix a b filter f = map (\x -> if f x then x else 0)  -- | Filter elements in the matrix. Filtered elements will be replaced by 0-ifilter :: (Int -> Int -> Double -> Bool) -> Matrix -> Matrix+ifilter :: I.Elem a b => (Int -> Int -> a -> Bool) -> Matrix a b -> Matrix a b ifilter f = imap (\r c x -> if f r c x then x else 0)  -- | Reduce matrix using user provided function applied to each element.-fold :: (a -> Double -> a) -> a -> Matrix -> a-fold f a Matrix{..} = VS.foldl (\a x -> f a (I.cast x)) a m_vals+fold :: I.Elem a b => (c -> a -> c) -> c -> Matrix a b -> c+fold f a (Matrix _ _ vals) = VS.foldl (\a x -> f a (I.cast x)) a vals  -- | Reduce matrix using user provided function applied to each element. This is strict version of 'fold'-fold' :: (a -> Double -> a) -> a -> Matrix -> a-fold' f a Matrix{..} = VS.foldl' (\a x -> f a (I.cast x)) a m_vals+fold' :: I.Elem a b => (c -> a -> c) -> c -> Matrix a b -> c+fold' f a (Matrix _ _ vals) = VS.foldl' (\a x -> f a (I.cast x)) a vals  -- | Reduce matrix using user provided function applied to each element and it's index-ifold :: (Int -> Int -> a -> Double -> a) -> a -> Matrix -> a-ifold f a Matrix{..} = VS.ifoldl (\a n x -> let (c,r) = divMod n m_rows in f r c a (I.cast x)) a m_vals+ifold :: I.Elem a b => (Int -> Int -> c -> a -> c) -> c -> Matrix a b -> c+ifold f a (Matrix rows _ vals) = VS.ifoldl (\a n x -> let (c,r) = divMod n rows in f r c a (I.cast x)) a vals  -- | Reduce matrix using user provided function applied to each element and it's index. This is strict version of 'ifold'-ifold' :: (Int -> Int -> a -> Double -> a) -> a -> Matrix -> a-ifold' f a Matrix{..} = VS.ifoldl' (\a n x -> let (c,r) = divMod n m_rows in f r c a (I.cast x)) a m_vals+ifold' :: I.Elem a b => (Int -> Int -> c -> a -> c) -> c -> Matrix a b -> c+ifold' f a (Matrix rows _ vals) = VS.ifoldl' (\a n x -> let (c,r) = divMod n rows in f r c a (I.cast x)) a vals +-- | Reduce matrix using user provided function applied to each element.+fold1 :: I.Elem a b => (a -> a -> a) -> Matrix a b -> a+fold1 f = foldl1 f . P.map I.cast . VS.toList . vals++-- | Reduce matrix using user provided function applied to each element. This is strict version of 'fold'+fold1' :: I.Elem a b => (a -> a -> a) -> Matrix a b -> a+fold1' f = foldl1' f . P.map I.cast . VS.toList . vals+ -- | Diagonal of the matrix-diagonal :: Matrix -> Matrix-diagonal = _unop (\(rows, cols) -> (min rows cols, 1)) I.c_diagonal+diagonal :: I.Elem a b => Matrix a b -> Matrix a b+diagonal = _unop (\(rows, cols) -> (min rows cols, 1)) I.diagonal  {- | Inverse of the matrix  For small fixed sizes up to 4x4, this method uses cofactors. In the general case, this method uses PartialPivLU decomposition -}-inverse :: Matrix -> Matrix-inverse m@Matrix{..}-    | m_rows == m_cols = _unop id I.c_inverse m+inverse :: I.Elem a b => Matrix a b -> Matrix a b+inverse m+    | square m = _unop id I.inverse m     | otherwise = error "Matrix.inverse: non-square matrix"  -- | Adjoint of the matrix-adjoint :: Matrix -> Matrix-adjoint = _unop swap I.c_adjoint+adjoint :: I.Elem a b => Matrix a b -> Matrix a b+adjoint = _unop swap I.adjoint  -- | Transpose of the matrix-transpose :: Matrix -> Matrix-transpose = _unop swap I.c_transpose+transpose :: I.Elem a b => Matrix a b -> Matrix a b+transpose = _unop swap I.transpose  -- | Conjugate of the matrix-conjugate :: Matrix -> Matrix-conjugate = _unop id I.c_conjugate+conjugate :: I.Elem a b => Matrix a b -> Matrix a b+conjugate = _unop id I.conjugate  -- | Nomalize the matrix by deviding it on its 'norm'-normalize :: Matrix -> Matrix-normalize Matrix{..} = I.performIO $ do-    vals <- VS.thaw m_vals+normalize :: I.Elem a b => Matrix a b -> Matrix a b+normalize (Matrix rows cols vals) = I.performIO $ do+    vals <- VS.thaw vals     VSM.unsafeWith vals $ \p ->-        I.call $ I.c_normalize p (I.cast m_rows) (I.cast m_cols)-    Matrix m_rows m_cols <$> VS.unsafeFreeze vals+        I.call $ I.normalize p (I.cast rows) (I.cast cols)+    Matrix rows cols <$> VS.unsafeFreeze vals + -- | Apply a destructive operation to a matrix. The operation will be performed in place if it is safe to do so and will modify a copy of the matrix otherwise.-modify :: (forall s. M.MMatrix s -> ST s ()) -> Matrix -> Matrix-modify f m@Matrix{..} = m { m_vals = VS.modify (f . M.MMatrix m_rows m_cols) m_vals } where+modify :: I.Elem a b => (forall s. M.MMatrix a b s -> ST s ()) -> Matrix a b -> Matrix a b+modify f (Matrix rows cols vals) = Matrix rows cols (VS.modify (f . M.MMatrix rows cols) vals)  -- | Yield an immutable copy of the mutable matrix-freeze :: PrimMonad m => M.MMatrix (PrimState m) -> m Matrix-freeze M.MMatrix{..} = VS.freeze mm_vals >>= return . Matrix mm_rows mm_cols+freeze :: I.Elem a b => PrimMonad m => M.MMatrix a b (PrimState m) -> m (Matrix a b)+freeze (M.MMatrix mrows mcols mvals) = VS.freeze mvals >>= return . Matrix mrows mcols  -- | Yield a mutable copy of the immutable matrix-thaw :: PrimMonad m => Matrix -> m (M.MMatrix (PrimState m))-thaw Matrix{..} = VS.thaw m_vals >>= return . M.MMatrix m_rows m_cols+thaw :: I.Elem a b => PrimMonad m => Matrix a b -> m (M.MMatrix a b (PrimState m))+thaw (Matrix rows cols vals) = VS.thaw vals >>= return . M.MMatrix rows cols  -- | Unsafe convert a mutable matrix to an immutable one without copying. The mutable matrix may not be used after this operation.-unsafeFreeze :: PrimMonad m => M.MMatrix (PrimState m) -> m Matrix-unsafeFreeze M.MMatrix{..} = VS.unsafeFreeze mm_vals >>= return . Matrix mm_rows mm_cols+unsafeFreeze :: I.Elem a b => PrimMonad m => M.MMatrix a b (PrimState m) -> m (Matrix a b)+unsafeFreeze (M.MMatrix mrows mcols mvals) = VS.unsafeFreeze mvals >>= return . Matrix mrows mcols  -- | Unsafely convert an immutable matrix to a mutable one without copying. The immutable matrix may not be used after this operation.-unsafeThaw :: PrimMonad m => Matrix -> m (M.MMatrix (PrimState m))-unsafeThaw Matrix{..} = VS.unsafeThaw m_vals >>= return . M.MMatrix m_rows m_cols+unsafeThaw :: I.Elem a b => PrimMonad m => Matrix a b -> m (M.MMatrix a b (PrimState m))+unsafeThaw (Matrix rows cols vals) = VS.unsafeThaw vals >>= return . M.MMatrix rows cols  -- | Pass a pointer to the matrix's data to the IO action. The data may not be modified through the pointer.-unsafeWith :: Matrix -> (Ptr CDouble -> CInt -> CInt -> IO a) -> IO a-unsafeWith m@Matrix{..} f+unsafeWith :: I.Elem a b => Matrix a b -> (Ptr b -> CInt -> CInt -> IO c) -> IO c+unsafeWith m@(Matrix rows cols vals) f     | not (valid m) = fail "matrix layout is invalid"-    | otherwise = VS.unsafeWith m_vals $ \p -> f p (I.cast m_rows) (I.cast m_cols)+    | otherwise = VS.unsafeWith vals $ \p -> f p (I.cast rows) (I.cast cols) -_prop :: (Ptr CDouble -> CInt -> CInt -> IO CDouble) -> Matrix -> Double-_prop f m = I.cast $ I.performIO $ unsafeWith m f+_prop :: I.Elem a b => (Ptr b -> Ptr b -> CInt -> CInt -> IO CString) -> Matrix a b -> a+_prop f m = I.cast $ I.performIO $ alloca $ \p -> do+    I.call $ unsafeWith m (f p)+    peek p -_binop :: ((Int, Int) -> (Int, Int) -> (Int, Int)) -> (Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString) -> Matrix -> Matrix -> Matrix+_binop :: I.Elem a b => ((Int, Int) -> (Int, Int) -> (Int, Int)) -> (Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString) -> Matrix a b -> Matrix a b -> Matrix a b _binop f g m1 m2 = I.performIO $ do     m0 <- uncurry M.new $ f (dims m1) (dims m2)     M.unsafeWith m0 $ \vals0 rows0 cols0 ->@@ -470,7 +506,7 @@                     vals2 rows2 cols2     unsafeFreeze m0 -_unop :: ((Int,Int) -> (Int,Int)) -> (Ptr CDouble -> CInt -> CInt -> Ptr CDouble -> CInt -> CInt -> IO CString) -> Matrix -> Matrix+_unop :: I.Elem a b => ((Int,Int) -> (Int,Int)) -> (Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString) -> Matrix a b -> Matrix a b _unop f g m1 = I.performIO $ do     m0 <- uncurry M.new $ f (dims m1)     M.unsafeWith m0 $ \vals0 rows0 cols0 ->

@@ -1,11 +1,17 @@-{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE RecordWildCards, ScopedTypeVariables #-} module Data.Eigen.Matrix.Mutable (     MMatrix(..),+    MMatrixXf,+    MMatrixXd,+    MMatrixXcf,+    MMatrixXcd,     IOMatrix,     STMatrix,     -- * Construction     new,     replicate,+    -- * Consistency check+    valid,     -- * Accessing individual elements     read,     write,@@ -22,46 +28,56 @@ import Control.Monad.Primitive import Foreign.Ptr import Foreign.C.Types+import Data.Complex import Text.Printf import qualified Data.Vector.Storable.Mutable as VSM import qualified Data.Eigen.Internal as I  -- | Mutable matrix. You can modify elements-data MMatrix s = MMatrix {+data MMatrix a b s = MMatrix {     mm_rows :: Int,     mm_cols :: Int,-    mm_vals :: VSM.MVector s CDouble+    mm_vals :: VSM.MVector s b } -type IOMatrix = MMatrix RealWorld-type STMatrix s = MMatrix s+-- | Alias for single precision mutable matrix+type MMatrixXf = MMatrix Float CFloat+-- | Alias for double precision mutable matrix+type MMatrixXd = MMatrix Double CDouble+-- | Alias for single previsiom mutable matrix of complex numbers+type MMatrixXcf = MMatrix (Complex Float) (I.CComplex CFloat)+-- | Alias for double prevision mutable matrix of complex numbers+type MMatrixXcd = MMatrix (Complex Double) (I.CComplex CDouble) +type IOMatrix a b = MMatrix a b RealWorld+type STMatrix a b s = MMatrix a b s+ -- | Verify matrix dimensions and memory layout-valid :: MMatrix s -> Bool+valid :: I.Elem a b => MMatrix a b s -> Bool valid MMatrix{..} = mm_rows >= 0 && mm_cols >= 0 && VSM.length mm_vals == mm_rows * mm_cols  -- | Create a mutable matrix of the given size and fill it with 0 as an initial value.-new :: PrimMonad m => Int -> Int -> m (MMatrix (PrimState m))+new :: forall a b m . (PrimMonad m, I.Elem a b) => Int -> Int -> m (MMatrix a b (PrimState m)) new rows cols = do-    vals <- VSM.replicate (rows * cols) 0+    vals <- VSM.replicate (rows * cols) (I.cast (0::a))     return $ MMatrix rows cols vals  -- | Create a mutable matrix of the given size and fill it with as an initial value.-replicate :: PrimMonad m => Int -> Int -> m (MMatrix (PrimState m))+replicate :: forall a b m . (PrimMonad m, I.Elem a b) => Int -> Int -> m (MMatrix a b (PrimState m)) replicate rows cols = do-    vals <- VSM.replicate (rows * cols) 0+    vals <- VSM.replicate (rows * cols) (I.cast (0::a))     return $ MMatrix rows cols vals  -- | Set all elements of the matrix to the given value-set :: PrimMonad m => (MMatrix (PrimState m)) -> Double -> m ()+set :: (PrimMonad m, I.Elem a b) => (MMatrix a b (PrimState m)) -> a -> m () set MMatrix{..} val = VSM.set mm_vals (I.cast val)  -- | Copy a matrix. The two matrices must have the same length and may not overlap.-copy :: PrimMonad m => (MMatrix (PrimState m)) -> (MMatrix (PrimState m)) -> m ()+copy :: (PrimMonad m, I.Elem a b) => (MMatrix a b (PrimState m)) -> (MMatrix a b (PrimState m)) -> m () copy m1 m2 = VSM.copy (mm_vals m1) (mm_vals m2)  -- | Yield the element at the given position.-read :: PrimMonad m => MMatrix (PrimState m) -> Int -> Int -> m Double+read :: (PrimMonad m, I.Elem a b) => MMatrix a b (PrimState m) -> Int -> Int -> m a read mm@MMatrix{..} row col     | not (valid mm) = fail "MMatrix.read: matrix layout is invalid"     | row < 0 || row >= mm_rows = fail $ printf "MMatrix.read: row %d is out of bounds [0..%d)" row mm_rows@@ -69,7 +85,7 @@     | otherwise = unsafeRead mm row col  -- | Replace the element at the given position.-write :: PrimMonad m => MMatrix (PrimState m) -> Int -> Int -> Double -> m ()+write :: (PrimMonad m, I.Elem a b) => MMatrix a b (PrimState m) -> Int -> Int -> a -> m () write mm@MMatrix{..} row col val     | not (valid mm) = fail "MMatrix.write: matrix layout is invalid"     | row < 0 || row >= mm_rows = fail $ printf "MMatrix.write: row %d is out of bounds [0..%d)" row mm_rows@@ -77,15 +93,15 @@     | otherwise = unsafeWrite mm row col val  -- | Yield the element at the given position. No bounds checks are performed.-unsafeRead :: PrimMonad m => MMatrix (PrimState m) -> Int -> Int -> m Double+unsafeRead :: (PrimMonad m, I.Elem a b) => MMatrix a b (PrimState m) -> Int -> Int -> m a unsafeRead MMatrix{..} row col = VSM.unsafeRead mm_vals (col * mm_rows + row) >>= \val -> return (I.cast val)  -- | Replace the element at the given position. No bounds checks are performed.-unsafeWrite :: PrimMonad m => MMatrix (PrimState m) -> Int -> Int -> Double -> m ()+unsafeWrite :: (PrimMonad m, I.Elem a b) => MMatrix a b (PrimState m) -> Int -> Int -> a -> m () unsafeWrite MMatrix{..} row col val = VSM.unsafeWrite mm_vals (col * mm_rows + row) (I.cast val)  -- | Pass a pointer to the matrix's data to the IO action. Modifying data through the pointer is unsafe if the matrix could have been frozen before the modification.-unsafeWith :: IOMatrix -> (Ptr CDouble -> CInt -> CInt -> IO a) -> IO a+unsafeWith :: I.Elem a b => IOMatrix a b -> (Ptr b -> CInt -> CInt -> IO c) -> IO c unsafeWith mm@MMatrix{..} f     | not (valid mm) = fail "mutable matrix layout is invalid"     | otherwise = VSM.unsafeWith mm_vals $ \p -> f p (I.cast mm_rows) (I.cast mm_cols)

@@ -2,14 +2,14 @@ Some Eigen's algorithms can exploit the multiple cores present in your hardware. To this end, it is enough to enable OpenMP on your compiler, for instance: GCC: -fopenmp. You can control the number of thread that will be used using either the OpenMP API or Eiegn's API using the following priority:-@-OMP_NUM_THREADS=n ./my_program-setNbThreads n-@++1. OMP_NUM_THREADS=n ./my_program+2. setNbThreads n+ Unless setNbThreads has been called, Eigen uses the number of threads specified by OpenMP. You can restore this bahavior by calling @setNbThreads n@ -Currently, the following algorithms can make use of multi-threading: general matrix - matrix products PartialPivLU+Currently, the following algorithms can make use of multi-threading: general matrix - matrix products PartialPivLU. -}  module Data.Eigen.Parallel where

@@ -20,190 +20,200 @@     throw eigen_assert_exception(os.str()); } -typedef Map< Matrix<double,Dynamic,Dynamic> > MapMatrix;-typedef Map< Matrix<double,1,Dynamic> > MapVector;-typedef Map< Vector3d > MapVector3d;-typedef Map< Vector4d > MapVector4d;--extern "C" {--const char* eigen_add(-    double* data, int rows, int cols,-    const double* data1, int rows1, int cols1,-    const double* data2, int rows2, int cols2)-{-    GUARD_START-    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1) + MapMatrix(data2, rows2, cols2);-    GUARD_END-}--const char* eigen_sub(-    double* data, int rows, int cols,-    const double* data1, int rows1, int cols1,-    const double* data2, int rows2, int cols2)-{-    GUARD_START-    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1) - MapMatrix(data2, rows2, cols2);-    GUARD_END-}--const char* eigen_mul(-    double* data, int rows, int cols,-    const double* data1, int rows1, int cols1,-    const double* data2, int rows2, int cols2)-{-    GUARD_START-    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1) * MapMatrix(data2, rows2, cols2);-    GUARD_END-}--double eigen_norm(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).norm();  }-double eigen_squaredNorm(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).squaredNorm(); }-double eigen_blueNorm(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).squaredNorm(); }-double eigen_hypotNorm(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).hypotNorm(); }-double eigen_sum(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).sum(); }-double eigen_prod(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).prod(); }-double eigen_mean(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).mean(); }-double eigen_trace(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).trace(); }-double eigen_determinant(const double* data, int rows, int cols) { return MapMatrix(data, rows, cols).determinant(); }+typedef float T0;+typedef double T1;+typedef std::complex<float> T2;+typedef std::complex<double> T3; -const char* eigen_inverse(double* data, int rows, int cols, const double* data1, int rows1, int cols1)-{-    GUARD_START-    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1).inverse();-    GUARD_END+template <class T>+Map< Matrix<T,Dynamic,Dynamic> > matrix(void* p, int r, int c) {+    return Map< Matrix<T,Dynamic,Dynamic> >((T*)p, r, c); } -const char* eigen_adjoint(double* data, int rows, int cols, const double* data1, int rows1, int cols1)-{-    GUARD_START-    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1).adjoint();-    GUARD_END+template <class T>+Map< Matrix<T,Dynamic,Dynamic> > matrix(const void* p, int r, int c) {+    return Map< Matrix<T,Dynamic,Dynamic> >((const T*)p, r, c); } -const char* eigen_conjugate(double* data, int rows, int cols, const double* data1, int rows1, int cols1)-{-    GUARD_START-    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1).conjugate();-    GUARD_END+#define BINOP(name,op) \+extern "C" const char* eigen_##name(\+    int code,\+    void* p, int r, int c,\+    const void* p1, int r1, int c1,\+    const void* p2, int r2, int c2)\+{\+    GUARD_START\+    switch (code) {\+        case 0: matrix<T0>(p,r,c) = matrix<T0>(p1,r1,c1) op matrix<T0>(p2,r2,c2); break;\+        case 1: matrix<T1>(p,r,c) = matrix<T1>(p1,r1,c1) op matrix<T1>(p2,r2,c2); break;\+        case 2: matrix<T2>(p,r,c) = matrix<T2>(p1,r1,c1) op matrix<T2>(p2,r2,c2); break;\+        case 3: matrix<T3>(p,r,c) = matrix<T3>(p1,r1,c1) op matrix<T3>(p2,r2,c2); break;\+    }\+    GUARD_END\ } -const char* eigen_diagonal(double* data, int rows, int cols, const double* data1, int rows1, int cols1)-{-    GUARD_START-    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1).diagonal();-    GUARD_END-}+BINOP(add,+);+BINOP(sub,-);+BINOP(mul,*); -const char* eigen_transpose(double* data, int rows, int cols, const double* data1, int rows1, int cols1)-{-    GUARD_START-    MapMatrix(data, rows, cols) = MapMatrix(data1, rows1, cols1).transpose();-    GUARD_END-}+#define PROP(name) \+extern "C" const char* eigen_##name(int code, void* q, const void* p, int r, int c) {\+        GUARD_START\+        switch (code) {\+            case 0: *(T0*)q = matrix<T0>(p,r,c).name(); break;\+            case 1: *(T1*)q = matrix<T1>(p,r,c).name(); break;\+            case 2: *(T2*)q = matrix<T2>(p,r,c).name(); break;\+            case 3: *(T3*)q = matrix<T3>(p,r,c).name(); break;\+        }\+        GUARD_END\+    } -const char* eigen_dot(double* retval, -    const double* data1, int size1,-    const double* data2, int size2)-{-    GUARD_START-    *retval = MapVector(data1, size1).dot(MapVector(data2, size2));-    GUARD_END-}+PROP(norm);+PROP(squaredNorm);+PROP(blueNorm);+PROP(hypotNorm);+PROP(sum);+PROP(prod);+PROP(mean);+PROP(trace);+PROP(determinant); -const char* eigen_cross(-    double* data,-    const double* data1,-    const double* data2)-{-    GUARD_START-    MapVector3d retval(data);-    retval = MapVector3d(data1).cross(MapVector3d(data2));-    GUARD_END-}+#define UNOP(name) \+extern "C" const char* eigen_##name(int code, void* p, int r, int c, const void* p1, int r1, int c1) {\+        GUARD_START\+        switch (code) {\+            case 0: matrix<T0>(p,r,c) = matrix<T0>(p1,r1,c1).name(); break;\+            case 1: matrix<T1>(p,r,c) = matrix<T1>(p1,r1,c1).name(); break;\+            case 2: matrix<T2>(p,r,c) = matrix<T2>(p1,r1,c1).name(); break;\+            case 3: matrix<T3>(p,r,c) = matrix<T3>(p1,r1,c1).name(); break;\+        }\+        GUARD_END\+    } -const char* eigen_cross3(-    double* data,-    const double* data1,-    const double* data2)-{-    GUARD_START-    MapVector4d retval(data);-    retval = MapVector4d(data1).cross3(MapVector4d(data2));-    GUARD_END-}+UNOP(inverse);+UNOP(adjoint);+UNOP(conjugate);+UNOP(diagonal);+UNOP(transpose); -const char* eigen_normalize(double* data, int rows, int cols)+extern "C" const char* eigen_normalize(int code, void* p, int r, int c) {     GUARD_START-    MapMatrix(data, rows, cols).normalize();+    switch (code) {+        case 0: matrix<T0>(p,r,c).normalize(); break;+        case 1: matrix<T1>(p,r,c).normalize(); break;+        case 2: matrix<T2>(p,r,c).normalize(); break;+        case 3: matrix<T3>(p,r,c).normalize(); break;+    }     GUARD_END } -const char* eigen_random(double* data, int rows, int cols)+extern "C" const char* eigen_random(int code, void* p, int r, int c) {     GUARD_START-    MapMatrix(data, rows, cols) = MatrixXd::Random(rows, cols);+    switch (code) {+        case 0: matrix<T0>(p,r,c) = MatrixXf::Random(r,c); break;+        case 1: matrix<T1>(p,r,c) = MatrixXd::Random(r,c); break;+        case 2: matrix<T2>(p,r,c) = MatrixXcf::Random(r,c); break;+        case 3: matrix<T3>(p,r,c) = MatrixXcd::Random(r,c); break;+    }     GUARD_END } -const char* eigen_rank(Decomposition d, int* r, const double* data, int rows, int cols) {-    GUARD_START-    MapMatrix A(data, rows, cols);+template <class T>+const char* rank(Decomposition d, int* v, const void* p, int r, int c) {+    typedef Map< Matrix<T,Dynamic,Dynamic> > MapMatrix;+    MapMatrix A((const T*)p,r,c);     switch (d) {         case ::FullPivLU:-            *r = A.fullPivLu().rank();+            *v = A.fullPivLu().rank();             break;         case ::ColPivHouseholderQR:-            *r = A.colPivHouseholderQr().rank();+            *v = A.colPivHouseholderQr().rank();             break;         case ::FullPivHouseholderQR:-            *r = A.fullPivHouseholderQr().rank();+            *v = A.fullPivHouseholderQr().rank();             break;         case ::JacobiSVD:-            *r = A.jacobiSvd(ComputeThinU | ComputeThinV).rank();+            *v = A.jacobiSvd(ComputeThinU | ComputeThinV).rank();             break;         default:             return strdup("Selected decomposition doesn't support rank revealing.");     }-    GUARD_END+    return 0; } -const char* eigen_kernel(Decomposition d, double** data0, int* rows0, int* cols0, const double* data1, int rows1, int cols1) {+extern "C" const char* eigen_rank(int code, Decomposition d, int* v, const void* p, int r, int c) {     GUARD_START+    switch (code) {+        case 0: return rank<T0>(d,v,p,r,c);+        case 1: return rank<T1>(d,v,p,r,c);+        case 2: return rank<T2>(d,v,p,r,c);+        case 3: return rank<T3>(d,v,p,r,c);+    }+    GUARD_END+}++template <class T>+const char* kernel(Decomposition d, void** p0, int* r0, int* c0, const void* p1, int r1, int c1) {+    typedef Map< Matrix<T,Dynamic,Dynamic> > MapMatrix;     if (d != ::FullPivLU)         return strdup("Selected decomposition doesn't support kernel revealing.");-    MapMatrix A(data1, rows1, cols1);-    MatrixXd B = A.fullPivLu().kernel();-    *rows0 = B.rows();-    *cols0 = B.cols();-    *data0 = (double*)malloc(*rows0 * *cols0 * sizeof(double));-    MapMatrix(*data0, *rows0, *cols0) = B;-    GUARD_END+    MapMatrix A((const T*)p1,r1,c1);+    Matrix<T,Dynamic,Dynamic> B = A.fullPivLu().kernel();+    *r0 = B.rows();+    *c0 = B.cols();+    *p0 = malloc(*r0 * *c0 * sizeof(T));+    MapMatrix((T*)*p0, *r0, *c0) = B;+    return 0; } -const char* eigen_image(Decomposition d, double** data0, int* rows0, int* cols0, const double* data1, int rows1, int cols1) {+extern "C" const char* eigen_kernel(int code, Decomposition d, void** p0, int* r0, int* c0, const void* p1, int r1, int c1) {     GUARD_START+    switch (code) {+        case 0: return kernel<T0>(d,p0,r0,c0,p1,r1,c1);+        case 1: return kernel<T1>(d,p0,r0,c0,p1,r1,c1);+        case 2: return kernel<T2>(d,p0,r0,c0,p1,r1,c1);+        case 3: return kernel<T3>(d,p0,r0,c0,p1,r1,c1);+    }+    GUARD_END+}++template <class T>+const char* image(Decomposition d, void** p0, int* r0, int* c0, const void* p1, int r1, int c1) {+    typedef Map< Matrix<T,Dynamic,Dynamic> > MapMatrix;     if (d != ::FullPivLU)         return strdup("Selected decomposition doesn't support image revealing.");-    MapMatrix A(data1, rows1, cols1);-    MatrixXd B = A.fullPivLu().image(A);-    *rows0 = B.rows();-    *cols0 = B.cols();-    *data0 = (double*)malloc(*rows0 * *cols0 * sizeof(double));-    MapMatrix(*data0, *rows0, *cols0) = B;+    MapMatrix A((const T*)p1,r1,c1);+    Matrix<T,Dynamic,Dynamic> B = A.fullPivLu().image(A);+    *r0 = B.rows();+    *c0 = B.cols();+    *p0 = malloc(*r0 * *c0 * sizeof(T));+    MapMatrix((T*)*p0, *r0, *c0) = B;+    return 0;+}++extern "C" const char* eigen_image(int code, Decomposition d, void** p0, int* r0, int* c0, const void* p1, int r1, int c1) {+    GUARD_START+    switch (code) {+        case 0: return image<T0>(d,p0,r0,c0,p1,r1,c1);+        case 1: return image<T1>(d,p0,r0,c0,p1,r1,c1);+        case 2: return image<T2>(d,p0,r0,c0,p1,r1,c1);+        case 3: return image<T3>(d,p0,r0,c0,p1,r1,c1);+    }     GUARD_END } -const char* eigen_solve(Decomposition d,-    double* px, int rx, int cx,-    const double* pa, int ra, int ca,-    const double* pb, int rb, int cb)+template <class T>+const char* solve(Decomposition d,+    void* px, int rx, int cx,+    const void* pa, int ra, int ca,+    const void* pb, int rb, int cb) {-    GUARD_START-    MapMatrix x(px, rx, cx);-    MapMatrix A(pa, ra, ca);-    MapMatrix b(pb, rb, cb);+    typedef Map< Matrix<T,Dynamic,Dynamic> > MapMatrix;+    MapMatrix x((T*)px, rx, cx);+    MapMatrix A((const T*)pa, ra, ca);+    MapMatrix b((const T*)pb, rb, cb);     switch (d) {         case ::PartialPivLU:             x = A.partialPivLu().solve(b);@@ -230,35 +240,64 @@             x = A.jacobiSvd(ComputeThinU | ComputeThinV).solve(b);             break;     }+    return 0;+}++extern "C" const char* eigen_solve(int code, Decomposition d,+    void* px, int rx, int cx,+    const void* pa, int ra, int ca,+    const void* pb, int rb, int cb)+{+    GUARD_START+    switch (code) {+        case 0: return solve<T0>(d,px,rx,cx,pa,ra,ca,pb,rb,cb);+        case 1: return solve<T1>(d,px,rx,cx,pa,ra,ca,pb,rb,cb);+        case 2: return solve<T2>(d,px,rx,cx,pa,ra,ca,pb,rb,cb);+        case 3: return solve<T3>(d,px,rx,cx,pa,ra,ca,pb,rb,cb);+    }     GUARD_END } +template <class T>+const char* relativeError(void* e,+    const void* px, int rx, int cx,+    const void* pa, int ra, int ca,+    const void* pb, int rb, int cb)+{+    typedef Map< Matrix<T,Dynamic,Dynamic> > MapMatrix;+    MapMatrix x((const T*)px, rx, cx);+    MapMatrix A((const T*)pa, ra, ca);+    MapMatrix b((const T*)pb, rb, cb);+    *(T*)e = (A*x - b).norm() / b.norm();+    return 0;+} -const char* eigen_relativeError(double* e,-    const double* px, int rx, int cx,-    const double* pa, int ra, int ca,-    const double* pb, int rb, int cb)++extern "C" const char* eigen_relativeError(int code, void* e,+    const void* px, int rx, int cx,+    const void* pa, int ra, int ca,+    const void* pb, int rb, int cb) {     GUARD_START-    MapMatrix x(px, rx, cx);-    MapMatrix A(pa, ra, ca);-    MapMatrix b(pb, rb, cb);-    *e = (A*x - b).norm() / b.norm();+    switch (code) {+        case 0: return relativeError<T0>(e,px,rx,cx,pa,ra,ca,pb,rb,cb);+        case 1: return relativeError<T1>(e,px,rx,cx,pa,ra,ca,pb,rb,cb);+        case 2: return relativeError<T2>(e,px,rx,cx,pa,ra,ca,pb,rb,cb);+        case 3: return relativeError<T3>(e,px,rx,cx,pa,ra,ca,pb,rb,cb);+    }     GUARD_END }  -bool eigen_initParallel() {+extern "C" bool eigen_initParallel() {     initParallel();     return true; } -void eigen_setNbThreads(int n) {+extern "C" void eigen_setNbThreads(int n) {     setNbThreads(n); } -int eigen_getNbThreads() {+extern "C" int eigen_getNbThreads() {     return nbThreads();-}- }

@@ -12,28 +12,25 @@  extern "C" { -const char* eigen_add(double*, int, int, const double*, int, int, const double*, int, int);-const char* eigen_sub(double*, int, int, const double*, int, int, const double*, int, int);-const char* eigen_mul(double*, int, int, const double*, int, int, const double*, int, int);-double eigen_norm(const double*, int, int);-double eigen_squaredNorm(const double*, int, int);-double eigen_blueNorm(const double*, int, int);-double eigen_hypotNorm(const double*, int, int);-double eigen_sum(const double*, int, int);-double eigen_prod(const double*, int, int);-double eigen_mean(const double*, int, int);-double eigen_trace(const double*, int, int);-double eigen_determinant(const double*, int, int);-const char* eigen_dot(double*, const double*, int, const double*, int);-const char* eigen_cross(double*, const double*, const double*);-const char* eigen_cross3(double*, const double*, const double*);-const char* eigen_diagonal(double*, int, int, const double*, int, int);-const char* eigen_transpose(double*, int, int, const double*, int, int);-const char* eigen_normalize(double*, int, int);-const char* eigen_random(double*, int, int);-const char* eigen_inverse(double*, int, int, const double*, int, int);-const char* eigen_conjugate(double*, int, int, const double*, int, int);-const char* eigen_adjoint(double*, int, int, const double*, int, int);+const char* eigen_add(int, void*, int, int, const void*, int, int, const void*, int, int);+const char* eigen_sub(int, void*, int, int, const void*, int, int, const void*, int, int);+const char* eigen_mul(int, void*, int, int, const void*, int, int, const void*, int, int);+const char* eigen_norm(int, void*, const void*, int, int);+const char* eigen_squaredNorm(int, void*, const void*, int, int);+const char* eigen_blueNorm(int, void*, const void*, int, int);+const char* eigen_hypotNorm(int, void*, const void*, int, int);+const char* eigen_sum(int, void*, const void*, int, int);+const char* eigen_prod(int, void*, const void*, int, int);+const char* eigen_mean(int, void*, const void*, int, int);+const char* eigen_trace(int, void*, const void*, int, int);+const char* eigen_determinant(int, void*, const void*, int, int);+const char* eigen_diagonal(int, void*, int, int, const void*, int, int);+const char* eigen_transpose(int, void*, int, int, const void*, int, int);+const char* eigen_normalize(int, void*, int, int);+const char* eigen_random(int, void*, int, int);+const char* eigen_inverse(int, void*, int, int, const void*, int, int);+const char* eigen_conjugate(int, void*, int, int, const void*, int, int);+const char* eigen_adjoint(int, void*, int, int, const void*, int, int); bool eigen_initParallel(); int eigen_getNbThreads(); void eigen_setNbThreads(int);@@ -53,23 +50,23 @@ 	// GeneralizedSelfAdjointEigenSolver }; -const char* eigen_rank(Decomposition d, int*, const double*, int, int);-const char* eigen_kernel(Decomposition d, double**, int*, int*, const double*, int, int);-const char* eigen_image(Decomposition d, double**, int*, int*, const double*, int, int);+const char* eigen_rank(int, Decomposition d, int*, const void*, int, int);+const char* eigen_kernel(int, Decomposition d, void**, int*, int*, const void*, int, int);+const char* eigen_image(int, Decomposition d, void**, int*, int*, const void*, int, int); -const char* eigen_solve(Decomposition d, // Ax=b -	double*, int, int, // x-	const double*, int, int, // A-	const double*, int, int); // b+const char* eigen_solve(int, Decomposition d, // Ax=b +	void*, int, int, // x+	const void*, int, int, // A+	const void*, int, int); // b  -const char* eigen_relativeError(double* e, -	const double*, int, int, // x-	const double*, int, int, // A-	const double*, int, int); // b;+const char* eigen_relativeError(int, void* e, +	const void*, int, int, // x+	const void*, int, int, // A+	const void*, int, int); // b;   } // end extern "C"  #define GUARD_START try { assert(inited); do {-#define GUARD_END } while(false); return 0; } catch (const eigen_assert_exception& ex) { return strdup(ex.what()); }+#define GUARD_END } while(false); return 0; } catch (const std::exception& ex) { return strdup(ex.what()); }

@@ -1,5 +1,5 @@ name:           eigen-version:        1.2.5+version:        2.0.0 homepage:       https://github.com/osidorkin/haskell-eigen synopsis:       Eigen C++ library (linear algebra: matrices, vectors, numerical solvers). description:    This module provides Haskell binding for Eigen C++ library.@@ -13,7 +13,7 @@                 Virtually any software may use Eigen. For example, closed-source software may use Eigen without having to disclose its own source code. Many proprietary and closed-source software projects are using Eigen right now, as well as many BSD-licensed projects.                 .                 Some parts of documentation strings are based or replicate original Eigen documentation which is available at <http://eigen.tuxfamily.org/dox/>.-category:       Data, Math, Algebra, Statistics, Algorithms+category:       Data, Math, Algebra, Statistics, Algorithms, Numeric license:        BSD3 license-file:   LICENSE copyright:      (c) 2013-2015, Oleg Sidorkin

@@ -2,13 +2,13 @@ import Data.Eigen.LA  main = do-    let a = fromList [[1,2,5],[2,1,4],[3,0,3]]+    let a = fromList [[1,2,5],[2,1,4],[3,0,3]] :: MatrixXf     putStrLn "Here is the matrix A:"     print a      putStrLn "The rank of A is:"     print $ rank FullPivLU a-   +     putStrLn "Here is a matrix whose columns form a basis of the null-space of A:"     print $ kernel FullPivLU a 

@@ -6,6 +6,7 @@  main = do     let+        a :: MatrixXd         a = fromList [             [1,3.02, 6.89],             [1,2.01, 5.39],@@ -13,7 +14,7 @@             [1,2.09, 5.55],             [1,2.58, 6.32]] -        b = fromList $ map return [-4.32,-3.79,-4.01,-3.86,-4.10]+        b = fromList $ Prelude.map return [-4.32,-3.79,-4.01,-3.86,-4.10]      print a     print b@@ -29,11 +30,12 @@     putStrLn "\n-2.34666 - 0.25349 x1 - 0.174965 x2"     putStrLn "done" -    print $ identity 4+    print $ (identity 4 :: MatrixXd)     print $ M.normalize a     print $ M.transpose a      let+        a :: MatrixXd         a = M.fromList [[0.68,  0.597,  -0.33],[-0.211,  0.823,  0.536],[ 0.566, -0.605, -0.444]]         b = M.inverse a     print a

@@ -3,6 +3,7 @@  main = do     let+        a :: MatrixXd         a = fromList [[1,2,3], [4,5,6], [7,8,10]]         b = fromList [[3],[3],[4]]         x = solve ColPivHouseholderQR a b