packages feed

matrix-sized 0.0.2 → 0.0.3

raw patch · 13 files changed

+357/−94 lines, 13 filesdep +matrix-sizeddep +tastydep +tasty-quickcheckdep −binarydep −bytestringdep −deepseqdep ~base

Dependencies added: matrix-sized, tasty, tasty-quickcheck

Dependencies removed: binary, bytestring, deepseq

Dependency ranges changed: base

Files

README.md view
@@ -18,18 +18,39 @@ -------  ```haskell-let mat = D.matrix [ [1,0,3]-                   , [0,5,6]-                   , [0,0,0] ] :: Matrix 3 3 Double-    mat' = D.convertAny mat :: SparseMatrix 3 3 Double+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE TypeOperators #-} -print mat-print mat'+import Data.Matrix.Static.LinearAlgebra+import qualified Data.Matrix.Static.Generic as G+import qualified Data.Matrix.Static.Dense as D+import qualified Data.Matrix.Static.Sparse as S+import Data.Singletons.Prelude hiding ((@@))+import Data.Singletons.TypeLits+import Data.Complex+import System.Random+import Control.Monad+import Data.Type.Equality -print $ eigs (sing :: Sing 1) mat == eigs (sing :: Sing 1) mat'+f :: (SingI n, (2 <= n - 2) ~ 'True)+  => Matrix n n Double -> Matrix n n (Complex Double)+f m = let (d, v) = eigs (sing :: Sing 2) m+      in v @@ S.diag d @@ G.transpose v -print $ cholesky mat+main :: IO ()+main = do+    n <- randomRIO (2, 6)+    vals <- replicateM (n*n) $ randomRIO (-100,100) :: IO [Double] -print $ mat %*% mat %*% mat-print $ mat' %*% mat' %*% mat+    withSomeSing (fromIntegral n) $ \sn@(SNat :: Sing n) ->+        let s0 = SNat :: Sing 2+            sn2 = sn %- s0+        in case testEquality (s0 %<= sn2) STrue of+            Just Refl -> do+                let mat = G.fromList vals :: Matrix n n Double+                print $ f mat+            Nothing -> error $ "Requiring Matrix size >= 4, but got: " <> show n ```
cbits/eigen-basic.cpp view
@@ -70,8 +70,9 @@     const void* val, const void* outer, const void* inner, int r2, int c2, int s,     const void* p1, int r1, int c1), (p,r,c,val,outer,inner,r2,c2,s,p1,r1,c1)); + template <class T>-RET ss_mul( void** v, void** o, void** i, int r, int c, int* s,+RET ss_mul( void** v, void* o, void** i, int r, int c, int* s,     const void* v1, const void* o1, const void* i1, int r1, int c1, int s1,     const void* v2, const void* o2, const void* i2, int r2, int c2, int s2) {@@ -79,24 +80,96 @@     MapSparseMatrix a(r1, c1, s1, (int*)o1, (int*)i1, (T*)v1);     MapSparseMatrix b(r2, c2, s2, (int*)o2, (int*)i2, (T*)v2);     SparseMatrix<T> M = (a * b).pruned();++    memcpy(o, M.outerIndexPtr(), (c+1) * sizeof(int));+     *s = M.nonZeros();-    T* p1 = (T*) malloc(*s * sizeof(T));-    memcpy(p1, M.valuePtr(), *s * sizeof(T));+    T* p1 = (T*) malloc((*s) * sizeof(T));+    memcpy(p1, M.valuePtr(), (*s) * sizeof(T));     *v = p1;-    int* p2 = (int*) malloc(*s * sizeof(int));-    memcpy(p2, M.innerIndexPtr(), *s * sizeof(T));+    int* p2 = (int*) malloc((*s) * sizeof(int));+    memcpy(p2, M.innerIndexPtr(), (*s) * sizeof(int));     *i = p2;-    int* p3 = (int*) malloc(c * sizeof(int));-    memcpy(p3, M.outerIndexPtr(), c * sizeof(int));-    *o = p3;     return 0; } API(ss_mul, (int code,-    void** v, void** o, void** i, int r, int c, int* s,+    void** v, void* o, void** i, int r, int c, int* s,     const void* v1, const void* o1, const void* i1, int r1, int c1, int s1,     const void* v2, const void* o2, const void* i2, int r2, int c2, int s2),     (v,o,i,r,c,s,v1,o1,i1,r1,c1,s1,v2,o2,i2,r2,c2,s2)); ++template <class T>+RET ss_plus( void** v, void* o, void** i, int r, int c, int* s,+    const void* v1, const void* o1, const void* i1, int r1, int c1, int s1,+    const void* v2, const void* o2, const void* i2, int r2, int c2, int s2)+{+    typedef Map<SparseMatrix<T> > MapSparseMatrix;+    MapSparseMatrix a(r1, c1, s1, (int*)o1, (int*)i1, (T*)v1);+    MapSparseMatrix b(r2, c2, s2, (int*)o2, (int*)i2, (T*)v2);+    SparseMatrix<T> M = a + b;++    memcpy(o, M.outerIndexPtr(), (c+1) * sizeof(int));++    *s = M.nonZeros();+    T* p1 = (T*) malloc((*s) * sizeof(T));+    memcpy(p1, M.valuePtr(), (*s) * sizeof(T));+    *v = p1;+    int* p2 = (int*) malloc((*s) * sizeof(int));+    memcpy(p2, M.innerIndexPtr(), (*s) * sizeof(int));+    *i = p2;+    return 0;+}+API(ss_plus, (int code,+    void** v, void* o, void** i, int r, int c, int* s,+    const void* v1, const void* o1, const void* i1, int r1, int c1, int s1,+    const void* v2, const void* o2, const void* i2, int r2, int c2, int s2),+    (v,o,i,r,c,s,v1,o1,i1,r1,c1,s1,v2,o2,i2,r2,c2,s2));++template <class T>+RET ss_cmul( void** v, void* o, void** i, int r, int c, int* s,+    const void* v1, const void* o1, const void* i1, int r1, int c1, int s1,+    const void* v2, const void* o2, const void* i2, int r2, int c2, int s2)+{+    typedef Map<SparseMatrix<T> > MapSparseMatrix;+    MapSparseMatrix a(r1, c1, s1, (int*)o1, (int*)i1, (T*)v1);+    MapSparseMatrix b(r2, c2, s2, (int*)o2, (int*)i2, (T*)v2);++    SparseMatrix<T> M = a.cwiseProduct(b);++    memcpy(o, M.outerIndexPtr(), (c+1) * sizeof(int));++    *s = M.nonZeros();+    T* p1 = (T*) malloc((*s) * sizeof(T));+    memcpy(p1, M.valuePtr(), (*s) * sizeof(T));+    *v = p1;+    int* p2 = (int*) malloc((*s) * sizeof(int));+    memcpy(p2, M.innerIndexPtr(), (*s) * sizeof(int));+    *i = p2;+    return 0;+}+API(ss_cmul, (int code,+    void** v, void* o, void** i, int r, int c, int* s,+    const void* v1, const void* o1, const void* i1, int r1, int c1, int s1,+    const void* v2, const void* o2, const void* i2, int r2, int c2, int s2),+    (v,o,i,r,c,s,v1,o1,i1,r1,c1,s1,v2,o2,i2,r2,c2,s2));++template <class T>+RET sd_plus( void* p, int r, int c,+    const void* val, const void* outer, const void* inner,+    int r2, int c2, int s,+    const void* p1, int r1, int c1)+{+    typedef Map< Matrix<T,Dynamic,Dynamic> > MapMatrix;+    MapMatrix x((T*)p, r, c);+    x = matrix<T>(p1,r1,c1);+    x += smatrix<T>(val, outer, inner, r2, c2, s);+    return 0;+}+API(sd_plus, (int code,+    void* p, int r, int c,+    const void* val, const void* outer, const void* inner, int r2, int c2, int s,+    const void* p1, int r1, int c1), (p,r,c,val,outer,inner,r2,c2,s,p1,r1,c1));  #define UNOP(name) \ extern "C" RET __attribute__((noinline)) eigen_##name(int code, void* p, int r, int c, const void* p1, int r1, int c1) {\
cbits/eigen-solver.cpp view
@@ -18,7 +18,9 @@     MapComplexMatrix V((T3*)v, n, k);      DenseGenMatProd<double> op(M);-    GenEigsSolver< double, LARGEST_MAGN, DenseGenMatProd<double> > eigs(&op, k, 2*k+1);+    int ncv = 2 * k + 1;+    ncv = (ncv <= n) ? ncv : n;+    GenEigsSolver< double, LARGEST_MAGN, DenseGenMatProd<double> > eigs(&op, k, ncv);     eigs.init();     int nconv = eigs.compute();     if(eigs.info() == 0)@@ -42,7 +44,9 @@     MapComplexMatrix V((T3*)v, n, k);      SparseGenMatProd<double> op(M);-    GenEigsSolver< double, LARGEST_MAGN, SparseGenMatProd<double> > eigs(&op, k, 2*k+1);+    int ncv = 2 * k + 1;+    ncv = (ncv <= n) ? ncv : n;+    GenEigsSolver< double, LARGEST_MAGN, SparseGenMatProd<double> > eigs(&op, k, ncv);     eigs.init();     int nconv = eigs.compute();     if(eigs.info() == 0)@@ -62,3 +66,18 @@ } API(cholesky, (int code,     void* px, const void* pa, int n), (px,pa,n));+++/*+template <class T>+RET bdcsvd(void* px, int r, int c+{+    typedef Map< Matrix<T,Dynamic,Dynamic> > MapMatrix;+    MapMatrix x((T*)px, n, n);+    MapMatrix A((T*)pa, n, n);+    x = A.llt().matrixL();+    return 0;+}+API(cholesky, (int code,+    void* px, const void* pa, int n), (px,pa,n));+*/
matrix-sized.cabal view
@@ -4,10 +4,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: 6910822a984416806867bac316c8f9685ade7b6a7672d4437c67fd4770977af1+-- hash: ef7b97a49f99defe3cd58f2c5a3411fe8c958f0ad40f4b3f0bd763720b643db8  name:           matrix-sized-version:        0.0.2+version:        0.0.3 synopsis:       Haskell matrix library with interface to C++ linear algebra libraries. description:    A Haskell implementation of matrices with statically known sizes. The library also comes with the bindings to high performance C++ linear algebra libraries such as Eigen and Spectra. category:       Math@@ -375,6 +375,11 @@   type: git   location: https://github.com/kaizhang/matrix-sized +flag parallel+  description: Enable multithreading+  manual: True+  default: False+ library   exposed-modules:       Data.Matrix.Static.Dense@@ -403,10 +408,27 @@       stdc++   build-depends:       base >=4.10 && <5-    , binary-    , bytestring-    , deepseq     , primitive >=0.6.4.0     , singletons+    , vector >=0.11+  if flag(parallel)+    cc-options: -fopenmp+    cxx-options: -fopenmp -std=c++11+    ld-options: -fopenmp+  default-language: Haskell2010++test-suite test+  type: exitcode-stdio-1.0+  main-is: Main.hs+  hs-source-dirs:+      tests+  ghc-options: -Wall -threaded+  build-depends:+      base+    , matrix-sized+    , primitive >=0.6.4.0+    , singletons+    , tasty+    , tasty-quickcheck     , vector >=0.11   default-language: Haskell2010
src/Data/Matrix/Static/Dense.hs view
@@ -41,8 +41,7 @@     , C.fromRows     , C.fromColumns     , C.unsafeFromVector--    , diag+    , replicate     , diagRect      -- * Conversions@@ -91,10 +90,9 @@     , C.create     ) where -import           Control.DeepSeq                   hiding (force) import           Control.Monad                     (liftM) import qualified Data.Vector.Generic               as G-import Prelude hiding (mapM, mapM_, zipWith, map, sequence, sequence_, zip, unzip, zipWith3, zip3, unzip3)+import Prelude hiding (replicate, mapM, mapM_, zipWith, map, sequence, sequence_, zip, unzip, zipWith3, zip3, unzip3) import GHC.TypeLits (type (<=)) import Data.Singletons import Data.Tuple (swap)@@ -136,9 +134,6 @@         recip = C.map recip         fromRational = undefined -instance NFData (v a) => NFData (Matrix r c v a) where-    rnf (Matrix vec) = rnf vec- instance G.Vector v a => C.Matrix Matrix v a where     -- | O(1) Return the size of matrix.     dim :: forall r c. Matrix r c v a -> (Int, Int)@@ -199,20 +194,14 @@     sequence_ (Matrix vec) = G.sequence_ vec     {-# INLINE sequence_ #-} -{---- | O(m*n) Create an identity matrix-ident :: (Num a, G.Vector v a)-      => Matrix n n v a-ident = diagRect 0 $ replicate 1-{-# INLINE ident #-}--}---- | O(m*n) Create a square matrix with given diagonal, other entries default to 0-diag :: (Num a, G.Vector v a, SingI n)-     => Matrix n 1 v a       -- ^ diagonal-     -> Matrix n n v a-diag = diagRect 0-{-# INLINE diag #-}+-- | O(m*n) Create a constant matrix.+replicate :: forall r c v a. (G.Vector v a, SingI r, SingI c)+          => a -> Matrix r c v a+replicate = C.unsafeFromVector . G.replicate (r*c)+  where+    r = fromIntegral $ fromSing (sing :: Sing r)+    c = fromIntegral $ fromSing (sing :: Sing c)+{-# INLINE replicate #-}  -- | O(m*n) Create a rectangular matrix with default values and given diagonal diagRect :: (G.Vector v a, SingI r, SingI c, n <= r, n <= c)
src/Data/Matrix/Static/Dense/Mutable.hs view
@@ -18,7 +18,6 @@    , C.replicate    ) where -import           Control.DeepSeq import qualified Data.Vector.Generic.Mutable as GM import           Prelude                     hiding (read, replicate) import Data.Singletons@@ -29,9 +28,6 @@ -- | Column-major mutable matrix. data MMatrix :: C.MMatrixKind where     MMatrix :: (SingI r, SingI c) => v s a -> MMatrix r c v s a--instance (NFData (v s a)) => NFData (MMatrix r c v s a) where-    rnf (MMatrix vec) = rnf vec  instance GM.MVector v a => C.MMatrix MMatrix v a where     dim :: forall r c s. MMatrix r c v s a -> (Int, Int)
src/Data/Matrix/Static/Generic.hs view
@@ -41,6 +41,7 @@ import Text.Printf (printf) import Prelude hiding (map, mapM, mapM_, sequence, sequence_) import qualified Data.List as L+import Data.Tuple (swap) import Data.Kind (Type) import GHC.TypeLits (Nat, type (<=)) import Data.Singletons (SingI, Sing, fromSing, sing, withSomeSing)@@ -62,7 +63,7 @@     -- | Convert matrix to vector in column order.     -- Default algorithm is O((m*n) * O(unsafeIndex)).     flatten :: mat r c v a -> v a-    flatten mat = G.generate (r*c) $ \i -> unsafeIndex mat (i `divMod` r)+    flatten mat = G.generate (r*c) $ \i -> unsafeIndex mat (swap $ i `divMod` r)       where         (r,c) = dim mat     {-# INLINE flatten #-}@@ -90,7 +91,7 @@      transpose :: (SingI r, SingI c) => mat r c v a -> mat c r v a     transpose mat = unsafeFromVector $ G.generate (r*c) $ \x ->-        unsafeIndex mat $ x `divMod` r+        unsafeIndex mat $ x `divMod` c       where        (r, c) = dim mat     {-# INLINE transpose #-}
src/Data/Matrix/Static/Internal.hs view
@@ -4,6 +4,9 @@     , c_ds_mul     , c_sd_mul     , c_ss_mul+    , c_ss_cmul+    , c_sd_plus+    , c_ss_plus     , c_inverse     , c_cholesky     , c_eigs@@ -39,13 +42,33 @@           -> Ptr a -> CInt -> CInt           -> IO CString +foreign import ccall "eigen_sd_plus"+    c_sd_plus :: CInt+          -> Ptr a -> CInt -> CInt+          -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt+          -> Ptr a -> CInt -> CInt+          -> IO CString+ foreign import ccall "eigen_ss_mul"     c_ss_mul :: CInt-          -> Ptr (Ptr a) -> Ptr (Ptr CInt) -> Ptr (Ptr CInt) -> CInt -> CInt -> Ptr CInt+          -> Ptr (Ptr a) -> Ptr CInt -> Ptr (Ptr CInt) -> CInt -> CInt -> Ptr CInt           -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt           -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt           -> IO CString +foreign import ccall "eigen_ss_cmul"+    c_ss_cmul :: CInt+          -> Ptr (Ptr a) -> Ptr CInt -> Ptr (Ptr CInt) -> CInt -> CInt -> Ptr CInt+          -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt+          -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt+          -> IO CString++foreign import ccall "eigen_ss_plus"+    c_ss_plus :: CInt+          -> Ptr (Ptr a) -> Ptr CInt -> Ptr (Ptr CInt) -> CInt -> CInt -> Ptr CInt+          -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt+          -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt+          -> IO CString               foreign import ccall "eigen_inverse"     c_inverse :: CInt
src/Data/Matrix/Static/LinearAlgebra.hs view
@@ -11,15 +11,15 @@ module Data.Matrix.Static.LinearAlgebra     ( Arithmetic(..)     , Factorization(..)-    , inverse     , module Data.Matrix.Static.LinearAlgebra.Types     ) where  import qualified Data.Vector.Storable as VS import System.IO.Unsafe (unsafePerformIO) import Data.Complex (Complex)-import Data.Singletons-import GHC.TypeLits (type (<=), type (-))+import Data.Singletons.Prelude hiding ((@@), type (==))+import Data.Type.Bool (If)+import Data.Type.Equality (type (==))  import qualified Data.Matrix.Static.Dense as D import qualified Data.Matrix.Static.Sparse as S@@ -28,35 +28,68 @@ import qualified Data.Matrix.Static.Internal as Internal import Data.Matrix.Static.LinearAlgebra.Types -class Arithmetic (mat1 :: C.MatrixKind)-                 (mat2 :: C.MatrixKind)-                 (mat3 :: C.MatrixKind) |-                 mat1 mat2 -> mat3 where-    (%*%) :: (Numeric a, SingI n, SingI m)-          => mat1 n p VS.Vector a-          -> mat2 p m VS.Vector a+class Arithmetic (mat1 :: C.MatrixKind) (mat2 :: C.MatrixKind) where+    -- | Matrix multiplication+    (@@) :: ( Numeric a, SingI n, SingI m+            , If (mat1 == mat2) mat1 D.Matrix ~ mat3 )+         => mat1 n p VS.Vector a+         -> mat2 p m VS.Vector a+         -> mat3 n m VS.Vector a+    infixr 8 @@++    (%+%) :: ( Numeric a, SingI n, SingI m+             , If (mat1 == mat2) mat1 D.Matrix ~ mat3 )+          => mat1 n m VS.Vector a+          -> mat2 n m VS.Vector a           -> mat3 n m VS.Vector a+    infixr 8 %+%++    (%-%) :: ( Numeric a, SingI n, SingI m+             , If (mat1 == mat2) mat1 D.Matrix ~ mat3 )+          => mat1 n m VS.Vector a+          -> mat2 n m VS.Vector a+          -> mat3 n m VS.Vector a+    infixr 8 %-%++    (%*%) :: ( Numeric a, SingI n, SingI m+             , If (mat1 == mat2) mat1 S.SparseMatrix ~ mat3 )+          => mat1 n m VS.Vector a+          -> mat2 n m VS.Vector a+          -> mat3 n m VS.Vector a     infixr 8 %*% -instance Arithmetic D.Matrix D.Matrix D.Matrix where-    (%*%) = withFun2 Internal.c_dd_mul+instance Arithmetic D.Matrix D.Matrix where+    (@@) = withFun2 Internal.c_dd_mul+    (%+%) = (+)+    (%-%) = (-)+    (%*%) = (*) -instance Arithmetic D.Matrix S.SparseMatrix D.Matrix where-    (%*%) = withDS Internal.c_ds_mul+instance Arithmetic D.Matrix S.SparseMatrix where+    (@@) = withDS Internal.c_ds_mul+    (%+%) = flip (%+%)+    (%-%) a b = a %+% C.map negate b+    (%*%) = undefined -instance Arithmetic S.SparseMatrix D.Matrix D.Matrix where-    (%*%) = withSD Internal.c_sd_mul+instance Arithmetic S.SparseMatrix D.Matrix where+    (@@) = withSD Internal.c_sd_mul+    (%+%) = withSD Internal.c_sd_plus+    (%-%) a b = a %+% C.map negate b+    (%*%) = undefined -instance Arithmetic S.SparseMatrix S.SparseMatrix S.SparseMatrix where-    (%*%) = withSS Internal.c_ss_mul+instance Arithmetic S.SparseMatrix S.SparseMatrix where+    (@@) = withSS Internal.c_ss_mul+    (%+%) = withSS Internal.c_ss_plus+    (%-%) a b = a %+% C.map negate b+    (%*%) = withSS Internal.c_ss_cmul -inverse :: (SingI n, Numeric a) => Matrix n n a -> Matrix n n a-inverse = withFun1 Internal.c_inverse  class Factorization mat where+    -- | Matrix inverse+    inverse :: (SingI n, Numeric a) => mat n n VS.Vector a -> mat n n VS.Vector a+     -- | Eigenvalues (not ordered) and     -- eigenvectors (as columns) of a general square matrix.-    eigs :: (SingI k, SingI n, k <= n - 2)+    eigs :: (SingI k, SingI n, (k <= n - 2) ~ 'True)          => Sing k          -> mat n n VS.Vector Double          -> (Matrix k 1 (Complex Double), Matrix n k (Complex Double))@@ -65,6 +98,9 @@     cholesky :: (Numeric a, SingI n) => mat n n VS.Vector a -> mat n n VS.Vector a  instance Factorization D.Matrix where++    inverse = withFun1 Internal.c_inverse+     eigs s mat = unsafePerformIO $ do         m1 <- CM.new         m2 <- CM.new@@ -82,6 +118,8 @@     {-# INLINE cholesky #-}  instance Factorization S.SparseMatrix where+    inverse = undefined+     eigs s mat = unsafePerformIO $ do         m1 <- CM.new         m2 <- CM.new
src/Data/Matrix/Static/LinearAlgebra/Types.hs view
@@ -38,7 +38,7 @@ import qualified Data.Matrix.Static.Generic.Mutable as CM import qualified Data.Matrix.Static.Generic as C -class (S.Zero a, Storable a) => Numeric a where+class (S.Zero a, Storable a, Num a) => Numeric a where     foreignType :: a -> CInt  instance Numeric Float where foreignType _ = 0@@ -127,16 +127,19 @@ {-# INLINE withSD #-}  mkSparseMatrix :: forall r c a. (Storable a, SingI r, SingI c)-    => (Ptr (Ptr a) -> Ptr (Ptr CInt) -> Ptr (Ptr CInt) -> IO Int)+    => (Ptr (Ptr a) -> Ptr (Ptr CInt) -> Ptr CInt -> IO Int)     -> IO (SparseMatrix r c a)-mkSparseMatrix f = alloca $ \ppv -> alloca $ \ppi -> alloca $ \ppo -> do-    n <- f ppv ppi ppo-    pv <- peek ppv >>= newForeignPtr finalizerFree-    pinner <- peek ppi >>= newForeignPtr finalizerFree-    pouter <- peek ppo >>= newForeignPtr finalizerFree+mkSparseMatrix f = do+    outer' <- VSM.new $ c + 1+    (n, pv, pinner) <- VSM.unsafeWith outer' $ \pouter -> alloca $ \ppv -> alloca $ \ppi -> do+        n <- f ppv ppi pouter+        pv <- peek ppv >>= newForeignPtr finalizerFree+        pinner <- peek ppi >>= newForeignPtr finalizerFree+        return (n, pv, pinner)+    outer <- VS.unsafeFreeze outer'     return $ S.SparseMatrix (VS.unsafeFromForeignPtr0 pv n)         (VS.unsafeFromForeignPtr0 pinner n)-        (VS.unsafeFromForeignPtr0 pouter $ c + 1)+        outer   where     c = fromIntegral $ fromSing (sing :: Sing c) {-# INLINE mkSparseMatrix #-}@@ -144,7 +147,7 @@ withSS :: forall r1 c1 r2 c2 r3 c3 a.             (SingI r3, SingI c3, Numeric a)        => ( CInt-         -> Ptr (Ptr a) -> Ptr (Ptr CInt) -> Ptr (Ptr CInt) -> CInt -> CInt -> Ptr CInt+         -> Ptr (Ptr a) -> Ptr CInt -> Ptr (Ptr CInt) -> CInt -> CInt -> Ptr CInt          -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt          -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt          -> IO CString )
src/Data/Matrix/Static/Sparse.hs view
@@ -34,10 +34,10 @@      -- * Construction     , C.empty+    , fromTriplet     , C.fromVector     , C.fromList     , C.unsafeFromVector-     , diag     , diagRect @@ -49,12 +49,13 @@     , C.convertAny    ) where -import           Control.DeepSeq import qualified Data.Vector.Generic as G+import qualified Data.Vector.Generic.Mutable as GM import qualified Data.Vector.Storable as S import qualified Data.Vector.Storable.Mutable as SM import Data.Singletons import Control.Monad+import Control.Monad.ST (runST) import           Data.Bits                         (shiftR) import Text.Printf (printf) import GHC.TypeLits (type (<=))@@ -106,9 +107,22 @@         (r,c) = C.dim mat         vals = unlines $ map (unwords . map show . G.toList) $ C.toRows mat -instance (NFData (v a)) => NFData (SparseMatrix r c v a) where-    rnf (SparseMatrix vec inner outer) = rnf vec+instance (SingI r, SingI c, G.Vector v a, Zero a, Num a) =>+    Num (SparseMatrix r c v a) where+        m1 + m2 = undefined+        m1 - m2 = undefined+        m1 * m2 = undefined+        negate = C.map negate+        abs = C.map abs+        signum = undefined+        fromInteger = undefined +instance (SingI r, SingI c, G.Vector v a, Zero a, Fractional a) =>+    Fractional (SparseMatrix r c v a) where+        m1 / m2 = undefined+        recip = C.map recip+        fromRational = undefined+ instance (G.Vector v a, Zero a) => C.Matrix SparseMatrix v a where     -- | O(1) Return the size of matrix.     dim :: forall r c. SparseMatrix r c v a -> (Int, Int)@@ -151,6 +165,8 @@         c = fromIntegral $ fromSing (sing :: Sing c)     {-# INLINE unsafeFromVector #-} +    transpose (SparseMatrix val inner outer) = undefined +     thaw = undefined     {-# INLINE thaw #-} @@ -167,7 +183,38 @@     imap = undefined     {-# INLINE map #-} --- | O(m*n) Create a square matrix with given diagonal.+-- | O(n) Create matrix from triplet. row and column indices *are not* assumed to be ordered+-- duplicate entries are carried over to the CSR represention+fromTriplet :: forall t r c v a. (Traversable t, G.Vector v a, SingI r, SingI c)+            => t (Int, Int, a) -> SparseMatrix r c v a+fromTriplet triplets = SparseMatrix val inner outer+  where+    outer = S.scanl (+) 0 $ S.create $ do+        vec <- SM.replicate c 0+        _ <- flip mapM triplets $ \(_, j, _) -> +            SM.modify vec (+1) j+        return vec+    (val, inner) = runST $ do+        outer' <- S.thaw outer+        val' <- GM.new nnz+        inner' <- SM.new nnz+        _ <- flip mapM triplets $ \(i, j, v) -> do+            idx <- fromIntegral <$> SM.read outer' j+            GM.write val' idx v+            SM.write inner' idx $ fromIntegral i+            SM.modify outer' (+1) j+        (,) <$> G.unsafeFreeze val' <*> S.unsafeFreeze inner'+    nnz = length triplets+    c = fromIntegral $ fromSing (sing :: Sing c)+{-# INLINE fromTriplet #-}++{-+toTriplet :: (G.Vector v1 (Int, Int, a), G.Vector v2 a, SingI r, SingI c)+          => SparseMatrix r c v2 a -> v1 (Int, Int, a)+toTriplet mat = +-}++-- | O(m*n) Create a rectangular matrix with default values and given diagonal diag :: (G.Vector v a, Zero a, SingI n)      => D.Matrix n 1 v a       -- ^ diagonal      -> SparseMatrix n n v a
src/Data/Matrix/Static/Sparse/Mutable.hs view
@@ -13,12 +13,10 @@      MSparseMatrix(..)    ) where -import           Control.DeepSeq import qualified Data.Vector.Generic.Mutable as GM import qualified Data.Vector.Storable as S import           Prelude                     hiding (read, replicate) import Data.Singletons-import Control.Monad.Primitive     (PrimMonad, PrimState)  import qualified Data.Matrix.Static.Generic.Mutable as C @@ -33,9 +31,6 @@                                       -- (resp. row) the index of the first                                       -- non-zero in the previous two arrays.                   -> MSparseMatrix r c v s a--instance (NFData (v s a)) => NFData (MSparseMatrix r c v s a) where-    rnf (MSparseMatrix vec inner outer) = rnf vec  instance GM.MVector v a => C.MMatrix MSparseMatrix v a where     dim :: forall r c s. MSparseMatrix r c v s a -> (Int, Int)
+ tests/Main.hs view
@@ -0,0 +1,36 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE FlexibleInstances #-}++module Main where++import Test.Tasty+import Data.Matrix.Static.LinearAlgebra+import qualified Data.Matrix.Static.Dense as D+import qualified Data.Matrix.Static.Sparse as S+import Data.Singletons+import Data.Vector.Storable (Storable)+import GHC.TypeNats (KnownNat)+import Control.Monad.IO.Class (liftIO)++import Test.Tasty.QuickCheck++instance (Arbitrary a, KnownNat m, KnownNat n) => Arbitrary (Matrix m n a) where+    arbitrary = D.fromList <$> vector (m*n)+      where+        m = fromIntegral $ fromSing (sing :: Sing m)+        n = fromIntegral $ fromSing (sing :: Sing n)+    shrink _v = []++{-+propTranspose :: Matrix 50 100 Double -> Bool+propTranspose m = D.transpose (D.transpose m) == m && +    D.convertAny (S.transpose $ S.transpose (D.convertAny m)) == m+-}++main :: IO ()+main = defaultMain $ testGroup "Main"+    [ testProperty "" square+    ]+