packages feed

matrix-sized 0.0.1 → 0.0.2

raw patch · 21 files changed

+1464/−1364 lines, 21 files

Files

README.md view
@@ -6,3 +6,30 @@ - The dimensions of matrices are statically typed.  - Provides bindings to high performance C++ linear algebra libraries such Eigen and Spectra.++Following GHC extensions may be needed:++- ScopedTypeVariables+- RankNTypes+- TypeFamilies+- DataKinds++Example+-------++```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++print mat+print mat'++print $ eigs (sing :: Sing 1) mat == eigs (sing :: Sing 1) mat'++print $ cholesky mat++print $ mat %*% mat %*% mat+print $ mat' %*% mat' %*% mat+```
cbits/eigen-basic.cpp view
@@ -16,7 +16,13 @@ }  template <class T>-Map< SparseMatrix<T> > smatrix(const void* val, const void* outer,+Map< SparseMatrix<T> > smatrix(void* val, void* outer,+    void* inner, int r, int c, int s) {+    return Map< SparseMatrix<T> >(r, c, s, (int*)outer, (int*)inner, (T*)val);+}++template <class T>+Map<SparseMatrix<T> > smatrix(const void* val, const void* outer,     const void* inner, int r, int c, int s) {     return Map< SparseMatrix<T> >(r, c, s, (int*)outer, (int*)inner, (T*)val); }@@ -64,22 +70,32 @@     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) {-    smatrix<T>(v,o,i,r,c,s) =-        (smatrix<T>(v1,o1,i1,r1,c1,s1) * smatrix<T>(v2,o2,i2,r2,c2,s2)).pruned();+    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).pruned();+    *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(T));+    *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));-*/   #define UNOP(name) \
matrix-sized.cabal view
@@ -4,10 +4,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: 32a85fe252d7c7ac4e14b5bdda92aa2b5b67b902bee14e17f98b64df94f58492+-- hash: 6910822a984416806867bac316c8f9685ade7b6a7672d4437c67fd4770977af1  name:           matrix-sized-version:        0.0.1+version:        0.0.2 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@@ -377,16 +377,16 @@  library   exposed-modules:-      Data.Matrix.Dense-      Data.Matrix.Dense.Mutable-      Data.Matrix.Sparse-      Data.Matrix.Sparse.Mutable-      Data.Matrix.LinearAlgebra-      Data.Matrix.LinearAlgebra.Types+      Data.Matrix.Static.Dense+      Data.Matrix.Static.Dense.Mutable+      Data.Matrix.Static.Sparse+      Data.Matrix.Static.Sparse.Mutable+      Data.Matrix.Static.LinearAlgebra+      Data.Matrix.Static.LinearAlgebra.Types+      Data.Matrix.Static.Generic+      Data.Matrix.Static.Generic.Mutable   other-modules:-      Data.Matrix.Internal.Class-      Data.Matrix.Internal.Class.Mutable-      Data.Matrix.Internal.LinearAlgebra+      Data.Matrix.Static.Internal   hs-source-dirs:       src   ghc-options: -Wall
− src/Data/Matrix/Dense.hs 
@@ -1,389 +0,0 @@-{-# LANGUAGE BangPatterns          #-}-{-# LANGUAGE DeriveGeneric         #-}-{-# LANGUAGE FlexibleContexts      #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE DataKinds #-}-{-# LANGUAGE RankNTypes #-}-{-# LANGUAGE GADTs #-}-{-# LANGUAGE InstanceSigs #-}-{-# LANGUAGE FlexibleInstances     #-}-{-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE TypeFamilies          #-}-{-# LANGUAGE TypeOperators #-}-module Data.Matrix.Dense-    (-    -- * Immutable Matrix-      Matrix(..)--    -- * Accessors-    -- ** length information-    , C.dim-    , C.rows-    , C.cols--    -- ** Query-    , (C.!)-    , C.takeRow-    , C.takeColumn-    , C.takeDiag--    -- ** Unsafe Query-    , C.unsafeIndex-    , C.unsafeTakeRow-    , C.unsafeTakeColumn--    -- * Construction-    , C.empty-    , C.matrix-    , C.fromVector-    , C.fromList-    , fromRows-    , fromColumns-    , C.unsafeFromVector--    , diag-    , diagRect--    -- * Conversions-    , C.flatten-    , C.toRows-    , C.toColumns-    , C.toList--    -- * Conversion between Different matrix types-    , convert-    , C.convertAny--    , transpose--    , C.mapM-    , C.imapM--    -- * Zipping-    , zip-    , zip3-    , zipWith-    , zipWith3-    , zipWith4-    , zipWith5-    , zipWith6-    , izipWith-    , izipWith3--    -- * Monadic Zipping-    , zipWithM-    , zipWithM_--    -- * Unzipping-    , unzip-    , unzip3--    , generate--    -- * Mutable matrix-    , C.thaw-    , C.unsafeThaw-    , C.freeze-    , C.unsafeFreeze-    , 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 GHC.TypeLits (type (<=))-import Data.Singletons-import Data.Tuple (swap)-import qualified Data.List as L-import Text.Printf (printf)--import           Data.Matrix.Dense.Mutable (MMatrix (..))-import qualified Data.Matrix.Dense.Mutable as DM-import qualified Data.Matrix.Internal.Class as C--type instance C.Mutable Matrix = MMatrix---- | Column-major matrix-data Matrix :: C.MatrixKind where-    Matrix :: (SingI r, SingI c) => v a -> Matrix r c v a--instance (G.Vector v a, Show a) => Show (Matrix r c v a) where-    show mat = printf "(%d x %d)\n%s" r c vals-      where-        (r,c) = C.dim mat-        vals = unlines $ L.map (unwords . L.map show . G.toList) $ C.toRows mat--instance (G.Vector v a, Eq (v a)) => Eq (Matrix r c v a) where-    (==) (Matrix v1) (Matrix v2) = v1 == v2--instance (SingI r, SingI c, G.Vector v a, Num a) =>-    Num (Matrix r c v a) where-        m1 + m2 = zipWith (+) m1 m2-        m1 - m2 = zipWith (-) m1 m2-        m1 * m2 = zipWith (*) m1 m2-        negate = C.map negate-        abs = C.map abs-        signum = undefined-        fromInteger = undefined--instance (SingI r, SingI c, G.Vector v a, Fractional a) =>-    Fractional (Matrix r c v a) where-        m1 / m2 = zipWith (/) m1 m2-        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)-    dim (Matrix _) = (r,c)-      where-        r = fromIntegral $ fromSing (sing :: Sing r)-        c = fromIntegral $ fromSing (sing :: Sing c)-    {-# INLINE dim #-}--    -- | O(1) Unsafe indexing without bound check.-    unsafeIndex mat@(Matrix vec) (i,j) = vec `G.unsafeIndex` idx-      where-        idx = i + j * fst (C.dim mat)-    {-# INLINE unsafeIndex #-}--    -- | O(1) Create matrix from vector.-    unsafeFromVector = Matrix-    {-# INLINE unsafeFromVector #-}--    -- | O(1) Extract a row.-    unsafeTakeColumn mat@(Matrix vec) i = G.slice (i*r) r vec-      where-        (r, _) = C.dim mat-    {-# INLINE unsafeTakeColumn #-}--    -- | Create a vector by concatenating columns.-    flatten (Matrix vec) = vec-    {-# INLINE flatten #-}--    thaw (Matrix v) = MMatrix <$> G.thaw v-    {-# INLINE thaw #-}--    unsafeThaw (Matrix v) = MMatrix <$> G.unsafeThaw v-    {-# INLINE unsafeThaw #-}--    freeze (MMatrix v) = Matrix <$> G.freeze v-    {-# INLINE freeze #-}--    unsafeFreeze (MMatrix v) = Matrix <$> G.unsafeFreeze v-    {-# INLINE unsafeFreeze #-}--    map f (Matrix vec) = Matrix $ G.map f vec-    {-# INLINE map #-}--    imap f m@(Matrix vec) = Matrix $ G.imap g vec-      where-        g i = f (toIndex (C.rows m) i)-    {-# INLINE imap #-}--    imapM_ f m@(Matrix vec) = G.imapM_ g vec-      where-        g i = f (toIndex (C.rows m) i)-    {-# INLINE imapM_ #-}--    sequence (Matrix vec) = Matrix <$> G.sequence vec-    {-# INLINE sequence #-}--    sequence_ (Matrix vec) = G.sequence_ vec-    {-# INLINE sequence_ #-}-----reshape :: G.Vector v a => Matrix v a -> (Int, Int) -> Matrix v a---- | O(m*n) Create matrix from rows-fromRows :: (G.Vector v a, SingI r, SingI c) => [v a] -> Matrix r c v a-fromRows = transpose . C.unsafeFromVector . G.concat-{-# INLINE fromRows #-}---- | O(m*n) Create matrix from columns-fromColumns :: (G.Vector v a, SingI r, SingI c)-            => [v a] -> Matrix r c v a-fromColumns = C.fromVector . G.concat-{-# INLINE fromColumns #-}---- | O(m*n) Matrix transpose-transpose :: G.Vector v a => Matrix m n v a -> Matrix n m v a-transpose mat@(Matrix vec) = C.unsafeFromVector $ G.generate (r*c) f-  where-    (r, c) = C.dim mat-    f i = vec G.! (i `mod` c * r + i `div` c)-{-# INLINE transpose #-}--{---- | 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 rectangular matrix with default values and given diagonal-diagRect :: (G.Vector v a, SingI r, SingI c, n <= r, n <= c)-         => a                    -- ^ default value-         -> Matrix n 1 v a       -- ^ diagonal-         -> Matrix r c v a-diagRect z0 d = C.create $ do-    mat <- DM.replicate z0-    C.imapM_ (DM.unsafeWrite mat) d-    return mat-{-# INLINE diagRect #-}--zipWith :: (G.Vector v a, G.Vector v b, G.Vector v c-           , SingI n, SingI m )-        => (a -> b -> c)-        -> Matrix n m v a -> Matrix n m v b -> Matrix n m v c-zipWith f m1 m2 = C.unsafeFromVector $ G.zipWith f (C.flatten m1) $ C.flatten m2-{-# INLINE zipWith #-}--zipWith3 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d-            , SingI n, SingI m )-         => (a -> b -> c -> d)-         -> Matrix n m v a -> Matrix n m v b -> Matrix n m v c-         -> Matrix n m v d-zipWith3 f m1 m2 m3 = C.unsafeFromVector $-    G.zipWith3 f (C.flatten m1) (C.flatten m2) $ C.flatten m3-{-# INLINE zipWith3 #-}--zipWith4 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d, G.Vector v e-            , SingI n, SingI m )-         => (a -> b -> c -> d -> e)-         -> Matrix n m v a-         -> Matrix n m v b-         -> Matrix n m v c-         -> Matrix n m v d-         -> Matrix n m v e-zipWith4 f m1 m2 m3 m4 = C.unsafeFromVector $-    G.zipWith4 f (C.flatten m1) (C.flatten m2) (C.flatten m3) $ C.flatten m4-{-# INLINE zipWith4 #-}--zipWith5 :: ( G.Vector v a, G.Vector v b, G.Vector v c,G.Vector v d-            , G.Vector v e, G.Vector v f-            , SingI n, SingI m )-         => (a -> b -> c -> d -> e -> f)-         -> Matrix n m v a-         -> Matrix n m v b-         -> Matrix n m v c-         -> Matrix n m v d-         -> Matrix n m v e-         -> Matrix n m v f-zipWith5 f m1 m2 m3 m4 m5 = C.unsafeFromVector $-    G.zipWith5 f (C.flatten m1) (C.flatten m2)-    (C.flatten m3) (C.flatten m4) $ C.flatten m5-{-# INLINE zipWith5 #-}--zipWith6 :: ( G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d-            , G.Vector v e, G.Vector v f, G.Vector v g-            , SingI n, SingI m )-         => (a -> b -> c -> d -> e -> f -> g)-         -> Matrix n m v a-         -> Matrix n m v b-         -> Matrix n m v c-         -> Matrix n m v d-         -> Matrix n m v e-         -> Matrix n m v f-         -> Matrix n m v g-zipWith6 f m1 m2 m3 m4 m5 m6 = C.unsafeFromVector $-    G.zipWith6 f (C.flatten m1) (C.flatten m2) (C.flatten m3)-    (C.flatten m4) (C.flatten m5) $ C.flatten m6-{-# INLINE zipWith6 #-}--izipWith :: ( G.Vector v a, G.Vector v b, G.Vector v c-            , SingI n, SingI m )-         => ((Int, Int) -> a -> b -> c)-         -> Matrix n m v a -> Matrix n m v b -> Matrix n m v c-izipWith f m1 m2 = C.unsafeFromVector $-    G.izipWith g (C.flatten m1) $ C.flatten m2-  where-    g i = f (toIndex (C.rows m1) i)-{-# INLINE izipWith #-}--izipWith3 :: ( G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d-             , SingI n, SingI m )-          => ((Int, Int) -> a -> b -> c -> d)-          -> Matrix n m v a -> Matrix n m v b -> Matrix n m v c-          -> Matrix n m v d-izipWith3 f m1 m2 m3 = C.unsafeFromVector $ G.izipWith3 g-    (C.flatten m1) (C.flatten m2) $ C.flatten m3-  where-    g i = f (toIndex (C.rows m1) i)-{-# INLINE izipWith3 #-}--zip :: (SingI n, SingI m, G.Vector v a, G.Vector v b, G.Vector v (a,b))-    => Matrix n m v a -> Matrix n m v b -> Matrix n m v (a,b)-zip m1 m2 = C.unsafeFromVector $ G.zip (C.flatten m1) $ C.flatten m2-{-# INLINE zip #-}--zip3 :: (SingI n, SingI m, G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v (a,b,c))-     => Matrix n m v a-     -> Matrix n m v b-     -> Matrix n m v c-     -> Matrix n m v (a,b,c)-zip3 m1 m2 m3 = C.unsafeFromVector $-    G.zip3 (C.flatten m1) (C.flatten m2) $ C.flatten m3-{-# INLINE zip3 #-}--zipWithM :: ( G.Vector v a, G.Vector v b, G.Vector v c-            , Monad monad, SingI n, SingI m )-         => (a -> b -> monad c)-         -> Matrix n m v a -> Matrix n m v b -> monad (Matrix n m v c)-zipWithM f m1 m2 = liftM C.unsafeFromVector $-    G.zipWithM f (C.flatten m1) $ C.flatten m2-{-# INLINE zipWithM #-}--zipWithM_ :: (G.Vector v a, G.Vector v b, G.Vector v c, Monad monad)-          => (a -> b -> monad c)-          -> Matrix n m v a -> Matrix n m v b -> monad ()-zipWithM_ f m1 m2 = G.zipWithM_ f (C.flatten m1) $ C.flatten m2-{-# INLINE zipWithM_ #-}--unzip :: ( G.Vector v a, G.Vector v b, G.Vector v (a,b)-         , SingI n, SingI m )-      => Matrix n m v (a,b) -> (Matrix n m v a, Matrix n m v b )-unzip m = (C.unsafeFromVector v1, C.unsafeFromVector v2)-  where-    (v1, v2) = G.unzip $ C.flatten m-{-# INLINE unzip #-}--unzip3 :: ( G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v (a,b,c)-          , SingI n, SingI m )-       => Matrix n m v (a, b, c)-       -> (Matrix n m v a, Matrix n m v b, Matrix n m v c)-unzip3 m = (C.unsafeFromVector v1, C.unsafeFromVector v2, C.unsafeFromVector v3)-  where-    (v1, v2, v3) = G.unzip3 $ C.flatten m-{-# INLINE unzip3 #-}--generate :: forall r c v a. (G.Vector v a, SingI r, SingI c)-         => ((Int, Int) -> a) -> Matrix r c v a-generate f = C.unsafeFromVector . G.generate (r*c) $ \i -> f (i `divMod` r)-  where-    r = fromIntegral $ fromSing (sing :: Sing r)-    c = fromIntegral $ fromSing (sing :: Sing c)-{-# INLINE generate #-}---- | O(m*n) Convert different matrix type-convert :: (G.Vector v a, G.Vector w a) => Matrix r c v a -> Matrix r c w a-convert (Matrix vec) = Matrix $ G.convert vec-{-# INLINE convert #-}----- Helper-toIndex :: Int -> Int -> (Int, Int)-toIndex r i = swap $ i `divMod` r-{-# INLINE toIndex #-}
− src/Data/Matrix/Dense/Mutable.hs 
@@ -1,77 +0,0 @@-{-# LANGUAGE FlexibleInstances     #-}-{-# LANGUAGE FlexibleContexts      #-}-{-# LANGUAGE DataKinds #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE Rank2Types #-}-{-# LANGUAGE KindSignatures #-}-{-# LANGUAGE InstanceSigs #-}-{-# LANGUAGE GADTs #-}-{-# LANGUAGE MultiParamTypeClasses #-}-module Data.Matrix.Dense.Mutable-   ( -- * Mutable Matrix-     MMatrix(..)-   , C.dim-   , C.unsafeWrite-   , C.unsafeRead-   , C.new-   , C.replicate-   ) where--import           Control.DeepSeq-import qualified Data.Vector.Generic.Mutable as GM-import           Prelude                     hiding (read, replicate)-import Data.Singletons-import Control.Monad.Primitive     (PrimMonad, PrimState)--import qualified Data.Matrix.Internal.Class.Mutable as C---- | 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)-    dim (MMatrix _) = (r,c)-      where-        r = fromIntegral $ fromSing (sing :: Sing r)-        c = fromIntegral $ fromSing (sing :: Sing c)-    {-# INLINE dim #-}--    unsafeRead mat@(MMatrix v) (i,j) = GM.unsafeRead v idx-      where-        (r, _) = C.dim mat-        idx = i + j * r-    {-# INLINE unsafeRead #-}--    unsafeWrite mat@(MMatrix v) (i,j) = GM.unsafeWrite v idx-      where -        (r, _) = C.dim mat-        idx = i + j * r-    {-# INLINE unsafeWrite #-}--    new :: forall r c s. (SingI r, SingI c, PrimMonad s)-        => s (MMatrix r c v (PrimState s) a)-    new = MMatrix <$> GM.new (r*c)-      where-        r = fromIntegral $ fromSing (sing :: Sing r)-        c = fromIntegral $ fromSing (sing :: Sing c)-    {-# INLINE new #-}--    replicate :: forall r c s. (SingI r, SingI c, PrimMonad s)-              => a -> s (MMatrix r c v (PrimState s) a)-    replicate x = MMatrix <$> GM.replicate (r*c) x-      where-        r = fromIntegral $ fromSing (sing :: Sing r)-        c = fromIntegral $ fromSing (sing :: Sing c)-    {-# INLINE replicate #-}--{--takeColumn :: GM.MVector v a => MMatrix v m a -> Int -> v m a-takeColumn (MMatrix _ c tda offset vec) i = GM.slice i' c vec-  where-    i' = offset + i * tda-{-# INLINE takeColumn #-}--}
− src/Data/Matrix/Internal/Class.hs 
@@ -1,212 +0,0 @@-{-# LANGUAGE FlexibleContexts      #-}-{-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE DataKinds #-}-{-# LANGUAGE TypeFamilyDependencies #-}-{-# LANGUAGE Rank2Types            #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TypeFamilies          #-}-{-# LANGUAGE TypeOperators #-}-module Data.Matrix.Internal.Class-    ( Mutable-    , Matrix(..)-    , MatrixKind--    -- * Derived mothods-    , rows-    , cols-    , (!)-    , takeColumn-    , takeRow-    , toRows-    , toColumns-    , empty-    , matrix-    , fromVector-    , fromList-    , toList-    , create-    , convertAny-    , mapM-    , imapM-    ) where--import           Control.Monad.Primitive     (PrimMonad, PrimState)-import           Control.Monad.ST            (ST, runST)-import qualified Data.Vector.Generic         as G-import Text.Printf (printf)-import Prelude hiding (map, mapM, mapM_, sequence, sequence_)-import qualified Data.List as L-import Data.Kind (Type)-import GHC.TypeLits (Nat, type (<=))-import Data.Singletons (SingI, Sing, fromSing, sing)--import Data.Matrix.Internal.Class.Mutable (MMatrix, MMatrixKind)--type MatrixKind = Nat -> Nat -> (Type -> Type) -> Type -> Type--type family Mutable (mat :: MatrixKind) = (mmat :: MMatrixKind) | mmat -> mat--class (MMatrix (Mutable mat) (G.Mutable v) a, G.Vector v a) => Matrix (mat :: MatrixKind) v a where-    dim :: mat r c v a -> (Int, Int)--    unsafeIndex :: mat r c v a -> (Int, Int) -> a--    unsafeFromVector :: (SingI r, SingI c) => v a -> mat r c v a--    -- | 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)-      where-        (r,c) = dim mat-    {-# INLINE flatten #-}--    -- | Extract a row. Default algorithm is O(n * O(unsafeIndex)).-    unsafeTakeRow :: mat r c v a -> Int -> v a-    unsafeTakeRow mat i = G.generate c $ \j -> unsafeIndex mat (i,j)-      where-        (_,c) = dim mat-    {-# INLINE unsafeTakeRow #-}--    -- | Extract a column. Default algorithm is O(m * O(unsafeIndex)).-    unsafeTakeColumn :: mat r c v a -> Int -> v a-    unsafeTakeColumn mat j = G.generate r $ \i -> unsafeIndex mat (i,j)-      where-        (r,_) = dim mat-    {-# INLINE unsafeTakeColumn #-}--    -- | Extract the diagonal. Default algorithm is O(min(m,n) * O(unsafeIndex)).-    takeDiag :: mat r c v a -> v a-    takeDiag mat = G.generate n $ \i -> unsafeIndex mat (i,i)-      where-        n = uncurry min . dim $ mat-    {-# INLINE takeDiag #-}--    thaw :: PrimMonad s-         => mat r c v a-         -> s ((Mutable mat) r c (G.Mutable v) (PrimState s) a)--    unsafeThaw :: PrimMonad s-               => mat r c v a-               -> s ((Mutable mat) r c (G.Mutable v) (PrimState s) a)--    freeze :: PrimMonad s-           => (Mutable mat) r c (G.Mutable v) (PrimState s) a-           -> s (mat r c v a)--    unsafeFreeze :: PrimMonad s-                 => (Mutable mat) r c (G.Mutable v) (PrimState s) a-                 -> s (mat r c v a)--    map :: G.Vector v b => (a -> b) -> mat r c v a -> mat r c v b-    imap :: G.Vector v b => ((Int, Int) -> a -> b) -> mat r c v a -> mat r c v b-    imapM_ :: (Monad monad, Matrix mat v a)-           => ((Int, Int) -> a -> monad b) -> mat r c v a -> monad ()-    sequence :: (G.Vector v (monad a), Monad monad)-             => mat r c v (monad a) -> monad (mat r c v a)-    sequence_ :: (G.Vector v (monad a), Monad monad) => mat r c v (monad a) -> monad ()---- | Derived methods---- | Return the number of rows-rows :: Matrix m v a => m r c v a -> Int-rows = fst . dim-{-# INLINE rows #-}---- | Return the number of columns-cols :: Matrix m v a => m r c v a -> Int-cols = snd . dim-{-# INLINE cols #-}---- | Indexing-(!) :: forall m r c v a i j. (Matrix m v a, i <= r, j <= c)-    => m r c v a -> (Sing i, Sing j) -> a-(!) m (si, sj) = unsafeIndex m (i,j)-  where-    i = fromIntegral $ fromSing si-    j = fromIntegral $ fromSing sj-{-# INLINE (!) #-}---- | Construct matrix from a vector containg columns.-fromVector :: forall m r c v a. (SingI r, SingI c, Matrix m v a)-           => v a -> m r c v a-fromVector vec | r*c /= n = error errMsg-               | otherwise = unsafeFromVector vec-  where-    errMsg = printf "fromVector: incorrect length (%d * %d != %d)" r c n-    n = G.length vec-    r = fromIntegral $ fromSing (sing :: Sing r)-    c = fromIntegral $ fromSing (sing :: Sing c)-{-# INLINE fromVector #-}--matrix :: (SingI r, SingI c, Matrix m v a)-       => [[a]] -> m r c v a-matrix = fromList . concat . L.transpose-{-# INLINE matrix #-}---- | Construct matrix from a list containg columns.-fromList :: (SingI r, SingI c, Matrix m v a)-         => [a] -> m r c v a-fromList = fromVector . G.fromList-{-# INLINE fromList #-}---- | O(m*n) Create a list by concatenating columns-toList :: Matrix m v a => m r c v a -> [a]-toList = G.toList . flatten-{-# INLINE toList #-}--empty :: Matrix m v a => m 0 0 v a-empty = unsafeFromVector G.empty-{-# INLINE empty #-}--create :: Matrix m v a-       => (forall s . ST s ((Mutable m) r c (G.Mutable v) s a)) -> m r c v a-create m = runST $ unsafeFreeze =<< m-{-# INLINE create #-}---- | O(m*n) Convert to any type of matrix.-convertAny :: (Matrix m1 v1 a, Matrix m2 v2 a, SingI r, SingI c)-           => m1 r c v1 a -> m2 r c v2 a-convertAny = unsafeFromVector . G.convert . flatten-{-# INLINE convertAny #-}---- | Extract a row.-takeRow :: forall m r c v a i. (i <= r, SingI i, Matrix m v a)-        => m r c v a -> Sing i -> v a-takeRow mat _ = unsafeTakeRow mat i-  where-    i = fromIntegral $ fromSing (sing :: Sing i)-{-# INLINE takeRow #-}---- | O(m) Return the rows-toRows :: Matrix m v a => m r c v a -> [v a]-toRows mat = L.map (unsafeTakeRow mat) [0..r-1]-  where-    (r,_) = dim mat-{-# INLINE toRows #-}---- | Extract a row.-takeColumn :: forall m r c v a j. (j <= c, SingI j, Matrix m v a)-           => m r c v a -> Sing j -> v a-takeColumn mat _ = unsafeTakeColumn mat j-  where-    j = fromIntegral $ fromSing (sing :: Sing j)-{-# INLINE takeColumn #-}---- | O(m*n) Return the columns-toColumns :: Matrix m v a => m r c v a -> [v a]-toColumns mat = L.map (unsafeTakeColumn mat) [0..c-1]-  where-    (_,c) = dim mat-{-# INLINE toColumns #-}--mapM :: (G.Vector v (monad b), Monad monad, Matrix mat v a, Matrix mat v b)-     => (a -> monad b) -> mat r c v a -> monad (mat r c v b)-mapM f = sequence . map f-{-# INLINE mapM #-}--imapM :: (G.Vector v (monad b), Monad monad, Matrix mat v a, Matrix mat v b)-      => ((Int, Int) -> a -> monad b)-      -> mat r c v a -> monad (mat r c v b)-imapM f = sequence . imap f-{-# INLINE imapM #-}
− src/Data/Matrix/Internal/Class/Mutable.hs 
@@ -1,51 +0,0 @@-{-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE DataKinds #-}-{-# LANGUAGE TypeFamilies          #-}--module Data.Matrix.Internal.Class.Mutable-    ( MMatrix(..)-    , MMatrixKind-    ) where--import           Control.Monad.Primitive     (PrimMonad, PrimState)-import qualified Data.Vector.Generic.Mutable as GM-import           Prelude                     hiding (read)-import Data.Kind (Type)-import GHC.TypeLits (Nat)-import Data.Singletons (SingI)--type MMatrixKind = Nat -> Nat -> (Type -> Type -> Type) -> Type -> Type -> Type--class GM.MVector v a => MMatrix (mat :: MMatrixKind) v a where-    dim :: mat r c v s a -> (Int, Int)--    unsafeRead :: PrimMonad s => mat r c v (PrimState s) a -> (Int, Int) -> s a--    unsafeWrite :: PrimMonad s => mat r c v (PrimState s) a -> (Int, Int) -> a -> s ()--    -- | Create a mutable matrix without initialization-    new :: (SingI r, SingI c, PrimMonad s) => s (mat r c v (PrimState s) a)--    replicate :: (SingI r, SingI c, PrimMonad s) => a -> s (mat r c v (PrimState s) a)--    {-# MINIMAL dim, unsafeRead, unsafeWrite, new, replicate #-}--{--write :: (PrimMonad s, MMatrix m v a)-      => m v (PrimState s) a -> (Int, Int) -> a -> s ()-write mat (i,j)-    | i < 0 || i >= r || j < 0 || j >= c = error "write: Index out of bounds"-    | otherwise = unsafeWrite mat (i,j)-  where-    (r,c) = dim mat-{-# INLINE write #-}--read :: (PrimMonad s, MMatrix m v a)-     => m v (PrimState s) a -> (Int, Int) -> s a-read mat (i,j)-    | i <0 || i >= r || j < 0 || j >= c = error "read: Index out of bounds"-    | otherwise = unsafeRead mat (i,j)-  where-    (r,c) = dim mat-{-# INLINE read #-}--}
− src/Data/Matrix/Internal/LinearAlgebra.hs 
@@ -1,59 +0,0 @@-{-# LANGUAGE ForeignFunctionInterface #-}-module Data.Matrix.Internal.LinearAlgebra-    ( c_dd_mul-    , c_ds_mul-    , c_sd_mul-    , c_inverse-    , c_cholesky-    , c_eigs-    , c_seigs-    ) where--import Data.Complex (Complex)-import Foreign-import Foreign.C.Types-import Foreign.C.String------------------------------------------------------------------------------------ Arithmetic---------------------------------------------------------------------------------foreign import ccall "eigen_dd_mul"-    c_dd_mul :: CInt-          -> Ptr a -> CInt -> CInt-          -> Ptr a -> CInt -> CInt-          -> Ptr a -> CInt -> CInt-          -> IO CString--foreign import ccall "eigen_ds_mul"-    c_ds_mul :: CInt-          -> Ptr a -> CInt -> CInt-          -> Ptr a -> CInt -> CInt-          -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt-          -> IO CString--foreign import ccall "eigen_sd_mul"-    c_sd_mul :: CInt-          -> Ptr a -> CInt -> CInt-          -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt-          -> Ptr a -> CInt -> CInt-          -> IO CString--             -foreign import ccall "eigen_inverse"-    c_inverse :: CInt-              -> Ptr a -> CInt -> CInt-              -> Ptr a -> CInt -> CInt-              -> IO CString--foreign import ccall "eigen_cholesky"-    c_cholesky :: CInt-               -> Ptr a -> Ptr a ->  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_seigs"-    c_seigs :: CInt -> Ptr (Complex Double) -> Ptr (Complex Double)-            -> Ptr Double -> Ptr CInt -> Ptr CInt-            -> CInt -> CInt -> IO CString
− src/Data/Matrix/LinearAlgebra.hs 
@@ -1,99 +0,0 @@-{-# LANGUAGE FlexibleInstances     #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE RankNTypes #-}-{-# LANGUAGE DataKinds #-}-{-# LANGUAGE KindSignatures #-}-{-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE FunctionalDependencies #-}-{-# LANGUAGE GADTs #-}-{-# LANGUAGE ExplicitNamespaces #-}-{-# LANGUAGE TypeOperators #-}-module Data.Matrix.LinearAlgebra-    ( Arithmetic(..)-    , Factorization(..)-    , inverse-    , module Data.Matrix.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 qualified Data.Matrix.Dense as D-import qualified Data.Matrix.Sparse as S-import qualified Data.Matrix.Internal.Class.Mutable as CM-import qualified Data.Matrix.Internal.Class as C-import qualified Data.Matrix.Internal.LinearAlgebra as Internal--import Data.Matrix.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-          -> 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 S.SparseMatrix D.Matrix where-    (%*%) = withDS Internal.c_ds_mul--instance Arithmetic S.SparseMatrix D.Matrix D.Matrix where-    (%*%) = withSD Internal.c_sd_mul--instance Arithmetic S.SparseMatrix S.SparseMatrix S.SparseMatrix where-    (%*%) = undefined--inverse :: (SingI n, Numeric a) => Matrix n n a -> Matrix n n a-inverse = withFun1 Internal.c_inverse--class Factorization mat where-    -- | Eigenvalues (not ordered) and-    -- eigenvectors (as columns) of a general square matrix.-    eigs :: (SingI k, SingI n, k <= n - 2)-         => Sing k-         -> mat n n VS.Vector Double-         -> (Matrix k 1 (Complex Double), Matrix n k (Complex Double))--    -- | Cholesky decomposition-    cholesky :: (Numeric a, SingI n) => mat n n VS.Vector a -> mat n n VS.Vector a--instance Factorization D.Matrix where-    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')-      where-        k = fromIntegral $ fromSing s-    {-# INLINE eigs #-}--    cholesky mat = flip withFun1 mat $-        \code p1 c1 _ p2 _ _ -> Internal.c_cholesky code p1 p2 c1-    {-# INLINE cholesky #-}--instance Factorization S.SparseMatrix where-    eigs 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_seigs 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 eigs #-}--    cholesky = undefined
− src/Data/Matrix/LinearAlgebra/Types.hs 
@@ -1,188 +0,0 @@-{-# LANGUAGE FlexibleInstances     #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE RankNTypes #-}-{-# LANGUAGE DataKinds #-}-module Data.Matrix.LinearAlgebra.Types-    ( Numeric(..)-    , Matrix-    , MMatrix-    , SparseMatrix-    , withFun1-    , withFun2-    , withDS-    , withSD-    , 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 Data.Int--import qualified Data.Matrix.Dense as D-import qualified Data.Matrix.Dense.Mutable as DM-import qualified Data.Matrix.Sparse as S-import qualified Data.Matrix.Internal.Class.Mutable as CM-import qualified Data.Matrix.Internal.Class as C-import qualified Data.Matrix.Internal.LinearAlgebra as Internal--class (S.Zero a, Storable a) => Numeric a where-    foreignType :: a -> CInt--instance Numeric Float where foreignType _ = 0-instance Numeric Double where foreignType _ =1-instance Numeric (Complex Float) where foreignType _ = 2-instance Numeric (Complex Double) where foreignType _ = 3--type Matrix r c a = D.Matrix r c Vector a-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 #-}--{--withSS :: forall r1 c1 r2 c2 r3 c3 a.-            (SingI r3, SingI c3, Numeric a)-       => ( CInt-         -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> 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 m2 m1 = unsafePerformIO $ do-    m0 <- CM.new-    _ <- 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 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/Sparse.hs 
@@ -1,196 +0,0 @@-{-# LANGUAGE FlexibleInstances     #-}-{-# LANGUAGE FlexibleContexts      #-}-{-# LANGUAGE DataKinds #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE Rank2Types #-}-{-# LANGUAGE KindSignatures #-}-{-# LANGUAGE InstanceSigs #-}-{-# LANGUAGE GADTs #-}-{-# LANGUAGE TypeOperators #-}-{-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE BangPatterns #-}-{-# LANGUAGE TypeFamilies #-}-module Data.Matrix.Sparse-   ( -- * Sparse matrix-     SparseMatrix(..)--   , Zero(..)--    -- * Accessors-    -- ** length information-    , C.dim-    , C.rows-    , C.cols--    -- ** Query-    , (C.!)-    , C.takeDiag--    -- ** Unsafe Query-    , C.unsafeIndex-    , C.unsafeTakeRow-    , C.unsafeTakeColumn--    -- * Construction-    , C.empty-    , C.fromVector-    , C.fromList-    , C.unsafeFromVector--    , diag-    , diagRect--    -- * Conversions-    , C.flatten-    , C.toList--    -- * Different matrix types-    , C.convertAny-   ) where--import           Control.DeepSeq-import qualified Data.Vector.Generic as G-import qualified Data.Vector.Storable as S-import qualified Data.Vector.Storable.Mutable as SM-import Data.Singletons-import Control.Monad-import           Data.Bits                         (shiftR)-import Text.Printf (printf)-import GHC.TypeLits (type (<=))-import Foreign.C.Types-import Data.Complex--import qualified Data.Matrix.Dense as D-import qualified Data.Matrix.Internal.Class as C-import Data.Matrix.Sparse.Mutable--type instance C.Mutable SparseMatrix = MSparseMatrix--class Eq a => Zero a where-    zero :: a--instance Zero Int where-    zero = 0--instance Zero Float where-    zero = 0.0--instance Zero Double where-    zero = 0.0--instance Zero (Complex Float) where-    zero = 0--instance Zero (Complex Double) where-    zero = 0--instance Eq a => Zero ([] a) where-    zero = []---- | Column-major mutable matrix.-data SparseMatrix :: C.MatrixKind where-    SparseMatrix :: (SingI r, SingI c)-                 => !(v a)           -- ^ Values: stores the coefficient values-                                     -- of the non-zeros.-                 -> !(S.Vector CInt)  -- ^ InnerIndices: stores the row-                                     -- (resp. column) indices of the non-zeros.-                 -> !(S.Vector CInt)  -- ^ OuterStarts: stores for each column-                                     -- (resp. row) the index of the first-                                     -- non-zero in the previous two arrays.-                 -> SparseMatrix r c v a--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 (NFData (v a)) => NFData (SparseMatrix r c v a) where-    rnf (SparseMatrix vec inner outer) = rnf vec--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)-    dim (SparseMatrix _ _ _) = (r,c)-      where-        r = fromIntegral $ fromSing (sing :: Sing r)-        c = fromIntegral $ fromSing (sing :: Sing c)-    {-# INLINE dim #-}--    -- | O(1) Unsafe indexing without bound check.-    unsafeIndex (SparseMatrix vec inner outer) (i,j) = -        case binarySearchByBounds inner (fromIntegral i) r0 r1 of-            Nothing -> zero-            Just k -> vec `G.unsafeIndex` k-      where-        r0 = fromIntegral $ outer `S.unsafeIndex` j-        r1 = fromIntegral $ outer `S.unsafeIndex` (j+1) - 1-    {-# INLINE unsafeIndex #-}--    -- | 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-      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-        r = fromIntegral $ fromSing (sing :: Sing r)-        c = fromIntegral $ fromSing (sing :: Sing c)-    {-# INLINE unsafeFromVector #-}--    thaw = undefined-    {-# INLINE thaw #-}--    unsafeThaw = undefined-    {-# INLINE unsafeThaw #-}--    freeze = undefined-    {-# INLINE freeze #-}--    unsafeFreeze = undefined-    {-# INLINE unsafeFreeze #-}--    map f (SparseMatrix vec inner outer) = SparseMatrix (G.map f vec) inner outer-    imap = undefined-    {-# INLINE map #-}---- | O(m*n) Create a square matrix with given diagonal.-diag :: (G.Vector v a, Zero a, SingI n)-     => D.Matrix n 1 v a       -- ^ diagonal-     -> SparseMatrix n n v a-diag = diagRect-{-# INLINE diag #-}---- | O(m*n) Create a rectangular matrix with default values and given diagonal-diagRect :: (G.Vector v a, Zero a, SingI r, SingI c, n <= r, n <= c)-         => D.Matrix n 1 v a       -- ^ diagonal-         -> SparseMatrix r c v a-diagRect d = SparseMatrix (C.flatten d) (S.enumFromN 0 n) (S.enumFromN 0 $ n + 1)-  where-    n = C.rows d-{-# INLINE diagRect #-}--binarySearchByBounds :: S.Vector CInt -> CInt -> Int -> Int -> Maybe Int-binarySearchByBounds vec x = loop-  where-    loop !l !u-        | l > u = Nothing-        | x == x' = Just k-        | x < x' = loop l (k-1)-        | otherwise = loop (k+1) u-      where-        k = (u+l) `shiftR` 1-        x' = vec `S.unsafeIndex` k-{-# INLINE binarySearchByBounds #-}
− src/Data/Matrix/Sparse/Mutable.hs 
@@ -1,75 +0,0 @@-{-# LANGUAGE FlexibleInstances     #-}-{-# LANGUAGE FlexibleContexts      #-}-{-# LANGUAGE DataKinds #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE Rank2Types #-}-{-# LANGUAGE KindSignatures #-}-{-# LANGUAGE InstanceSigs #-}-{-# LANGUAGE GADTs #-}-{-# LANGUAGE MultiParamTypeClasses #-}-module Data.Matrix.Sparse.Mutable-   ( -- * Mutable sparse matrix-     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.Internal.Class.Mutable as C---- | Column-major mutable matrix.-data MSparseMatrix :: C.MMatrixKind where-    MSparseMatrix :: (SingI r, SingI c)-                  => !(v s a)         -- ^ Values: stores the coefficient values-                                      -- of the non-zeros.-                  -> !(S.Vector Int)  -- ^ InnerIndices: stores the row-                                      -- (resp. column) indices of the non-zeros.-                  -> !(S.Vector Int)  -- ^ OuterStarts: stores for each column-                                      -- (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)-    dim (MSparseMatrix _ _ _) = (r,c)-      where-        r = fromIntegral $ fromSing (sing :: Sing r)-        c = fromIntegral $ fromSing (sing :: Sing c)-    {-# INLINE dim #-}--{--    unsafeRead mat@(MMatrix v) (i,j) = GM.unsafeRead v idx-      where-        (r, _) = C.dim mat-        idx = i + j * r-    {-# INLINE unsafeRead #-}--    unsafeWrite mat@(MMatrix v) (i,j) = GM.unsafeWrite v idx-      where -        (r, _) = C.dim mat-        idx = i + j * r-    {-# INLINE unsafeWrite #-}--    new :: forall r c s. (SingI r, SingI c, PrimMonad s)-        => s (MMatrix r c v (PrimState s) a)-    new = MMatrix <$> GM.new (r*c)-      where-        r = fromIntegral $ fromSing (sing :: Sing r)-        c = fromIntegral $ fromSing (sing :: Sing c)-    {-# INLINE new #-}--    replicate :: forall r c s. (SingI r, SingI c, PrimMonad s)-              => a -> s (MMatrix r c v (PrimState s) a)-    replicate x = MMatrix <$> GM.replicate (r*c) x-      where-        r = fromIntegral $ fromSing (sing :: Sing r)-        c = fromIntegral $ fromSing (sing :: Sing c)-    {-# INLINE replicate #-}-    -}
+ src/Data/Matrix/Static/Dense.hs view
@@ -0,0 +1,370 @@+{-# LANGUAGE BangPatterns          #-}+{-# LANGUAGE DeriveGeneric         #-}+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE StrictData #-}+{-# LANGUAGE InstanceSigs #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE TypeFamilies          #-}+{-# LANGUAGE TypeOperators #-}+module Data.Matrix.Static.Dense+    (+    -- * Immutable Matrix+      Matrix(..)++    -- * Accessors+    -- ** length information+    , C.dim+    , C.rows+    , C.cols++    -- ** Query+    , (C.!)+    , C.takeRow+    , C.takeColumn+    , C.takeDiag++    -- ** Unsafe Query+    , C.unsafeIndex+    , C.unsafeTakeRow+    , C.unsafeTakeColumn++    -- * Construction+    , C.empty+    , C.matrix+    , C.fromVector+    , C.fromList+    , C.fromRows+    , C.fromColumns+    , C.unsafeFromVector++    , diag+    , diagRect++    -- * Conversions+    , C.flatten+    , C.toRows+    , C.toColumns+    , C.toList++    -- * Conversion between Different matrix types+    , convert+    , C.convertAny++    , C.transpose++    , C.map+    , C.mapM+    , C.imap+    , C.imapM++    -- * Zipping+    , zip+    , zip3+    , zipWith+    , zipWith3+    , zipWith4+    , zipWith5+    , zipWith6+    , izipWith+    , izipWith3++    -- * Monadic Zipping+    , zipWithM+    , zipWithM_++    -- * Unzipping+    , unzip+    , unzip3++    , generate++    -- * Mutable matrix+    , C.thaw+    , C.unsafeThaw+    , C.freeze+    , C.unsafeFreeze+    , 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 GHC.TypeLits (type (<=))+import Data.Singletons+import Data.Tuple (swap)+import qualified Data.List as L+import Text.Printf (printf)++import           Data.Matrix.Static.Dense.Mutable (MMatrix (..))+import qualified Data.Matrix.Static.Dense.Mutable as DM+import qualified Data.Matrix.Static.Generic as C++type instance C.Mutable Matrix = MMatrix++-- | Column-major matrix+data Matrix :: C.MatrixKind where+    Matrix :: (SingI r, SingI c) => v a -> Matrix r c v a++instance (G.Vector v a, Show a) => Show (Matrix r c v a) where+    show mat = printf "(%d x %d)\n%s" r c vals+      where+        (r,c) = C.dim mat+        vals = unlines $ L.map (unwords . L.map show . G.toList) $ C.toRows mat++instance (G.Vector v a, Eq (v a)) => Eq (Matrix r c v a) where+    (==) (Matrix v1) (Matrix v2) = v1 == v2++instance (SingI r, SingI c, G.Vector v a, Num a) =>+    Num (Matrix r c v a) where+        m1 + m2 = zipWith (+) m1 m2+        m1 - m2 = zipWith (-) m1 m2+        m1 * m2 = zipWith (*) m1 m2+        negate = C.map negate+        abs = C.map abs+        signum = undefined+        fromInteger = undefined++instance (SingI r, SingI c, G.Vector v a, Fractional a) =>+    Fractional (Matrix r c v a) where+        m1 / m2 = zipWith (/) m1 m2+        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)+    dim (Matrix _) = (r,c)+      where+        r = fromIntegral $ fromSing (sing :: Sing r)+        c = fromIntegral $ fromSing (sing :: Sing c)+    {-# INLINE dim #-}++    -- | O(1) Unsafe indexing without bound check.+    unsafeIndex mat@(Matrix vec) (i,j) = vec `G.unsafeIndex` idx+      where+        idx = i + j * fst (C.dim mat)+    {-# INLINE unsafeIndex #-}++    -- | O(1) Create matrix from vector.+    unsafeFromVector = Matrix+    {-# INLINE unsafeFromVector #-}++    -- | O(1) Extract a row.+    unsafeTakeColumn mat@(Matrix vec) i = G.slice (i*r) r vec+      where+        (r, _) = C.dim mat+    {-# INLINE unsafeTakeColumn #-}++    -- | Create a vector by concatenating columns.+    flatten (Matrix vec) = vec+    {-# INLINE flatten #-}++    thaw (Matrix v) = MMatrix <$> G.thaw v+    {-# INLINE thaw #-}++    unsafeThaw (Matrix v) = MMatrix <$> G.unsafeThaw v+    {-# INLINE unsafeThaw #-}++    freeze (MMatrix v) = Matrix <$> G.freeze v+    {-# INLINE freeze #-}++    unsafeFreeze (MMatrix v) = Matrix <$> G.unsafeFreeze v+    {-# INLINE unsafeFreeze #-}++    map f (Matrix vec) = Matrix $ G.map f vec+    {-# INLINE map #-}++    imap f m@(Matrix vec) = Matrix $ G.imap g vec+      where+        g i = f (toIndex (C.rows m) i)+    {-# INLINE imap #-}++    imapM_ f m@(Matrix vec) = G.imapM_ g vec+      where+        g i = f (toIndex (C.rows m) i)+    {-# INLINE imapM_ #-}++    sequence (Matrix vec) = Matrix <$> G.sequence vec+    {-# INLINE sequence #-}++    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 rectangular matrix with default values and given diagonal+diagRect :: (G.Vector v a, SingI r, SingI c, n <= r, n <= c)+         => a                    -- ^ default value+         -> Matrix n 1 v a       -- ^ diagonal+         -> Matrix r c v a+diagRect z0 d = C.create $ do+    mat <- DM.replicate z0+    C.imapM_ (DM.unsafeWrite mat) d+    return mat+{-# INLINE diagRect #-}++zipWith :: (G.Vector v a, G.Vector v b, G.Vector v c+           , SingI n, SingI m )+        => (a -> b -> c)+        -> Matrix n m v a -> Matrix n m v b -> Matrix n m v c+zipWith f m1 m2 = C.unsafeFromVector $ G.zipWith f (C.flatten m1) $ C.flatten m2+{-# INLINE zipWith #-}++zipWith3 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d+            , SingI n, SingI m )+         => (a -> b -> c -> d)+         -> Matrix n m v a -> Matrix n m v b -> Matrix n m v c+         -> Matrix n m v d+zipWith3 f m1 m2 m3 = C.unsafeFromVector $+    G.zipWith3 f (C.flatten m1) (C.flatten m2) $ C.flatten m3+{-# INLINE zipWith3 #-}++zipWith4 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d, G.Vector v e+            , SingI n, SingI m )+         => (a -> b -> c -> d -> e)+         -> Matrix n m v a+         -> Matrix n m v b+         -> Matrix n m v c+         -> Matrix n m v d+         -> Matrix n m v e+zipWith4 f m1 m2 m3 m4 = C.unsafeFromVector $+    G.zipWith4 f (C.flatten m1) (C.flatten m2) (C.flatten m3) $ C.flatten m4+{-# INLINE zipWith4 #-}++zipWith5 :: ( G.Vector v a, G.Vector v b, G.Vector v c,G.Vector v d+            , G.Vector v e, G.Vector v f+            , SingI n, SingI m )+         => (a -> b -> c -> d -> e -> f)+         -> Matrix n m v a+         -> Matrix n m v b+         -> Matrix n m v c+         -> Matrix n m v d+         -> Matrix n m v e+         -> Matrix n m v f+zipWith5 f m1 m2 m3 m4 m5 = C.unsafeFromVector $+    G.zipWith5 f (C.flatten m1) (C.flatten m2)+    (C.flatten m3) (C.flatten m4) $ C.flatten m5+{-# INLINE zipWith5 #-}++zipWith6 :: ( G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d+            , G.Vector v e, G.Vector v f, G.Vector v g+            , SingI n, SingI m )+         => (a -> b -> c -> d -> e -> f -> g)+         -> Matrix n m v a+         -> Matrix n m v b+         -> Matrix n m v c+         -> Matrix n m v d+         -> Matrix n m v e+         -> Matrix n m v f+         -> Matrix n m v g+zipWith6 f m1 m2 m3 m4 m5 m6 = C.unsafeFromVector $+    G.zipWith6 f (C.flatten m1) (C.flatten m2) (C.flatten m3)+    (C.flatten m4) (C.flatten m5) $ C.flatten m6+{-# INLINE zipWith6 #-}++izipWith :: ( G.Vector v a, G.Vector v b, G.Vector v c+            , SingI n, SingI m )+         => ((Int, Int) -> a -> b -> c)+         -> Matrix n m v a -> Matrix n m v b -> Matrix n m v c+izipWith f m1 m2 = C.unsafeFromVector $+    G.izipWith g (C.flatten m1) $ C.flatten m2+  where+    g i = f (toIndex (C.rows m1) i)+{-# INLINE izipWith #-}++izipWith3 :: ( G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d+             , SingI n, SingI m )+          => ((Int, Int) -> a -> b -> c -> d)+          -> Matrix n m v a -> Matrix n m v b -> Matrix n m v c+          -> Matrix n m v d+izipWith3 f m1 m2 m3 = C.unsafeFromVector $ G.izipWith3 g+    (C.flatten m1) (C.flatten m2) $ C.flatten m3+  where+    g i = f (toIndex (C.rows m1) i)+{-# INLINE izipWith3 #-}++zip :: (SingI n, SingI m, G.Vector v a, G.Vector v b, G.Vector v (a,b))+    => Matrix n m v a -> Matrix n m v b -> Matrix n m v (a,b)+zip m1 m2 = C.unsafeFromVector $ G.zip (C.flatten m1) $ C.flatten m2+{-# INLINE zip #-}++zip3 :: (SingI n, SingI m, G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v (a,b,c))+     => Matrix n m v a+     -> Matrix n m v b+     -> Matrix n m v c+     -> Matrix n m v (a,b,c)+zip3 m1 m2 m3 = C.unsafeFromVector $+    G.zip3 (C.flatten m1) (C.flatten m2) $ C.flatten m3+{-# INLINE zip3 #-}++zipWithM :: ( G.Vector v a, G.Vector v b, G.Vector v c+            , Monad monad, SingI n, SingI m )+         => (a -> b -> monad c)+         -> Matrix n m v a -> Matrix n m v b -> monad (Matrix n m v c)+zipWithM f m1 m2 = liftM C.unsafeFromVector $+    G.zipWithM f (C.flatten m1) $ C.flatten m2+{-# INLINE zipWithM #-}++zipWithM_ :: (G.Vector v a, G.Vector v b, G.Vector v c, Monad monad)+          => (a -> b -> monad c)+          -> Matrix n m v a -> Matrix n m v b -> monad ()+zipWithM_ f m1 m2 = G.zipWithM_ f (C.flatten m1) $ C.flatten m2+{-# INLINE zipWithM_ #-}++unzip :: ( G.Vector v a, G.Vector v b, G.Vector v (a,b)+         , SingI n, SingI m )+      => Matrix n m v (a,b) -> (Matrix n m v a, Matrix n m v b )+unzip m = (C.unsafeFromVector v1, C.unsafeFromVector v2)+  where+    (v1, v2) = G.unzip $ C.flatten m+{-# INLINE unzip #-}++unzip3 :: ( G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v (a,b,c)+          , SingI n, SingI m )+       => Matrix n m v (a, b, c)+       -> (Matrix n m v a, Matrix n m v b, Matrix n m v c)+unzip3 m = (C.unsafeFromVector v1, C.unsafeFromVector v2, C.unsafeFromVector v3)+  where+    (v1, v2, v3) = G.unzip3 $ C.flatten m+{-# INLINE unzip3 #-}++generate :: forall r c v a. (G.Vector v a, SingI r, SingI c)+         => ((Int, Int) -> a) -> Matrix r c v a+generate f = C.unsafeFromVector . G.generate (r*c) $ \i -> f (i `divMod` r)+  where+    r = fromIntegral $ fromSing (sing :: Sing r)+    c = fromIntegral $ fromSing (sing :: Sing c)+{-# INLINE generate #-}++-- | O(m*n) Convert different matrix type+convert :: (G.Vector v a, G.Vector w a) => Matrix r c v a -> Matrix r c w a+convert (Matrix vec) = Matrix $ G.convert vec+{-# INLINE convert #-}+++-- Helper+toIndex :: Int -> Int -> (Int, Int)+toIndex r i = swap $ i `divMod` r+{-# INLINE toIndex #-}
+ src/Data/Matrix/Static/Dense/Mutable.hs view
@@ -0,0 +1,78 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE InstanceSigs #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE StrictData #-}+{-# LANGUAGE MultiParamTypeClasses #-}+module Data.Matrix.Static.Dense.Mutable+   ( -- * Mutable Matrix+     MMatrix(..)+   , C.dim+   , C.unsafeWrite+   , C.unsafeRead+   , C.new+   , C.replicate+   ) where++import           Control.DeepSeq+import qualified Data.Vector.Generic.Mutable as GM+import           Prelude                     hiding (read, replicate)+import Data.Singletons+import Control.Monad.Primitive     (PrimMonad, PrimState)++import qualified Data.Matrix.Static.Generic.Mutable as C++-- | 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)+    dim (MMatrix _) = (r,c)+      where+        r = fromIntegral $ fromSing (sing :: Sing r)+        c = fromIntegral $ fromSing (sing :: Sing c)+    {-# INLINE dim #-}++    unsafeRead mat@(MMatrix v) (i,j) = GM.unsafeRead v idx+      where+        (r, _) = C.dim mat+        idx = i + j * r+    {-# INLINE unsafeRead #-}++    unsafeWrite mat@(MMatrix v) (i,j) = GM.unsafeWrite v idx+      where +        (r, _) = C.dim mat+        idx = i + j * r+    {-# INLINE unsafeWrite #-}++    new :: forall r c s. (SingI r, SingI c, PrimMonad s)+        => s (MMatrix r c v (PrimState s) a)+    new = MMatrix <$> GM.new (r*c)+      where+        r = fromIntegral $ fromSing (sing :: Sing r)+        c = fromIntegral $ fromSing (sing :: Sing c)+    {-# INLINE new #-}++    replicate :: forall r c s. (SingI r, SingI c, PrimMonad s)+              => a -> s (MMatrix r c v (PrimState s) a)+    replicate x = MMatrix <$> GM.replicate (r*c) x+      where+        r = fromIntegral $ fromSing (sing :: Sing r)+        c = fromIntegral $ fromSing (sing :: Sing c)+    {-# INLINE replicate #-}++{-+takeColumn :: GM.MVector v a => MMatrix v m a -> Int -> v m a+takeColumn (MMatrix _ c tda offset vec) i = GM.slice i' c vec+  where+    i' = offset + i * tda+{-# INLINE takeColumn #-}+-}
+ src/Data/Matrix/Static/Generic.hs view
@@ -0,0 +1,263 @@+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE TypeFamilyDependencies #-}+{-# LANGUAGE Rank2Types            #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies          #-}+{-# LANGUAGE TypeOperators #-}+module Data.Matrix.Static.Generic+    ( Mutable+    , Matrix(..)+    , MatrixKind++    -- * Derived mothods+    , rows+    , cols+    , (!)+    , takeColumn+    , takeRow+    , toRows+    , toColumns+    , empty+    , matrix+    , withMatrix+    , fromRows+    , withRows+    , fromColumns+    , withColumns+    , fromVector+    , fromList+    , toList+    , create+    , convertAny+    , mapM+    , imapM+    ) where++import           Control.Monad.Primitive     (PrimMonad, PrimState)+import           Control.Monad.ST            (ST, runST)+import qualified Data.Vector.Generic         as G+import Text.Printf (printf)+import Prelude hiding (map, mapM, mapM_, sequence, sequence_)+import qualified Data.List as L+import Data.Kind (Type)+import GHC.TypeLits (Nat, type (<=))+import Data.Singletons (SingI, Sing, fromSing, sing, withSomeSing)+import Data.Singletons.TypeLits++import Data.Matrix.Static.Generic.Mutable (MMatrix, MMatrixKind)++type MatrixKind = Nat -> Nat -> (Type -> Type) -> Type -> Type++type family Mutable (mat :: MatrixKind) = (mmat :: MMatrixKind) | mmat -> mat++class (MMatrix (Mutable mat) (G.Mutable v) a, G.Vector v a) => Matrix (mat :: MatrixKind) v a where+    dim :: mat r c v a -> (Int, Int)++    unsafeIndex :: mat r c v a -> (Int, Int) -> a++    unsafeFromVector :: (SingI r, SingI c) => v a -> mat r c v a++    -- | 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)+      where+        (r,c) = dim mat+    {-# INLINE flatten #-}++    -- | Extract a row. Default algorithm is O(n * O(unsafeIndex)).+    unsafeTakeRow :: mat r c v a -> Int -> v a+    unsafeTakeRow mat i = G.generate c $ \j -> unsafeIndex mat (i,j)+      where+        (_,c) = dim mat+    {-# INLINE unsafeTakeRow #-}++    -- | Extract a column. Default algorithm is O(m * O(unsafeIndex)).+    unsafeTakeColumn :: mat r c v a -> Int -> v a+    unsafeTakeColumn mat j = G.generate r $ \i -> unsafeIndex mat (i,j)+      where+        (r,_) = dim mat+    {-# INLINE unsafeTakeColumn #-}++    -- | Extract the diagonal. Default algorithm is O(min(m,n) * O(unsafeIndex)).+    takeDiag :: mat r c v a -> v a+    takeDiag mat = G.generate n $ \i -> unsafeIndex mat (i,i)+      where+        n = uncurry min . dim $ mat+    {-# INLINE takeDiag #-}++    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+      where+       (r, c) = dim mat+    {-# INLINE transpose #-}++    thaw :: PrimMonad s+         => mat r c v a+         -> s ((Mutable mat) r c (G.Mutable v) (PrimState s) a)++    unsafeThaw :: PrimMonad s+               => mat r c v a+               -> s ((Mutable mat) r c (G.Mutable v) (PrimState s) a)++    freeze :: PrimMonad s+           => (Mutable mat) r c (G.Mutable v) (PrimState s) a+           -> s (mat r c v a)++    unsafeFreeze :: PrimMonad s+                 => (Mutable mat) r c (G.Mutable v) (PrimState s) a+                 -> s (mat r c v a)++    map :: G.Vector v b => (a -> b) -> mat r c v a -> mat r c v b+    imap :: G.Vector v b => ((Int, Int) -> a -> b) -> mat r c v a -> mat r c v b+    imapM_ :: (Monad monad, Matrix mat v a)+           => ((Int, Int) -> a -> monad b) -> mat r c v a -> monad ()+    sequence :: (G.Vector v (monad a), Monad monad)+             => mat r c v (monad a) -> monad (mat r c v a)+    sequence_ :: (G.Vector v (monad a), Monad monad) => mat r c v (monad a) -> monad ()++-- | Derived methods++-- | Return the number of rows+rows :: Matrix m v a => m r c v a -> Int+rows = fst . dim+{-# INLINE rows #-}++-- | Return the number of columns+cols :: Matrix m v a => m r c v a -> Int+cols = snd . dim+{-# INLINE cols #-}++-- | Indexing+(!) :: forall m r c v a i j. (Matrix m v a, i <= r, j <= c)+    => m r c v a -> (Sing i, Sing j) -> a+(!) m (si, sj) = unsafeIndex m (i,j)+  where+    i = fromIntegral $ fromSing si+    j = fromIntegral $ fromSing sj+{-# INLINE (!) #-}++-- | Construct matrix from a vector containg columns.+fromVector :: forall m r c v a. (SingI r, SingI c, Matrix m v a)+           => v a -> m r c v a+fromVector vec | r*c /= n = error errMsg+               | otherwise = unsafeFromVector vec+  where+    errMsg = printf "fromVector: incorrect length (%d * %d != %d)" r c n+    n = G.length vec+    r = fromIntegral $ fromSing (sing :: Sing r)+    c = fromIntegral $ fromSing (sing :: Sing c)+{-# INLINE fromVector #-}++matrix :: (SingI r, SingI c, Matrix m v a)+       => [[a]] -> m r c v a+matrix = fromList . concat . L.transpose+{-# INLINE matrix #-}++withMatrix :: forall mat v a b. Matrix mat v a+           => [[a]] -> (forall r c. mat r c v a -> b) -> b+withMatrix xs f = withSomeSing n $ \(SNat :: Sing n) -> +    withSomeSing m $ \(SNat :: Sing m) -> f (matrix xs :: mat n m v a)+  where+    n = fromIntegral $ length xs+    m = fromIntegral $ length $ head xs+{-# INLINE withMatrix #-}++-- | Construct matrix from a list containg columns.+fromList :: (SingI r, SingI c, Matrix m v a)+         => [a] -> m r c v a+fromList = fromVector . G.fromList+{-# INLINE fromList #-}++-- | O(m*n) Create matrix from rows+fromRows :: (Matrix m v a, SingI r, SingI c) => [v a] -> m r c v a+fromRows = transpose . fromColumns+{-# INLINE fromRows #-}++withRows :: forall mat v a b. Matrix mat v a+          => [v a] -> (forall r c. mat r c v a -> b) -> b+withRows xs f = withSomeSing n $ \(SNat :: Sing n) -> +    withSomeSing m $ \(SNat :: Sing m) -> f (fromRows xs :: mat n m v a)+  where+    n = fromIntegral $ length xs+    m = fromIntegral $ G.length $ head xs+{-# INLINE withRows #-}++-- | O(m*n) Create matrix from columns+fromColumns :: (Matrix m v a, SingI r, SingI c)+            => [v a] -> m r c v a+fromColumns = fromVector . G.concat+{-# INLINE fromColumns #-}++withColumns :: forall mat v a b. Matrix mat v a+            => [v a] -> (forall r c. mat r c v a -> b) -> b+withColumns xs f = withSomeSing n $ \(SNat :: Sing n) -> +    withSomeSing m $ \(SNat :: Sing m) -> f (fromRows xs :: mat n m v a)+  where+    m = fromIntegral $ length xs+    n = fromIntegral $ G.length $ head xs+{-# INLINE withColumns #-}++-- | O(m*n) Create a list by concatenating columns+toList :: Matrix m v a => m r c v a -> [a]+toList = G.toList . flatten+{-# INLINE toList #-}++empty :: Matrix m v a => m 0 0 v a+empty = unsafeFromVector G.empty+{-# INLINE empty #-}++create :: Matrix m v a+       => (forall s . ST s ((Mutable m) r c (G.Mutable v) s a)) -> m r c v a+create m = runST $ unsafeFreeze =<< m+{-# INLINE create #-}++-- | O(m*n) Convert to any type of matrix.+convertAny :: (Matrix m1 v1 a, Matrix m2 v2 a, SingI r, SingI c)+           => m1 r c v1 a -> m2 r c v2 a+convertAny = unsafeFromVector . G.convert . flatten+{-# INLINE convertAny #-}++-- | Extract a row.+takeRow :: forall m r c v a i. (i <= r, SingI i, Matrix m v a)+        => m r c v a -> Sing i -> v a+takeRow mat _ = unsafeTakeRow mat i+  where+    i = fromIntegral $ fromSing (sing :: Sing i)+{-# INLINE takeRow #-}++-- | O(m) Return the rows+toRows :: Matrix m v a => m r c v a -> [v a]+toRows mat = L.map (unsafeTakeRow mat) [0..r-1]+  where+    (r,_) = dim mat+{-# INLINE toRows #-}++-- | Extract a row.+takeColumn :: forall m r c v a j. (j <= c, SingI j, Matrix m v a)+           => m r c v a -> Sing j -> v a+takeColumn mat _ = unsafeTakeColumn mat j+  where+    j = fromIntegral $ fromSing (sing :: Sing j)+{-# INLINE takeColumn #-}++-- | O(m*n) Return the columns+toColumns :: Matrix m v a => m r c v a -> [v a]+toColumns mat = L.map (unsafeTakeColumn mat) [0..c-1]+  where+    (_,c) = dim mat+{-# INLINE toColumns #-}++mapM :: (G.Vector v (monad b), Monad monad, Matrix mat v a, Matrix mat v b)+     => (a -> monad b) -> mat r c v a -> monad (mat r c v b)+mapM f = sequence . map f+{-# INLINE mapM #-}++imapM :: (G.Vector v (monad b), Monad monad, Matrix mat v a, Matrix mat v b)+      => ((Int, Int) -> a -> monad b)+      -> mat r c v a -> monad (mat r c v b)+imapM f = sequence . imap f+{-# INLINE imapM #-}
+ src/Data/Matrix/Static/Generic/Mutable.hs view
@@ -0,0 +1,51 @@+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE TypeFamilies          #-}++module Data.Matrix.Static.Generic.Mutable+    ( MMatrix(..)+    , MMatrixKind+    ) where++import           Control.Monad.Primitive     (PrimMonad, PrimState)+import qualified Data.Vector.Generic.Mutable as GM+import           Prelude                     hiding (read)+import Data.Kind (Type)+import GHC.TypeLits (Nat)+import Data.Singletons (SingI)++type MMatrixKind = Nat -> Nat -> (Type -> Type -> Type) -> Type -> Type -> Type++class GM.MVector v a => MMatrix (mat :: MMatrixKind) v a where+    dim :: mat r c v s a -> (Int, Int)++    unsafeRead :: PrimMonad s => mat r c v (PrimState s) a -> (Int, Int) -> s a++    unsafeWrite :: PrimMonad s => mat r c v (PrimState s) a -> (Int, Int) -> a -> s ()++    -- | Create a mutable matrix without initialization+    new :: (SingI r, SingI c, PrimMonad s) => s (mat r c v (PrimState s) a)++    replicate :: (SingI r, SingI c, PrimMonad s) => a -> s (mat r c v (PrimState s) a)++    {-# MINIMAL dim, unsafeRead, unsafeWrite, new, replicate #-}++{-+write :: (PrimMonad s, MMatrix m v a)+      => m v (PrimState s) a -> (Int, Int) -> a -> s ()+write mat (i,j)+    | i < 0 || i >= r || j < 0 || j >= c = error "write: Index out of bounds"+    | otherwise = unsafeWrite mat (i,j)+  where+    (r,c) = dim mat+{-# INLINE write #-}++read :: (PrimMonad s, MMatrix m v a)+     => m v (PrimState s) a -> (Int, Int) -> s a+read mat (i,j)+    | i <0 || i >= r || j < 0 || j >= c = error "read: Index out of bounds"+    | otherwise = unsafeRead mat (i,j)+  where+    (r,c) = dim mat+{-# INLINE read #-}+-}
+ src/Data/Matrix/Static/Internal.hs view
@@ -0,0 +1,67 @@+{-# LANGUAGE ForeignFunctionInterface #-}+module Data.Matrix.Static.Internal+    ( c_dd_mul+    , c_ds_mul+    , c_sd_mul+    , c_ss_mul+    , c_inverse+    , c_cholesky+    , c_eigs+    , c_seigs+    ) where++import Data.Complex (Complex)+import Foreign+import Foreign.C.Types+import Foreign.C.String++-------------------------------------------------------------------------------+-- Arithmetic+-------------------------------------------------------------------------------+foreign import ccall "eigen_dd_mul"+    c_dd_mul :: CInt+          -> Ptr a -> CInt -> CInt+          -> Ptr a -> CInt -> CInt+          -> Ptr a -> CInt -> CInt+          -> IO CString++foreign import ccall "eigen_ds_mul"+    c_ds_mul :: CInt+          -> Ptr a -> CInt -> CInt+          -> Ptr a -> CInt -> CInt+          -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt+          -> IO CString++foreign import ccall "eigen_sd_mul"+    c_sd_mul :: 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 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+              -> Ptr a -> CInt -> CInt+              -> Ptr a -> CInt -> CInt+              -> IO CString++foreign import ccall "eigen_cholesky"+    c_cholesky :: CInt+               -> Ptr a -> Ptr a ->  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_seigs"+    c_seigs :: CInt -> Ptr (Complex Double) -> Ptr (Complex Double)+            -> Ptr Double -> Ptr CInt -> Ptr CInt+            -> CInt -> CInt -> IO CString
+ src/Data/Matrix/Static/LinearAlgebra.hs view
@@ -0,0 +1,98 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE ExplicitNamespaces #-}+{-# LANGUAGE TypeOperators #-}+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 qualified Data.Matrix.Static.Dense as D+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 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+          -> 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 S.SparseMatrix D.Matrix where+    (%*%) = withDS Internal.c_ds_mul++instance Arithmetic S.SparseMatrix D.Matrix D.Matrix where+    (%*%) = withSD Internal.c_sd_mul++instance Arithmetic S.SparseMatrix S.SparseMatrix S.SparseMatrix where+    (%*%) = withSS Internal.c_ss_mul++inverse :: (SingI n, Numeric a) => Matrix n n a -> Matrix n n a+inverse = withFun1 Internal.c_inverse++class Factorization mat where+    -- | Eigenvalues (not ordered) and+    -- eigenvectors (as columns) of a general square matrix.+    eigs :: (SingI k, SingI n, k <= n - 2)+         => Sing k+         -> mat n n VS.Vector Double+         -> (Matrix k 1 (Complex Double), Matrix n k (Complex Double))++    -- | Cholesky decomposition+    cholesky :: (Numeric a, SingI n) => mat n n VS.Vector a -> mat n n VS.Vector a++instance Factorization D.Matrix where+    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')+      where+        k = fromIntegral $ fromSing s+    {-# INLINE eigs #-}++    cholesky mat = flip withFun1 mat $+        \code p1 c1 _ p2 _ _ -> Internal.c_cholesky code p1 p2 c1+    {-# INLINE cholesky #-}++instance Factorization S.SparseMatrix where+    eigs 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_seigs 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 eigs #-}++    cholesky = undefined
+ src/Data/Matrix/Static/LinearAlgebra/Types.hs view
@@ -0,0 +1,203 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE DataKinds #-}+{-# OPTIONS_GHC -Wno-simplifiable-class-constraints #-}+module Data.Matrix.Static.LinearAlgebra.Types+    ( Numeric(..)+    , 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) => Numeric a where+    foreignType :: a -> CInt++instance Numeric Float where foreignType _ = 0+instance Numeric Double where foreignType _ =1+instance Numeric (Complex Float) where foreignType _ = 2+instance Numeric (Complex Double) where foreignType _ = 3++type Matrix r c a = D.Matrix r c Vector a+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 (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+    return $ S.SparseMatrix (VS.unsafeFromForeignPtr0 pv n)+        (VS.unsafeFromForeignPtr0 pinner n)+        (VS.unsafeFromForeignPtr0 pouter $ c + 1)+  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 (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
@@ -0,0 +1,197 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE InstanceSigs #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE StrictData #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE TypeFamilies #-}+module Data.Matrix.Static.Sparse+   ( -- * Sparse matrix+     SparseMatrix(..)++   , Zero(..)++    -- * Accessors+    -- ** length information+    , C.dim+    , C.rows+    , C.cols++    -- ** Query+    , (C.!)+    , C.takeDiag++    -- ** Unsafe Query+    , C.unsafeIndex+    , C.unsafeTakeRow+    , C.unsafeTakeColumn++    -- * Construction+    , C.empty+    , C.fromVector+    , C.fromList+    , C.unsafeFromVector++    , diag+    , diagRect++    -- * Conversions+    , C.flatten+    , C.toList++    -- * Different matrix types+    , C.convertAny+   ) where++import           Control.DeepSeq+import qualified Data.Vector.Generic as G+import qualified Data.Vector.Storable as S+import qualified Data.Vector.Storable.Mutable as SM+import Data.Singletons+import Control.Monad+import           Data.Bits                         (shiftR)+import Text.Printf (printf)+import GHC.TypeLits (type (<=))+import Foreign.C.Types+import Data.Complex++import qualified Data.Matrix.Static.Dense as D+import qualified Data.Matrix.Static.Generic as C+import Data.Matrix.Static.Sparse.Mutable++type instance C.Mutable SparseMatrix = MSparseMatrix++class Eq a => Zero a where+    zero :: a++instance Zero Int where+    zero = 0++instance Zero Float where+    zero = 0.0++instance Zero Double where+    zero = 0.0++instance Zero (Complex Float) where+    zero = 0++instance Zero (Complex Double) where+    zero = 0++instance Eq a => Zero ([] a) where+    zero = []++-- | Column-major mutable matrix.+data SparseMatrix :: C.MatrixKind where+    SparseMatrix :: (SingI r, SingI c)+                 => (v a)           -- ^ Values: stores the coefficient values+                                     -- of the non-zeros.+                 -> (S.Vector CInt)  -- ^ InnerIndices: stores the row+                                     -- (resp. column) indices of the non-zeros.+                 -> (S.Vector CInt)  -- ^ OuterStarts: stores for each column+                                     -- (resp. row) the index of the first+                                     -- non-zero in the previous two arrays.+                 -> SparseMatrix r c v a++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 (NFData (v a)) => NFData (SparseMatrix r c v a) where+    rnf (SparseMatrix vec inner outer) = rnf vec++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)+    dim (SparseMatrix _ _ _) = (r,c)+      where+        r = fromIntegral $ fromSing (sing :: Sing r)+        c = fromIntegral $ fromSing (sing :: Sing c)+    {-# INLINE dim #-}++    -- | O(1) Unsafe indexing without bound check.+    unsafeIndex (SparseMatrix vec inner outer) (i,j) = +        case binarySearchByBounds inner (fromIntegral i) r0 r1 of+            Nothing -> zero+            Just k -> vec `G.unsafeIndex` k+      where+        r0 = fromIntegral $ outer `S.unsafeIndex` j+        r1 = fromIntegral $ outer `S.unsafeIndex` (j+1) - 1+    {-# INLINE unsafeIndex #-}++    -- | 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+      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+        r = fromIntegral $ fromSing (sing :: Sing r)+        c = fromIntegral $ fromSing (sing :: Sing c)+    {-# INLINE unsafeFromVector #-}++    thaw = undefined+    {-# INLINE thaw #-}++    unsafeThaw = undefined+    {-# INLINE unsafeThaw #-}++    freeze = undefined+    {-# INLINE freeze #-}++    unsafeFreeze = undefined+    {-# INLINE unsafeFreeze #-}++    map f (SparseMatrix vec inner outer) = SparseMatrix (G.map f vec) inner outer+    imap = undefined+    {-# INLINE map #-}++-- | O(m*n) Create a square matrix with given diagonal.+diag :: (G.Vector v a, Zero a, SingI n)+     => D.Matrix n 1 v a       -- ^ diagonal+     -> SparseMatrix n n v a+diag = diagRect+{-# INLINE diag #-}++-- | O(m*n) Create a rectangular matrix with default values and given diagonal+diagRect :: (G.Vector v a, Zero a, SingI r, SingI c, n <= r, n <= c)+         => D.Matrix n 1 v a       -- ^ diagonal+         -> SparseMatrix r c v a+diagRect d = SparseMatrix (C.flatten d) (S.enumFromN 0 n) (S.enumFromN 0 $ n + 1)+  where+    n = C.rows d+{-# INLINE diagRect #-}++binarySearchByBounds :: S.Vector CInt -> CInt -> Int -> Int -> Maybe Int+binarySearchByBounds vec x = loop+  where+    loop !l !u+        | l > u = Nothing+        | x == x' = Just k+        | x < x' = loop l (k-1)+        | otherwise = loop (k+1) u+      where+        k = (u+l) `shiftR` 1+        x' = vec `S.unsafeIndex` k+{-# INLINE binarySearchByBounds #-}
+ src/Data/Matrix/Static/Sparse/Mutable.hs view
@@ -0,0 +1,76 @@+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE InstanceSigs #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE StrictData #-}+{-# LANGUAGE MultiParamTypeClasses #-}+module Data.Matrix.Static.Sparse.Mutable+   ( -- * Mutable sparse matrix+     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++-- | Column-major mutable matrix.+data MSparseMatrix :: C.MMatrixKind where+    MSparseMatrix :: (SingI r, SingI c)+                  => (v s a)         -- ^ Values: stores the coefficient values+                                      -- of the non-zeros.+                  -> (S.Vector Int)  -- ^ InnerIndices: stores the row+                                      -- (resp. column) indices of the non-zeros.+                  -> (S.Vector Int)  -- ^ OuterStarts: stores for each column+                                      -- (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)+    dim (MSparseMatrix _ _ _) = (r,c)+      where+        r = fromIntegral $ fromSing (sing :: Sing r)+        c = fromIntegral $ fromSing (sing :: Sing c)+    {-# INLINE dim #-}++{-+    unsafeRead mat@(MMatrix v) (i,j) = GM.unsafeRead v idx+      where+        (r, _) = C.dim mat+        idx = i + j * r+    {-# INLINE unsafeRead #-}++    unsafeWrite mat@(MMatrix v) (i,j) = GM.unsafeWrite v idx+      where +        (r, _) = C.dim mat+        idx = i + j * r+    {-# INLINE unsafeWrite #-}++    new :: forall r c s. (SingI r, SingI c, PrimMonad s)+        => s (MMatrix r c v (PrimState s) a)+    new = MMatrix <$> GM.new (r*c)+      where+        r = fromIntegral $ fromSing (sing :: Sing r)+        c = fromIntegral $ fromSing (sing :: Sing c)+    {-# INLINE new #-}++    replicate :: forall r c s. (SingI r, SingI c, PrimMonad s)+              => a -> s (MMatrix r c v (PrimState s) a)+    replicate x = MMatrix <$> GM.replicate (r*c) x+      where+        r = fromIntegral $ fromSing (sing :: Sing r)+        c = fromIntegral $ fromSing (sing :: Sing c)+    {-# INLINE replicate #-}+    -}