packages feed

eigen 2.0.1 → 2.1.0

raw patch · 9 files changed

+839/−251 lines, 9 filesdep +bytestring

Dependencies added: bytestring

Files

− Data/Eigen/Internal.hs 
@@ -1,147 +0,0 @@-{-# 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, 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-instance Elem (Complex Float) (CComplex CFloat) where-instance Elem (Complex Double) (CComplex CDouble) where--class Cast a b where-    cast :: a -> b--data CComplex a = CComplex !a !a--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 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)---performIO :: IO a -> a-performIO = unsafeDupablePerformIO--foreign import ccall "eigen-proxy.h free" c_freeString :: CString -> IO ()--call :: IO CString -> IO ()-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      :: CInt -> Ptr a -> CInt -> CInt -> IO CString-foreign import ccall "eigen-proxy.h eigen_identity"    c_identity    :: 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 -> 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-identity    :: 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))-identity    = c_identity    (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))
+ Data/Eigen/Internal.hsc view
@@ -0,0 +1,187 @@+{-# LANGUAGE MultiParamTypeClasses, ForeignFunctionInterface, ScopedTypeVariables, FunctionalDependencies, FlexibleInstances, EmptyDataDecls, CPP #-}+module Data.Eigen.Internal where++import Foreign.Ptr+import Foreign.ForeignPtr+import Foreign.Storable+import Foreign.C.Types+import Foreign.C.String+import Control.Monad+import Control.Applicative+import System.IO.Unsafe+import Data.Complex+import Data.IORef+import Data.Bits+import qualified Data.Vector.Storable as VS+import qualified Data.ByteString as BS+import qualified Data.ByteString.Lazy as BSL+import qualified Data.ByteString.Internal as BSI++class (Num 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+instance Elem (Complex Float) (CComplex CFloat) where+instance Elem (Complex Double) (CComplex CDouble) where++class Cast a b where+	cast :: a -> b++-- | Complex number for FFI with the same memory layout as std::complex\<T\>+data CComplex a = CComplex !a !a++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++data CTriplet a = CTriplet !CInt !CInt !a++instance Storable a => Storable (CTriplet a) where+	sizeOf _ = sizeOf (undefined :: a) + sizeOf (undefined :: CInt) * 2+	alignment _ = alignment (undefined :: CInt)+	poke p (CTriplet row col val) = do+		pokeElemOff (castPtr p) 0 row+		pokeElemOff (castPtr p) 1 col+		pokeByteOff p (sizeOf (undefined :: CInt) * 2) val+	peek p = CTriplet+		<$> peekElemOff (castPtr p) 0+		<*> peekElemOff (castPtr p) 1+		<*> peekByteOff p (sizeOf (undefined :: CInt) * 2)++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)++intSize :: Int+intSize = sizeOf (undefined :: CInt)++encodeInt :: CInt -> BS.ByteString+encodeInt x = BSI.unsafeCreate (sizeOf x) $ (`poke` x) . castPtr++decodeInt :: BS.ByteString -> CInt+decodeInt (BSI.PS fp fo fs)+	| fs == sizeOf x = x+	| otherwise = error "decodeInt: wrong buffer size"+	where x = performIO $ withForeignPtr fp $ peek . (`plusPtr` fo)++data CSparseMatrix a b+type CSparseMatrixPtr a b = Ptr (CSparseMatrix a b)++performIO :: IO a -> a+performIO = unsafeDupablePerformIO++plusForeignPtr :: ForeignPtr a -> Int -> ForeignPtr b+plusForeignPtr fp fo = castForeignPtr fp' where+	vs :: VS.Vector CChar+	vs = VS.unsafeFromForeignPtr (castForeignPtr fp) fo 0+	(fp', _) = VS.unsafeToForeignPtr0 vs++foreign import ccall "eigen-proxy.h free" c_freeString :: CString -> IO ()++call :: IO CString -> IO ()+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++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++magicCode :: Code a => a -> CInt+magicCode x = code x `xor` 0x45696730++#let api name, args = "foreign import ccall \"eigen_%s\" c_%s :: CInt -> %s\n%s :: forall b . Code b => %s\n%s = c_%s (code (undefined :: b))", #name, #name, args, #name, args, #name, #name++#api random, 		"Ptr b -> CInt -> CInt -> IO CString"+#api identity, 		"Ptr b -> CInt -> CInt -> IO CString"+#api add, 			"Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"+#api sub, 			"Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"+#api mul, 			"Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"+#api diagonal, 		"Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"+#api transpose, 	"Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"+#api inverse, 		"Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"+#api adjoint, 		"Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"+#api conjugate, 	"Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"+#api normalize, 	"Ptr b -> CInt -> CInt -> IO CString"+#api sum, 			"Ptr b -> Ptr b -> CInt -> CInt -> IO CString"+#api prod, 			"Ptr b -> Ptr b -> CInt -> CInt -> IO CString"+#api mean, 			"Ptr b -> Ptr b -> CInt -> CInt -> IO CString"+#api norm, 			"Ptr b -> Ptr b -> CInt -> CInt -> IO CString"+#api trace, 		"Ptr b -> Ptr b -> CInt -> CInt -> IO CString"+#api squaredNorm, 	"Ptr b -> Ptr b -> CInt -> CInt -> IO CString"+#api blueNorm, 		"Ptr b -> Ptr b -> CInt -> CInt -> IO CString"+#api hypotNorm, 	"Ptr b -> Ptr b -> CInt -> CInt -> IO CString"+#api determinant, 	"Ptr b -> Ptr b -> CInt -> CInt -> IO CString"+#api rank, 			"CInt -> Ptr CInt -> Ptr b -> CInt -> CInt -> IO CString"+#api image, 		"CInt -> Ptr (Ptr b) -> Ptr CInt -> Ptr CInt -> Ptr b -> CInt -> CInt -> IO CString"+#api kernel, 		"CInt -> Ptr (Ptr b) -> Ptr CInt -> Ptr CInt -> Ptr b -> CInt -> CInt -> IO CString"+#api solve, 		"CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"+#api relativeError, "Ptr b -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"++#let api2 name, args = "foreign import ccall \"eigen_%s\" c_%s :: CInt -> %s\n%s :: forall a b . Code b => %s\n%s = c_%s (code (undefined :: b))", #name, #name, args, #name, args, #name, #name++#api2 sparse_fromList, 		"CInt -> CInt -> Ptr (CTriplet b) -> CInt -> Ptr (CSparseMatrixPtr a b) -> IO CString"+#api2 sparse_toList, 		"CSparseMatrixPtr a b -> Ptr (CTriplet b) -> CInt -> IO CString"+#api2 sparse_free, 			"CSparseMatrixPtr a b -> IO CString"+#api2 sparse_compress, 		"CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"+#api2 sparse_uncompress, 	"CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"+#api2 sparse_isCompressed, 	"CSparseMatrixPtr a b -> Ptr CInt -> IO CString"+#api2 sparse_transpose, 	"CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"+#api2 sparse_adjoint, 		"CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"+#api2 sparse_pruned, 		"CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"+#api2 sparse_prunedRef, 	"CSparseMatrixPtr a b -> Ptr b -> Ptr (CSparseMatrixPtr a b) -> IO CString"+#api2 sparse_scale, 		"CSparseMatrixPtr a b -> Ptr b -> Ptr (CSparseMatrixPtr a b) -> IO CString"+#api2 sparse_lowerTriangle, "CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"+#api2 sparse_upperTriangle, "CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"+#api2 sparse_nonZeros, 		"CSparseMatrixPtr a b -> Ptr CInt -> IO CString"+#api2 sparse_innerSize, 	"CSparseMatrixPtr a b -> Ptr CInt -> IO CString"+#api2 sparse_outerSize, 	"CSparseMatrixPtr a b -> Ptr CInt -> IO CString"+#api2 sparse_coeff, 		"CSparseMatrixPtr a b -> CInt -> CInt -> Ptr b -> IO CString"+#api2 sparse_cols, 			"CSparseMatrixPtr a b -> Ptr CInt -> IO CString"+#api2 sparse_rows, 			"CSparseMatrixPtr a b -> Ptr CInt -> IO CString"+#api2 sparse_norm, 			"CSparseMatrixPtr a b -> Ptr b -> IO CString"+#api2 sparse_squaredNorm, 	"CSparseMatrixPtr a b -> Ptr b -> IO CString"+#api2 sparse_blueNorm, 		"CSparseMatrixPtr a b -> Ptr b -> IO CString"+#api2 sparse_add, 			"CSparseMatrixPtr a b -> CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"+#api2 sparse_sub, 			"CSparseMatrixPtr a b -> CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"+#api2 sparse_mul, 			"CSparseMatrixPtr a b -> CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"+#api2 sparse_block,			"CSparseMatrixPtr a b -> CInt -> CInt -> CInt -> CInt -> Ptr (CSparseMatrixPtr a b) -> IO CString"+#api2 sparse_fromMatrix,	"Ptr b -> CInt -> CInt -> Ptr (CSparseMatrixPtr a b) -> IO CString"+#api2 sparse_toMatrix,		"CSparseMatrixPtr a b -> Ptr b -> CInt -> CInt -> IO CString"++openStream :: BSL.ByteString -> IO (IORef BSL.ByteString)+openStream = newIORef++readStream :: IORef BSL.ByteString -> Int -> IO BS.ByteString+readStream ref size = readIORef ref >>= \a ->+	let (b,c) = BSL.splitAt (fromIntegral size) a+	in if BSL.length b /= fromIntegral size+		then error "readStream: stream exhausted"+		else do+			writeIORef ref c+			return . BS.concat . BSL.toChunks $ b++closeStream :: IORef BSL.ByteString -> IO ()+closeStream ref = BSL.null <$> readIORef ref >>= (`unless` fail "closeStream: stream underrun")++readInt :: IORef BSL.ByteString -> IO CInt+readInt st = decodeInt <$> readStream st intSize
Data/Eigen/LA.hs view
@@ -16,6 +16,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@@ -66,7 +67,7 @@ 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]] :: MatrixXd     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)
Data/Eigen/Matrix.hs view
@@ -1,11 +1,14 @@ {-# LANGUAGE RecordWildCards, MultiParamTypeClasses, Rank2Types, ScopedTypeVariables, FunctionalDependencies, GADTs #-} module Data.Eigen.Matrix (     -- * Matrix type+    -- | Matrix aliases follows Eigen naming convention     Matrix(..),     MatrixXf,     MatrixXd,     MatrixXcf,     MatrixXcd,+    I.Elem,+    I.CComplex,     valid,     -- * Matrix conversions     fromList,@@ -74,8 +77,12 @@     conjugate,     normalize,     modify,+    convert,     upperTriangle,     lowerTriangle,+    -- * Matrix serialization+    encode,+    decode,     -- * Mutable matrices     thaw,     freeze,@@ -86,11 +93,12 @@ ) where  import qualified Prelude as P+import qualified Data.List as L import Prelude hiding (null, sum, all, any, map, filter)-import Data.List (intercalate, foldl1') import Data.Tuple import Data.Complex hiding (conjugate) import Foreign.Ptr+import Foreign.ForeignPtr import Foreign.C.Types import Foreign.C.String import Foreign.Storable@@ -104,10 +112,12 @@ import qualified Data.Vector.Storable.Mutable as VSM import qualified Data.Eigen.Internal as I import qualified Data.Eigen.Matrix.Mutable as M+import qualified Data.ByteString.Lazy as BSL+import qualified Data.ByteString.Internal as BSI  -- | 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 a b where +data Matrix a b where     Matrix :: I.Elem a b => !Int -> !Int -> !(VS.Vector b) -> Matrix a b  -- | Alias for single precision matrix@@ -123,7 +133,7 @@ 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"]+        "\n", L.intercalate "\n" $ P.map (L.intercalate "\t" . P.map show) $ toList m, "\n"]   -- | Shortcuts for basic matrix math@@ -134,6 +144,7 @@     fromInteger = constant 1 1 . fromInteger     signum = map signum     abs = map abs+    negate = map negate  -- | Empty 0x0 matrix empty :: I.Elem a b => Matrix a b@@ -250,7 +261,7 @@ fromList :: I.Elem a b => [[a]] -> Matrix a b fromList list = Matrix rows cols vals where     rows = length list-    cols = maximum $ P.map length list+    cols = L.foldl' max 0 $ P.map length list     vals = VS.create $ do         vm <- VSM.replicate (rows * cols) (I.cast (0 `asTypeOf` (head (head list))))         forM_ (zip [0..] list) $ \(row, vals) ->@@ -349,7 +360,7 @@  Here is an example how to implement scalar matrix multiplication: ->>> let a = fromList [[1,2],[3,4]]+>>> let a = fromList [[1,2],[3,4]] :: MatrixXf  >>> a Matrix 2x2@@ -367,9 +378,9 @@  {- | Apply a given function to each element of the matrix. -Here is an example how getting upper triangular matrix can be implemented:+Here is an example how upper triangular matrix can be implemented: ->>> let a = fromList [[1,2,3],[4,5,6],[7,8,9]]+>>> let a = fromList [[1,2,3],[4,5,6],[7,8,9]] :: MatrixXf  >>> a Matrix 3x3@@ -426,7 +437,7 @@  -- | 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+fold1' f = L.foldl1' f . P.map I.cast . VS.toList . vals  -- | Diagonal of the matrix diagonal :: I.Elem a b => Matrix a b -> Matrix a b@@ -461,11 +472,14 @@         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 :: 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) +-- | Convert matrix to different type using user provided element converter+convert :: (I.Elem a b, I.Elem c d) => (a -> c) -> Matrix a b -> Matrix c d+convert f (Matrix rows cols vals) = Matrix rows cols $ VS.map (I.cast . f . I.cast) vals+ -- | Yield an immutable copy of the mutable matrix 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@@ -485,8 +499,33 @@ -- | Pass a pointer to the matrix's data to the IO action. The data may not be modified through the pointer. 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"+    | not (valid m) = fail "Matrix.unsafeWith: matrix layout is invalid"     | otherwise = VS.unsafeWith vals $ \p -> f p (I.cast rows) (I.cast cols)++-- | Encode the matrix as a lazy byte string+encode :: I.Elem a b => Matrix a b -> BSL.ByteString+encode m@(Matrix rows cols vals)+    | valid m = BSL.fromChunks [+            I.encodeInt (I.magicCode $ VS.head vals),+            I.encodeInt (I.cast rows),+            I.encodeInt (I.cast cols),+            let (fp, fs) = VS.unsafeToForeignPtr0 vals in BSI.PS (castForeignPtr fp) 0 (fs * sizeOf (VS.head vals))]+    | otherwise = error "Matrix.encode: matrix layout is invalid"++-- | Decode matrix from the lazy byte string+decode :: I.Elem a b => BSL.ByteString -> Matrix a b+decode st = Matrix rows cols vals where+    (rows, cols, vals) = I.performIO $ do+        st <- I.openStream st+        code <- I.readInt st+        when (code /= I.magicCode (VS.head vals)) $+            fail "Matrix.decode: wrong matrix type"+        rows <- I.cast <$> I.readInt st+        cols <- I.cast <$> I.readInt st+        BSI.PS fp fo _ <- I.readStream st (rows * cols * sizeOf (VS.head vals))+        I.closeStream st+        return (rows, cols, VS.unsafeFromForeignPtr0 (I.plusForeignPtr fp fo) (rows * cols))+  _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
Data/Eigen/Matrix/Mutable.hs view
@@ -20,6 +20,7 @@     -- * Modifying matrices     set,     copy,+    unsafeCopy,     -- * Raw pointers     unsafeWith ) where@@ -57,24 +58,27 @@ 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 :: 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) (I.cast (0::a))-    return $ MMatrix rows cols vals+new :: (PrimMonad m, I.Elem a b) => Int -> Int -> m (MMatrix a b (PrimState m))+new rows cols = replicate rows cols 0  -- | Create a mutable matrix of the given size and fill it with as an initial value.-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) (I.cast (0::a))+replicate :: (PrimMonad m, I.Elem a b) => Int -> Int -> a -> m (MMatrix a b (PrimState m))+replicate rows cols val = do+    vals <- VSM.replicate (rows * cols) (I.cast val)     return $ MMatrix rows cols vals  -- | Set all elements of the matrix to the given value 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 a matrix. The two matrices must have the same size and may not overlap. 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)+copy m1 m2+    | not (valid m1) = fail "MMatrix.copy: lhs matrix layout is invalid"+    | not (valid m2) = fail "MMatrix.copy: rhs matrix layout is invalid"+    | mm_rows m1 /= mm_rows m2 = fail "MMatrix.copy: matrices have different number of cols"+    | mm_cols m1 /= mm_cols m2 = fail "MMatrix.copy: matrices have different number of rows"+    | otherwise = VSM.copy (mm_vals m1) (mm_vals m2)  -- | Yield the element at the given position. read :: (PrimMonad m, I.Elem a b) => MMatrix a b (PrimState m) -> Int -> Int -> m a@@ -91,6 +95,10 @@     | row < 0 || row >= mm_rows = fail $ printf "MMatrix.write: row %d is out of bounds [0..%d)" row mm_rows     | col < 0 || col >= mm_cols = fail $ printf "MMatrix.write: col %d is out of bounds [0..%d)" col mm_cols     | otherwise = unsafeWrite mm row col val++-- | Copy a matrix. The two matrices must have the same size and may not overlap however no bounds check performaned to it may SEGFAULT for incorrect input.+unsafeCopy :: (PrimMonad m, I.Elem a b) => (MMatrix a b (PrimState m)) -> (MMatrix a b (PrimState m)) -> m ()+unsafeCopy m1 m2 = VSM.unsafeCopy (mm_vals m1) (mm_vals m2)  -- | Yield the element at the given position. No bounds checks are performed. unsafeRead :: (PrimMonad m, I.Elem a b) => MMatrix a b (PrimState m) -> Int -> Int -> m a
+ Data/Eigen/SparseMatrix.hs view
@@ -0,0 +1,339 @@+{-# LANGUAGE GADTs, ScopedTypeVariables #-}+module Data.Eigen.SparseMatrix (+    -- * SparseMatrix type+    -- | SparseMatrix aliases follows Eigen naming convention+    SparseMatrix,+    SparseMatrixXf,+    SparseMatrixXd,+    SparseMatrixXcf,+    SparseMatrixXcd,+    -- * Accessing matrix data+    cols,+    rows,+    coeff,+    (!),+    -- * Matrix conversions+    fromList,+    toList,+    fromDenseList,+    toDenseList,+    fromMatrix,+    toMatrix,+    -- * Matrix properties+    norm,+    squaredNorm,+    blueNorm,+    block,+    --compress,+    --uncompress,+    --compressed,+    nonZeros,+    --innerSize,+    --outerSize,+    -- * Basic matrix algebra+    add,+    sub,+    mul,+    -- * Matrix transformations+    pruned,+    scale,+    transpose,+    adjoint,+    --lowerTriangle,+    --upperTriangle,+    -- * Matrix serialization+    encode,+    decode,+) where++import qualified Prelude as P+import Prelude hiding (map)+import qualified Data.List as L+import Data.Complex+import Foreign.C.Types+import Foreign.C.String+import Foreign.Storable+import Foreign.Ptr+import Foreign.ForeignPtr+import Foreign.Marshal.Alloc+import Control.Monad+import Control.Applicative+import qualified Data.Eigen.Matrix as M+import qualified Data.Eigen.Matrix.Mutable as MM+import qualified Foreign.Concurrent as FC+import qualified Data.Eigen.Internal as I+import qualified Data.Vector.Storable as VS+import qualified Data.Vector.Storable.Mutable as VSM+import qualified Data.ByteString as BS+import qualified Data.ByteString.Lazy as BSL+import qualified Data.ByteString.Internal as BSI++{-| A versatible sparse matrix representation.++This class implements a more versatile variants of the common compressed row/column storage format.+Each colmun's (resp. row) non zeros are stored as a pair of value with associated row (resp. colmiun) index.+All the non zeros are stored in a single large buffer.+Unlike the compressed format, there might be extra space inbetween the nonzeros of two successive colmuns+(resp. rows) such that insertion of new non-zero can be done with limited memory reallocation and copies.++A call to the function 'makeCompressed' turns the matrix into the standard compressed format compatible with many library.++Implementation deails of SparseMatrix are intentionally hidden behind ForeignPtr bacause Eigen doesn't provide mapping over plain data for sparse matricies.++For more infomration please see Eigen documentation page: <http://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html>+-}++data SparseMatrix a b where+    SparseMatrix :: I.Elem a b => !(ForeignPtr (I.CSparseMatrix a b)) -> SparseMatrix a b++-- | Alias for single precision sparse matrix+type SparseMatrixXf = SparseMatrix Float CFloat+-- | Alias for double precision sparse matrix+type SparseMatrixXd = SparseMatrix Double CDouble+-- | Alias for single previsiom sparse matrix of complex numbers+type SparseMatrixXcf = SparseMatrix (Complex Float) (I.CComplex CFloat)+-- | Alias for double prevision sparse matrix of complex numbers+type SparseMatrixXcd = SparseMatrix (Complex Double) (I.CComplex CDouble)++-- | Pretty prints the sparse matrix+instance (I.Elem a b, Show a) => Show (SparseMatrix a b) where+    show m = concat [+        "SparseMatrix ", show (rows m), "x", show (cols m),+        "\n", L.intercalate "\n" $ P.map (L.intercalate "\t" . P.map show) $ toDenseList m, "\n"]++-- | Shortcuts for basic matrix math+instance I.Elem a b => Num (SparseMatrix a b) where+    (*) = mul+    (+) = add+    (-) = sub+    fromInteger x = fromList 1 1 [(0,0,fromInteger x)]+    signum = map signum+    abs = map abs+    negate = map negate++-- | Not exposed, For internal use donly+map :: I.Elem a b => (a -> a) -> SparseMatrix a b -> SparseMatrix a b+map f m = fromList (rows m) (cols m) . P.map (\(r,c,v) -> (r,c,f v)) . toList $ m++mk :: I.Elem a b => Ptr (I.CSparseMatrix a b) -> IO (SparseMatrix a b)+mk p = SparseMatrix <$> FC.newForeignPtr p (I.call $ I.sparse_free p)++-- | Number of rows for the sparse matrix+rows :: I.Elem a b => SparseMatrix a b -> Int+rows = _unop I.sparse_rows (return . I.cast)++-- | Number of columns for the sparse matrix+cols :: I.Elem a b => SparseMatrix a b -> Int+cols = _unop I.sparse_cols (return . I.cast)++-- | Matrix coefficient at given row and col+coeff :: I.Elem a b => Int -> Int -> SparseMatrix a b -> a+coeff row col (SparseMatrix fp) = I.performIO $ withForeignPtr fp $ \p -> alloca $ \pq -> do+    I.call $ I.sparse_coeff p (I.cast row) (I.cast col) pq+    I.cast <$> peek pq++-- | Matrix coefficient at given row and col+(!) :: I.Elem a b => SparseMatrix a b -> (Int, Int) -> a+(!) m (row, col) = coeff row col m++{-| 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 :: I.Elem a b => SparseMatrix a b -> a+norm = _unop I.sparse_norm (return . I.cast)++-- | 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 :: I.Elem a b => SparseMatrix a b -> a+squaredNorm = _unop I.sparse_squaredNorm (return . I.cast)++-- | 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 :: I.Elem a b => SparseMatrix a b -> a+blueNorm = _unop I.sparse_blueNorm (return . I.cast)++-- | Extract rectangular block from sparse matrix defined by startRow startCol blockRows blockCols+block :: I.Elem a b => Int -> Int -> Int -> Int -> SparseMatrix a b -> SparseMatrix a b+block row col rows cols = _unop (\p pq -> I.sparse_block p (I.cast row) (I.cast col) (I.cast rows) (I.cast cols) pq) mk++-- | Number of non-zeros elements in the sparse matrix+nonZeros :: I.Elem a b => SparseMatrix a b -> Int+nonZeros = _unop I.sparse_nonZeros (return . I.cast)++{-++-- | not exposed currently+compress :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b+compress = _unop I.sparse_compress mk++-- | not exposed currently+uncompress :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b+uncompress = _unop I.sparse_uncompress mk++-- | not exposed currently+compressed :: I.Elem a b => SparseMatrix a b -> Bool+compressed = _unop I.sparse_isCompressed (return . (/=0))+++-- | not exposed currently+innerSize :: I.Elem a b => SparseMatrix a b -> Int+innerSize = _unop I.sparse_innerSize (return . I.cast)++-- | not exposed currently+outerSize :: I.Elem a b => SparseMatrix a b -> Int+outerSize = _unop I.sparse_outerSize (return . I.cast)++-}++-- | Suppresses all nonzeros which are much smaller than reference under the tolerence epsilon+pruned :: I.Elem a b => a -> SparseMatrix a b -> SparseMatrix a b+pruned r = _unop (\p pq -> alloca $ \pr -> poke pr (I.cast r) >> I.sparse_prunedRef p pr pq) mk++-- | Multiply matrix on a given scalar+scale :: I.Elem a b => a -> SparseMatrix a b -> SparseMatrix a b+scale x = _unop (\p pq -> alloca $ \px -> poke px (I.cast x) >> I.sparse_scale p px pq) mk++-- | Transpose of the sparse matrix+transpose :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b+transpose = _unop I.sparse_transpose mk++-- | Adjoint of the sparse matrix+adjoint :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b+adjoint = _unop I.sparse_adjoint mk++{-+-- | not exposed currently+lowerTriangle :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b+lowerTriangle = _unop I.sparse_lowerTriangle mk++-- | not exposed currently+upperTriangle :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b+upperTriangle = _unop I.sparse_upperTriangle mk+-}++-- | Adding two sparse matrices by adding the corresponding entries together. You can use @(+)@ function as well.+add :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b -> SparseMatrix a b+add = _binop I.sparse_add mk++-- | Subtracting two sparse matrices by subtracting the corresponding entries together. You can use @(-)@ function as well.+sub :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b -> SparseMatrix a b+sub = _binop I.sparse_sub mk++-- | Matrix multiplication. You can use @(*)@ function as well.+mul :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b -> SparseMatrix a b+mul = _binop I.sparse_mul mk++-- | Construct sparse matrix of given size from the list of triplets (row, col, val)+fromList :: I.Elem a b => Int -> Int -> [(Int, Int, a)] -> SparseMatrix a b+fromList rows cols tris = I.performIO $ VS.unsafeWith vs $ \p -> alloca $ \pq -> do+    I.call $ I.sparse_fromList (I.cast rows) (I.cast cols) p (I.cast $ VS.length vs) pq+    peek pq >>= mk+    where vs = VS.fromList $ P.map (\(row,col,val) -> I.CTriplet (I.cast row) (I.cast col) (I.cast val)) tris++-- | Convert sparse matrix to the list of triplets (row, col, val). Compressed elements will not be included+toList :: I.Elem a b => SparseMatrix a b -> [(Int, Int, a)]+toList m@(SparseMatrix fp) = I.performIO $ do+    let size = nonZeros m+    tris <- VSM.new size+    withForeignPtr fp $ \p ->+        VSM.unsafeWith tris $ \q ->+            I.call $ I.sparse_toList p q (I.cast size)+    let f (I.CTriplet row col val) = (I.cast row, I.cast col, I.cast val)+    P.map f . VS.toList <$> VS.unsafeFreeze tris++-- | Construct sparse matrix of two-dimensional list of values. Matrix dimensions will be detected automatically. Zero values will be compressed.+fromDenseList :: (I.Elem a b, Eq a) => [[a]] -> SparseMatrix a b+fromDenseList list = fromList rows cols $ do+    (row, vals) <- zip [0..] list+    (col, val) <- zip [0..] vals+    guard $ val /= 0+    return (row, col, val)+    where+        rows = length list+        cols = L.foldl' max 0 $ P.map length list++-- | Convert sparse matrix to (rows X cols) dense list of values+toDenseList :: I.Elem a b => SparseMatrix a b -> [[a]]+toDenseList m = [[coeff row col m | col <- [0 .. cols m - 1]] | row <- [0 .. rows m - 1]]++-- | Construct sparse matrix from dense matrix. Zero elements will be compressed+fromMatrix :: I.Elem a b => M.Matrix a b -> SparseMatrix a b+fromMatrix m1 = I.performIO $ alloca $ \pm0 ->+    M.unsafeWith m1 $ \vals rows cols -> do+        I.call $ I.sparse_fromMatrix vals rows cols pm0+        peek pm0 >>= mk++-- | Construct dense matrix from sparse matrix+toMatrix :: I.Elem a b => SparseMatrix a b -> M.Matrix a b+toMatrix m1@(SparseMatrix fp) = I.performIO $ do+    m0 <- MM.new (rows m1) (cols m1)+    MM.unsafeWith m0 $ \vals rows cols ->+        withForeignPtr fp $ \pm1 ->+            I.call $ I.sparse_toMatrix pm1 vals rows cols+    M.unsafeFreeze m0++-- | Encode the sparse matrix as a lazy byte string+encode :: I.Elem a b => SparseMatrix a b -> BSL.ByteString+encode m@(SparseMatrix fp) = I.performIO $ do+    let size = nonZeros m+    tris <- VSM.new size+    withForeignPtr fp $ \p ->+        VSM.unsafeWith tris $ \q ->+            I.call $ I.sparse_toList p q (I.cast size)+    tris <- VS.unsafeFreeze tris+    let +        tri@(I.CTriplet _ _ val) = VS.head tris++    return $ BSL.fromChunks [+        encodeInt (I.magicCode val),+        encodeInt (I.cast $ rows m),+        encodeInt (I.cast $ cols m),+        encodeInt (I.cast $ size),+        let (fp, fs) = VS.unsafeToForeignPtr0 tris in BSI.PS (castForeignPtr fp) 0 (fs * sizeOf tri)]++    where+        encodeInt :: CInt -> BS.ByteString+        encodeInt x = BSI.unsafeCreate (sizeOf x) $ (`poke` x) . castPtr+        ++-- | Decode sparse matrix from the lazy byte string+decode :: forall a b . I.Elem a b => BSL.ByteString -> SparseMatrix a b+decode st = I.performIO $ do+    let+        val = undefined :: b+        tri = undefined :: I.CTriplet b+        triSize = sizeOf tri++    st <- I.openStream st++    code <- I.decodeInt <$> I.readStream st I.intSize+    when (code /= I.magicCode val) $+        fail "SparseMatrix.decode: wrong matrix type"+    +    rows <- I.readInt st+    cols <- I.readInt st+    size <- I.readInt st+    +    BSI.PS fp fo _ <- I.readStream st (I.cast size * triSize)++    I.closeStream st++    let tris = VS.unsafeFromForeignPtr0 (I.plusForeignPtr fp fo) (I.cast size)+    +    VS.unsafeWith tris $ \p -> alloca $ \pq -> do+        I.call $ I.sparse_fromList rows cols p size pq+        peek pq >>= mk++_unop :: Storable c => (I.CSparseMatrixPtr a b -> Ptr c -> IO CString) -> (c -> IO d) -> SparseMatrix a b -> d+_unop f g (SparseMatrix fp) = I.performIO $+    withForeignPtr fp $ \p ->+        alloca $ \pq -> do+            I.call (f p pq)+            peek pq >>= g++_binop :: Storable c => (I.CSparseMatrixPtr a b -> I.CSparseMatrixPtr a b -> Ptr c -> IO CString) -> (c -> IO d) -> SparseMatrix a b -> SparseMatrix a b -> d+_binop f g (SparseMatrix fp1) (SparseMatrix fp2) = I.performIO $+    withForeignPtr fp1 $ \p1 ->+        withForeignPtr fp2 $ \p2 ->+            alloca $ \pq -> do+                I.call (f p1 p2 pq)+                peek pq >>= g
cbits/eigen-proxy.cpp view
@@ -35,29 +35,41 @@     return Map< Matrix<T,Dynamic,Dynamic> >((const T*)p, r, c); } -#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)\-{\+#define RET const char* ++#define API(name,args,call) \+extern "C" RET eigen_##name args {\     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;\+        case 0: return name<T0>call;\+        case 1: return name<T1>call;\+        case 2: return name<T2>call;\+        case 3: return name<T3>call;\     }\     GUARD_END\ } ++#define BINOP(name,op) \+template <class T>\+RET name(void* p, int r, int c,\+    const void* p1, int r1, int c1,\+    const void* p2, int r2, int c2)\+{\+    matrix<T>(p,r,c) = matrix<T>(p1,r1,c1) op matrix<T>(p2,r2,c2);\+    return 0;\+}\+API(name, (int code,\+    void* p, int r, int c,\+    const void* p1, int r1, int c1,\+    const void* p2, int r2, int c2), (p,r,c,p1,r1,c1,p2,r2,c2));+ BINOP(add,+); BINOP(sub,-); BINOP(mul,*);  #define PROP(name) \-extern "C" const char* eigen_##name(int code, void* q, const void* p, int r, int c) {\+extern "C" RET 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;\@@ -79,7 +91,7 @@ PROP(determinant);  #define UNOP(name) \-extern "C" const char* eigen_##name(int code, void* p, int r, int c, const void* p1, int r1, int c1) {\+extern "C" RET 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;\@@ -96,7 +108,7 @@ UNOP(diagonal); UNOP(transpose); -extern "C" const char* eigen_normalize(int code, void* p, int r, int c)+extern "C" RET eigen_normalize(int code, void* p, int r, int c) {     GUARD_START     switch (code) {@@ -108,7 +120,7 @@     GUARD_END } -extern "C" const char* eigen_random(int code, void* p, int r, int c)+extern "C" RET eigen_random(int code, void* p, int r, int c) {     GUARD_START     switch (code) {@@ -120,7 +132,7 @@     GUARD_END } -extern "C" const char* eigen_identity(int code, void* p, int r, int c)+extern "C" RET eigen_identity(int code, void* p, int r, int c) {     GUARD_START     switch (code) {@@ -133,7 +145,7 @@ }  template <class T>-const char* rank(Decomposition d, int* v, const void* p, int r, int c) {+RET 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) {@@ -154,20 +166,10 @@     }     return 0; }--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-}+API(rank, (int code, Decomposition d, int* v, const void* p, int r, int c), (d,v,p,r,c));  template <class T>-const char* kernel(Decomposition d, void** p0, int* r0, int* c0, const void* p1, int r1, int c1) {+RET 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.");@@ -179,20 +181,10 @@     MapMatrix((T*)*p0, *r0, *c0) = B;     return 0; }--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-}+API(kernel, (int code, Decomposition d, void** p0, int* r0, int* c0, const void* p1, int r1, int c1), (d,p0,r0,c0,p1,r1,c1));  template <class T>-const char* image(Decomposition d, void** p0, int* r0, int* c0, const void* p1, int r1, int c1) {+RET 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.");@@ -204,20 +196,10 @@     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-}+API(image, (int code, Decomposition d, void** p0, int* r0, int* c0, const void* p1, int r1, int c1), (d,p0,r0,c0,p1,r1,c1));  template <class T>-const char* solve(Decomposition d,+RET solve(Decomposition d,     void* px, int rx, int cx,     const void* pa, int ra, int ca,     const void* pb, int rb, int cb)@@ -254,24 +236,13 @@     }     return 0; }--extern "C" const char* eigen_solve(int code, Decomposition d,+API(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-}+    const void* pb, int rb, int cb), (d,px,rx,cx,pa,ra,ca,pb,rb,cb));  template <class T>-const char* relativeError(void* e,+RET relativeError(void* e,     const void* px, int rx, int cx,     const void* pa, int ra, int ca,     const void* pb, int rb, int cb)@@ -283,23 +254,10 @@     *(T*)e = (A*x - b).norm() / b.norm();     return 0; }---extern "C" const char* eigen_relativeError(int code, void* e,+API(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-    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-}-+    const void* pb, int rb, int cb), (e,px,rx,cx,pa,ra,ca,pb,rb,cb));  extern "C" bool eigen_initParallel() {     initParallel();@@ -313,3 +271,167 @@ extern "C" int eigen_getNbThreads() {     return nbThreads(); }++template <class T>+RET sparse_fromList(int rows, int cols, void* data, int size, void** pr) {+    typedef SparseMatrix<T> M;+    typedef Triplet<T> E;+    std::auto_ptr<M> a(new M(rows, cols));+    a->setFromTriplets((E*)data, (E*)data + size);+    *(M**)pr = a.release();+    return 0;+}+API(sparse_fromList, (int code, int rows, int cols, void* data, int size, void** pr), (rows,cols,data,size,pr));++template <class  T>+RET sparse_toList(void* p, void* q, int s) {+    int n = 0;+    typedef SparseMatrix<T> M;+    typedef Triplet<T> E;+    M* m = (M*)p;+    for (int k = 0; k < m->outerSize(); ++k) {+        for (typename M::InnerIterator i(*m, k); i; ++i) {+            if (n >= s)+                return strdup("sparse_toList: buffer overrun detected");+            ((E*)q)[n++] = E(i.row(), i.col(), i.value());+        }+    }+    return n < s ? strdup("sparse_toList: buffer underrun detected") : 0;+}+API(sparse_toList, (int code, void* p, void* q, int s), (p,q,s));++template <class T>+RET sparse_free(void* p) {+    delete (SparseMatrix<T>*)p;+    return 0;+}+API(sparse_free, (int code, void* p), (p));++#define SPARSE_UNOP_INPLACE(name)\+template <class T>\+RET sparse_##name(void* p, void** pr) {\+    typedef SparseMatrix<T> M;\+    std::auto_ptr<M> a(new M(*(M*)p));\+    a->name();\+    *(M**)pr = a.release();\+    return 0;\+}\+API(sparse_##name, (int code, void* p, void** pr), (p, pr));++SPARSE_UNOP_INPLACE(makeCompressed);+SPARSE_UNOP_INPLACE(uncompress);++#define SPARSE_UNOP(name)\+template <class T>\+RET sparse_##name(void* p, void** pr) {\+    typedef SparseMatrix<T> M;\+    *(M**)pr = new M(((M*)p)->name());\+    return 0;\+}\+API(sparse_##name, (int code, void* p, void** pr), (p, pr));++SPARSE_UNOP(adjoint);+SPARSE_UNOP(transpose);++template <class T>+RET sparse_pruned(void* p, void** pr) {+    typedef SparseMatrix<T> M;+    std::auto_ptr<M> a(new M(*(M*)p));+    a->prune(T(1));+    *(M**)pr = a.release();+    return 0;+}+API(sparse_pruned, (int code, void* p, void** pr), (p, pr));++template <class T>+RET sparse_prunedRef(void* p, void* q, void** pr) {+    typedef SparseMatrix<T> M;+    std::auto_ptr<M> a(new M(*(M*)p));+    a->prune(*(T*)q);+    *(M**)pr = a.release();+    return 0;+}+API(sparse_prunedRef, (int code, void* p, void* q, void** pr), (p, q, pr));++template <class T>+RET sparse_scale(void* p, void* q, void** pr) {+    typedef SparseMatrix<T> M;+    *(M**)pr = new M(*(T*)q * *(M*)p);+    return 0;+}+API(sparse_scale, (int code, void* p, void* q, void** pr), (p, q, pr));++template <class T>+RET sparse_coeff(void* p, int row, int col, void* pr) {+    *(T*)pr = ((SparseMatrix<T>*)p)->coeff(row, col);+    return 0;+}+API(sparse_coeff, (int code, void* p, int row, int col, void* pr), (p, row, col, pr));++#define SPARSE_PROP(name,type)\+template <class T>\+RET sparse_##name(void* p, void* pr) {\+    *(type*)pr = ((SparseMatrix<T>*)p)->name();\+    return 0;\+}\+API(sparse_##name, (int code, void* p, void* pr), (p, pr));++SPARSE_PROP(cols, int);+SPARSE_PROP(rows, int);+SPARSE_PROP(innerSize, int);+SPARSE_PROP(outerSize, int);+SPARSE_PROP(nonZeros, int);+SPARSE_PROP(isCompressed, int);+SPARSE_PROP(norm, T);+SPARSE_PROP(squaredNorm, T);+SPARSE_PROP(blueNorm, T);++#define SPARSE_BINOP(name,op)\+template <class T>\+RET sparse_##name(void* p, void* q, void** pr) {\+    typedef SparseMatrix<T> M;\+    *(M**)pr = new M(*(M*)p op *(M*)q);\+    return 0;\+}\+API(sparse_##name, (int code, void* p, void* q, void** pr), (p, q, pr));++SPARSE_BINOP(add, +);+SPARSE_BINOP(sub, -);+SPARSE_BINOP(mul, *);++template <class T>+RET sparse_block(void* p, int row, int col, int rows, int cols, void** pr) {+    typedef SparseMatrix<T> M;+    *(M**)pr = new M(((M*)p)->block(row,col,rows,cols));+    return 0;+}+API(sparse_block, (int code, void* p, int row, int col, int rows, int cols, void** pr), (p, row, col, rows, cols, pr));++template <class T> bool isZero(T x) { return x == 0; }+template <class T> bool isZero(std::complex<T> x) { return x.real() == 0 && x.imag() == 0; }++template <class T>+RET sparse_fromMatrix(void* p, int rows, int cols, void** pq) {+    typedef SparseMatrix<T> M;+    typedef Map< Matrix<T,Dynamic,Dynamic> > MapMatrix;+    MapMatrix src((const T*)p, rows, cols);+    std::auto_ptr<M> dst(new M(rows, cols));+    for (int row = 0; row < rows; ++row) {+        for (int col = 0; col < cols; ++col) {+            T val = src.coeff(row,col);+            if (!isZero(val))+                dst->insert(row, col) = val;+        }+    }+    *(M**)pq = dst.release();+    return 0;+}+API(sparse_fromMatrix, (int code, void* p, int rows, int cols, void** pq), (p,rows,cols,pq));++template <class T>+RET sparse_toMatrix(void* p, void* q, int rows, int cols) {+    typedef Map< Matrix<T,Dynamic,Dynamic> > MapMatrix;+    MapMatrix((T*)q, rows, cols) = *(SparseMatrix<T>*)p;+    return 0;+}+API(sparse_toMatrix, (int code, void* p, void* q, int rows, int cols), (p,q,rows,cols));
cbits/eigen-proxy.h view
@@ -7,6 +7,7 @@     if (!(x)) eigen_assert_fail(#x, __PRETTY_FUNCTION__, __FILE__, __LINE__);\ } while(false) #include <Eigen/Core>+#include <Eigen/SparseCore>  using namespace Eigen; @@ -35,6 +36,42 @@ bool eigen_initParallel(); int eigen_getNbThreads(); void eigen_setNbThreads(int);++#define SPARSE_PROP(name) const char* eigen_sparse_##name(int, void*, void*)+SPARSE_PROP(cols);+SPARSE_PROP(rows);+SPARSE_PROP(innerSize);+SPARSE_PROP(outerSize);+SPARSE_PROP(nonZeros);+SPARSE_PROP(isCompressed);+SPARSE_PROP(norm);+SPARSE_PROP(squaredNorm);+SPARSE_PROP(blueNorm);+#undef SPARSE_PROP++#define SPARSE_BINOP(name) const char* eigen_sparse_##name(int, void*, void*, void**);+SPARSE_BINOP(add);+SPARSE_BINOP(sub);+SPARSE_BINOP(mul);+SPARSE_BINOP(prunedRef);+SPARSE_BINOP(scale);+#undef SPARSE_BINOP++#define SPARSE_UNOP(name) const char* eigen_sparse_##name(int, void*, void**);+SPARSE_UNOP(makeCompressed);+SPARSE_UNOP(uncompress);+SPARSE_UNOP(adjont);+SPARSE_UNOP(transponse);+SPARSE_UNOP(pruned);+#undef SPARSE_UNOP++const char* eigen_sparse_fromList(int, int, int, void*, int, void**);+const char* eigen_sparse_toList(int, void*, void*, int);+const char* eigen_sparse_free(int, void*);+const char* eigen_sparse_coeff(int, void*, int, int, void*);+const char* eigen_sparse_block(int, void*, int, int, int, int, void**);+const char* eigen_sparse_fromMatrix(int, void*, int, int, void**);+const char* eigen_sparse_toMatrix(int, void*, void*, int, int);  enum Decomposition { 	PartialPivLU,
eigen.cabal view
@@ -1,5 +1,5 @@ name:           eigen-version:        2.0.1+version:        2.1.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.@@ -364,6 +364,7 @@     exposed-modules:    Data.Eigen.LA                         Data.Eigen.Matrix                         Data.Eigen.Matrix.Mutable+                        Data.Eigen.SparseMatrix                         Data.Eigen.Parallel      other-modules:      Data.Eigen.Internal@@ -371,7 +372,8 @@     ghc-options:        -Wall -fno-warn-name-shadowing     build-depends:      base >= 3 && < 5,                         vector >= 0.5 && < 0.11,-                        primitive >= 0.1 && < 0.6+                        primitive >= 0.1 && < 0.6,+                        bytestring      include-dirs:       eigen3, cbits     c-sources:          cbits/eigen-proxy.cpp