diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -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
+```
diff --git a/cbits/eigen-basic.cpp b/cbits/eigen-basic.cpp
--- a/cbits/eigen-basic.cpp
+++ b/cbits/eigen-basic.cpp
@@ -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) \
diff --git a/matrix-sized.cabal b/matrix-sized.cabal
--- a/matrix-sized.cabal
+++ b/matrix-sized.cabal
@@ -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
diff --git a/src/Data/Matrix/Dense.hs b/src/Data/Matrix/Dense.hs
deleted file mode 100644
--- a/src/Data/Matrix/Dense.hs
+++ /dev/null
@@ -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 #-}
diff --git a/src/Data/Matrix/Dense/Mutable.hs b/src/Data/Matrix/Dense/Mutable.hs
deleted file mode 100644
--- a/src/Data/Matrix/Dense/Mutable.hs
+++ /dev/null
@@ -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 #-}
--}
diff --git a/src/Data/Matrix/Internal/Class.hs b/src/Data/Matrix/Internal/Class.hs
deleted file mode 100644
--- a/src/Data/Matrix/Internal/Class.hs
+++ /dev/null
@@ -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 #-}
diff --git a/src/Data/Matrix/Internal/Class/Mutable.hs b/src/Data/Matrix/Internal/Class/Mutable.hs
deleted file mode 100644
--- a/src/Data/Matrix/Internal/Class/Mutable.hs
+++ /dev/null
@@ -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 #-}
--}
diff --git a/src/Data/Matrix/Internal/LinearAlgebra.hs b/src/Data/Matrix/Internal/LinearAlgebra.hs
deleted file mode 100644
--- a/src/Data/Matrix/Internal/LinearAlgebra.hs
+++ /dev/null
@@ -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
diff --git a/src/Data/Matrix/LinearAlgebra.hs b/src/Data/Matrix/LinearAlgebra.hs
deleted file mode 100644
--- a/src/Data/Matrix/LinearAlgebra.hs
+++ /dev/null
@@ -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
diff --git a/src/Data/Matrix/LinearAlgebra/Types.hs b/src/Data/Matrix/LinearAlgebra/Types.hs
deleted file mode 100644
--- a/src/Data/Matrix/LinearAlgebra/Types.hs
+++ /dev/null
@@ -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 #-}
diff --git a/src/Data/Matrix/Sparse.hs b/src/Data/Matrix/Sparse.hs
deleted file mode 100644
--- a/src/Data/Matrix/Sparse.hs
+++ /dev/null
@@ -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 #-}
diff --git a/src/Data/Matrix/Sparse/Mutable.hs b/src/Data/Matrix/Sparse/Mutable.hs
deleted file mode 100644
--- a/src/Data/Matrix/Sparse/Mutable.hs
+++ /dev/null
@@ -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 #-}
-    -}
diff --git a/src/Data/Matrix/Static/Dense.hs b/src/Data/Matrix/Static/Dense.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Static/Dense.hs
@@ -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 #-}
diff --git a/src/Data/Matrix/Static/Dense/Mutable.hs b/src/Data/Matrix/Static/Dense/Mutable.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Static/Dense/Mutable.hs
@@ -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 #-}
+-}
diff --git a/src/Data/Matrix/Static/Generic.hs b/src/Data/Matrix/Static/Generic.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Static/Generic.hs
@@ -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 #-}
diff --git a/src/Data/Matrix/Static/Generic/Mutable.hs b/src/Data/Matrix/Static/Generic/Mutable.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Static/Generic/Mutable.hs
@@ -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 #-}
+-}
diff --git a/src/Data/Matrix/Static/Internal.hs b/src/Data/Matrix/Static/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Static/Internal.hs
@@ -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
diff --git a/src/Data/Matrix/Static/LinearAlgebra.hs b/src/Data/Matrix/Static/LinearAlgebra.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Static/LinearAlgebra.hs
@@ -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
diff --git a/src/Data/Matrix/Static/LinearAlgebra/Types.hs b/src/Data/Matrix/Static/LinearAlgebra/Types.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Static/LinearAlgebra/Types.hs
@@ -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 #-}
diff --git a/src/Data/Matrix/Static/Sparse.hs b/src/Data/Matrix/Static/Sparse.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Static/Sparse.hs
@@ -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 #-}
diff --git a/src/Data/Matrix/Static/Sparse/Mutable.hs b/src/Data/Matrix/Static/Sparse/Mutable.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Static/Sparse/Mutable.hs
@@ -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 #-}
+    -}