packages feed

matrix-sized 0.0.3 → 0.0.4

raw patch · 15 files changed

+797/−298 lines, 15 filesdep +conduitdep +data-ordlistdep +ieee754

Dependencies added: conduit, data-ordlist, ieee754

Files

README.md view
@@ -37,7 +37,7 @@  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+f m = let (d, v) = eigS (sing :: Sing 2) m       in v @@ S.diag d @@ G.transpose v  main :: IO ()
cbits/eigen-runtime.h view
@@ -39,6 +39,18 @@ 	GUARD_END\ } +#define API2(name,args,call) \+extern "C" RET eigen_##name args {\+	GUARD_START\+	switch (code) {\+		case 0: return name<T0, T0>call;\+		case 1: return name<T1, T1>call;\+		case 2: return name<T2, T0>call;\+		case 3: return name<T3, T1>call;\+	}\+	GUARD_END\+}+ typedef float T0; typedef double T1; typedef std::complex<float> T2;
cbits/eigen-solver.cpp view
@@ -1,4 +1,5 @@ #include <Spectra/GenEigsSolver.h>+#include <Spectra/SymEigsSolver.h> #include "eigen-runtime.h" #include <Eigen/Sparse> #include <Spectra/MatOp/SparseGenMatProd.h>@@ -6,6 +7,20 @@ using namespace Spectra; using namespace Eigen; +extern "C" RET eigen_eig( +    void* d, void* v, +    const void* p, int n)+{+    typedef Map< Matrix<T1,Dynamic,Dynamic> > MapMatrix;+    typedef Map< Matrix<T3,Dynamic,Dynamic> > MapComplexMatrix;+    MapMatrix M((T1*)p, n, n);+    MapComplexMatrix D((T3*)d, n, 1);+    MapComplexMatrix V((T3*)v, n, n);+    EigenSolver<MatrixXd> es(M);+    D = es.eigenvalues();+    V = es.eigenvectors();+}+ extern "C" RET spectral_eigs(      int k,     void* d, void* v, @@ -29,6 +44,28 @@     return 0; } +extern "C" RET spectral_eigsh( +    int k,+    void* d, void* v, +    const void* p, int n)+{+    typedef Map< Matrix<T1,Dynamic,Dynamic> > MapMatrix;+    MapMatrix M((T1*)p, n, n);+    MapMatrix D((T1*)d, k, 1);+    MapMatrix V((T1*)v, n, k);++    DenseGenMatProd<double> op(M);+    int ncv = 2 * k;+    ncv = (ncv <= n) ? ncv : n;+    SymEigsSolver< double, LARGEST_MAGN, DenseGenMatProd<double> > eigsh(&op, k, ncv);+    eigsh.init();+    int nconv = eigsh.compute();+    if(eigsh.info() == 0)+        D = eigsh.eigenvalues();+        V = eigsh.eigenvectors();+    return 0;+}+ extern "C" RET spectral_seigs(      int k,     void* d, void* v,@@ -55,6 +92,33 @@     return 0; } +extern "C" RET spectral_seigsh( +    int k,+    void* d, void* v,+    const void* values,+    const void* outerIndexPtr,+    const void* innerIndices,+    int n, int s)+{+    typedef Map< Matrix<T1,Dynamic,Dynamic> > MapMatrix;+    typedef Map<const SparseMatrix<T1> > MapSparseMatrix;+    MapSparseMatrix M(n, n, s, (int*)outerIndexPtr, (int*)innerIndices, (T1*)values);+    MapMatrix D((T1*)d, k, 1);+    MapMatrix V((T1*)v, n, k);++    SparseGenMatProd<double> op(M);+    int ncv = 2 * k;+    ncv = (ncv <= n) ? ncv : n;+    SymEigsSolver< double, LARGEST_MAGN, SparseGenMatProd<double> > eigsh(&op, k, ncv);+    eigsh.init();+    int nconv = eigsh.compute();+    if(eigsh.info() == 0)+        D = eigsh.eigenvalues();+        V = eigsh.eigenvectors();+    return 0;+}++ template <class T> RET cholesky(void* px, const void* pa, int n) {@@ -68,16 +132,23 @@     void* px, const void* pa, int n), (px,pa,n));  -/*-template <class T>-RET bdcsvd(void* px, int r, int c+template <class T, class TT>+RET bdcsvd(+    void* pu, void* ps, void* pv, +    const void* px, int r, int c) {+    int m = r < c ? r : c;     typedef Map< Matrix<T,Dynamic,Dynamic> > MapMatrix;-    MapMatrix x((T*)px, n, n);-    MapMatrix A((T*)pa, n, n);-    x = A.llt().matrixL();+    typedef Map< Matrix<TT,Dynamic,Dynamic> > MapMatrix2;+    MapMatrix A((T*)px, r, c);+    MapMatrix U((T*)pu, r, m);+    MapMatrix2 s((TT*)ps, m, 1);+    MapMatrix V((T*)pv, c, m);+    BDCSVD< Matrix<T,Dynamic,Dynamic>> svd(A, ComputeThinU|ComputeThinV);+    U = svd.matrixU();+    V = svd.matrixV();+    s = svd.singularValues();     return 0; }-API(cholesky, (int code,-    void* px, const void* pa, int n), (px,pa,n));-*/+API2(bdcsvd, (int code,+    void* pu, void* ps, void* pv, const void* px, int r, int c), (pu,ps,pv,px,r,c));
matrix-sized.cabal view
@@ -4,10 +4,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: ef7b97a49f99defe3cd58f2c5a3411fe8c958f0ad40f4b3f0bd763720b643db8+-- hash: a94165a01f28a1a01026b95b46341900f92789d9a5b17aec1fe9ab414b560a26  name:           matrix-sized-version:        0.0.3+version:        0.0.4 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@@ -392,6 +392,7 @@       Data.Matrix.Static.Generic.Mutable   other-modules:       Data.Matrix.Static.Internal+      Data.Matrix.Static.LinearAlgebra.Internal   hs-source-dirs:       src   ghc-options: -Wall@@ -408,6 +409,7 @@       stdc++   build-depends:       base >=4.10 && <5+    , conduit     , primitive >=0.6.4.0     , singletons     , vector >=0.11@@ -420,11 +422,18 @@ test-suite test   type: exitcode-stdio-1.0   main-is: Main.hs+  other-modules:+      Test.Base+      Test.LinearAlgebra+      Test.Utils   hs-source-dirs:       tests   ghc-options: -Wall -threaded   build-depends:       base+    , conduit+    , data-ordlist+    , ieee754     , matrix-sized     , primitive >=0.6.4.0     , singletons
src/Data/Matrix/Static/Dense.hs view
@@ -42,6 +42,7 @@     , C.fromColumns     , C.unsafeFromVector     , replicate+    , diag     , diagRect      -- * Conversions@@ -88,11 +89,17 @@     , C.freeze     , C.unsafeFreeze     , C.create++    , sum+    , all+    , any     ) where  import           Control.Monad                     (liftM) import qualified Data.Vector.Generic               as G-import Prelude hiding (replicate, 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, sum, all, any ) import GHC.TypeLits (type (<=)) import Data.Singletons import Data.Tuple (swap)@@ -203,6 +210,17 @@     c = fromIntegral $ fromSing (sing :: Sing c) {-# INLINE replicate #-} +-- | O(m*n) Create a square matrix with default values and given diagonal+diag :: (G.Vector v a, SingI n)+     => a                    -- ^ default value+     -> Matrix n 1 v a       -- ^ diagonal+     -> Matrix n n v a+diag z0 d = C.create $ do+    mat <- DM.replicate z0+    C.imapM_ (DM.unsafeWrite mat) d+    return mat+{-# INLINE diag #-}+ -- | 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)          => a                    -- ^ default value@@ -352,8 +370,21 @@ convert (Matrix vec) = Matrix $ G.convert vec {-# INLINE convert #-} +sum :: (Num a, G.Vector v a) => Matrix r c v a -> a+sum (Matrix vec) = G.sum vec+{-# INLINE sum #-} +all :: G.Vector v a => (a -> Bool) -> Matrix r c v a -> Bool+all f (Matrix vec) = G.all f vec+{-# INLINE all #-}++any :: G.Vector v a => (a -> Bool) -> Matrix r c v a -> Bool+any f (Matrix vec) = G.any f vec+{-# INLINE any #-}+ -- Helper-toIndex :: Int -> Int -> (Int, Int)+toIndex :: Int   -- ^ Number of rows+        -> Int   -- ^ 1-d index+        -> (Int, Int)   -- ^ 2-d index toIndex r i = swap $ i `divMod` r {-# INLINE toIndex #-}
src/Data/Matrix/Static/Internal.hs view
@@ -9,8 +9,12 @@     , c_ss_plus     , c_inverse     , c_cholesky+    , c_eig     , c_eigs+    , c_eigsh     , c_seigs+    , c_seigsh+    , c_bdcsvd     ) where  import Data.Complex (Complex)@@ -80,11 +84,27 @@     c_cholesky :: CInt                -> Ptr a -> Ptr a ->  CInt -> IO CString +foreign import ccall "eigen_eig"+    c_eig :: Ptr (Complex Double) -> Ptr (Complex Double)+          -> Ptr Double -> CInt -> IO CString+ foreign import ccall "spectral_eigs"     c_eigs :: CInt -> Ptr (Complex Double)            -> Ptr (Complex Double) -> Ptr Double -> CInt -> IO CString +foreign import ccall "spectral_eigsh"+    c_eigsh :: CInt -> Ptr Double -> Ptr Double -> Ptr Double -> CInt -> IO CString+ foreign import ccall "spectral_seigs"     c_seigs :: CInt -> Ptr (Complex Double) -> Ptr (Complex Double)             -> Ptr Double -> Ptr CInt -> Ptr CInt             -> CInt -> CInt -> IO CString++foreign import ccall "spectral_seigsh"+    c_seigsh :: CInt -> Ptr Double -> Ptr Double+             -> Ptr Double -> Ptr CInt -> Ptr CInt+             -> CInt -> CInt -> IO CString++foreign import ccall "eigen_bdcsvd"+    c_bdcsvd :: CInt -> Ptr a -> Ptr b -> Ptr a+             -> Ptr a -> CInt -> CInt -> IO CString
src/Data/Matrix/Static/LinearAlgebra.hs view
@@ -6,15 +6,26 @@ {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE FunctionalDependencies #-} {-# LANGUAGE GADTs #-}+{-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE ExplicitNamespaces #-}+{-# LANGUAGE TypeFamilies #-} {-# LANGUAGE TypeOperators #-}+{-# LANGUAGE InstanceSigs #-} module Data.Matrix.Static.LinearAlgebra-    ( Arithmetic(..)+    ( module Data.Matrix.Static.LinearAlgebra.Types+    , Arithmetic(..)     , Factorization(..)-    , module Data.Matrix.Static.LinearAlgebra.Types+    , LinearAlgebra(..)++      -- * Dense matrix operation+    , inverse+    , eig+    , svd+    , cond     ) where  import qualified Data.Vector.Storable as VS+import Data.Vector.Storable (Vector) import System.IO.Unsafe (unsafePerformIO) import Data.Complex (Complex) import Data.Singletons.Prelude hiding ((@@), type (==))@@ -27,35 +38,39 @@ import qualified Data.Matrix.Static.Generic as C import qualified Data.Matrix.Static.Internal as Internal import Data.Matrix.Static.LinearAlgebra.Types+import Data.Matrix.Static.LinearAlgebra.Internal  class Arithmetic (mat1 :: C.MatrixKind) (mat2 :: C.MatrixKind) where-    -- | Matrix multiplication+    -- | Matrix multiplication between different types of matrices.     (@@) :: ( 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+         => mat1 n p Vector a+         -> mat2 p m Vector a+         -> mat3 n m Vector a     infixr 8 @@ +    -- | Element-wise addition between different types of matrices.     (%+%) :: ( 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+          => mat1 n m Vector a+          -> mat2 n m Vector a+          -> mat3 n m Vector a     infixr 8 %+% +    -- | Element-wise substraction between different types of matrices.     (%-%) :: ( 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+          => mat1 n m Vector a+          -> mat2 n m Vector a+          -> mat3 n m Vector a     infixr 8 %-% +    -- | Element-wise multiplication between different types of matrices.     (%*%) :: ( 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+          => mat1 n m Vector a+          -> mat2 n m Vector a+          -> mat3 n m Vector a     infixr 8 %*%  instance Arithmetic D.Matrix D.Matrix where@@ -82,45 +97,69 @@     (%-%) a b = a %+% C.map negate b     (%*%) = withSS Internal.c_ss_cmul +class LinearAlgebra (mat :: C.MatrixKind) where+    ident :: (Numeric a, SingI n) => mat n n Vector a -class Factorization mat where-    -- | Matrix inverse-    inverse :: (SingI n, Numeric a) => mat n n VS.Vector a -> mat n n VS.Vector a+instance LinearAlgebra D.Matrix where+    ident = D.diag 0 $ D.replicate 1 -    -- | Eigenvalues (not ordered) and+instance LinearAlgebra S.SparseMatrix where+    ident = S.diag $ D.replicate 1++class Factorization mat where+    -- | Eigenvalues (from largest to smallest) and     -- eigenvectors (as columns) of a general square matrix.-    eigs :: (SingI k, SingI n, (k <= n - 2) ~ 'True)+    eigS :: (SingI k, SingI n, (k <= n - 2) ~ 'True)          => Sing k-         -> mat n n VS.Vector Double+         -> mat n n Vector Double          -> (Matrix k 1 (Complex Double), Matrix n k (Complex Double)) +    -- | Eigenvalues (from largest to smallest) and+    -- eigenvectors (as columns) of a symmetric square matrix.+    eigSH :: (SingI k, SingI n, (k <= n - 1) ~ 'True)+          => Sing k+          -> mat n n Vector Double+          -> (Matrix k 1 Double, Matrix n k Double)+     -- | Cholesky decomposition-    cholesky :: (Numeric a, SingI n) => mat n n VS.Vector a -> mat n n VS.Vector a+    cholesky :: (Numeric a, SingI n) => mat n n Vector a -> mat n n Vector a -instance Factorization D.Matrix where -    inverse = withFun1 Internal.c_inverse+instance Factorization D.Matrix where+    eigS s mat+        | D.all (==0) mat = ( D.replicate 0, D.replicate 1)+        | otherwise = unsafePerformIO $ do+            m1 <- CM.new+            m2 <- CM.new+            _ <- unsafeWith' m1 $ \v1 _ _ -> unsafeWith' m2 $ \v2 _ _ ->+                unsafeWith mat $ \v n _ -> Internal.c_eigs k v1 v2 v n+            m1' <- C.unsafeFreeze m1+            m2' <- C.unsafeFreeze m2+            return (m1', m2')+      where+        k = fromIntegral $ fromSing s+    {-# INLINE eigS #-} -    eigs s mat = unsafePerformIO $ do-        m1 <- CM.new-        m2 <- CM.new-        _ <- unsafeWith' m1 $ \v1 _ _ -> unsafeWith' m2 $ \v2 _ _ -> do-            unsafeWith mat $ \v n _ -> Internal.c_eigs k v1 v2 v n-        m1' <- C.unsafeFreeze m1-        m2' <- C.unsafeFreeze m2-        return (m1', m2')+    eigSH s mat+        | D.all (==0) mat = (D.replicate 0, D.replicate 1)+        | otherwise = unsafePerformIO $ do+            m1 <- CM.new+            m2 <- CM.new+            _ <- unsafeWith' m1 $ \v1 _ _ -> unsafeWith' m2 $ \v2 _ _ ->+                unsafeWith mat $ \v n _ -> Internal.c_eigsh k v1 v2 v n+            m1' <- C.unsafeFreeze m1+            m2' <- C.unsafeFreeze m2+            return (m1', m2')       where         k = fromIntegral $ fromSing s-    {-# INLINE eigs #-}+    {-# INLINE eigSH #-}      cholesky mat = flip withFun1 mat $         \code p1 c1 _ p2 _ _ -> Internal.c_cholesky code p1 p2 c1     {-# INLINE cholesky #-}  instance Factorization S.SparseMatrix where-    inverse = undefined--    eigs s mat = unsafePerformIO $ do+    eigS s mat = unsafePerformIO $ do         m1 <- CM.new         m2 <- CM.new         _ <- unsafeWith' m1 $ \v1 _ _ -> unsafeWith' m2 $ \v2 _ _ ->@@ -131,6 +170,72 @@         return (m1', m2')       where         k = fromIntegral $ fromSing s-    {-# INLINE eigs #-}+    {-# INLINE eigS #-} +    eigSH s mat = unsafePerformIO $ do+        m1 <- CM.new+        m2 <- CM.new+        _ <- unsafeWith' m1 $ \v1 _ _ -> unsafeWith' m2 $ \v2 _ _ ->+            unsafeWithS mat $ \pv pin po n _ size ->+                Internal.c_seigsh k v1 v2 pv po pin n size+        m1' <- C.unsafeFreeze m1+        m2' <- C.unsafeFreeze m2+        return (m1', m2')+      where+        k = fromIntegral $ fromSing s+    {-# INLINE eigSH #-}+     cholesky = undefined++type family R a where+    R Float = Float+    R Double = Double+    R (Complex Double) = Double+    R (Complex Float) = Float+    +-- | The inverse of a dense matrix.+inverse :: (SingI n, Numeric a) => Matrix n n a -> Matrix n n a+inverse = withFun1 Internal.c_inverse+{-# INLINE inverse #-}++-- | Compute the full eigendecomposition for dense matrix.+eig :: forall n . SingI n+    => Matrix n n Double+    -> (Matrix n 1 (Complex Double), Matrix n n (Complex Double))+eig mat = unsafePerformIO $ do+    m1 <- CM.new+    m2 <- CM.new+    _ <- unsafeWith' m1 $ \v1 _ _ -> unsafeWith' m2 $ \v2 _ _ ->+        unsafeWith mat $ \v n _ -> Internal.c_eig v1 v2 v n+    m1' <- C.unsafeFreeze m1+    m2' <- C.unsafeFreeze m2+    return (m1', m2')+{-# INLINE eig #-}++-- | Compute the full singular value decomposition for dense matrix.+svd :: forall n p a m. (Numeric (R a), Numeric a, SingI n, SingI p, SingI m, m ~ Min n p)+    => Matrix n p a+    -> (Matrix n m a, Matrix m 1 (R a), Matrix p m a)+svd mat = unsafePerformIO $ do+    mu <- CM.new+    ms <- CM.new+    mv <- CM.new+    checkResult $ unsafeWith' mu $ \pu _ _ -> unsafeWith' ms $ \ps _ _ ->+        unsafeWith' mv $ \pv _ _ -> unsafeWith mat $ \px r c ->+            Internal.c_bdcsvd (foreignType (undefined :: a))+                pu ps pv px r c+    u <- C.unsafeFreeze mu+    s <- C.unsafeFreeze ms+    v <- C.unsafeFreeze mv+    return (u, s, v)+{-# INLINE svd #-}++-- | Condition number.+cond :: ( Numeric a, Numeric (R a), Ord (R a), Fractional (R a)+        , SingI n, SingI m, SingI (Min n m))+     => Matrix n m a -> R a+cond mat = VS.maximum val / VS.minimum val+  where+    val = VS.filter (/=0) $ D.flatten s+    (_,s,_) = svd mat+{-# INLINE cond #-}
+ src/Data/Matrix/Static/LinearAlgebra/Internal.hs view
@@ -0,0 +1,189 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE DataKinds #-}+{-# OPTIONS_GHC -Wno-simplifiable-class-constraints #-}+module Data.Matrix.Static.LinearAlgebra.Internal+    ( withFun1+    , withFun2+    , withDS+    , withSD+    , withSS+    , checkResult+    , unsafeWith+    , unsafeWith'+    , unsafeWithS+    ) where++import Data.Vector.Storable (Storable)+import qualified Data.Vector.Storable as VS+import qualified Data.Vector.Storable.Mutable as VSM+import System.IO.Unsafe (unsafePerformIO)+import Control.Monad (when)+import Control.Monad.ST (RealWorld)+import Data.Singletons+import Foreign+import Foreign.C.Types+import Foreign.C.String++import qualified Data.Matrix.Static.Dense as D+import qualified Data.Matrix.Static.Dense.Mutable as DM+import qualified Data.Matrix.Static.Sparse as S+import qualified Data.Matrix.Static.Generic.Mutable as CM+import qualified Data.Matrix.Static.Generic as C+import Data.Matrix.Static.LinearAlgebra.Types++withFun1 :: forall r1 c1 r2 c2 a. (SingI r2, SingI c2, Numeric a)+         => (CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString)+         -> Matrix r1 c1 a -> Matrix r2 c2 a+withFun1 f m1 = unsafePerformIO $ do+    m0 <- CM.new+    checkResult $ unsafeWith' m0 $ \vals0 rows0 cols0 ->+        unsafeWith m1 $ \vals1 rows1 cols1 -> f (foreignType (undefined :: a))+            vals0 rows0 cols0+            vals1 rows1 cols1+    C.unsafeFreeze m0+{-# INLINE withFun1 #-}++withFun2 :: forall r1 c1 r2 c2 r3 c3 a.+            (SingI r3, SingI c3, Numeric a)+         => ( CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt+           -> Ptr a -> CInt -> CInt -> IO CString )+         -> Matrix r1 c1 a+         -> Matrix r2 c2 a+         -> Matrix r3 c3 a+withFun2 f m1 m2 = unsafePerformIO $ do+    m0 <- CM.new+    checkResult $ unsafeWith' m0 $ \vals0 rows0 cols0 ->+        unsafeWith m1 $ \vals1 rows1 cols1 ->+            unsafeWith m2 $ \vals2 rows2 cols2 ->+                f (foreignType (undefined :: a))+                    vals0 rows0 cols0+                    vals1 rows1 cols1+                    vals2 rows2 cols2+    C.unsafeFreeze m0+{-# INLINE withFun2 #-}++withDS :: forall r1 c1 r2 c2 r3 c3 a.+            (SingI r3, SingI c3, Numeric a)+       => ( CInt+         -> Ptr a -> CInt -> CInt+         -> Ptr a -> CInt -> CInt+         -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt+         -> IO CString )+       -> Matrix r1 c1 a+       -> SparseMatrix r2 c2 a+       -> Matrix r3 c3 a+withDS f m1 m2 = unsafePerformIO $ do+    m0 <- CM.new+    checkResult $ unsafeWith' m0 $ \v0 r0 c0 ->+        unsafeWith m1 $ \v1 r1 c1 ->+            unsafeWithS m2 $ \v2 inner outer r2 c2 s ->+                f (foreignType (undefined :: a))+                    v0 r0 c0+                    v1 r1 c1+                    v2 outer inner r2 c2 s+    C.unsafeFreeze m0+{-# INLINE withDS #-}++withSD :: forall r1 c1 r2 c2 r3 c3 a.+            (SingI r3, SingI c3, Numeric a)+       => ( CInt+         -> Ptr a -> CInt -> CInt+         -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt+         -> Ptr a -> CInt -> CInt+         -> IO CString )+       -> SparseMatrix r2 c2 a+       -> Matrix r1 c1 a+       -> Matrix r3 c3 a+withSD f m2 m1 = unsafePerformIO $ do+    m0 <- CM.new+    checkResult $ unsafeWith' m0 $ \v0 r0 c0 ->+        unsafeWith m1 $ \v1 r1 c1 ->+            unsafeWithS m2 $ \v2 inner outer r2 c2 s ->+                f (foreignType (undefined :: a))+                    v0 r0 c0+                    v2 outer inner r2 c2 s+                    v1 r1 c1+    C.unsafeFreeze m0+{-# INLINE withSD #-}++mkSparseMatrix :: forall r c a. (Storable a, SingI r, SingI c)+    => (Ptr (Ptr a) -> Ptr (Ptr CInt) -> Ptr CInt -> IO Int)+    -> IO (SparseMatrix r c a)+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)+        outer+  where+    c = fromIntegral $ fromSing (sing :: Sing c)+{-# INLINE mkSparseMatrix #-}++withSS :: forall r1 c1 r2 c2 r3 c3 a.+            (SingI r3, SingI c3, Numeric a)+       => ( 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 )+       -> SparseMatrix r1 c1 a+       -> SparseMatrix r2 c2 a+       -> SparseMatrix r3 c3 a+withSS f m1 m2 = unsafePerformIO $ mkSparseMatrix $ \v0 inner0 outer0 ->+    alloca $ \pn -> unsafeWithS m1 $ \v1 inner1 outer1 r1 c1 s1 ->+        unsafeWithS m2 $ \v2 inner2 outer2 r2 c2 s2 -> do+            checkResult $ f (foreignType (undefined :: a))+                v0 outer0 inner0 r c pn+                v1 outer1 inner1 r1 c1 s1+                v2 outer2 inner2 r2 c2 s2+            fromIntegral <$> peek pn+  where+    r = fromIntegral $ fromSing (sing :: Sing r3)+    c = fromIntegral $ fromSing (sing :: Sing c3)+{-# INLINE withSS #-}++checkResult :: IO CString -> IO ()+checkResult func = func >>= \c_str -> when (c_str /= nullPtr) $+    peekCString c_str >>= \str -> error str+{-# INLINE checkResult #-}++-------------------------------------------------------------------------------+-- Raw pointers+-------------------------------------------------------------------------------++-- | Pass a pointer to the matrix's data to the IO action.+-- The data may not be modified through the pointer.+unsafeWith :: Storable a => Matrix n m a -> (Ptr a -> CInt -> CInt -> IO b) -> IO b+unsafeWith mat@(D.Matrix vec) f = VS.unsafeWith vec $ \p ->+    f p (fromIntegral r) $ fromIntegral c +  where+    (r,c) = C.dim mat+{-# INLINE unsafeWith #-}++unsafeWith' :: Storable a => MMatrix n m RealWorld a -> (Ptr a -> CInt -> CInt -> IO b) -> IO b+unsafeWith' mat@(DM.MMatrix vec) f = VSM.unsafeWith vec $ \p ->+    f p (fromIntegral r) $ fromIntegral c+  where+    (r,c) = CM.dim mat+{-# INLINE unsafeWith' #-}++-- | Pass a pointer to the matrix's data to the IO action.+-- The data may not be modified through the pointer.+unsafeWithS :: (Storable a, S.Zero a)+            => SparseMatrix n m a+            -> (Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt -> IO b)+            -> IO b+unsafeWithS mat@(S.SparseMatrix val inner outer) f = VS.unsafeWith val $ \pval ->+    VS.unsafeWith inner $ \pinner -> VS.unsafeWith outer $ \pouter ->+        f pval pinner pouter (fromIntegral r) (fromIntegral c) (fromIntegral $ VS.length val)+  where+    (r,c) = C.dim mat+{-# INLINE unsafeWithS #-}
src/Data/Matrix/Static/LinearAlgebra/Types.hs view
@@ -9,39 +9,22 @@     , Matrix     , MMatrix     , SparseMatrix-    , withFun1-    , withFun2-    , withDS-    , withSD-    , withSS-    , unsafeWith-    , unsafeWith'-    , unsafeWithS     ) where  import Data.Vector.Storable (Vector, Storable)-import qualified Data.Vector.Storable as VS-import qualified Data.Vector.Storable.Mutable as VSM import Data.Vector.Storable.Mutable (MVector)-import System.IO.Unsafe (unsafePerformIO)-import Control.Monad (when) import Data.Complex (Complex)-import Control.Monad.ST (RealWorld)-import Data.Singletons-import Foreign import Foreign.C.Types-import Foreign.C.String  import qualified Data.Matrix.Static.Dense as D import qualified Data.Matrix.Static.Dense.Mutable as DM import qualified Data.Matrix.Static.Sparse as S-import qualified Data.Matrix.Static.Generic.Mutable as CM-import qualified Data.Matrix.Static.Generic as C  class (S.Zero a, Storable a, Num a) => Numeric a where     foreignType :: a -> CInt  instance Numeric Float where foreignType _ = 0+instance Numeric CFloat where foreignType _ = 0 instance Numeric Double where foreignType _ =1 instance Numeric (Complex Float) where foreignType _ = 2 instance Numeric (Complex Double) where foreignType _ = 3@@ -50,157 +33,3 @@ type MMatrix r c s a = DM.MMatrix r c MVector s a  type SparseMatrix r c a = S.SparseMatrix r c Vector a--withFun1 :: forall r1 c1 r2 c2 a. (SingI r2, SingI c2, Numeric a)-         => (CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString)-         -> Matrix r1 c1 a -> Matrix r2 c2 a-withFun1 f m1 = unsafePerformIO $ do-    m0 <- CM.new-    checkResult $ unsafeWith' m0 $ \vals0 rows0 cols0 ->-        unsafeWith m1 $ \vals1 rows1 cols1 -> f (foreignType (undefined :: a))-            vals0 rows0 cols0-            vals1 rows1 cols1-    C.unsafeFreeze m0-{-# INLINE withFun1 #-}--withFun2 :: forall r1 c1 r2 c2 r3 c3 a.-            (SingI r3, SingI c3, Numeric a)-         => ( CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt-           -> Ptr a -> CInt -> CInt -> IO CString )-         -> Matrix r1 c1 a-         -> Matrix r2 c2 a-         -> Matrix r3 c3 a-withFun2 f m1 m2 = unsafePerformIO $ do-    m0 <- CM.new-    checkResult $ unsafeWith' m0 $ \vals0 rows0 cols0 ->-        unsafeWith m1 $ \vals1 rows1 cols1 ->-            unsafeWith m2 $ \vals2 rows2 cols2 ->-                f (foreignType (undefined :: a))-                    vals0 rows0 cols0-                    vals1 rows1 cols1-                    vals2 rows2 cols2-    C.unsafeFreeze m0-{-# INLINE withFun2 #-}--withDS :: forall r1 c1 r2 c2 r3 c3 a.-            (SingI r3, SingI c3, Numeric a)-       => ( CInt-         -> Ptr a -> CInt -> CInt-         -> Ptr a -> CInt -> CInt-         -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt-         -> IO CString )-       -> Matrix r1 c1 a-       -> SparseMatrix r2 c2 a-       -> Matrix r3 c3 a-withDS f m1 m2 = unsafePerformIO $ do-    m0 <- CM.new-    checkResult $ unsafeWith' m0 $ \v0 r0 c0 ->-        unsafeWith m1 $ \v1 r1 c1 ->-            unsafeWithS m2 $ \v2 inner outer r2 c2 s ->-                f (foreignType (undefined :: a))-                    v0 r0 c0-                    v1 r1 c1-                    v2 outer inner r2 c2 s-    C.unsafeFreeze m0-{-# INLINE withDS #-}--withSD :: forall r1 c1 r2 c2 r3 c3 a.-            (SingI r3, SingI c3, Numeric a)-       => ( CInt-         -> Ptr a -> CInt -> CInt-         -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt-         -> Ptr a -> CInt -> CInt-         -> IO CString )-       -> SparseMatrix r2 c2 a-       -> Matrix r1 c1 a-       -> Matrix r3 c3 a-withSD f m2 m1 = unsafePerformIO $ do-    m0 <- CM.new-    checkResult $ unsafeWith' m0 $ \v0 r0 c0 ->-        unsafeWith m1 $ \v1 r1 c1 ->-            unsafeWithS m2 $ \v2 inner outer r2 c2 s ->-                f (foreignType (undefined :: a))-                    v0 r0 c0-                    v2 outer inner r2 c2 s-                    v1 r1 c1-    C.unsafeFreeze m0-{-# INLINE withSD #-}--mkSparseMatrix :: forall r c a. (Storable a, SingI r, SingI c)-    => (Ptr (Ptr a) -> Ptr (Ptr CInt) -> Ptr CInt -> IO Int)-    -> IO (SparseMatrix r c a)-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)-        outer-  where-    c = fromIntegral $ fromSing (sing :: Sing c)-{-# INLINE mkSparseMatrix #-}--withSS :: forall r1 c1 r2 c2 r3 c3 a.-            (SingI r3, SingI c3, Numeric a)-       => ( 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 )-       -> SparseMatrix r1 c1 a-       -> SparseMatrix r2 c2 a-       -> SparseMatrix r3 c3 a-withSS f m1 m2 = unsafePerformIO $ mkSparseMatrix $ \v0 inner0 outer0 ->-    alloca $ \pn -> unsafeWithS m1 $ \v1 inner1 outer1 r1 c1 s1 ->-        unsafeWithS m2 $ \v2 inner2 outer2 r2 c2 s2 -> do-            checkResult $ f (foreignType (undefined :: a))-                v0 outer0 inner0 r c pn-                v1 outer1 inner1 r1 c1 s1-                v2 outer2 inner2 r2 c2 s2-            fromIntegral <$> peek pn-  where-    r = fromIntegral $ fromSing (sing :: Sing r3)-    c = fromIntegral $ fromSing (sing :: Sing c3)-{-# INLINE withSS #-}--checkResult :: IO CString -> IO ()-checkResult func = func >>= \c_str -> when (c_str /= nullPtr) $-    peekCString c_str >>= \str -> error str-{-# INLINE checkResult #-}------------------------------------------------------------------------------------ Raw pointers------------------------------------------------------------------------------------ | Pass a pointer to the matrix's data to the IO action.--- The data may not be modified through the pointer.-unsafeWith :: Storable a => Matrix n m a -> (Ptr a -> CInt -> CInt -> IO b) -> IO b-unsafeWith mat@(D.Matrix vec) f = VS.unsafeWith vec $ \p ->-    f p (fromIntegral r) $ fromIntegral c -  where-    (r,c) = C.dim mat-{-# INLINE unsafeWith #-}--unsafeWith' :: Storable a => MMatrix n m RealWorld a -> (Ptr a -> CInt -> CInt -> IO b) -> IO b-unsafeWith' mat@(DM.MMatrix vec) f = VSM.unsafeWith vec $ \p ->-    f p (fromIntegral r) $ fromIntegral c-  where-    (r,c) = CM.dim mat-{-# INLINE unsafeWith' #-}---- | Pass a pointer to the matrix's data to the IO action.--- The data may not be modified through the pointer.-unsafeWithS :: (Storable a, S.Zero a)-            => SparseMatrix n m a-            -> (Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt -> IO b)-            -> IO b-unsafeWithS mat@(S.SparseMatrix val inner outer) f = VS.unsafeWith val $ \pval ->-    VS.unsafeWith inner $ \pinner -> VS.unsafeWith outer $ \pouter ->-        f pval pinner pouter (fromIntegral r) (fromIntegral c) (fromIntegral $ VS.length val)-  where-    (r,c) = C.dim mat-{-# INLINE unsafeWithS #-}
src/Data/Matrix/Static/Sparse.hs view
@@ -35,6 +35,8 @@     -- * Construction     , C.empty     , fromTriplet+    , fromTripletC+    , toTriplet     , C.fromVector     , C.fromList     , C.unsafeFromVector@@ -54,10 +56,12 @@ 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 Conduit+import Data.Conduit.Internal (zipSinks)+import Data.Tuple (swap) import GHC.TypeLits (type (<=)) import Foreign.C.Types import Data.Complex@@ -77,6 +81,9 @@ instance Zero Float where     zero = 0.0 +instance Zero CFloat where+    zero = 0.0+ instance Zero Double where     zero = 0.0 @@ -101,28 +108,16 @@                                      -- non-zero in the previous two arrays.                  -> SparseMatrix r c v a +instance (G.Vector v a, Eq (v a)) => Eq (SparseMatrix r c v a) where+    (==) (SparseMatrix a b c) (SparseMatrix a' b' c') =+        a == a' && b == b' && c == c'+ instance (G.Vector v a, Zero a, Show a) => Show (SparseMatrix r c v a) where     show mat = printf "(%d x %d)\n%s" r c vals       where         (r,c) = C.dim mat         vals = unlines $ map (unwords . map show . G.toList) $ C.toRows mat -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)@@ -145,27 +140,18 @@     -- | O(1) Create matrix from vector containing columns.     unsafeFromVector :: forall r c. (G.Vector v a, SingI r, SingI c)            => v a -> SparseMatrix r c v a-    unsafeFromVector vec = SparseMatrix-        (G.generate n (G.unsafeIndex vec . S.unsafeIndex nz))-        inner outer+    unsafeFromVector vec = fromTriplet vec'       where-        inner = S.map fromIntegral $ S.map (`mod` c) nz-        outer = S.create $ do-            v <- SM.replicate (c+1) 0-            S.forM_ nz $ \x -> do-                let i = x `div` r-                SM.unsafeModify v succ (i+1)-            forM_ [1..c] $ \i -> do-                x <- SM.unsafeRead v (i-1)-                SM.unsafeModify v (+x) i-            return v-        nz = S.filter (\i -> vec `G.unsafeIndex` i /= zero) $ S.enumFromN 0 (r*c)-        n = S.length nz+        vec' = map (\((a,b),c) -> (a,b,c)) $ filter ((/=zero) . snd) $+            zipWith (\i x -> (toIndex i, x)) [0..] $ G.toList vec+        toIndex i = swap $ i `divMod` r         r = fromIntegral $ fromSing (sing :: Sing r)-        c = fromIntegral $ fromSing (sing :: Sing c)     {-# INLINE unsafeFromVector #-} -    transpose (SparseMatrix val inner outer) = undefined +    transpose mat = runIdentity $ fromTripletC source+      where+        source = toTriplet mat .| mapC (\(i,j,x) -> (j,i,x))+    {-# INLINE transpose #-}      thaw = undefined     {-# INLINE thaw #-}@@ -180,9 +166,15 @@     {-# INLINE unsafeFreeze #-}      map f (SparseMatrix vec inner outer) = SparseMatrix (G.map f vec) inner outer-    imap = undefined     {-# INLINE map #-} +    imap f mat@(SparseMatrix _ inner outer) = SparseMatrix vec' inner outer+      where+        vec' = runST $ runConduit $ toTriplet mat .| mapC g .| sinkVector+        g (i,j,x) = f (i,j) x+    {-# INLINE imap #-}+    + -- | 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)@@ -192,28 +184,71 @@     outer = S.scanl (+) 0 $ S.create $ do         vec <- SM.replicate c 0         _ <- flip mapM triplets $ \(_, j, _) -> -            SM.modify vec (+1) j+            SM.unsafeModify 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+            idx <- fromIntegral <$> SM.unsafeRead outer' j+            GM.unsafeWrite val' idx v+            SM.unsafeWrite inner' idx $ fromIntegral i+            SM.unsafeModify 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(n) Create matrix from triplet. row and column indices *are not* assumed to be ordered+-- duplicate entries are carried over to the CSR represention+fromTripletC :: forall m r c v a. (Monad m, G.Vector v a, SingI r, SingI c)+             => ConduitT () (Int, Int, a) m ()+             -> m (SparseMatrix r c v a)+fromTripletC triplets = do+    (nnz, outer) <- runConduit $ triplets .| zipSinks lengthC sinkOuter+    (val, inner, _) <- runConduit $ triplets .| sinkValInner nnz (clone outer)+    return $ SparseMatrix val inner outer+  where+    sinkOuter = S.scanl (+) 0 <$> foldlC f (S.replicate c 0)+      where+        f vec (_, j, _) = S.modify (\v -> SM.unsafeModify v (+1) j) vec+    sinkValInner nnz outer0 = foldlC f (val0, inner0, outer0)+      where+        val0 = G.create $ GM.new nnz+        inner0 = S.create $ SM.new nnz+        f (val, inner, outer) (i, j, v) = (val', inner', outer')+          where+            idx = fromIntegral $ outer `S.unsafeIndex` j+            val' = G.create $ do+                vec <- G.unsafeThaw val+                GM.unsafeWrite vec idx v+                return vec+            inner' = S.create $ do+                vec <- S.unsafeThaw inner+                SM.unsafeWrite vec idx $ fromIntegral i+                return vec+            outer' = S.create $ do+                vec <- S.unsafeThaw outer+                SM.unsafeModify vec (+1) j+                return vec+    c = fromIntegral $ fromSing (sing :: Sing c)+    clone x = S.create $ S.thaw x+{-# INLINE fromTripletC #-} +toTriplet :: (Monad m, G.Vector v a, SingI r, SingI c)+          => SparseMatrix r c v a -> ConduitT i (Int, Int, a) m ()+toTriplet (SparseMatrix val inner outer) =+    G.ifoldM_ go (fromIntegral $ G.head outer) outer+  where+    go start curC end = do+        enumFromToC start (end'-1) .| mapC f+        return end'+      where+        end' = fromIntegral end+        f i = (fromIntegral $ inner `G.unsafeIndex` i, fromIntegral curC - 1, val `G.unsafeIndex` i)+{-# INLINE toTriplet #-}+ -- | 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@@ -242,3 +277,12 @@         k = (u+l) `shiftR` 1         x' = vec `S.unsafeIndex` k {-# INLINE binarySearchByBounds #-}+++-------------------------------------------------------------------------------+-- Helper+-------------------------------------------------------------------------------++--getIndex :: Int -> (Int, Int)+--getIndex = +--{-# INLINE getIndex #-}
src/Data/Matrix/Static/Sparse/Mutable.hs view
@@ -16,6 +16,7 @@ import qualified Data.Vector.Generic.Mutable as GM import qualified Data.Vector.Storable as S import           Prelude                     hiding (read, replicate)+import Foreign.C.Types import Data.Singletons  import qualified Data.Matrix.Static.Generic.Mutable as C@@ -25,16 +26,18 @@     MSparseMatrix :: (SingI r, SingI c)                   => (v s a)         -- ^ Values: stores the coefficient values                                       -- of the non-zeros.-                  -> (S.Vector Int)  -- ^ InnerIndices: stores the row+                  -> (S.Vector CInt)  -- ^ InnerIndices: stores the row                                       -- (resp. column) indices of the non-zeros.-                  -> (S.Vector Int)  -- ^ OuterStarts: stores for each column+                  -> (S.Vector CInt)  -- ^ OuterStarts: stores for each column                                       -- (resp. row) the index of the first                                       -- non-zero in the previous two arrays.+                  -> (S.Vector CInt)  -- InnerNNZs: stores the number of non-zeros+                                      -- of each column (resp. row).                    -> MSparseMatrix r c v s a  instance GM.MVector v a => C.MMatrix MSparseMatrix v a where     dim :: forall r c s. MSparseMatrix r c v s a -> (Int, Int)-    dim (MSparseMatrix _ _ _) = (r,c)+    dim (MSparseMatrix _ _ _ _) = (r,c)       where         r = fromIntegral $ fromSing (sing :: Sing r)         c = fromIntegral $ fromSing (sing :: Sing c)
tests/Main.hs view
@@ -1,36 +1,13 @@-{-# 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--}+import Test.Base+import Test.LinearAlgebra  main :: IO () main = defaultMain $ testGroup "Main"-    [ testProperty "" square+    [ base+    , linearAlgebra     ]
+ tests/Test/Base.hs view
@@ -0,0 +1,61 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE IncoherentInstances #-}++module Test.Base (base) where++import Test.Tasty+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 hiding ((@@))+import Data.Singletons.Prelude (Min)+import Data.Vector (Vector)+import Control.Monad.IO.Class (liftIO)+import Test.Tasty.QuickCheck+import Conduit++import Test.Utils++base :: TestTree+base = testGroup "Base"+    [ pConversion+    , pTranspose+    ]++pConversion = testGroup "Conversion"+    [ testProperty "Dense: fromVector . flatten" t1+    , testProperty "Sparse: fromVector . flatten" t2+    , testProperty "Sparse -- fromTriplet . toTriplet" tTri+    , testProperty "Sparse -- fromTripletC . toTriplet" tTriC+    , testProperty "Sparse -- dense" t3 ]+  where+    t1 :: D.Matrix 80 60 Vector Int -> Bool+    t1 mat = (D.fromVector $ D.flatten mat) == mat+    t2 :: S.SparseMatrix 80 60 Vector Int -> Bool+    t2 mat = mat == S.fromVector (S.flatten mat)+    t3 :: D.Matrix 80 60 Vector Int -> Bool+    t3 mat = (D.fromVector $ S.flatten mat') == mat+      where+        mat' = S.fromVector $ D.flatten mat :: S.SparseMatrix 80 60 Vector Int+    tTri :: S.SparseMatrix 80 60 Vector Int -> Bool+    tTri mat = S.fromTriplet xs == mat+      where+        xs = runIdentity $ runConduit $ S.toTriplet mat .| sinkList+    tTriC :: S.SparseMatrix 80 60 Vector Int -> Bool+    tTriC mat = mat == runIdentity (S.fromTripletC (S.toTriplet mat))++pTranspose = testGroup "Transpose"+    [ testProperty "Dense" t1+    , testProperty "Sparse" tSp+    ]+  where+    t1 :: D.Matrix 80 40 Vector Int -> Bool+    t1 mat = D.transpose (D.transpose mat) == mat+    tSp :: S.SparseMatrix 80 40 Vector Int -> Bool+    tSp mat = G.transpose (G.transpose mat) == mat
+ tests/Test/LinearAlgebra.hs view
@@ -0,0 +1,83 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE IncoherentInstances #-}++module Test.LinearAlgebra (linearAlgebra) where++import Test.Tasty+import Data.Matrix.Static.LinearAlgebra+import qualified Data.Matrix.Static.Dense as D+import qualified Data.Matrix.Static.Generic as G+import qualified Data.Matrix.Static.Sparse as S+import Data.Singletons hiding ((@@))+import Data.Singletons.Prelude (Min)+import Data.Complex+import Data.Vector.Storable (Storable)+import GHC.TypeNats (KnownNat)+import Control.Monad.IO.Class (liftIO)+import Test.Tasty.QuickCheck++import Test.Utils++linearAlgebra :: TestTree+linearAlgebra = testGroup "Linear algebra"+    [ svdTest+    , eigenTest+    ]++svdTest = testGroup "SVD"+    [ testProperty "SVD (Float)" svd1+    , testProperty "SVD (Double)" svd2+    ]+  where+    svd1 :: Matrix 50 30 Float -> Bool+    svd1 m = m' ~= m+      where+        m' = u @@ S.diag d @@ D.transpose v+        (u,d,v) = svd m+    svd2 :: Matrix 50 30 Double -> Bool+    svd2 m = m' ~= m+      where+        m' = u @@ S.diag d @@ D.transpose v+        (u,d,v) = svd m++eigenTest = testGroup "Eigendecomposition"+    [ testProperty "Full" eigen1+    , testProperty "Partial dense" eigen2+    , testProperty "Partial symmetric dense" eigen3+    , testProperty "Partial symmetric sparse" eigen4+    ]+  where+    eigen1 :: Matrix 100 100 Double -> Bool+    eigen1 m = (m' @@ v) ~= (v @@ S.diag d)+      where+        m' = D.map (\x -> mkPolar x 0) m+        (d, v)= eig m+    eigen2 :: Matrix 10 10 Double -> Bool+    eigen2 m = m' @@ v ~= v @@ S.diag d+      where+        m' = D.map (\x -> mkPolar x 0) m+        (d, v)= eigS (sing :: Sing 8) m+    eigen3 :: Matrix 100 100 Double -> Bool+    eigen3 raw = m @@ v ~= v @@ S.diag d+      where+        m = raw %+% D.transpose raw+        (d, v)= eigSH (sing :: Sing 99) m+    eigen4 :: SparseMatrix 100 100 Double -> Bool+    eigen4 raw = m @@ v ~= v @@ S.diag d+      where+        m = raw %+% D.transpose raw+        (d, v)= eigSH (sing :: Sing 99) m++++{-+propTranspose :: Matrix 50 100 Double -> Bool+propTranspose m = D.transpose (D.transpose m) == m && +    D.convertAny (S.transpose $ S.transpose (D.convertAny m)) == m+-}
+ tests/Test/Utils.hs view
@@ -0,0 +1,65 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE IncoherentInstances #-}++module Test.Utils where++import Test.Tasty.QuickCheck+import Data.Complex+import Data.Matrix.Static.LinearAlgebra+import Data.Singletons+import qualified Data.Matrix.Static.Dense as D+import qualified Data.Matrix.Static.Sparse as S+import Data.Vector.Storable (Storable)+import Data.List.Ordered+import Data.Ord+import qualified Data.Vector.Generic as G+import Data.AEq+import Control.Monad++class Approx a where+    (~=) :: a -> a -> Bool+    infixl 4 ~=++instance Approx Double where+    a ~= b = abs (a - b) < 1e-10++instance Approx Float where+    a ~= b = abs (a - b) < 1e-3++instance Approx (Complex Double) where+    a ~= b = r1 ~= r2 && i1 ~= i2+      where+        (r1, i1) = polar a+        (r2, i2) = polar b++instance (SingI m, SingI n, Storable a, Approx a) => Approx (Matrix m n a) where+    m1 ~= m2 = D.all id $ D.zipWith (~=) m1 m2++instance (G.Vector v a, Arbitrary a, SingI m, SingI n)+    => Arbitrary (D.Matrix m n v a) where+        arbitrary = D.fromList <$> vector (m*n)+          where+            m = fromIntegral $ fromSing (sing :: Sing m)+            n = fromIntegral $ fromSing (sing :: Sing n)+        shrink _v = []++instance (G.Vector v a, Arbitrary a, SingI m, SingI n, S.Zero a)+    => Arbitrary (S.SparseMatrix m n v a) where+        arbitrary = do+            vals <- filter (/=S.zero) <$> vector p+            xs <- fmap (nubSortBy (comparing fst)) $ forM vals $ \v -> do+                i <- choose (0, m-1)+                j <- choose (0, n-1)+                return ((i,j),v)+            return $ S.fromTriplet $ map (\((a,b),c) -> (a,b,c)) xs+          where+            p = (m * n) `div` 10+            m = fromIntegral $ fromSing (sing :: Sing m)+            n = fromIntegral $ fromSing (sing :: Sing n)+        shrink _v = []