diff --git a/matrices.cabal b/matrices.cabal
--- a/matrices.cabal
+++ b/matrices.cabal
@@ -1,8 +1,5 @@
--- Initial matrices.cabal generated by cabal init.  For further
--- documentation, see http://haskell.org/cabal/users-guide/
-
 name:                matrices
-version:             0.4.5
+version:             0.5.0
 synopsis:            native matrix based on vector
 description:         Pure Haskell matrix library, supporting creating, indexing,
                      and modifying dense/sparse matrices.
@@ -10,7 +7,7 @@
 license-file:        LICENSE
 author:              Kai Zhang
 maintainer:          kai@kzhang.org
-copyright:           (c) 2015-2017 Kai Zhang
+copyright:           (c) 2015-2018 Kai Zhang
 category:            Data
 build-type:          Simple
 cabal-version:       >=1.10
@@ -25,11 +22,11 @@
     Data.Matrix.Unboxed.Mutable
     Data.Matrix.Generic
     Data.Matrix.Generic.Mutable
-    Data.Matrix.Dense.Generic
-    Data.Matrix.Dense.Generic.Mutable
     Data.Matrix.Sparse.Generic
-    Data.Matrix.Symmetric
-    Data.Matrix.Symmetric.Mutable
+    Data.Matrix.Symmetric.Generic
+    Data.Matrix.Symmetric.Generic.Mutable
+    Data.Matrix.Class
+    Data.Matrix.Class.Mutable
 
   ghc-options: -Wall -funbox-strict-fields
 
diff --git a/src/Data/Matrix.hs b/src/Data/Matrix.hs
--- a/src/Data/Matrix.hs
+++ b/src/Data/Matrix.hs
@@ -1,10 +1,394 @@
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE ConstraintKinds #-}
+{-# LANGUAGE KindSignatures #-}
+
 module Data.Matrix
     ( Matrix
-    , module Data.Matrix.Dense.Generic
-    )where
 
-import           Data.Matrix.Dense.Generic hiding (Matrix)
-import qualified Data.Matrix.Dense.Generic as MG
-import qualified Data.Vector               as V
+    -- * Accessors
+    -- ** length information
+    , dim
+    , rows
+    , cols
 
-type Matrix = MG.Matrix V.Vector
+    -- ** Indexing
+    , unsafeIndex
+    , (!)
+    , takeRow
+    , takeColumn
+    , takeDiag
+
+    -- * Construction
+    , unsafeFromVector
+    , fromVector
+    , matrix
+    , fromList
+    , fromLists
+    , fromRows
+    , fromColumns
+    , empty
+
+    -- * Conversions
+    , flatten
+    , toRows
+    , toColumns
+    , toList
+    , toLists
+
+    , tr
+    , subMatrix
+    , ident
+    , diag
+    , diagRect
+    , fromBlocks
+    , isSymmetric
+    , force
+
+    , foldl
+
+    -- * Mapping
+    , map
+    , imap
+
+    -- * Monadic mapping
+    , mapM
+    , imapM
+    , mapM_
+    , imapM_
+    , forM
+    , forM_
+
+    -- * Zipping
+    , zipWith
+    , zipWith3
+    , zipWith4
+    , zipWith5
+    , zipWith6
+    , izipWith
+    , izipWith3
+    , izipWith4
+    , izipWith5
+    , izipWith6
+    , zip
+    , zip3
+    , zip4
+    , zip5
+    , zip6
+
+    -- * Monadic Zipping
+    , zipWithM
+    , zipWithM_
+
+    -- * Unzipping
+    , unzip
+    , unzip3
+    , unzip4
+    , unzip5
+    , unzip6
+
+    -- * Monadic sequencing
+    , sequence
+    , sequence_
+
+    , generate
+
+    -- * Mutable matrix
+    , thaw
+    , unsafeThaw
+    , freeze
+    , unsafeFreeze
+    , create
+    ) where
+
+import GHC.Exts (Constraint)
+import Prelude hiding (sequence, sequence_, mapM_, zip, zip, zip3, zipWith, zipWith3, foldl, unzip, map, mapM, unzip3)
+import           Control.Monad.Primitive     (PrimMonad, PrimState)
+import           Control.Monad.ST            (ST)
+import Data.Foldable (Foldable)
+import Data.Vector (Vector)
+
+import qualified Data.Matrix.Generic as MG
+import Data.Matrix.Mutable (MMatrix)
+
+type Matrix = MG.Matrix Vector
+type Context x = (() :: Constraint)
+
+dim :: Context a => Matrix a -> (Int, Int)
+dim = MG.dim
+
+rows :: Context a => Matrix a -> Int
+rows = MG.rows
+
+cols :: Context a => Matrix a -> Int
+cols = MG.cols
+
+unsafeIndex :: Context a => Matrix a -> (Int, Int) -> a
+unsafeIndex = MG.unsafeIndex
+
+(!) :: Context a => Matrix a -> (Int, Int) -> a
+(!) = (MG.!)
+
+takeRow :: Context a => Matrix a -> Int -> Vector a
+takeRow = MG.takeRow
+
+takeColumn :: Context a => Matrix a -> Int -> Vector a
+takeColumn = MG.takeColumn
+
+takeDiag :: Context a => Matrix a -> Vector a
+takeDiag = MG.takeDiag
+
+unsafeFromVector :: Context a => (Int, Int) -> Vector a -> Matrix a
+unsafeFromVector = MG.unsafeFromVector
+
+fromVector :: Context a => (Int, Int) -> Vector a -> Matrix a
+fromVector = MG.fromVector
+
+-- | O(m*n) Matrix construction
+matrix :: Context a => Int -> [a] -> Matrix a
+matrix = MG.matrix
+
+fromList :: Context a => (Int, Int) -> [a] -> Matrix a
+fromList = MG.fromList
+
+-- | O(m*n) Create matrix from list of lists, it doesn't check if the list of
+-- list is a valid matrix
+fromLists :: Context a => [[a]] -> Matrix a
+fromLists = MG.fromLists
+
+-- | O(m*n) Create matrix from rows
+fromRows :: Context a => [Vector a] -> Matrix a
+fromRows = MG.fromRows
+
+-- | O(m*n) Create matrix from columns
+fromColumns :: Context a => [Vector a] -> Matrix a
+fromColumns = MG.fromColumns
+
+empty :: Context a => Matrix a
+empty = MG.empty
+
+flatten :: Context a => Matrix a -> Vector a
+flatten = MG.flatten
+
+-- | O(m) Return the rows
+toRows :: Context a => Matrix a -> [Vector a]
+toRows = MG.toRows
+
+toColumns :: Context a => Matrix a -> [Vector a]
+toColumns = MG.toColumns
+
+-- | O(m*n) Create a list by concatenating rows
+toList :: Context a => Matrix a -> [a]
+toList = MG.toList
+
+-- | O(m*n) List of lists
+toLists :: Context a => Matrix a -> [[a]]
+toLists = MG.toLists
+
+-- | O(m*n) Matrix transpose
+tr :: Context a => Matrix a -> Matrix a
+tr = MG.tr
+
+-- | O(1) Extract sub matrix
+subMatrix :: Context a
+          => (Int, Int)  -- ^ upper left corner of the submatrix
+          -> (Int, Int)  -- ^ bottom right corner of the submatrix
+          -> Matrix a -> Matrix a
+subMatrix = MG.subMatrix
+
+-- | O(m*n) Create an identity matrix
+ident :: (Context a, Num a) => Int -> Matrix a
+ident = MG.ident
+
+-- | O(m*n) Create a square matrix with given diagonal, other entries default to 0
+diag :: (Context a, Num a, Foldable t)
+     => t a  -- ^ diagonal
+     -> Matrix a
+diag = MG.diag
+
+-- | O(m*n) Create a rectangular matrix with default values and given diagonal
+diagRect :: (Context a, Foldable t)
+         => a         -- ^ default value
+         -> (Int, Int)
+         -> t a       -- ^ diagonal
+         -> Matrix a
+diagRect = MG.diagRect
+
+fromBlocks :: Context a
+           => a    -- ^ default value
+           -> [[Matrix a]]
+           -> Matrix a
+fromBlocks = MG.fromBlocks
+
+isSymmetric :: (Context a, Eq a) => Matrix a -> Bool
+isSymmetric = MG.isSymmetric
+
+force :: Context a => Matrix a -> Matrix a
+force = MG.force
+
+foldl :: Context b => (a -> b -> a) -> a -> Matrix b -> a
+foldl = MG.foldl
+
+map :: (Context a, Context b) => (a -> b) -> Matrix a -> Matrix b
+map = MG.map
+
+imap :: (Context a, Context b) => ((Int, Int) -> a -> b) -> Matrix a -> Matrix b
+imap = MG.imap
+
+mapM :: (Context a, Context b, Monad m) => (a -> m b) -> Matrix a -> m (Matrix b)
+mapM = MG.mapM
+
+-- | O(m*n) Apply the monadic action to every element and its index,
+-- yielding a matrix of results.
+imapM :: (Context a, Context b, Monad m) => ((Int, Int) -> a -> m b) -> Matrix a -> m (Matrix b)
+imapM = MG.imapM
+
+mapM_ :: (Context a, Monad m) => (a -> m b) -> Matrix a -> m ()
+mapM_ = MG.mapM_
+
+-- | O(m*n) Apply the monadic action to every element and its index,
+-- ignoring the results.
+imapM_ :: (Context a, Monad m) => ((Int, Int) -> a -> m b) -> Matrix a -> m ()
+imapM_ = MG.imapM_
+
+forM :: (Context a, Context b, Monad m) => Matrix a -> (a -> m b) -> m (Matrix b)
+forM = MG.forM
+
+forM_ :: (Context a, Monad m) => Matrix a -> (a -> m b) -> m ()
+forM_ = MG.forM_
+
+zipWith :: ( Context a, Context b, Context c)
+        => (a -> b -> c) -> Matrix a -> Matrix b -> Matrix c
+zipWith = MG.zipWith
+
+zipWith3 :: ( Context a, Context b, Context c, Context d)
+         => (a -> b -> c -> d) -> Matrix a -> Matrix b -> Matrix c
+         -> Matrix d
+zipWith3 = MG.zipWith3
+
+zipWith4 :: ( Context a, Context b, Context c, Context d, Context e)
+         => (a -> b -> c -> d -> e) -> Matrix a -> Matrix b -> Matrix c
+         -> Matrix d -> Matrix e
+zipWith4 = MG.zipWith4
+
+zipWith5 :: ( Context a, Context b, Context c, Context d, Context e, Context f)
+         => (a -> b -> c -> d -> e -> f) -> Matrix a -> Matrix b
+         -> Matrix c -> Matrix d -> Matrix e -> Matrix f
+zipWith5 = MG.zipWith5
+
+zipWith6 :: ( Context a, Context b, Context c, Context d, Context e, Context f
+            , Context g )
+         => (a -> b -> c -> d -> e -> f -> g) -> Matrix a -> Matrix b
+         -> Matrix c -> Matrix d -> Matrix e -> Matrix f -> Matrix g
+zipWith6 = MG.zipWith6
+
+izipWith :: ( Context a, Context b, Context c)
+         => ((Int, Int) -> a -> b -> c) -> Matrix a -> Matrix b -> Matrix c
+izipWith = MG.izipWith
+
+izipWith3 :: ( Context a, Context b, Context c, Context d)
+          => ((Int, Int) -> a -> b -> c -> d) -> Matrix a -> Matrix b
+          -> Matrix c -> Matrix d
+izipWith3 = MG.izipWith3
+
+izipWith4 :: ( Context a, Context b, Context c, Context d, Context e)
+          => ((Int, Int) -> a -> b -> c -> d -> e) -> Matrix a -> Matrix b
+          -> Matrix c -> Matrix d -> Matrix e
+izipWith4 = MG.izipWith4
+
+izipWith5 :: ( Context a, Context b, Context c, Context d, Context e, Context f)
+          => ((Int, Int) -> a -> b -> c -> d -> e -> f) -> Matrix a
+          -> Matrix b -> Matrix c -> Matrix d -> Matrix e -> Matrix f
+izipWith5 = MG.izipWith5
+
+izipWith6 :: ( Context a, Context b, Context c, Context d, Context e, Context f
+             , Context g )
+          => ((Int, Int) -> a -> b -> c -> d -> e -> f -> g) -> Matrix a
+          -> Matrix b -> Matrix c -> Matrix d -> Matrix e -> Matrix f
+          -> Matrix g
+izipWith6 = MG.izipWith6
+
+zip :: ( Context a, Context b
+       , Context (a,b) )
+    => Matrix a -> Matrix b -> Matrix (a,b)
+zip = MG.zip
+
+zip3 :: ( Context a, Context b, Context c
+        , Context (a,b,c) )
+     => Matrix a -> Matrix b -> Matrix c -> Matrix (a,b,c)
+zip3 = MG.zip3
+
+zip4 :: ( Context a, Context b, Context c, Context d
+        , Context (a,b,c,d) )
+     => Matrix a -> Matrix b -> Matrix c -> Matrix d -> Matrix (a,b,c,d)
+zip4 = MG.zip4
+
+zip5 :: ( Context a, Context b, Context c, Context d, Context e
+        , Context (a,b,c,d,e) )
+     => Matrix a -> Matrix b -> Matrix c -> Matrix d -> Matrix e
+     -> Matrix (a,b,c,d,e)
+zip5 = MG.zip5
+
+zip6 :: ( Context a, Context b, Context c, Context d, Context e, Context f
+        , Context (a,b,c,d,e,f) )
+     => Matrix a -> Matrix b -> Matrix c -> Matrix d -> Matrix e
+     -> Matrix f -> Matrix (a,b,c,d,e,f)
+zip6 = MG.zip6
+
+zipWithM :: (Context a, Context b, Context c, Monad m)
+         => (a -> b -> m c) -> Matrix a -> Matrix b -> m (Matrix c)
+zipWithM = MG.zipWithM
+
+zipWithM_ :: (Context a, Context b, Monad m)
+          => (a -> b -> m c) -> Matrix a -> Matrix b -> m ()
+zipWithM_ = MG.zipWithM_
+
+unzip :: (Context a, Context b, Context (a,b))
+      => Matrix (a,b) -> (Matrix a, Matrix b )
+unzip = MG.unzip
+
+unzip3 :: ( Context a, Context b, Context c
+          , Context (a,b,c) )
+       => Matrix (a,b,c) -> (Matrix a, Matrix b, Matrix c)
+unzip3 = MG.unzip3
+
+unzip4 :: ( Context a, Context b, Context c, Context d
+          , Context (a,b,c,d) )
+       => Matrix (a,b,c,d) -> (Matrix a, Matrix b, Matrix c, Matrix d)
+unzip4 = MG.unzip4
+
+unzip5 :: ( Context a, Context b, Context c, Context d, Context e
+          , Context (a,b,c,d,e) )
+       => Matrix (a,b,c,d,e)
+       -> (Matrix a, Matrix b, Matrix c, Matrix d, Matrix e)
+unzip5 = MG.unzip5
+
+unzip6 :: ( Context a, Context b, Context c, Context d, Context e, Context f
+          , Context (a,b,c,d,e,f) )
+       => Matrix (a,b,c,d,e,f)
+       -> (Matrix a, Matrix b, Matrix c, Matrix d, Matrix e, Matrix f)
+unzip6 = MG.unzip6
+
+sequence :: Monad m => Matrix (m a) -> m (Matrix a)
+sequence = MG.sequence
+
+sequence_ :: Monad m
+          => Matrix (m a) -> m ()
+sequence_ = MG.sequence_
+
+generate :: Context a => (Int, Int) -> ((Int, Int) -> a) -> Matrix a
+generate = MG.generate
+
+thaw :: (Context a, PrimMonad s) => Matrix a -> s (MMatrix (PrimState s) a)
+thaw = MG.thaw
+
+unsafeThaw :: (Context a, PrimMonad s) => Matrix a -> s (MMatrix (PrimState s) a)
+unsafeThaw = MG.unsafeThaw
+
+freeze :: (Context a, PrimMonad s) => MMatrix (PrimState s) a -> s (Matrix a)
+freeze = MG.freeze
+
+unsafeFreeze :: (Context a, PrimMonad s) => MMatrix (PrimState s) a -> s (Matrix a)
+unsafeFreeze = MG.unsafeFreeze
+
+create :: Context a => (forall s . ST s (MMatrix s a)) -> Matrix a
+create = MG.create
diff --git a/src/Data/Matrix/Class.hs b/src/Data/Matrix/Class.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Class.hs
@@ -0,0 +1,198 @@
+{-# LANGUAGE FlexibleContexts      #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE Rank2Types            #-}
+{-# LANGUAGE TypeFamilies          #-}
+module Data.Matrix.Class
+    ( Mutable
+    , Matrix(..)
+
+    -- * Derived mothods
+    , rows
+    , cols
+    , (!)
+    , fromVector
+    , fromList
+    , empty
+    , toList
+    , fromLists
+    , matrix
+    , fromRows
+    , takeRow
+    , toRows
+    , takeColumn
+    , toColumns
+    , toLists
+    , create
+    ) where
+
+import           Control.Monad.Primitive     (PrimMonad, PrimState)
+import           Control.Monad.ST            (ST, runST)
+import qualified Data.Vector.Generic         as G
+import           Text.Printf
+
+import qualified Data.Matrix.Class.Mutable as MM
+
+type family Mutable (m :: (* -> *) -> * -> *) :: (* -> * -> *) -> * -> * -> *
+
+class (MM.MMatrix (Mutable m) (G.Mutable v) a, G.Vector v a) => Matrix m v a where
+    dim :: m v a -> (Int, Int)
+
+    unsafeIndex :: m v a -> (Int, Int) -> a
+
+    unsafeFromVector :: (Int, Int) -> v a -> m v a
+
+    -- | Default algorithm is O((m*n) * O(unsafeIndex)).
+    flatten :: m v a -> v a
+    flatten mat = G.generate (r*c) $ \i -> unsafeIndex mat (i `divMod` c)
+      where
+        (r,c) = dim mat
+    {-# INLINE flatten #-}
+
+    -- | Extract a row. Default algorithm is O(n * O(unsafeIndex)).
+    unsafeTakeRow :: m 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 :: m 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 :: m 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 => m v a -> s ((Mutable m) (G.Mutable v) (PrimState s) a)
+
+    unsafeThaw :: PrimMonad s
+               => m v a -> s ((Mutable m) (G.Mutable v) (PrimState s) a)
+
+    freeze :: PrimMonad s
+           => (Mutable m) (G.Mutable v) (PrimState s) a -> s (m v a)
+
+    unsafeFreeze :: PrimMonad s
+                 => (Mutable m) (G.Mutable v) (PrimState s) a -> s (m v a)
+
+    {-# MINIMAL dim, unsafeIndex, unsafeFromVector, thaw, unsafeThaw, freeze, unsafeFreeze #-}
+
+-- | Derived methods
+
+-- | Return the number of rows
+rows :: Matrix m v a => m v a -> Int
+rows = fst . dim
+{-# INLINE rows #-}
+
+-- | Return the number of columns
+cols :: Matrix m v a => m v a -> Int
+cols = snd . dim
+{-# INLINE cols #-}
+
+-- | Indexing
+(!) :: Matrix m v a => m v a -> (Int, Int) -> a
+(!) mat (i,j) | i < 0 || i >= r || j < 0 || j >= c =
+                error "Index out of bounds"
+              | otherwise = unsafeIndex mat (i,j)
+  where
+    (r,c) = dim mat
+{-# INLINE (!) #-}
+
+-- | O(m*n) Create a list by concatenating rows
+toList :: Matrix m v a => m v a -> [a]
+toList = G.toList . flatten
+{-# INLINE toList #-}
+
+empty :: Matrix m v a => m v a
+empty = fromVector (0,0) G.empty
+{-# INLINE empty #-}
+
+fromVector :: Matrix m v a => (Int, Int) -> v a -> m v a
+fromVector (r,c) vec | r*c /= n = error errMsg
+                     | otherwise = unsafeFromVector (r,c) vec
+  where
+    errMsg = printf "fromVector: incorrect length (%d * %d != %d)" r c n
+    n = G.length vec
+{-# INLINE fromVector #-}
+
+fromList :: Matrix m v a => (Int, Int) -> [a] -> m v a
+fromList (r,c) = fromVector (r,c) . G.fromList
+{-# INLINE fromList #-}
+
+-- | O(m*n) Matrix construction
+matrix :: Matrix m v a
+       => Int  -- ^ number of columns
+       -> [a]  -- ^ row list
+       -> m v a
+matrix ncol xs | n `mod` ncol /= 0 = error "incorrect length"
+               | otherwise = unsafeFromVector (nrow,ncol) vec
+  where
+    vec = G.fromList xs
+    nrow = n `div` ncol
+    n = G.length vec
+{-# INLINE matrix #-}
+
+-- | O(m*n) Create matrix from list of lists, it doesn't check if the list of
+-- list is a valid matrix
+fromLists :: Matrix m v a => [[a]] -> m v a
+fromLists xs | null xs = empty
+             | otherwise = fromVector (r,c) . G.fromList . concat $ xs
+  where
+    r = length xs
+    c = length . head $ xs
+{-# INLINE fromLists #-}
+
+-- | O(m*n) Create matrix from rows
+fromRows :: Matrix m v a => [v a] -> m v a
+fromRows xs | null xs = empty
+            | otherwise = fromVector (r,c) . G.concat $ xs
+  where
+    r = length xs
+    c = G.length . head $ xs
+{-# INLINE fromRows #-}
+
+-- | Extract a row.
+takeRow :: Matrix m v a => m v a -> Int -> v a
+takeRow mat i | i < 0 || i >= r =
+                error $ printf "index out of bounds: (%d,%d)" i r
+              | otherwise = unsafeTakeRow mat i
+  where
+    (r,_) = dim mat
+{-# INLINE takeRow #-}
+
+-- | O(m) Return the rows
+toRows :: Matrix m v a => m v a -> [v a]
+toRows mat = map (unsafeTakeRow mat) [0..r-1]
+  where
+    (r,_) = dim mat
+{-# INLINE toRows #-}
+
+-- | Extract a row.
+takeColumn :: Matrix m v a => m v a -> Int -> v a
+takeColumn mat j | j < 0 || j >= c =
+                   error $ printf "index out of bounds: (%d,%d)" j c
+                 | otherwise = unsafeTakeColumn mat j
+  where
+    (_,c) = dim mat
+{-# INLINE takeColumn #-}
+
+-- | O(m*n) Return the columns
+toColumns :: Matrix m v a => m v a -> [v a]
+toColumns mat = map (unsafeTakeColumn mat) [0..c-1]
+  where
+    (_,c) = dim mat
+{-# INLINE toColumns #-}
+
+-- | O(m*n) List of lists
+toLists :: Matrix m v a => m v a -> [[a]]
+toLists = map G.toList . toRows
+{-# INLINE toLists #-}
+
+create :: Matrix m v a => (forall s . ST s ((Mutable m) (G.Mutable v) s a)) -> m v a
+create m = runST $ unsafeFreeze =<< m
+{-# INLINE create #-}
diff --git a/src/Data/Matrix/Class/Mutable.hs b/src/Data/Matrix/Class/Mutable.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Class/Mutable.hs
@@ -0,0 +1,46 @@
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE TypeFamilies          #-}
+
+module Data.Matrix.Class.Mutable
+    ( MMatrix(..)
+    , write
+    , read
+    ) where
+
+import           Control.Monad.Primitive     (PrimMonad, PrimState)
+import qualified Data.Vector.Generic.Mutable as GM
+import           Prelude                     hiding (read)
+
+class GM.MVector v a => MMatrix m v a where
+    dim ::  m v s a -> (Int, Int)
+
+    unsafeRead :: PrimMonad s => m v (PrimState s) a -> (Int, Int) -> s a
+
+    unsafeWrite :: PrimMonad s => m v (PrimState s) a -> (Int, Int) -> a -> s ()
+
+    -- | Create a mutable matrix without initialization
+    new :: PrimMonad s => (Int, Int) -> s (m v (PrimState s) a)
+
+    replicate :: PrimMonad s => (Int, Int) -> a -> s (m v (PrimState s) a)
+
+    {-# MINIMAL dim, unsafeRead, unsafeWrite, new, replicate #-}
+
+-- | Derived methods
+
+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/Dense/Generic.hs b/src/Data/Matrix/Dense/Generic.hs
deleted file mode 100644
--- a/src/Data/Matrix/Dense/Generic.hs
+++ /dev/null
@@ -1,603 +0,0 @@
-{-# LANGUAGE BangPatterns          #-}
-{-# LANGUAGE DeriveGeneric         #-}
-{-# LANGUAGE FlexibleContexts      #-}
-{-# LANGUAGE FlexibleInstances     #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE TypeFamilies          #-}
-module Data.Matrix.Dense.Generic
-    (
-    -- * Immutable Matrix
-      Matrix(..)
-
-    -- * Accessors
-    -- ** length information
-    , MG.dim
-    , MG.rows
-    , MG.cols
-
-    -- ** Indexing
-    , MG.unsafeIndex
-    , (MG.!)
-    , MG.takeRow
-    , MG.takeColumn
-    , MG.takeDiag
-
-    -- * Construction
-    , MG.unsafeFromVector
-    , MG.fromVector
-    , MG.matrix
-    , MG.fromLists
-    , MG.fromRows
-    , fromColumns
-    , MG.empty
-
-    -- * Conversions
-    , MG.flatten
-    , MG.toRows
-    , MG.toColumns
-    , MG.toList
-    , MG.toLists
-
-    -- * Different matrix types
-    , convert
-
-    , tr
-    , subMatrix
-    , ident
-    , diag
-    , diagRect
-    , fromBlocks
-    , isSymmetric
-    , force
-
-    , Data.Matrix.Dense.Generic.foldl
-
-    -- * Mapping
-    , Data.Matrix.Dense.Generic.map
-    , imap
-
-    -- * Monadic mapping
-    , mapM
-    , imapM
-    , mapM_
-    , imapM_
-    , forM
-    , forM_
-
-    -- * Zipping
-    , Data.Matrix.Dense.Generic.zipWith
-    , Data.Matrix.Dense.Generic.zipWith3
-    , zipWith4
-    , zipWith5
-    , zipWith6
-    , izipWith
-    , izipWith3
-    , izipWith4
-    , izipWith5
-    , izipWith6
-    , Data.Matrix.Dense.Generic.zip
-    , Data.Matrix.Dense.Generic.zip3
-    , zip4
-    , zip5
-    , zip6
-
-    -- * Monadic Zipping
-    , zipWithM
-    , zipWithM_
-
-    -- * Unzipping
-    , Data.Matrix.Dense.Generic.unzip
-    , Data.Matrix.Dense.Generic.unzip3
-    , unzip4
-    , unzip5
-    , unzip6
-
-    -- * Monadic sequencing
-    , Data.Matrix.Dense.Generic.sequence
-    , Data.Matrix.Dense.Generic.sequence_
-
-    , generate
-
-    -- * Mutable matrix
-    , MG.thaw
-    , MG.unsafeThaw
-    , MG.freeze
-    , MG.unsafeFreeze
-    , MG.create
-    ) where
-
-import           Control.Arrow                     ((&&&), (***))
-import           Control.DeepSeq                   hiding (force)
-import           Control.Monad                     (foldM, foldM_, liftM)
-import qualified Data.Foldable                     as F
-import qualified Data.Vector.Generic               as G
-import qualified Data.Vector.Generic.Mutable       as GM
-import           Prelude                           hiding (mapM, mapM_)
-
-import           Data.Matrix.Dense.Generic.Mutable (MMatrix (..))
-import qualified Data.Matrix.Generic               as MG
-import           GHC.Generics                      (Generic)
-
-type instance MG.Mutable Matrix = MMatrix
-
--- | Row-major matrix supporting efficient slice
-data Matrix v a = Matrix !Int    -- number of rows
-                         !Int    -- number of cols
-                         !Int    -- physical row dimension
-                         !Int    -- offset
-                         !(v a)  -- flat matrix
-    deriving (Show, Read, Eq, Generic)
-
-instance NFData (v a) => NFData (Matrix v a) where
- rnf (Matrix _ _ _ _ vec) = rnf vec
-
-instance G.Vector v a => MG.Matrix Matrix v a where
-    -- | O(1) Return the size of matrix.
-    dim (Matrix r c _ _ _) = (r,c)
-    {-# INLINE dim #-}
-
-    -- | O(1) Unsafe indexing without bound check.
-    unsafeIndex (Matrix _ _ tda offset vec) (i,j) = vec `G.unsafeIndex` idx
-      where
-        idx = offset + i * tda + j
-    {-# INLINE unsafeIndex #-}
-
-    -- | O(1) Create matrix from vector.
-    unsafeFromVector (r,c) = Matrix r c c 0
-    {-# INLINE unsafeFromVector #-}
-
-    -- | O(1) Extract a row.
-    unsafeTakeRow (Matrix _ c tda offset vec) i = G.slice i' c vec
-      where
-        i' = offset + i * tda
-    {-# INLINE unsafeTakeRow #-}
-
-    -- | Create a vector by concatenating rows.
-    flatten (Matrix r c tda offset vec)
-        | c == tda = G.slice offset (r*c) vec
-        | otherwise = G.generate (r*c) $ \i ->
-            vec `G.unsafeIndex` (offset + (i `div` c) * tda + (i `mod` c))
-    {-# INLINE flatten #-}
-
-    thaw (Matrix r c tda offset v) = MMatrix r c tda offset `liftM` G.thaw v
-    {-# INLINE thaw #-}
-
-    unsafeThaw (Matrix r c tda offset v) = MMatrix r c tda offset `liftM` G.unsafeThaw v
-    {-# INLINE unsafeThaw #-}
-
-    freeze (MMatrix r c tda offset v) = Matrix r c tda offset `liftM` G.freeze v
-    {-# INLINE freeze #-}
-
-    unsafeFreeze (MMatrix r c tda offset v) = Matrix r c tda offset `liftM` G.unsafeFreeze v
-    {-# INLINE unsafeFreeze #-}
-
---reshape :: G.Vector v a => Matrix v a -> (Int, Int) -> Matrix v a
-
--- | O(m*n) Create matrix from columns
-fromColumns :: G.Vector v a => [v a] -> Matrix v a
-fromColumns = tr . MG.fromRows
-{-# INLINE fromColumns #-}
-
----- | construct upper triangular matrix from vector
---upperTriangular :: (Num a, G.Vector v a) => Int -> v a -> Matrix v a
---upperTriangular n vec =
-
--- | O(m*n) Convert different matrix type
-convert :: (G.Vector v a, G.Vector w a) => Matrix v a -> Matrix w a
-convert (Matrix r c tda offset vec) = Matrix r c tda offset . G.convert $ vec
-{-# INLINE convert #-}
-
--- | O(1) Extract sub matrix
-subMatrix :: G.Vector v a
-          => (Int, Int)  -- ^ upper left corner of the submatrix
-          -> (Int, Int)  -- ^ bottom right corner of the submatrix
-          -> Matrix v a -> Matrix v a
-subMatrix (i,j) (i',j') (Matrix _ n tda offset vec)
-    | m' <= 0 || n' <= 0 = MG.empty
-    | otherwise = Matrix m' n' tda offset' vec
-  where
-    m' = i' - i + 1
-    n' = j' - j + 1
-    offset' = offset + i * n + j
-{-# INLINE subMatrix #-}
-
--- | O(m*n) Matrix transpose
-tr :: G.Vector v a => Matrix v a -> Matrix v a
-tr (Matrix r c tda offset vec) = MG.fromVector (c,r) $ G.generate (r*c) f
-  where
-    f i = vec G.! (offset + i `mod` r * tda + i `div` r)
-{-# INLINE tr #-}
-
--- | O(m*n) Create an identity matrix
-ident :: (Num a, G.Vector v a) => Int -> Matrix v a
-ident n = diagRect 0 (n,n) $ replicate n 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, F.Foldable t)
-     => t a  -- ^ diagonal
-     -> Matrix v a
-diag d = diagRect 0 (n,n) d
-  where n = length . F.toList $ d
-{-# INLINE diag #-}
-
--- | O(m*n) Create a rectangular matrix with default values and given diagonal
-diagRect :: (G.Vector v a, F.Foldable t)
-         => a         -- ^ default value
-         -> (Int, Int)
-         -> t a       -- ^ diagonal
-         -> Matrix v a
-diagRect z0 (r,c) d = MG.fromVector (r,c) $ G.create $ GM.replicate n z0 >>= go d c
-  where
-    go xs c' v = F.foldlM f 0 xs >> return v
-      where
-        f !i x = GM.unsafeWrite v (i*(c'+1)) x >> return (i+1)
-    n = r * c
-{-# INLINE diagRect #-}
-
-fromBlocks :: G.Vector v a
-           => a               -- ^ default value
-           -> [[Matrix v a]]
-           -> Matrix v a
-fromBlocks d ms = MG.fromVector (m,n) $ G.create $ GM.replicate (m*n) d >>= go n ms
-  where
-    go n' xss v = foldM_ f 0 xss >> return v
-      where
-        f !cr xs = do (r', _) <- foldM g (0, 0) xs
-                      return $ cr + r'
-          where
-            g (!maxR, !cc) x = do
-                let (r,c) = MG.dim x
-                    vec = MG.flatten x
-                    step i u = do
-                        GM.unsafeWrite v ((cr + i `div` c) * n' + i `mod` c + cc) u
-                        return (i+1)
-                G.foldM'_ step (0::Int) vec
-                return (max maxR r, cc + c)
-    -- figure out the dimension of the new matrix
-    (m, n) = (sum *** maximum) . Prelude.unzip . Prelude.map ((maximum *** sum) .
-                Prelude.unzip . Prelude.map (MG.rows &&& MG.cols)) $ ms
-{-# INLINE fromBlocks #-}
-
-isSymmetric :: (Eq a, G.Vector v a) => Matrix v a -> Bool
-isSymmetric m@(Matrix r c _ _ _) | r /= c = False
-                                 | otherwise = all f [0 .. r-1]
-  where
-    f i = all g [i + 1 .. c-1]
-      where g j = m MG.! (i,j) == m MG.! (j,i)
-{-# INLINE isSymmetric #-}
-
-force :: G.Vector v a => Matrix v a -> Matrix v a
-force m@(Matrix r c _ _ _) = MG.fromVector (r,c) . G.force . MG.flatten $ m
-{-# INLINE force #-}
-
-map :: (G.Vector v a, G.Vector v b) => (a -> b) -> Matrix v a -> Matrix v b
-map f m@(Matrix r c _ _ _) = MG.fromVector (r,c) $ G.map f . MG.flatten $ m
-{-# INLINE map #-}
-
-imap :: (G.Vector v a, G.Vector v b) => ((Int, Int) -> a -> b) -> Matrix v a -> Matrix v b
-imap f m@(Matrix r c _ _ _) = MG.fromVector (r,c) $ G.imap f' . MG.flatten $ m
-  where
-    f' i = f (i `div` c, i `mod` c)
-{-# INLINE imap #-}
-
-foldl :: G.Vector v b => (a -> b -> a) -> a -> Matrix v b -> a
-foldl f acc m = G.foldl f acc . MG.flatten $ m
-{-# INLINE foldl #-}
-
-mapM :: (G.Vector v a, G.Vector v b, Monad m)
-     => (a -> m b) -> Matrix v a -> m (Matrix v b)
-mapM f m@(Matrix r c _ _ _) = liftM (MG.fromVector (r,c)) $ G.mapM f $ MG.flatten m
-{-# INLINE mapM #-}
-
--- | O(m*n) Apply the monadic action to every element and its index,
--- yielding a matrix of results.
-imapM :: (G.Vector v a, G.Vector v b, Monad m)
-      => ((Int, Int) -> a -> m b) -> Matrix v a -> m (Matrix v b)
-imapM f m@(Matrix r c _ _ _) = fmap (MG.fromVector (r,c)) $ G.imapM f' $
-    MG.flatten m
-  where
-    f' i = f (i `div` c, i `mod` c)
-{-# INLINE imapM #-}
-
-mapM_ :: (G.Vector v a, Monad m) => (a -> m b) -> Matrix v a -> m ()
-mapM_ f = G.mapM_ f . MG.flatten
-{-# INLINE mapM_ #-}
-
--- | O(m*n) Apply the monadic action to every element and its index,
--- ignoring the results.
-imapM_ :: (G.Vector v a, Monad m)
-       => ((Int, Int) -> a -> m b) -> Matrix v a -> m ()
-imapM_ f m@(Matrix _ c _ _ _) = G.imapM_ f' $ MG.flatten m
-  where
-    f' i = f (i `div` c, i `mod` c)
-{-# INLINE imapM_ #-}
-
-forM :: (G.Vector v a, G.Vector v b, Monad m)
-     => Matrix v a -> (a -> m b) -> m (Matrix v b)
-forM = flip mapM
-{-# INLINE forM #-}
-
-forM_ :: (G.Vector v a, Monad m) => Matrix v a -> (a -> m b) -> m ()
-forM_ = flip mapM_
-{-# INLINE forM_ #-}
-
-zipWith :: (G.Vector v a, G.Vector v b, G.Vector v c)
-        => (a -> b -> c) -> Matrix v a -> Matrix v b -> Matrix v c
-zipWith f m1 m2
-    | MG.dim m1 /= MG.dim m2 = error "zipWith: Dimensions don't match."
-    | otherwise = MG.fromVector (MG.dim m1) $
-                  G.zipWith f (MG.flatten m1) $ MG.flatten m2
-{-# INLINE zipWith #-}
-
-zipWith3 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d)
-         => (a -> b -> c -> d) -> Matrix v a -> Matrix v b -> Matrix v c
-         -> Matrix v d
-zipWith3 f m1 m2 m3
-    | MG.dim m1 /= MG.dim m2 ||
-      MG.dim m2 /= MG.dim m3 = error "zipWith3: Dimensions don't match."
-    | otherwise = MG.fromVector (MG.dim m1) $
-                  G.zipWith3 f (MG.flatten m1) (MG.flatten m2) $ MG.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)
-         => (a -> b -> c -> d -> e) -> Matrix v a -> Matrix v b -> Matrix v c
-         -> Matrix v d -> Matrix v e
-zipWith4 f m1 m2 m3 m4
-    | MG.dim m1 /= MG.dim m2 ||
-      MG.dim m2 /= MG.dim m3 ||
-      MG.dim m3 /= MG.dim m4 = error "zipWith4: Dimensions don't match."
-    | otherwise = MG.fromVector (MG.dim m1) $
-                  G.zipWith4 f (MG.flatten m1) (MG.flatten m2)
-                  (MG.flatten m3) $ MG.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 )
-         => (a -> b -> c -> d -> e -> f) -> Matrix v a -> Matrix v b
-         -> Matrix v c -> Matrix v d -> Matrix v e -> Matrix v f
-zipWith5 f m1 m2 m3 m4 m5
-    | MG.dim m1 /= MG.dim m2 ||
-      MG.dim m2 /= MG.dim m3 ||
-      MG.dim m3 /= MG.dim m4 ||
-      MG.dim m4 /= MG.dim m5 = error "zipWith5: Dimensions don't match."
-    | otherwise = MG.fromVector (MG.dim m1) $
-                  G.zipWith5 f (MG.flatten m1) (MG.flatten m2)
-                  (MG.flatten m3) (MG.flatten m4) $ MG.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 )
-         => (a -> b -> c -> d -> e -> f -> g) -> Matrix v a -> Matrix v b
-         -> Matrix v c -> Matrix v d -> Matrix v e -> Matrix v f -> Matrix v g
-zipWith6 f m1 m2 m3 m4 m5 m6
-    | MG.dim m1 /= MG.dim m2 ||
-      MG.dim m2 /= MG.dim m3 ||
-      MG.dim m3 /= MG.dim m4 ||
-      MG.dim m4 /= MG.dim m5 ||
-      MG.dim m5 /= MG.dim m6 = error "zipWith6: Dimensions don't match."
-    | otherwise = MG.fromVector (MG.dim m1) $
-                  G.zipWith6 f (MG.flatten m1) (MG.flatten m2) (MG.flatten m3)
-                  (MG.flatten m4) (MG.flatten m5) $ MG.flatten m6
-{-# INLINE zipWith6 #-}
-
-izipWith :: (G.Vector v a, G.Vector v b, G.Vector v c)
-         => ((Int, Int) -> a -> b -> c) -> Matrix v a -> Matrix v b -> Matrix v c
-izipWith f m1 m2
-    | MG.dim m1 /= MG.dim m2 = error "izipWith: Dimensions don't match."
-    | otherwise = MG.fromVector (MG.dim m1) $
-                  G.izipWith f' (MG.flatten m1) $ MG.flatten m2
-  where
-    c = MG.cols m1
-    f' i = f (i `div` c, i `mod` c)
-{-# INLINE izipWith #-}
-
-izipWith3 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d)
-          => ((Int, Int) -> a -> b -> c -> d) -> Matrix v a -> Matrix v b
-          -> Matrix v c -> Matrix v d
-izipWith3 f m1 m2 m3
-    | MG.dim m1 /= MG.dim m2 ||
-      MG.dim m2 /= MG.dim m3 = error "izipWith3: Dimensions don't match."
-    | otherwise = MG.fromVector (MG.dim m1) $
-                  G.izipWith3 f' (MG.flatten m1) (MG.flatten m2) $ MG.flatten m3
-  where
-    c = MG.cols m1
-    f' i = f (i `div` c, i `mod` c)
-{-# INLINE izipWith3 #-}
-
-izipWith4 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d, G.Vector v e)
-          => ((Int, Int) -> a -> b -> c -> d -> e) -> Matrix v a -> Matrix v b
-          -> Matrix v c -> Matrix v d -> Matrix v e
-izipWith4 f m1 m2 m3 m4
-    | MG.dim m1 /= MG.dim m2 ||
-      MG.dim m2 /= MG.dim m3 ||
-      MG.dim m3 /= MG.dim m4 = error "izipWith4: Dimensions don't match."
-    | otherwise = MG.fromVector (MG.dim m1) $
-                  G.izipWith4 f' (MG.flatten m1) (MG.flatten m2)
-                  (MG.flatten m3) $ MG.flatten m4
-  where
-    c = MG.cols m1
-    f' i = f (i `div` c, i `mod` c)
-{-# INLINE izipWith4 #-}
-
-izipWith5 :: ( G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d
-             , G.Vector v e, G.Vector v f )
-          => ((Int, Int) -> a -> b -> c -> d -> e -> f) -> Matrix v a
-          -> Matrix v b -> Matrix v c -> Matrix v d -> Matrix v e -> Matrix v f
-izipWith5 f m1 m2 m3 m4 m5
-    | MG.dim m1 /= MG.dim m2 ||
-      MG.dim m2 /= MG.dim m3 ||
-      MG.dim m3 /= MG.dim m4 ||
-      MG.dim m4 /= MG.dim m5 = error "izipWith5: Dimensions don't match."
-    | otherwise = MG.fromVector (MG.dim m1) $
-                  G.izipWith5 f' (MG.flatten m1) (MG.flatten m2)
-                  (MG.flatten m3) (MG.flatten m4) $ MG.flatten m5
-  where
-    c = MG.cols m1
-    f' i = f (i `div` c, i `mod` c)
-{-# INLINE izipWith5 #-}
-
-izipWith6 :: ( 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 )
-          => ((Int, Int) -> a -> b -> c -> d -> e -> f -> g) -> Matrix v a
-          -> Matrix v b -> Matrix v c -> Matrix v d -> Matrix v e -> Matrix v f
-          -> Matrix v g
-izipWith6 f m1 m2 m3 m4 m5 m6
-    | MG.dim m1 /= MG.dim m2 ||
-      MG.dim m2 /= MG.dim m3 ||
-      MG.dim m3 /= MG.dim m4 ||
-      MG.dim m4 /= MG.dim m5 ||
-      MG.dim m5 /= MG.dim m6 = error "izipWith6: Dimensions don't match."
-    | otherwise = MG.fromVector (MG.dim m1) $
-                  G.izipWith6 f' (MG.flatten m1) (MG.flatten m2) (MG.flatten m3)
-                  (MG.flatten m4) (MG.flatten m5) $ MG.flatten m6
-  where
-    c = MG.cols m1
-    f' i = f (i `div` c, i `mod` c)
-{-# INLINE izipWith6 #-}
-
-
-zip :: (G.Vector v a, G.Vector v b, G.Vector v (a,b))
-    => Matrix v a -> Matrix v b -> Matrix v (a,b)
-zip m1 m2
-    | MG.dim m1 /= MG.dim m2 = error "zip: Dimensions don't match."
-    | otherwise = MG.fromVector (MG.dim m1) $
-                  G.zip (MG.flatten m1) $ MG.flatten m2
-{-# INLINE zip #-}
-
-zip3 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v (a,b,c))
-     => Matrix v a -> Matrix v b -> Matrix v c -> Matrix v (a,b,c)
-zip3 m1 m2 m3
-    | MG.dim m1 /= MG.dim m2 ||
-      MG.dim m2 /= MG.dim m3 = error "zip3: Dimensions don't match."
-    | otherwise = MG.fromVector (MG.dim m1) $
-                  G.zip3 (MG.flatten m1) (MG.flatten m2) $ MG.flatten m3
-{-# INLINE zip3 #-}
-
-zip4 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d, G.Vector v (a,b,c,d))
-     => Matrix v a -> Matrix v b -> Matrix v c -> Matrix v d -> Matrix v (a,b,c,d)
-zip4 m1 m2 m3 m4
-    | MG.dim m1 /= MG.dim m2 ||
-      MG.dim m2 /= MG.dim m3 ||
-      MG.dim m3 /= MG.dim m4 = error "zip4: Dimensions don't match."
-    | otherwise = MG.fromVector (MG.dim m1) $
-                  G.zip4 (MG.flatten m1) (MG.flatten m2)
-                  (MG.flatten m3) $ MG.flatten m4
-{-# INLINE zip4 #-}
-
-zip5 :: ( G.Vector v a, G.Vector v b, G.Vector v c
-        , G.Vector v d, G.Vector v e, G.Vector v (a,b,c,d,e) )
-     => Matrix v a -> Matrix v b -> Matrix v c -> Matrix v d -> Matrix v e
-     -> Matrix v (a,b,c,d,e)
-zip5 m1 m2 m3 m4 m5
-    | MG.dim m1 /= MG.dim m2 ||
-      MG.dim m2 /= MG.dim m3 ||
-      MG.dim m3 /= MG.dim m4 ||
-      MG.dim m4 /= MG.dim m5 = error "zip5: Dimensions don't match."
-    | otherwise = MG.fromVector (MG.dim m1) $
-                  G.zip5 (MG.flatten m1) (MG.flatten m2)
-                  (MG.flatten m3) (MG.flatten m4) $ MG.flatten m5
-{-# INLINE zip5 #-}
-
-zip6 :: ( 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 (a,b,c,d,e,f) )
-     => Matrix v a -> Matrix v b -> Matrix v c -> Matrix v d -> Matrix v e
-     -> Matrix v f -> Matrix v (a,b,c,d,e,f)
-zip6 m1 m2 m3 m4 m5 m6
-    | MG.dim m1 /= MG.dim m2 ||
-      MG.dim m2 /= MG.dim m3 ||
-      MG.dim m3 /= MG.dim m4 ||
-      MG.dim m4 /= MG.dim m5 ||
-      MG.dim m5 /= MG.dim m6 = error "zip6: Dimensions don't match."
-    | otherwise = MG.fromVector (MG.dim m1) $
-                  G.zip6 (MG.flatten m1) (MG.flatten m2) (MG.flatten m3)
-                  (MG.flatten m4) (MG.flatten m5) $ MG.flatten m6
-{-# INLINE zip6 #-}
-
-zipWithM :: (Monad m, G.Vector v a, G.Vector v b, G.Vector v c)
-         => (a -> b -> m c) -> Matrix v a -> Matrix v b -> m (Matrix v c)
-zipWithM f m1 m2
-    | MG.dim m1 /= MG.dim m2 = error "zipWithM: Dimensions don't match."
-    | otherwise = liftM (MG.fromVector $ MG.dim m1) $
-                  G.zipWithM f (MG.flatten m1) $ MG.flatten m2
-{-# INLINE zipWithM #-}
-
-zipWithM_ :: (Monad m, G.Vector v a, G.Vector v b)
-          => (a -> b -> m c) -> Matrix v a -> Matrix v b -> m ()
-zipWithM_ f m1 m2
-    | MG.dim m1 /= MG.dim m2 = error "zipWithM_: Dimensions don't match."
-    | otherwise = G.zipWithM_ f (MG.flatten m1) $ MG.flatten m2
-{-# INLINE zipWithM_ #-}
-
-unzip :: (G.Vector v a, G.Vector v b, G.Vector v (a,b))
-      => Matrix v (a,b) -> (Matrix v a, Matrix v b )
-unzip m = (MG.fromVector d v1, MG.fromVector d v2)
-  where
-    d = MG.dim m
-    (v1, v2) = G.unzip $ MG.flatten m
-{-# INLINE unzip #-}
-
-unzip3 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v (a,b,c))
-       => Matrix v (a,b, c) -> (Matrix v a, Matrix v b, Matrix v c)
-unzip3 m = (MG.fromVector d v1, MG.fromVector d v2, MG.fromVector d v3)
-  where
-    d = MG.dim m
-    (v1, v2, v3) = G.unzip3 $ MG.flatten m
-{-# INLINE unzip3 #-}
-
-unzip4 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d, G.Vector v (a,b,c,d))
-       => Matrix v (a,b,c,d) -> (Matrix v a, Matrix v b, Matrix v c, Matrix v d)
-unzip4 m = ( MG.fromVector d v1
-           , MG.fromVector d v2
-           , MG.fromVector d v3
-           , MG.fromVector d v4
-           )
-  where
-    d = MG.dim m
-    (v1, v2, v3, v4) = G.unzip4 $ MG.flatten m
-{-# INLINE unzip4 #-}
-
-unzip5 :: ( G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d
-          , G.Vector v e, G.Vector v (a,b,c,d,e) )
-       => Matrix v (a,b,c,d,e)
-       -> (Matrix v a, Matrix v b, Matrix v c, Matrix v d, Matrix v e)
-unzip5 m = ( MG.fromVector d v1
-           , MG.fromVector d v2
-           , MG.fromVector d v3
-           , MG.fromVector d v4
-           , MG.fromVector d v5
-           )
-  where
-    d = MG.dim m
-    (v1, v2, v3, v4, v5) = G.unzip5 $ MG.flatten m
-{-# INLINE unzip5 #-}
-
-unzip6 :: ( 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 (a,b,c,d,e,f) )
-       => Matrix v (a,b,c,d,e,f)
-       -> (Matrix v a, Matrix v b, Matrix v c, Matrix v d, Matrix v e, Matrix v f)
-unzip6 m = ( MG.fromVector d v1
-           , MG.fromVector d v2
-           , MG.fromVector d v3
-           , MG.fromVector d v4
-           , MG.fromVector d v5
-           , MG.fromVector d v6
-           )
-  where
-    d = MG.dim m
-    (v1, v2, v3, v4, v5, v6) = G.unzip6 $ MG.flatten m
-{-# INLINE unzip6 #-}
-
-sequence :: (G.Vector v a, G.Vector v (m a), Monad m)
-         => Matrix v (m a) -> m (Matrix v a)
-sequence (Matrix r c tda offset vec) = liftM (Matrix r c tda offset) . G.sequence $ vec
-{-# INLINE sequence #-}
-
-sequence_ :: (G.Vector v (m a), Monad m)
-          => Matrix v (m a) -> m ()
-sequence_ (Matrix _ _ _ _ vec) = G.sequence_ vec
-{-# INLINE sequence_ #-}
-
-generate :: G.Vector v a => (Int, Int) -> ((Int, Int) -> a) -> Matrix v a
-generate (r,c) f = MG.fromVector (r,c) . G.generate (r*c) $ \i -> f (i `div` c, i `mod` c)
-{-# INLINE generate #-}
diff --git a/src/Data/Matrix/Dense/Generic/Mutable.hs b/src/Data/Matrix/Dense/Generic/Mutable.hs
deleted file mode 100644
--- a/src/Data/Matrix/Dense/Generic/Mutable.hs
+++ /dev/null
@@ -1,52 +0,0 @@
-{-# LANGUAGE FlexibleInstances     #-}
-{-# LANGUAGE FlexibleContexts      #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-module Data.Matrix.Dense.Generic.Mutable
-   ( -- * Mutable Matrix
-     MMatrix(..)
-   , C.dim
-   , takeRow
-   , C.write
-   , C.unsafeWrite
-   , C.read
-   , C.unsafeRead
-   , C.new
-   , C.replicate
-   ) where
-
-import           Control.Monad               (liftM)
-import           Control.DeepSeq
-import qualified Data.Vector.Generic.Mutable as GM
-import           Prelude                     hiding (read, replicate)
-
-import qualified Data.Matrix.Generic.Mutable as C
-
--- | mutable matrix
-data MMatrix v s a = MMatrix !Int !Int !Int !Int !(v s a)
-
-instance (NFData (v s a)) => NFData (MMatrix v s a) where
- rnf (MMatrix _ _ _ _ vec) = rnf vec
-
-instance GM.MVector v a => C.MMatrix MMatrix v a where
-    dim (MMatrix r c _ _ _) = (r,c)
-    {-# INLINE dim #-}
-
-    unsafeRead (MMatrix _ _ tda offset v) (i,j) = GM.unsafeRead v idx
-      where idx = offset + i * tda + j
-    {-# INLINE unsafeRead #-}
-
-    unsafeWrite (MMatrix _ _ tda offset v) (i,j) = GM.unsafeWrite v idx
-      where idx = offset + i * tda + j
-    {-# INLINE unsafeWrite #-}
-
-    new (r,c) = MMatrix r c c 0 `liftM` GM.new (r*c)
-    {-# INLINE new #-}
-
-    replicate (r,c) x = MMatrix r c c 0 `liftM` GM.replicate (r*c) x
-    {-# INLINE replicate #-}
-
-takeRow :: GM.MVector v a => MMatrix v m a -> Int -> v m a
-takeRow (MMatrix _ c tda offset vec) i = GM.slice i' c vec
-  where
-    i' = offset + i * tda
-{-# INLINE takeRow #-}
diff --git a/src/Data/Matrix/Generic.hs b/src/Data/Matrix/Generic.hs
--- a/src/Data/Matrix/Generic.hs
+++ b/src/Data/Matrix/Generic.hs
@@ -1,198 +1,607 @@
+{-# LANGUAGE BangPatterns          #-}
+{-# LANGUAGE DeriveGeneric         #-}
 {-# LANGUAGE FlexibleContexts      #-}
+{-# LANGUAGE FlexibleInstances     #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE Rank2Types            #-}
 {-# LANGUAGE TypeFamilies          #-}
 module Data.Matrix.Generic
-    ( Mutable
-    , Matrix(..)
+    (
+    -- * Immutable Matrix
+      Matrix(..)
 
-    -- * Derived mothods
-    , rows
-    , cols
-    , (!)
-    , fromVector
-    , fromList
-    , empty
-    , toList
-    , fromLists
-    , matrix
-    , fromRows
-    , takeRow
-    , toRows
-    , takeColumn
-    , toColumns
-    , toLists
-    , create
+    -- * Accessors
+    -- ** length information
+    , MG.dim
+    , MG.rows
+    , MG.cols
+
+    -- ** Indexing
+    , MG.unsafeIndex
+    , (MG.!)
+    , MG.takeRow
+    , MG.takeColumn
+    , MG.takeDiag
+
+    -- * Construction
+    , MG.unsafeFromVector
+    , MG.fromVector
+    , MG.matrix
+    , MG.fromList
+    , MG.fromLists
+    , MG.fromRows
+    , fromColumns
+    , MG.empty
+
+    -- * Conversions
+    , MG.flatten
+    , MG.toRows
+    , MG.toColumns
+    , MG.toList
+    , MG.toLists
+
+    -- * Different matrix types
+    , convert
+
+    , tr
+    , subMatrix
+    , ident
+    , diag
+    , diagRect
+    , fromBlocks
+    , isSymmetric
+    , force
+
+    , Data.Matrix.Generic.foldl
+
+    -- * Mapping
+    , Data.Matrix.Generic.map
+    , imap
+
+    -- * Monadic mapping
+    , mapM
+    , imapM
+    , mapM_
+    , imapM_
+    , forM
+    , forM_
+
+    -- * Zipping
+    , Data.Matrix.Generic.zipWith
+    , Data.Matrix.Generic.zipWith3
+    , zipWith4
+    , zipWith5
+    , zipWith6
+    , izipWith
+    , izipWith3
+    , izipWith4
+    , izipWith5
+    , izipWith6
+    , Data.Matrix.Generic.zip
+    , Data.Matrix.Generic.zip3
+    , zip4
+    , zip5
+    , zip6
+
+    -- * Monadic Zipping
+    , zipWithM
+    , zipWithM_
+
+    -- * Unzipping
+    , Data.Matrix.Generic.unzip
+    , Data.Matrix.Generic.unzip3
+    , unzip4
+    , unzip5
+    , unzip6
+
+    -- * Monadic sequencing
+    , Data.Matrix.Generic.sequence
+    , Data.Matrix.Generic.sequence_
+
+    , generate
+
+    -- * Mutable matrix
+    , MG.thaw
+    , MG.unsafeThaw
+    , MG.freeze
+    , MG.unsafeFreeze
+    , MG.create
     ) where
 
-import           Control.Monad.Primitive     (PrimMonad, PrimState)
-import           Control.Monad.ST            (ST, runST)
-import qualified Data.Vector.Generic         as G
-import           Text.Printf
+import           Control.Arrow                     ((&&&), (***))
+import           Control.DeepSeq                   hiding (force)
+import           Control.Monad                     (foldM, foldM_, liftM)
+import qualified Data.Foldable                     as F
+import qualified Data.Vector.Generic               as G
+import qualified Data.Vector.Generic.Mutable       as GM
+import           Prelude                           hiding (mapM, mapM_)
 
-import qualified Data.Matrix.Generic.Mutable as MM
+import           Data.Matrix.Generic.Mutable (MMatrix (..))
+import qualified Data.Matrix.Class as MG
+import           GHC.Generics                      (Generic)
 
-type family Mutable (m :: (* -> *) -> * -> *) :: (* -> * -> *) -> * -> * -> *
+type instance MG.Mutable Matrix = MMatrix
 
-class (MM.MMatrix (Mutable m) (G.Mutable v) a, G.Vector v a) => Matrix m v a where
-    dim :: m v a -> (Int, Int)
+-- | Row-major matrix supporting efficient slice.
+data Matrix v a = Matrix !Int    -- number of rows
+                         !Int    -- number of cols
+                         !Int    -- physical row dimension
+                         !Int    -- offset
+                         !(v a)  -- flat matrix
+    deriving (Show, Read, Generic)
 
-    unsafeIndex :: m v a -> (Int, Int) -> a
+instance (G.Vector v a, Eq (v a)) => Eq (Matrix v a) where
+    (==) m1 m2 = MG.flatten m1 == MG.flatten m2
 
-    unsafeFromVector :: (Int, Int) -> v a -> m v a
+instance NFData (v a) => NFData (Matrix v a) where
+    rnf (Matrix _ _ _ _ vec) = rnf vec
 
-    -- | Default algorithm is O((m*n) * O(unsafeIndex)).
-    flatten :: m v a -> v a
-    flatten mat = G.generate (r*c) $ \i -> unsafeIndex mat (i `div` c, i `mod` c)
+instance G.Vector v a => MG.Matrix Matrix v a where
+    -- | O(1) Return the size of matrix.
+    dim (Matrix r c _ _ _) = (r,c)
+    {-# INLINE dim #-}
+
+    -- | O(1) Unsafe indexing without bound check.
+    unsafeIndex (Matrix _ _ tda offset vec) (i,j) = vec `G.unsafeIndex` idx
       where
-        (r,c) = dim mat
-    {-# INLINE flatten #-}
+        idx = offset + i * tda + j
+    {-# INLINE unsafeIndex #-}
 
-    -- | Extract a row. Default algorithm is O(n * O(unsafeIndex)).
-    unsafeTakeRow :: m v a -> Int -> v a
-    unsafeTakeRow mat i = G.generate c $ \j -> unsafeIndex mat (i,j)
+    -- | O(1) Create matrix from vector.
+    unsafeFromVector (r,c) = Matrix r c c 0
+    {-# INLINE unsafeFromVector #-}
+
+    -- | O(1) Extract a row.
+    unsafeTakeRow (Matrix _ c tda offset vec) i = G.slice i' c vec
       where
-        (_,c) = dim mat
+        i' = offset + i * tda
     {-# INLINE unsafeTakeRow #-}
 
-    -- | Extract a column. Default algorithm is O(m * O(unsafeIndex)).
-    unsafeTakeColumn :: m v a -> Int -> v a
-    unsafeTakeColumn mat j = G.generate r $ \i -> unsafeIndex mat (i,j)
+    -- | Create a vector by concatenating rows.
+    flatten (Matrix r c tda offset vec)
+        | c == tda = G.slice offset (r*c) vec
+        | otherwise = G.generate (r*c) $ \i ->
+            vec `G.unsafeIndex` (offset + (i `div` c) * tda + (i `mod` c))
+    {-# INLINE flatten #-}
+
+    thaw (Matrix r c tda offset v) = MMatrix r c tda offset `liftM` G.thaw v
+    {-# INLINE thaw #-}
+
+    unsafeThaw (Matrix r c tda offset v) = MMatrix r c tda offset `liftM` G.unsafeThaw v
+    {-# INLINE unsafeThaw #-}
+
+    freeze (MMatrix r c tda offset v) = Matrix r c tda offset `liftM` G.freeze v
+    {-# INLINE freeze #-}
+
+    unsafeFreeze (MMatrix r c tda offset v) = Matrix r c tda offset `liftM` G.unsafeFreeze v
+    {-# INLINE unsafeFreeze #-}
+
+--reshape :: G.Vector v a => Matrix v a -> (Int, Int) -> Matrix v a
+
+-- | O(m*n) Create matrix from columns
+fromColumns :: G.Vector v a => [v a] -> Matrix v a
+fromColumns = tr . MG.fromRows
+{-# INLINE fromColumns #-}
+
+---- | construct upper triangular matrix from vector
+--upperTriangular :: (Num a, G.Vector v a) => Int -> v a -> Matrix v a
+--upperTriangular n vec =
+
+-- | O(m*n) Convert different matrix type
+convert :: (G.Vector v a, G.Vector w a) => Matrix v a -> Matrix w a
+convert (Matrix r c tda offset vec) = Matrix r c tda offset . G.convert $ vec
+{-# INLINE convert #-}
+
+-- | O(1) Extract sub matrix
+subMatrix :: G.Vector v a
+          => (Int, Int)  -- ^ upper left corner of the submatrix
+          -> (Int, Int)  -- ^ bottom right corner of the submatrix
+          -> Matrix v a -> Matrix v a
+subMatrix (i,j) (i',j') (Matrix _ _ tda offset vec)
+    | m' <= 0 || n' <= 0 = MG.empty
+    | otherwise = Matrix m' n' tda offset' vec
+  where
+    m' = i' - i + 1
+    n' = j' - j + 1
+    offset' = offset + i * tda + j
+{-# INLINE subMatrix #-}
+
+-- | O(m*n) Matrix transpose
+tr :: G.Vector v a => Matrix v a -> Matrix v a
+tr (Matrix r c tda offset vec) = MG.fromVector (c,r) $ G.generate (r*c) f
+  where
+    f i = vec G.! (offset + i `mod` r * tda + i `div` r)
+{-# INLINE tr #-}
+
+-- | O(m*n) Create an identity matrix
+ident :: (Num a, G.Vector v a) => Int -> Matrix v a
+ident n = diagRect 0 (n,n) $ replicate n 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, F.Foldable t)
+     => t a  -- ^ diagonal
+     -> Matrix v a
+diag d = diagRect 0 (n,n) d
+  where n = length . F.toList $ d
+{-# INLINE diag #-}
+
+-- | O(m*n) Create a rectangular matrix with default values and given diagonal
+diagRect :: (G.Vector v a, F.Foldable t)
+         => a         -- ^ default value
+         -> (Int, Int)
+         -> t a       -- ^ diagonal
+         -> Matrix v a
+diagRect z0 (r,c) d = MG.fromVector (r,c) $ G.create $ GM.replicate n z0 >>= go d c
+  where
+    go xs c' v = F.foldlM f 0 xs >> return v
       where
-        (r,_) = dim mat
-    {-# INLINE unsafeTakeColumn #-}
+        f !i x = GM.unsafeWrite v (i*(c'+1)) x >> return (i+1)
+    n = r * c
+{-# INLINE diagRect #-}
 
-    -- | Extract the diagonal. Default algorithm is O(min(m,n) * O(unsafeIndex)).
-    takeDiag :: m v a -> v a
-    takeDiag mat = G.generate n $ \i -> unsafeIndex mat (i,i)
+fromBlocks :: G.Vector v a
+           => a               -- ^ default value
+           -> [[Matrix v a]]
+           -> Matrix v a
+fromBlocks d ms = MG.fromVector (m,n) $ G.create $ GM.replicate (m*n) d >>= go n ms
+  where
+    go n' xss v = foldM_ f 0 xss >> return v
       where
-        n = uncurry min . dim $ mat
-    {-# INLINE takeDiag #-}
+        f !cr xs = do (r', _) <- foldM g (0, 0) xs
+                      return $ cr + r'
+          where
+            g (!maxR, !cc) x = do
+                let (r,c) = MG.dim x
+                    vec = MG.flatten x
+                    step i u = do
+                        GM.unsafeWrite v ((cr + i `div` c) * n' + i `mod` c + cc) u
+                        return (i+1)
+                G.foldM'_ step (0::Int) vec
+                return (max maxR r, cc + c)
+    -- figure out the dimension of the new matrix
+    (m, n) = (sum *** maximum) . Prelude.unzip . Prelude.map ((maximum *** sum) .
+                Prelude.unzip . Prelude.map (MG.rows &&& MG.cols)) $ ms
+{-# INLINE fromBlocks #-}
 
-    thaw :: PrimMonad s => m v a -> s ((Mutable m) (G.Mutable v) (PrimState s) a)
+isSymmetric :: (Eq a, G.Vector v a) => Matrix v a -> Bool
+isSymmetric m@(Matrix r c _ _ _) | r /= c = False
+                                 | otherwise = all f [0 .. r-1]
+  where
+    f i = all g [i + 1 .. c-1]
+      where g j = m MG.! (i,j) == m MG.! (j,i)
+{-# INLINE isSymmetric #-}
 
-    unsafeThaw :: PrimMonad s
-               => m v a -> s ((Mutable m) (G.Mutable v) (PrimState s) a)
+force :: G.Vector v a => Matrix v a -> Matrix v a
+force m@(Matrix r c _ _ _) = MG.fromVector (r,c) . G.force . MG.flatten $ m
+{-# INLINE force #-}
 
-    freeze :: PrimMonad s
-           => (Mutable m) (G.Mutable v) (PrimState s) a -> s (m v a)
+map :: (G.Vector v a, G.Vector v b) => (a -> b) -> Matrix v a -> Matrix v b
+map f m@(Matrix r c _ _ _) = MG.fromVector (r,c) $ G.map f . MG.flatten $ m
+{-# INLINE map #-}
 
-    unsafeFreeze :: PrimMonad s
-                 => (Mutable m) (G.Mutable v) (PrimState s) a -> s (m v a)
+imap :: (G.Vector v a, G.Vector v b) => ((Int, Int) -> a -> b) -> Matrix v a -> Matrix v b
+imap f m@(Matrix r c _ _ _) = MG.fromVector (r,c) $ G.imap f' . MG.flatten $ m
+  where
+    f' i = f (i `divMod` c)
+{-# INLINE imap #-}
 
-    {-# MINIMAL dim, unsafeIndex, unsafeFromVector, thaw, unsafeThaw, freeze, unsafeFreeze #-}
+foldl :: G.Vector v b => (a -> b -> a) -> a -> Matrix v b -> a
+foldl f acc m = G.foldl f acc . MG.flatten $ m
+{-# INLINE foldl #-}
 
--- | Derived methods
+mapM :: (G.Vector v a, G.Vector v b, Monad m)
+     => (a -> m b) -> Matrix v a -> m (Matrix v b)
+mapM f m@(Matrix r c _ _ _) = liftM (MG.fromVector (r,c)) $ G.mapM f $ MG.flatten m
+{-# INLINE mapM #-}
 
--- | Return the number of rows
-rows :: Matrix m v a => m v a -> Int
-rows = fst . dim
-{-# INLINE rows #-}
+-- | O(m*n) Apply the monadic action to every element and its index,
+-- yielding a matrix of results.
+imapM :: (G.Vector v a, G.Vector v b, Monad m)
+      => ((Int, Int) -> a -> m b) -> Matrix v a -> m (Matrix v b)
+imapM f m@(Matrix r c _ _ _) = fmap (MG.fromVector (r,c)) $ G.imapM f' $
+    MG.flatten m
+  where
+    f' i = f (i `divMod` c)
+{-# INLINE imapM #-}
 
--- | Return the number of columns
-cols :: Matrix m v a => m v a -> Int
-cols = snd . dim
-{-# INLINE cols #-}
+mapM_ :: (G.Vector v a, Monad m) => (a -> m b) -> Matrix v a -> m ()
+mapM_ f = G.mapM_ f . MG.flatten
+{-# INLINE mapM_ #-}
 
--- | Indexing
-(!) :: Matrix m v a => m v a -> (Int, Int) -> a
-(!) mat (i,j) | i < 0 || i >= r || j < 0 || j >= c =
-                error "Index out of bounds"
-              | otherwise = unsafeIndex mat (i,j)
+-- | O(m*n) Apply the monadic action to every element and its index,
+-- ignoring the results.
+imapM_ :: (G.Vector v a, Monad m)
+       => ((Int, Int) -> a -> m b) -> Matrix v a -> m ()
+imapM_ f m@(Matrix _ c _ _ _) = G.imapM_ f' $ MG.flatten m
   where
-    (r,c) = dim mat
-{-# INLINE (!) #-}
+    f' i = f (i `divMod` c)
+{-# INLINE imapM_ #-}
 
--- | O(m*n) Create a list by concatenating rows
-toList :: Matrix m v a => m v a -> [a]
-toList = G.toList . flatten
-{-# INLINE toList #-}
+forM :: (G.Vector v a, G.Vector v b, Monad m)
+     => Matrix v a -> (a -> m b) -> m (Matrix v b)
+forM = flip mapM
+{-# INLINE forM #-}
 
-empty :: Matrix m v a => m v a
-empty = fromVector (0,0) G.empty
-{-# INLINE empty #-}
+forM_ :: (G.Vector v a, Monad m) => Matrix v a -> (a -> m b) -> m ()
+forM_ = flip mapM_
+{-# INLINE forM_ #-}
 
-fromVector :: Matrix m v a => (Int, Int) -> v a -> m v a
-fromVector (r,c) vec | r*c /= n = error errMsg
-                     | otherwise = unsafeFromVector (r,c) vec
+zipWith :: (G.Vector v a, G.Vector v b, G.Vector v c)
+        => (a -> b -> c) -> Matrix v a -> Matrix v b -> Matrix v c
+zipWith f m1 m2
+    | MG.dim m1 /= MG.dim m2 = error "zipWith: Dimensions don't match."
+    | otherwise = MG.fromVector (MG.dim m1) $
+                  G.zipWith f (MG.flatten m1) $ MG.flatten m2
+{-# INLINE zipWith #-}
+
+zipWith3 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d)
+         => (a -> b -> c -> d) -> Matrix v a -> Matrix v b -> Matrix v c
+         -> Matrix v d
+zipWith3 f m1 m2 m3
+    | MG.dim m1 /= MG.dim m2 ||
+      MG.dim m2 /= MG.dim m3 = error "zipWith3: Dimensions don't match."
+    | otherwise = MG.fromVector (MG.dim m1) $
+                  G.zipWith3 f (MG.flatten m1) (MG.flatten m2) $ MG.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)
+         => (a -> b -> c -> d -> e) -> Matrix v a -> Matrix v b -> Matrix v c
+         -> Matrix v d -> Matrix v e
+zipWith4 f m1 m2 m3 m4
+    | MG.dim m1 /= MG.dim m2 ||
+      MG.dim m2 /= MG.dim m3 ||
+      MG.dim m3 /= MG.dim m4 = error "zipWith4: Dimensions don't match."
+    | otherwise = MG.fromVector (MG.dim m1) $
+                  G.zipWith4 f (MG.flatten m1) (MG.flatten m2)
+                  (MG.flatten m3) $ MG.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 )
+         => (a -> b -> c -> d -> e -> f) -> Matrix v a -> Matrix v b
+         -> Matrix v c -> Matrix v d -> Matrix v e -> Matrix v f
+zipWith5 f m1 m2 m3 m4 m5
+    | MG.dim m1 /= MG.dim m2 ||
+      MG.dim m2 /= MG.dim m3 ||
+      MG.dim m3 /= MG.dim m4 ||
+      MG.dim m4 /= MG.dim m5 = error "zipWith5: Dimensions don't match."
+    | otherwise = MG.fromVector (MG.dim m1) $
+                  G.zipWith5 f (MG.flatten m1) (MG.flatten m2)
+                  (MG.flatten m3) (MG.flatten m4) $ MG.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 )
+         => (a -> b -> c -> d -> e -> f -> g) -> Matrix v a -> Matrix v b
+         -> Matrix v c -> Matrix v d -> Matrix v e -> Matrix v f -> Matrix v g
+zipWith6 f m1 m2 m3 m4 m5 m6
+    | MG.dim m1 /= MG.dim m2 ||
+      MG.dim m2 /= MG.dim m3 ||
+      MG.dim m3 /= MG.dim m4 ||
+      MG.dim m4 /= MG.dim m5 ||
+      MG.dim m5 /= MG.dim m6 = error "zipWith6: Dimensions don't match."
+    | otherwise = MG.fromVector (MG.dim m1) $
+                  G.zipWith6 f (MG.flatten m1) (MG.flatten m2) (MG.flatten m3)
+                  (MG.flatten m4) (MG.flatten m5) $ MG.flatten m6
+{-# INLINE zipWith6 #-}
+
+izipWith :: (G.Vector v a, G.Vector v b, G.Vector v c)
+         => ((Int, Int) -> a -> b -> c) -> Matrix v a -> Matrix v b -> Matrix v c
+izipWith f m1 m2
+    | MG.dim m1 /= MG.dim m2 = error "izipWith: Dimensions don't match."
+    | otherwise = MG.fromVector (MG.dim m1) $
+                  G.izipWith f' (MG.flatten m1) $ MG.flatten m2
   where
-    errMsg = printf "fromVector: incorrect length (%d * %d != %d)" r c n
-    n = G.length vec
-{-# INLINE fromVector #-}
+    c = MG.cols m1
+    f' i = f (i `divMod` c)
+{-# INLINE izipWith #-}
 
-fromList :: Matrix m v a => (Int, Int) -> [a] -> m v a
-fromList (r,c) = fromVector (r,c) . G.fromList
-{-# INLINE fromList #-}
+izipWith3 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d)
+          => ((Int, Int) -> a -> b -> c -> d) -> Matrix v a -> Matrix v b
+          -> Matrix v c -> Matrix v d
+izipWith3 f m1 m2 m3
+    | MG.dim m1 /= MG.dim m2 ||
+      MG.dim m2 /= MG.dim m3 = error "izipWith3: Dimensions don't match."
+    | otherwise = MG.fromVector (MG.dim m1) $
+                  G.izipWith3 f' (MG.flatten m1) (MG.flatten m2) $ MG.flatten m3
+  where
+    c = MG.cols m1
+    f' i = f (i `divMod` c)
+{-# INLINE izipWith3 #-}
 
--- | O(m*n) Matrix construction
-matrix :: Matrix m v a
-       => Int  -- ^ number of columns
-       -> [a]  -- ^ row list
-       -> m v a
-matrix ncol xs | n `mod` ncol /= 0 = error "incorrect length"
-               | otherwise = unsafeFromVector (nrow,ncol) vec
+izipWith4 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d, G.Vector v e)
+          => ((Int, Int) -> a -> b -> c -> d -> e) -> Matrix v a -> Matrix v b
+          -> Matrix v c -> Matrix v d -> Matrix v e
+izipWith4 f m1 m2 m3 m4
+    | MG.dim m1 /= MG.dim m2 ||
+      MG.dim m2 /= MG.dim m3 ||
+      MG.dim m3 /= MG.dim m4 = error "izipWith4: Dimensions don't match."
+    | otherwise = MG.fromVector (MG.dim m1) $
+                  G.izipWith4 f' (MG.flatten m1) (MG.flatten m2)
+                  (MG.flatten m3) $ MG.flatten m4
   where
-    vec = G.fromList xs
-    nrow = n `div` ncol
-    n = G.length vec
-{-# INLINE matrix #-}
+    c = MG.cols m1
+    f' i = f (i `divMod` c)
+{-# INLINE izipWith4 #-}
 
--- | O(m*n) Create matrix from list of lists, it doesn't check if the list of
--- list is a valid matrix
-fromLists :: Matrix m v a => [[a]] -> m v a
-fromLists xs | null xs = empty
-             | otherwise = fromVector (r,c) . G.fromList . concat $ xs
+izipWith5 :: ( G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d
+             , G.Vector v e, G.Vector v f )
+          => ((Int, Int) -> a -> b -> c -> d -> e -> f) -> Matrix v a
+          -> Matrix v b -> Matrix v c -> Matrix v d -> Matrix v e -> Matrix v f
+izipWith5 f m1 m2 m3 m4 m5
+    | MG.dim m1 /= MG.dim m2 ||
+      MG.dim m2 /= MG.dim m3 ||
+      MG.dim m3 /= MG.dim m4 ||
+      MG.dim m4 /= MG.dim m5 = error "izipWith5: Dimensions don't match."
+    | otherwise = MG.fromVector (MG.dim m1) $
+                  G.izipWith5 f' (MG.flatten m1) (MG.flatten m2)
+                  (MG.flatten m3) (MG.flatten m4) $ MG.flatten m5
   where
-    r = length xs
-    c = length . head $ xs
-{-# INLINE fromLists #-}
+    c = MG.cols m1
+    f' i = f (i `divMod` c)
+{-# INLINE izipWith5 #-}
 
--- | O(m*n) Create matrix from rows
-fromRows :: Matrix m v a => [v a] -> m v a
-fromRows xs | null xs = empty
-            | otherwise = fromVector (r,c) . G.concat $ xs
+izipWith6 :: ( 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 )
+          => ((Int, Int) -> a -> b -> c -> d -> e -> f -> g) -> Matrix v a
+          -> Matrix v b -> Matrix v c -> Matrix v d -> Matrix v e -> Matrix v f
+          -> Matrix v g
+izipWith6 f m1 m2 m3 m4 m5 m6
+    | MG.dim m1 /= MG.dim m2 ||
+      MG.dim m2 /= MG.dim m3 ||
+      MG.dim m3 /= MG.dim m4 ||
+      MG.dim m4 /= MG.dim m5 ||
+      MG.dim m5 /= MG.dim m6 = error "izipWith6: Dimensions don't match."
+    | otherwise = MG.fromVector (MG.dim m1) $
+                  G.izipWith6 f' (MG.flatten m1) (MG.flatten m2) (MG.flatten m3)
+                  (MG.flatten m4) (MG.flatten m5) $ MG.flatten m6
   where
-    r = length xs
-    c = G.length . head $ xs
-{-# INLINE fromRows #-}
+    c = MG.cols m1
+    f' i = f (i `divMod` c)
+{-# INLINE izipWith6 #-}
 
--- | Extract a row.
-takeRow :: Matrix m v a => m v a -> Int -> v a
-takeRow mat i | i < 0 || i >= r =
-                error $ printf "index out of bounds: (%d,%d)" i r
-              | otherwise = unsafeTakeRow mat i
+
+zip :: (G.Vector v a, G.Vector v b, G.Vector v (a,b))
+    => Matrix v a -> Matrix v b -> Matrix v (a,b)
+zip m1 m2
+    | MG.dim m1 /= MG.dim m2 = error "zip: Dimensions don't match."
+    | otherwise = MG.fromVector (MG.dim m1) $
+                  G.zip (MG.flatten m1) $ MG.flatten m2
+{-# INLINE zip #-}
+
+zip3 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v (a,b,c))
+     => Matrix v a -> Matrix v b -> Matrix v c -> Matrix v (a,b,c)
+zip3 m1 m2 m3
+    | MG.dim m1 /= MG.dim m2 ||
+      MG.dim m2 /= MG.dim m3 = error "zip3: Dimensions don't match."
+    | otherwise = MG.fromVector (MG.dim m1) $
+                  G.zip3 (MG.flatten m1) (MG.flatten m2) $ MG.flatten m3
+{-# INLINE zip3 #-}
+
+zip4 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d, G.Vector v (a,b,c,d))
+     => Matrix v a -> Matrix v b -> Matrix v c -> Matrix v d -> Matrix v (a,b,c,d)
+zip4 m1 m2 m3 m4
+    | MG.dim m1 /= MG.dim m2 ||
+      MG.dim m2 /= MG.dim m3 ||
+      MG.dim m3 /= MG.dim m4 = error "zip4: Dimensions don't match."
+    | otherwise = MG.fromVector (MG.dim m1) $
+                  G.zip4 (MG.flatten m1) (MG.flatten m2)
+                  (MG.flatten m3) $ MG.flatten m4
+{-# INLINE zip4 #-}
+
+zip5 :: ( G.Vector v a, G.Vector v b, G.Vector v c
+        , G.Vector v d, G.Vector v e, G.Vector v (a,b,c,d,e) )
+     => Matrix v a -> Matrix v b -> Matrix v c -> Matrix v d -> Matrix v e
+     -> Matrix v (a,b,c,d,e)
+zip5 m1 m2 m3 m4 m5
+    | MG.dim m1 /= MG.dim m2 ||
+      MG.dim m2 /= MG.dim m3 ||
+      MG.dim m3 /= MG.dim m4 ||
+      MG.dim m4 /= MG.dim m5 = error "zip5: Dimensions don't match."
+    | otherwise = MG.fromVector (MG.dim m1) $
+                  G.zip5 (MG.flatten m1) (MG.flatten m2)
+                  (MG.flatten m3) (MG.flatten m4) $ MG.flatten m5
+{-# INLINE zip5 #-}
+
+zip6 :: ( 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 (a,b,c,d,e,f) )
+     => Matrix v a -> Matrix v b -> Matrix v c -> Matrix v d -> Matrix v e
+     -> Matrix v f -> Matrix v (a,b,c,d,e,f)
+zip6 m1 m2 m3 m4 m5 m6
+    | MG.dim m1 /= MG.dim m2 ||
+      MG.dim m2 /= MG.dim m3 ||
+      MG.dim m3 /= MG.dim m4 ||
+      MG.dim m4 /= MG.dim m5 ||
+      MG.dim m5 /= MG.dim m6 = error "zip6: Dimensions don't match."
+    | otherwise = MG.fromVector (MG.dim m1) $
+                  G.zip6 (MG.flatten m1) (MG.flatten m2) (MG.flatten m3)
+                  (MG.flatten m4) (MG.flatten m5) $ MG.flatten m6
+{-# INLINE zip6 #-}
+
+zipWithM :: (Monad m, G.Vector v a, G.Vector v b, G.Vector v c)
+         => (a -> b -> m c) -> Matrix v a -> Matrix v b -> m (Matrix v c)
+zipWithM f m1 m2
+    | MG.dim m1 /= MG.dim m2 = error "zipWithM: Dimensions don't match."
+    | otherwise = liftM (MG.fromVector $ MG.dim m1) $
+                  G.zipWithM f (MG.flatten m1) $ MG.flatten m2
+{-# INLINE zipWithM #-}
+
+zipWithM_ :: (Monad m, G.Vector v a, G.Vector v b)
+          => (a -> b -> m c) -> Matrix v a -> Matrix v b -> m ()
+zipWithM_ f m1 m2
+    | MG.dim m1 /= MG.dim m2 = error "zipWithM_: Dimensions don't match."
+    | otherwise = G.zipWithM_ f (MG.flatten m1) $ MG.flatten m2
+{-# INLINE zipWithM_ #-}
+
+unzip :: (G.Vector v a, G.Vector v b, G.Vector v (a,b))
+      => Matrix v (a,b) -> (Matrix v a, Matrix v b )
+unzip m = (MG.fromVector d v1, MG.fromVector d v2)
   where
-    (r,_) = dim mat
-{-# INLINE takeRow #-}
+    d = MG.dim m
+    (v1, v2) = G.unzip $ MG.flatten m
+{-# INLINE unzip #-}
 
--- | O(m) Return the rows
-toRows :: Matrix m v a => m v a -> [v a]
-toRows mat = map (unsafeTakeRow mat) [0..r-1]
+unzip3 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v (a,b,c))
+       => Matrix v (a,b, c) -> (Matrix v a, Matrix v b, Matrix v c)
+unzip3 m = (MG.fromVector d v1, MG.fromVector d v2, MG.fromVector d v3)
   where
-    (r,_) = dim mat
-{-# INLINE toRows #-}
+    d = MG.dim m
+    (v1, v2, v3) = G.unzip3 $ MG.flatten m
+{-# INLINE unzip3 #-}
 
--- | Extract a row.
-takeColumn :: Matrix m v a => m v a -> Int -> v a
-takeColumn mat j | j < 0 || j >= c =
-                   error $ printf "index out of bounds: (%d,%d)" j c
-                 | otherwise = unsafeTakeColumn mat j
+unzip4 :: (G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d, G.Vector v (a,b,c,d))
+       => Matrix v (a,b,c,d) -> (Matrix v a, Matrix v b, Matrix v c, Matrix v d)
+unzip4 m = ( MG.fromVector d v1
+           , MG.fromVector d v2
+           , MG.fromVector d v3
+           , MG.fromVector d v4
+           )
   where
-    (_,c) = dim mat
-{-# INLINE takeColumn #-}
+    d = MG.dim m
+    (v1, v2, v3, v4) = G.unzip4 $ MG.flatten m
+{-# INLINE unzip4 #-}
 
--- | O(m*n) Return the columns
-toColumns :: Matrix m v a => m v a -> [v a]
-toColumns mat = map (unsafeTakeColumn mat) [0..c-1]
+unzip5 :: ( G.Vector v a, G.Vector v b, G.Vector v c, G.Vector v d
+          , G.Vector v e, G.Vector v (a,b,c,d,e) )
+       => Matrix v (a,b,c,d,e)
+       -> (Matrix v a, Matrix v b, Matrix v c, Matrix v d, Matrix v e)
+unzip5 m = ( MG.fromVector d v1
+           , MG.fromVector d v2
+           , MG.fromVector d v3
+           , MG.fromVector d v4
+           , MG.fromVector d v5
+           )
   where
-    (_,c) = dim mat
-{-# INLINE toColumns #-}
+    d = MG.dim m
+    (v1, v2, v3, v4, v5) = G.unzip5 $ MG.flatten m
+{-# INLINE unzip5 #-}
 
--- | O(m*n) List of lists
-toLists :: Matrix m v a => m v a -> [[a]]
-toLists = map G.toList . toRows
-{-# INLINE toLists #-}
+unzip6 :: ( 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 (a,b,c,d,e,f) )
+       => Matrix v (a,b,c,d,e,f)
+       -> (Matrix v a, Matrix v b, Matrix v c, Matrix v d, Matrix v e, Matrix v f)
+unzip6 m = ( MG.fromVector d v1
+           , MG.fromVector d v2
+           , MG.fromVector d v3
+           , MG.fromVector d v4
+           , MG.fromVector d v5
+           , MG.fromVector d v6
+           )
+  where
+    d = MG.dim m
+    (v1, v2, v3, v4, v5, v6) = G.unzip6 $ MG.flatten m
+{-# INLINE unzip6 #-}
 
-create :: Matrix m v a => (forall s . ST s ((Mutable m) (G.Mutable v) s a)) -> m v a
-create m = runST $ unsafeFreeze =<< m
-{-# INLINE create #-}
+sequence :: (G.Vector v a, G.Vector v (m a), Monad m)
+         => Matrix v (m a) -> m (Matrix v a)
+sequence (Matrix r c tda offset vec) = liftM (Matrix r c tda offset) . G.sequence $ vec
+{-# INLINE sequence #-}
+
+sequence_ :: (G.Vector v (m a), Monad m)
+          => Matrix v (m a) -> m ()
+sequence_ (Matrix _ _ _ _ vec) = G.sequence_ vec
+{-# INLINE sequence_ #-}
+
+generate :: G.Vector v a => (Int, Int) -> ((Int, Int) -> a) -> Matrix v a
+generate (r,c) f = MG.fromVector (r,c) . G.generate (r*c) $ \i -> f (i `divMod` c)
+{-# INLINE generate #-}
diff --git a/src/Data/Matrix/Generic/Mutable.hs b/src/Data/Matrix/Generic/Mutable.hs
--- a/src/Data/Matrix/Generic/Mutable.hs
+++ b/src/Data/Matrix/Generic/Mutable.hs
@@ -1,46 +1,52 @@
+{-# LANGUAGE FlexibleInstances     #-}
+{-# LANGUAGE FlexibleContexts      #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE TypeFamilies          #-}
-
 module Data.Matrix.Generic.Mutable
-    ( MMatrix(..)
-    , write
-    , read
-    ) where
+   ( -- * Mutable Matrix
+     MMatrix(..)
+   , C.dim
+   , takeRow
+   , C.write
+   , C.unsafeWrite
+   , C.read
+   , C.unsafeRead
+   , C.new
+   , C.replicate
+   ) where
 
-import           Control.Monad.Primitive     (PrimMonad, PrimState)
+import           Control.Monad               (liftM)
+import           Control.DeepSeq
 import qualified Data.Vector.Generic.Mutable as GM
-import           Prelude                     hiding (read)
+import           Prelude                     hiding (read, replicate)
 
-class GM.MVector v a => MMatrix m v a where
-    dim ::  m v s a -> (Int, Int)
+import qualified Data.Matrix.Class.Mutable as C
 
-    unsafeRead :: PrimMonad s => m v (PrimState s) a -> (Int, Int) -> s a
+-- | mutable matrix
+data MMatrix v s a = MMatrix !Int !Int !Int !Int !(v s a)
 
-    unsafeWrite :: PrimMonad s => m v (PrimState s) a -> (Int, Int) -> a -> s ()
+instance (NFData (v s a)) => NFData (MMatrix v s a) where
+ rnf (MMatrix _ _ _ _ vec) = rnf vec
 
-    -- | Create a mutable matrix without initialization
-    new :: PrimMonad s => (Int, Int) -> s (m v (PrimState s) a)
+instance GM.MVector v a => C.MMatrix MMatrix v a where
+    dim (MMatrix r c _ _ _) = (r,c)
+    {-# INLINE dim #-}
 
-    replicate :: PrimMonad s => (Int, Int) -> a -> s (m v (PrimState s) a)
+    unsafeRead (MMatrix _ _ tda offset v) (i,j) = GM.unsafeRead v idx
+      where idx = offset + i * tda + j
+    {-# INLINE unsafeRead #-}
 
-    {-# MINIMAL dim, unsafeRead, unsafeWrite, new, replicate #-}
+    unsafeWrite (MMatrix _ _ tda offset v) (i,j) = GM.unsafeWrite v idx
+      where idx = offset + i * tda + j
+    {-# INLINE unsafeWrite #-}
 
--- | Derived methods
+    new (r,c) = MMatrix r c c 0 `liftM` GM.new (r*c)
+    {-# INLINE new #-}
 
-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 #-}
+    replicate (r,c) x = MMatrix r c c 0 `liftM` GM.replicate (r*c) x
+    {-# INLINE replicate #-}
 
-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)
+takeRow :: GM.MVector v a => MMatrix v m a -> Int -> v m a
+takeRow (MMatrix _ c tda offset vec) i = GM.slice i' c vec
   where
-    (r,c) = dim mat
-{-# INLINE read #-}
+    i' = offset + i * tda
+{-# INLINE takeRow #-}
diff --git a/src/Data/Matrix/Mutable.hs b/src/Data/Matrix/Mutable.hs
--- a/src/Data/Matrix/Mutable.hs
+++ b/src/Data/Matrix/Mutable.hs
@@ -1,10 +1,49 @@
+{-# LANGUAGE ConstraintKinds #-}
+{-# LANGUAGE KindSignatures #-}
 module Data.Matrix.Mutable
-    ( MMatrix
-    , module Data.Matrix.Dense.Generic.Mutable
+    ( -- * Mutable Matrix
+      MMatrix
+    , dim
+    , takeRow
+    , write
+    , unsafeWrite
+    , read
+    , unsafeRead
+    , new
+    , replicate
     ) where
 
-import           Data.Matrix.Dense.Generic.Mutable hiding (MMatrix)
-import qualified Data.Matrix.Dense.Generic.Mutable as MG
-import qualified Data.Vector.Mutable               as VM
+import GHC.Exts (Constraint)
+import Prelude hiding (read, replicate)
+import           Control.Monad.Primitive     (PrimMonad, PrimState)
+import Data.Vector.Mutable (MVector)
 
-type MMatrix a = MG.MMatrix VM.MVector a
+import qualified Data.Matrix.Generic.Mutable as MG
+
+type MMatrix a = MG.MMatrix MVector a
+type Context x = (() :: Constraint)
+
+dim :: Context a => MMatrix s a -> (Int, Int)
+dim = MG.dim
+
+takeRow :: Context a => MMatrix s a -> Int -> MVector s a
+takeRow = MG.takeRow
+
+write :: Context a => PrimMonad s => MMatrix (PrimState s) a -> (Int, Int) -> a -> s ()
+write = MG.write
+
+unsafeWrite :: (Context a, PrimMonad s) => MMatrix (PrimState s) a -> (Int, Int) -> a -> s ()
+unsafeWrite = MG.unsafeWrite
+
+read :: (Context a, PrimMonad s) => MMatrix (PrimState s) a -> (Int, Int) -> s a
+read = MG.read
+
+unsafeRead :: (Context a, PrimMonad s) => MMatrix (PrimState s) a -> (Int, Int) -> s a
+unsafeRead = MG.unsafeRead
+
+-- | Create a mutable matrix without initialization
+new :: (Context a, PrimMonad s) => (Int, Int) -> s (MMatrix (PrimState s) a)
+new = MG.new
+
+replicate :: (Context a, PrimMonad s) => (Int, Int) -> a -> s (MMatrix (PrimState s) a)
+replicate = MG.replicate
diff --git a/src/Data/Matrix/Sparse/Generic.hs b/src/Data/Matrix/Sparse/Generic.hs
--- a/src/Data/Matrix/Sparse/Generic.hs
+++ b/src/Data/Matrix/Sparse/Generic.hs
@@ -39,7 +39,6 @@
     , MG.toLists
     ) where
 
-import           Control.Applicative               ((<$>))
 import           Control.Monad                     (foldM, forM_, when)
 import           Control.Monad.ST                  (runST)
 import           Data.Bits                         (shiftR)
@@ -49,8 +48,8 @@
 import           Text.Printf                       (printf)
 import           GHC.Generics          (Generic)
 
-import           Data.Matrix.Dense.Generic.Mutable (MMatrix)
-import qualified Data.Matrix.Generic               as MG
+import           Data.Matrix.Generic.Mutable (MMatrix)
+import qualified Data.Matrix.Class as MG
 
 class Eq a => Zero a where
     zero :: a
diff --git a/src/Data/Matrix/Storable.hs b/src/Data/Matrix/Storable.hs
--- a/src/Data/Matrix/Storable.hs
+++ b/src/Data/Matrix/Storable.hs
@@ -1,10 +1,384 @@
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE ConstraintKinds #-}
+{-# LANGUAGE KindSignatures #-}
+
 module Data.Matrix.Storable
     ( Matrix
-    , module Data.Matrix.Dense.Generic
+
+    -- * Accessors
+    -- ** length information
+    , dim
+    , rows
+    , cols
+
+    -- ** Indexing
+    , unsafeIndex
+    , (!)
+    , takeRow
+    , takeColumn
+    , takeDiag
+
+    -- * Construction
+    , unsafeFromVector
+    , fromVector
+    , matrix
+    , fromList
+    , fromLists
+    , fromRows
+    , fromColumns
+    , empty
+
+    -- * Conversions
+    , flatten
+    , toRows
+    , toColumns
+    , toList
+    , toLists
+
+    , tr
+    , subMatrix
+    , ident
+    , diag
+    , diagRect
+    , fromBlocks
+    , isSymmetric
+    , force
+
+    , foldl
+
+    -- * Mapping
+    , map
+    , imap
+
+    -- * Monadic mapping
+    , mapM
+    , imapM
+    , mapM_
+    , imapM_
+    , forM
+    , forM_
+
+    -- * Zipping
+    , zipWith
+    , zipWith3
+    , zipWith4
+    , zipWith5
+    , zipWith6
+    , izipWith
+    , izipWith3
+    , izipWith4
+    , izipWith5
+    , izipWith6
+    , zip
+    , zip3
+    , zip4
+    , zip5
+    , zip6
+
+    -- * Monadic Zipping
+    , zipWithM
+    , zipWithM_
+
+    -- * Unzipping
+    , unzip
+    , unzip3
+    , unzip4
+    , unzip5
+    , unzip6
+
+    , generate
+
+    -- * Mutable matrix
+    , thaw
+    , unsafeThaw
+    , freeze
+    , unsafeFreeze
+    , create
     ) where
 
-import           Data.Matrix.Dense.Generic hiding (Matrix)
-import qualified Data.Matrix.Dense.Generic as MG
-import qualified Data.Vector.Storable      as V
+import GHC.Exts (Constraint)
+import Prelude hiding (sequence, sequence_, mapM_, zip, zip, zip3, zipWith, zipWith3, foldl, unzip, map, mapM, unzip3)
+import           Control.Monad.Primitive     (PrimMonad, PrimState)
+import           Control.Monad.ST            (ST)
+import Data.Foldable (Foldable)
+import Data.Vector.Storable (Vector, Storable)
 
-type Matrix a = MG.Matrix V.Vector a
+import qualified Data.Matrix.Generic as MG
+import Data.Matrix.Storable.Mutable (MMatrix)
+
+type Matrix = MG.Matrix Vector
+type Context x = (Storable x :: Constraint)
+
+dim :: Context a => Matrix a -> (Int, Int)
+dim = MG.dim
+
+rows :: Context a => Matrix a -> Int
+rows = MG.rows
+
+cols :: Context a => Matrix a -> Int
+cols = MG.cols
+
+unsafeIndex :: Context a => Matrix a -> (Int, Int) -> a
+unsafeIndex = MG.unsafeIndex
+
+(!) :: Context a => Matrix a -> (Int, Int) -> a
+(!) = (MG.!)
+
+takeRow :: Context a => Matrix a -> Int -> Vector a
+takeRow = MG.takeRow
+
+takeColumn :: Context a => Matrix a -> Int -> Vector a
+takeColumn = MG.takeColumn
+
+takeDiag :: Context a => Matrix a -> Vector a
+takeDiag = MG.takeDiag
+
+unsafeFromVector :: Context a => (Int, Int) -> Vector a -> Matrix a
+unsafeFromVector = MG.unsafeFromVector
+
+fromVector :: Context a => (Int, Int) -> Vector a -> Matrix a
+fromVector = MG.fromVector
+
+-- | O(m*n) Matrix construction
+matrix :: Context a => Int -> [a] -> Matrix a
+matrix = MG.matrix
+
+fromList :: Context a => (Int, Int) -> [a] -> Matrix a
+fromList = MG.fromList
+
+-- | O(m*n) Create matrix from list of lists, it doesn't check if the list of
+-- list is a valid matrix
+fromLists :: Context a => [[a]] -> Matrix a
+fromLists = MG.fromLists
+
+-- | O(m*n) Create matrix from rows
+fromRows :: Context a => [Vector a] -> Matrix a
+fromRows = MG.fromRows
+
+-- | O(m*n) Create matrix from columns
+fromColumns :: Context a => [Vector a] -> Matrix a
+fromColumns = MG.fromColumns
+
+empty :: Context a => Matrix a
+empty = MG.empty
+
+flatten :: Context a => Matrix a -> Vector a
+flatten = MG.flatten
+
+-- | O(m) Return the rows
+toRows :: Context a => Matrix a -> [Vector a]
+toRows = MG.toRows
+
+toColumns :: Context a => Matrix a -> [Vector a]
+toColumns = MG.toColumns
+
+-- | O(m*n) Create a list by concatenating rows
+toList :: Context a => Matrix a -> [a]
+toList = MG.toList
+
+-- | O(m*n) List of lists
+toLists :: Context a => Matrix a -> [[a]]
+toLists = MG.toLists
+
+-- | O(m*n) Matrix transpose
+tr :: Context a => Matrix a -> Matrix a
+tr = MG.tr
+
+-- | O(1) Extract sub matrix
+subMatrix :: Context a
+          => (Int, Int)  -- ^ upper left corner of the submatrix
+          -> (Int, Int)  -- ^ bottom right corner of the submatrix
+          -> Matrix a -> Matrix a
+subMatrix = MG.subMatrix
+
+-- | O(m*n) Create an identity matrix
+ident :: (Context a, Num a) => Int -> Matrix a
+ident = MG.ident
+
+-- | O(m*n) Create a square matrix with given diagonal, other entries default to 0
+diag :: (Context a, Num a, Foldable t)
+     => t a  -- ^ diagonal
+     -> Matrix a
+diag = MG.diag
+
+-- | O(m*n) Create a rectangular matrix with default values and given diagonal
+diagRect :: (Context a, Foldable t)
+         => a         -- ^ default value
+         -> (Int, Int)
+         -> t a       -- ^ diagonal
+         -> Matrix a
+diagRect = MG.diagRect
+
+fromBlocks :: Context a
+           => a    -- ^ default value
+           -> [[Matrix a]]
+           -> Matrix a
+fromBlocks = MG.fromBlocks
+
+isSymmetric :: (Context a, Eq a) => Matrix a -> Bool
+isSymmetric = MG.isSymmetric
+
+force :: Context a => Matrix a -> Matrix a
+force = MG.force
+
+foldl :: Context b => (a -> b -> a) -> a -> Matrix b -> a
+foldl = MG.foldl
+
+map :: (Context a, Context b) => (a -> b) -> Matrix a -> Matrix b
+map = MG.map
+
+imap :: (Context a, Context b) => ((Int, Int) -> a -> b) -> Matrix a -> Matrix b
+imap = MG.imap
+
+mapM :: (Context a, Context b, Monad m) => (a -> m b) -> Matrix a -> m (Matrix b)
+mapM = MG.mapM
+
+-- | O(m*n) Apply the monadic action to every element and its index,
+-- yielding a matrix of results.
+imapM :: (Context a, Context b, Monad m) => ((Int, Int) -> a -> m b) -> Matrix a -> m (Matrix b)
+imapM = MG.imapM
+
+mapM_ :: (Context a, Monad m) => (a -> m b) -> Matrix a -> m ()
+mapM_ = MG.mapM_
+
+-- | O(m*n) Apply the monadic action to every element and its index,
+-- ignoring the results.
+imapM_ :: (Context a, Monad m) => ((Int, Int) -> a -> m b) -> Matrix a -> m ()
+imapM_ = MG.imapM_
+
+forM :: (Context a, Context b, Monad m) => Matrix a -> (a -> m b) -> m (Matrix b)
+forM = MG.forM
+
+forM_ :: (Context a, Monad m) => Matrix a -> (a -> m b) -> m ()
+forM_ = MG.forM_
+
+zipWith :: ( Context a, Context b, Context c)
+        => (a -> b -> c) -> Matrix a -> Matrix b -> Matrix c
+zipWith = MG.zipWith
+
+zipWith3 :: ( Context a, Context b, Context c, Context d)
+         => (a -> b -> c -> d) -> Matrix a -> Matrix b -> Matrix c
+         -> Matrix d
+zipWith3 = MG.zipWith3
+
+zipWith4 :: ( Context a, Context b, Context c, Context d, Context e)
+         => (a -> b -> c -> d -> e) -> Matrix a -> Matrix b -> Matrix c
+         -> Matrix d -> Matrix e
+zipWith4 = MG.zipWith4
+
+zipWith5 :: ( Context a, Context b, Context c, Context d, Context e, Context f)
+         => (a -> b -> c -> d -> e -> f) -> Matrix a -> Matrix b
+         -> Matrix c -> Matrix d -> Matrix e -> Matrix f
+zipWith5 = MG.zipWith5
+
+zipWith6 :: ( Context a, Context b, Context c, Context d, Context e, Context f
+            , Context g )
+         => (a -> b -> c -> d -> e -> f -> g) -> Matrix a -> Matrix b
+         -> Matrix c -> Matrix d -> Matrix e -> Matrix f -> Matrix g
+zipWith6 = MG.zipWith6
+
+izipWith :: ( Context a, Context b, Context c)
+         => ((Int, Int) -> a -> b -> c) -> Matrix a -> Matrix b -> Matrix c
+izipWith = MG.izipWith
+
+izipWith3 :: ( Context a, Context b, Context c, Context d)
+          => ((Int, Int) -> a -> b -> c -> d) -> Matrix a -> Matrix b
+          -> Matrix c -> Matrix d
+izipWith3 = MG.izipWith3
+
+izipWith4 :: ( Context a, Context b, Context c, Context d, Context e)
+          => ((Int, Int) -> a -> b -> c -> d -> e) -> Matrix a -> Matrix b
+          -> Matrix c -> Matrix d -> Matrix e
+izipWith4 = MG.izipWith4
+
+izipWith5 :: ( Context a, Context b, Context c, Context d, Context e, Context f)
+          => ((Int, Int) -> a -> b -> c -> d -> e -> f) -> Matrix a
+          -> Matrix b -> Matrix c -> Matrix d -> Matrix e -> Matrix f
+izipWith5 = MG.izipWith5
+
+izipWith6 :: ( Context a, Context b, Context c, Context d, Context e, Context f
+             , Context g )
+          => ((Int, Int) -> a -> b -> c -> d -> e -> f -> g) -> Matrix a
+          -> Matrix b -> Matrix c -> Matrix d -> Matrix e -> Matrix f
+          -> Matrix g
+izipWith6 = MG.izipWith6
+
+zip :: ( Context a, Context b
+       , Context (a,b) )
+    => Matrix a -> Matrix b -> Matrix (a,b)
+zip = MG.zip
+
+zip3 :: ( Context a, Context b, Context c
+        , Context (a,b,c) )
+     => Matrix a -> Matrix b -> Matrix c -> Matrix (a,b,c)
+zip3 = MG.zip3
+
+zip4 :: ( Context a, Context b, Context c, Context d
+        , Context (a,b,c,d) )
+     => Matrix a -> Matrix b -> Matrix c -> Matrix d -> Matrix (a,b,c,d)
+zip4 = MG.zip4
+
+zip5 :: ( Context a, Context b, Context c, Context d, Context e
+        , Context (a,b,c,d,e) )
+     => Matrix a -> Matrix b -> Matrix c -> Matrix d -> Matrix e
+     -> Matrix (a,b,c,d,e)
+zip5 = MG.zip5
+
+zip6 :: ( Context a, Context b, Context c, Context d, Context e, Context f
+        , Context (a,b,c,d,e,f) )
+     => Matrix a -> Matrix b -> Matrix c -> Matrix d -> Matrix e
+     -> Matrix f -> Matrix (a,b,c,d,e,f)
+zip6 = MG.zip6
+
+zipWithM :: (Context a, Context b, Context c, Monad m)
+         => (a -> b -> m c) -> Matrix a -> Matrix b -> m (Matrix c)
+zipWithM = MG.zipWithM
+
+zipWithM_ :: (Context a, Context b, Monad m)
+          => (a -> b -> m c) -> Matrix a -> Matrix b -> m ()
+zipWithM_ = MG.zipWithM_
+
+unzip :: (Context a, Context b, Context (a,b))
+      => Matrix (a,b) -> (Matrix a, Matrix b )
+unzip = MG.unzip
+
+unzip3 :: ( Context a, Context b, Context c
+          , Context (a,b,c) )
+       => Matrix (a,b,c) -> (Matrix a, Matrix b, Matrix c)
+unzip3 = MG.unzip3
+
+unzip4 :: ( Context a, Context b, Context c, Context d
+          , Context (a,b,c,d) )
+       => Matrix (a,b,c,d) -> (Matrix a, Matrix b, Matrix c, Matrix d)
+unzip4 = MG.unzip4
+
+unzip5 :: ( Context a, Context b, Context c, Context d, Context e
+          , Context (a,b,c,d,e) )
+       => Matrix (a,b,c,d,e)
+       -> (Matrix a, Matrix b, Matrix c, Matrix d, Matrix e)
+unzip5 = MG.unzip5
+
+unzip6 :: ( Context a, Context b, Context c, Context d, Context e, Context f
+          , Context (a,b,c,d,e,f) )
+       => Matrix (a,b,c,d,e,f)
+       -> (Matrix a, Matrix b, Matrix c, Matrix d, Matrix e, Matrix f)
+unzip6 = MG.unzip6
+
+generate :: Context a => (Int, Int) -> ((Int, Int) -> a) -> Matrix a
+generate = MG.generate
+
+thaw :: (Context a, PrimMonad s) => Matrix a -> s (MMatrix (PrimState s) a)
+thaw = MG.thaw
+
+unsafeThaw :: (Context a, PrimMonad s) => Matrix a -> s (MMatrix (PrimState s) a)
+unsafeThaw = MG.unsafeThaw
+
+freeze :: (Context a, PrimMonad s) => MMatrix (PrimState s) a -> s (Matrix a)
+freeze = MG.freeze
+
+unsafeFreeze :: (Context a, PrimMonad s) => MMatrix (PrimState s) a -> s (Matrix a)
+unsafeFreeze = MG.unsafeFreeze
+
+create :: Context a => (forall s . ST s (MMatrix s a)) -> Matrix a
+create = MG.create
diff --git a/src/Data/Matrix/Storable/Mutable.hs b/src/Data/Matrix/Storable/Mutable.hs
--- a/src/Data/Matrix/Storable/Mutable.hs
+++ b/src/Data/Matrix/Storable/Mutable.hs
@@ -1,10 +1,49 @@
+{-# LANGUAGE ConstraintKinds #-}
+{-# LANGUAGE KindSignatures #-}
 module Data.Matrix.Storable.Mutable
-    ( MMatrix
-    , module Data.Matrix.Dense.Generic.Mutable
+    ( -- * Mutable Matrix
+      MMatrix
+    , dim
+    , takeRow
+    , write
+    , unsafeWrite
+    , read
+    , unsafeRead
+    , new
+    , replicate
     ) where
 
-import           Data.Matrix.Dense.Generic.Mutable hiding (MMatrix)
-import qualified Data.Matrix.Dense.Generic.Mutable as MG
-import qualified Data.Vector.Storable.Mutable      as VM
+import GHC.Exts (Constraint)
+import Prelude hiding (read, replicate)
+import           Control.Monad.Primitive     (PrimMonad, PrimState)
+import Data.Vector.Storable.Mutable (MVector, Storable)
 
-type MMatrix a = MG.MMatrix VM.MVector a
+import qualified Data.Matrix.Generic.Mutable as MG
+
+type MMatrix a = MG.MMatrix MVector a
+type Context x = (Storable x :: Constraint)
+
+dim :: Context a => MMatrix s a -> (Int, Int)
+dim = MG.dim
+
+takeRow :: Context a => MMatrix s a -> Int -> MVector s a
+takeRow = MG.takeRow
+
+write :: Context a => PrimMonad s => MMatrix (PrimState s) a -> (Int, Int) -> a -> s ()
+write = MG.write
+
+unsafeWrite :: (Context a, PrimMonad s) => MMatrix (PrimState s) a -> (Int, Int) -> a -> s ()
+unsafeWrite = MG.unsafeWrite
+
+read :: (Context a, PrimMonad s) => MMatrix (PrimState s) a -> (Int, Int) -> s a
+read = MG.read
+
+unsafeRead :: (Context a, PrimMonad s) => MMatrix (PrimState s) a -> (Int, Int) -> s a
+unsafeRead = MG.unsafeRead
+
+-- | Create a mutable matrix without initialization
+new :: (Context a, PrimMonad s) => (Int, Int) -> s (MMatrix (PrimState s) a)
+new = MG.new
+
+replicate :: (Context a, PrimMonad s) => (Int, Int) -> a -> s (MMatrix (PrimState s) a)
+replicate = MG.replicate
diff --git a/src/Data/Matrix/Symmetric.hs b/src/Data/Matrix/Symmetric.hs
deleted file mode 100644
--- a/src/Data/Matrix/Symmetric.hs
+++ /dev/null
@@ -1,120 +0,0 @@
-{-# LANGUAGE BangPatterns          #-}
-{-# LANGUAGE FlexibleContexts      #-}
-{-# LANGUAGE FlexibleInstances     #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-{-# LANGUAGE TypeFamilies          #-}
-{-# LANGUAGE DeriveGeneric      #-}
-module Data.Matrix.Symmetric
-    ( SymMatrix(..)
-    , dim
-    , rows
-    , cols
-    , unsafeIndex
-    , (!)
-    , flatten
-    , unsafeFromVector
-    , fromVector
-    , takeRow
-    , thaw
-    , unsafeThaw
-    , freeze
-    , unsafeFreeze
-    , create
-    , Data.Matrix.Symmetric.map
-    , imap
-    , zip
-    , zipWith
-    ) where
-
-import           Control.Monad                 (liftM)
-import           Data.Bits                     (shiftR)
-import qualified Data.Vector.Generic           as G
-import           Prelude                       hiding (zip, zipWith)
-import           GHC.Generics          (Generic)
-
-import           Data.Matrix.Generic
-import           Data.Matrix.Symmetric.Mutable (SymMMatrix (..), new,
-                                                unsafeWrite)
-
-type instance Mutable SymMatrix = SymMMatrix
-
--- | Symmetric square matrix
-data SymMatrix v a = SymMatrix !Int !(v a)
-    deriving (Show, Read, Generic, Eq)
-
-
---------------------------------------------------------------------------------
--- Instances
---------------------------------------------------------------------------------
-
-instance G.Vector v a => Matrix SymMatrix v a where
-    dim (SymMatrix n _) = (n,n)
-    {-# INLINE dim #-}
-
-    unsafeIndex (SymMatrix n vec) (i,j) = vec `G.unsafeIndex` idx n i j
-    {-# INLINE unsafeIndex #-}
-
-    unsafeFromVector (r,c) vec | r /= c = error "columns /= rows"
-                               | otherwise = SymMatrix r . G.concat . Prelude.map f $ [0..r-1]
-      where
-        f i = G.slice (i*(c+1)) (c-i) vec
---        n = ((r+1)*r) `shiftR` 1
-    {-# INLINE unsafeFromVector #-}
-
-    thaw (SymMatrix n v) = SymMMatrix n `liftM` G.thaw v
-    {-# INLINE thaw #-}
-
-    unsafeThaw (SymMatrix n v) = SymMMatrix n `liftM` G.thaw v
-    {-# INLINE unsafeThaw #-}
-
-    freeze (SymMMatrix n v) = SymMatrix n `liftM` G.freeze v
-    {-# INLINE freeze #-}
-
-    unsafeFreeze (SymMMatrix n v) = SymMatrix n `liftM` G.unsafeFreeze v
-    {-# INLINE unsafeFreeze #-}
-
---------------------------------------------------------------------------------
-
-map :: (G.Vector v a, G.Vector v b) => (a -> b) -> SymMatrix v a -> SymMatrix v b
-map f (SymMatrix n vec) = SymMatrix n $ G.map f vec
-{-# INLINE map #-}
-
--- | Upper triangular imap, i.e., i <= j
-imap :: (G.Vector v a, G.Vector v b) => ((Int, Int) -> a -> b) -> SymMatrix v a -> SymMatrix v b
-imap f mat = create $ do
-    mat' <- new (n,n)
-    let loop m !i !j
-            | i >= n = return ()
-            | j >= n = loop m (i+1) (i+1)
-            | otherwise = unsafeWrite m (i,j) (f (i,j) x) >>
-                          loop m i (j+1)
-          where
-            x = unsafeIndex mat (i,j)
-    loop mat' 0 0
-    return mat'
-  where
-    n = rows mat
-{-# INLINE imap #-}
-
-zip :: (G.Vector v a, G.Vector v b, G.Vector v (a,b))
-    => SymMatrix v a -> SymMatrix v b -> SymMatrix v (a,b)
-zip (SymMatrix n1 v1) (SymMatrix n2 v2)
-    | n1 /= n2 = error "imcompatible size"
-    | otherwise = SymMatrix n1 $ G.zip v1 v2
-{-# INLINE zip #-}
-
-zipWith :: (G.Vector v a, G.Vector v b, G.Vector v c)
-        => (a -> b -> c) -> SymMatrix v a -> SymMatrix v b -> SymMatrix v c
-zipWith f (SymMatrix n1 v1) (SymMatrix n2 v2)
-    | n1 /= n2 = error "imcompatible size"
-    | otherwise = SymMatrix n1 . G.zipWith f v1 $ v2
-{-# INLINE zipWith #-}
-
-
--- helper
-
--- row major upper triangular indexing
-idx :: Int -> Int -> Int -> Int
-idx n i j | i <= j = (i * (2 * n - i - 1)) `shiftR` 1 + j
-          | otherwise = (j * (2 * n - j - 1)) `shiftR` 1 + i
-{-# INLINE idx #-}
diff --git a/src/Data/Matrix/Symmetric/Generic.hs b/src/Data/Matrix/Symmetric/Generic.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Symmetric/Generic.hs
@@ -0,0 +1,120 @@
+{-# LANGUAGE BangPatterns          #-}
+{-# LANGUAGE FlexibleContexts      #-}
+{-# LANGUAGE FlexibleInstances     #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE TypeFamilies          #-}
+{-# LANGUAGE DeriveGeneric      #-}
+module Data.Matrix.Symmetric.Generic
+    ( SymMatrix(..)
+    , dim
+    , rows
+    , cols
+    , unsafeIndex
+    , (!)
+    , flatten
+    , unsafeFromVector
+    , fromVector
+    , takeRow
+    , thaw
+    , unsafeThaw
+    , freeze
+    , unsafeFreeze
+    , create
+    , Data.Matrix.Symmetric.Generic.map
+    , imap
+    , zip
+    , zipWith
+    ) where
+
+import           Control.Monad                 (liftM)
+import           Data.Bits                     (shiftR)
+import qualified Data.Vector.Generic           as G
+import           Prelude                       hiding (zip, zipWith)
+import           GHC.Generics          (Generic)
+
+import           Data.Matrix.Class
+import           Data.Matrix.Symmetric.Generic.Mutable (SymMMatrix (..), new,
+                                                unsafeWrite)
+
+type instance Mutable SymMatrix = SymMMatrix
+
+-- | Symmetric square matrix
+data SymMatrix v a = SymMatrix !Int !(v a)
+    deriving (Show, Read, Generic, Eq)
+
+
+--------------------------------------------------------------------------------
+-- Instances
+--------------------------------------------------------------------------------
+
+instance G.Vector v a => Matrix SymMatrix v a where
+    dim (SymMatrix n _) = (n,n)
+    {-# INLINE dim #-}
+
+    unsafeIndex (SymMatrix n vec) (i,j) = vec `G.unsafeIndex` idx n i j
+    {-# INLINE unsafeIndex #-}
+
+    unsafeFromVector (r,c) vec | r /= c = error "columns /= rows"
+                               | otherwise = SymMatrix r . G.concat . Prelude.map f $ [0..r-1]
+      where
+        f i = G.slice (i*(c+1)) (c-i) vec
+--        n = ((r+1)*r) `shiftR` 1
+    {-# INLINE unsafeFromVector #-}
+
+    thaw (SymMatrix n v) = SymMMatrix n `liftM` G.thaw v
+    {-# INLINE thaw #-}
+
+    unsafeThaw (SymMatrix n v) = SymMMatrix n `liftM` G.thaw v
+    {-# INLINE unsafeThaw #-}
+
+    freeze (SymMMatrix n v) = SymMatrix n `liftM` G.freeze v
+    {-# INLINE freeze #-}
+
+    unsafeFreeze (SymMMatrix n v) = SymMatrix n `liftM` G.unsafeFreeze v
+    {-# INLINE unsafeFreeze #-}
+
+--------------------------------------------------------------------------------
+
+map :: (G.Vector v a, G.Vector v b) => (a -> b) -> SymMatrix v a -> SymMatrix v b
+map f (SymMatrix n vec) = SymMatrix n $ G.map f vec
+{-# INLINE map #-}
+
+-- | Upper triangular imap, i.e., i <= j
+imap :: (G.Vector v a, G.Vector v b) => ((Int, Int) -> a -> b) -> SymMatrix v a -> SymMatrix v b
+imap f mat = create $ do
+    mat' <- new (n,n)
+    let loop m !i !j
+            | i >= n = return ()
+            | j >= n = loop m (i+1) (i+1)
+            | otherwise = unsafeWrite m (i,j) (f (i,j) x) >>
+                          loop m i (j+1)
+          where
+            x = unsafeIndex mat (i,j)
+    loop mat' 0 0
+    return mat'
+  where
+    n = rows mat
+{-# INLINE imap #-}
+
+zip :: (G.Vector v a, G.Vector v b, G.Vector v (a,b))
+    => SymMatrix v a -> SymMatrix v b -> SymMatrix v (a,b)
+zip (SymMatrix n1 v1) (SymMatrix n2 v2)
+    | n1 /= n2 = error "imcompatible size"
+    | otherwise = SymMatrix n1 $ G.zip v1 v2
+{-# INLINE zip #-}
+
+zipWith :: (G.Vector v a, G.Vector v b, G.Vector v c)
+        => (a -> b -> c) -> SymMatrix v a -> SymMatrix v b -> SymMatrix v c
+zipWith f (SymMatrix n1 v1) (SymMatrix n2 v2)
+    | n1 /= n2 = error "imcompatible size"
+    | otherwise = SymMatrix n1 . G.zipWith f v1 $ v2
+{-# INLINE zipWith #-}
+
+
+-- helper
+
+-- row major upper triangular indexing
+idx :: Int -> Int -> Int -> Int
+idx n i j | i <= j = (i * (2 * n - i - 1)) `shiftR` 1 + j
+          | otherwise = (j * (2 * n - j - 1)) `shiftR` 1 + i
+{-# INLINE idx #-}
diff --git a/src/Data/Matrix/Symmetric/Generic/Mutable.hs b/src/Data/Matrix/Symmetric/Generic/Mutable.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Symmetric/Generic/Mutable.hs
@@ -0,0 +1,46 @@
+{-# LANGUAGE FlexibleInstances     #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+module Data.Matrix.Symmetric.Generic.Mutable
+   ( -- * Mutable Matrix
+     SymMMatrix(..)
+   , C.dim
+   , C.write
+   , C.unsafeWrite
+   , C.read
+   , C.unsafeRead
+   , C.new
+   , C.replicate
+   ) where
+
+import           Control.Monad               (liftM)
+import           Data.Bits                   (shiftR)
+import qualified Data.Vector.Generic.Mutable as GM
+import           Prelude                     hiding (read, replicate)
+
+import qualified Data.Matrix.Class.Mutable as C
+
+-- | mutable matrix
+data SymMMatrix v s a = SymMMatrix !Int !(v s a)
+
+instance GM.MVector v a => C.MMatrix SymMMatrix v a where
+    dim (SymMMatrix n _) = (n,n)
+    {-# INLINE dim #-}
+
+    unsafeRead (SymMMatrix n v) (i,j) = GM.unsafeRead v (idx n i j)
+    {-# INLINE unsafeRead #-}
+
+    unsafeWrite (SymMMatrix n v) (i,j) = GM.unsafeWrite v (idx n i j)
+    {-# INLINE unsafeWrite #-}
+
+    new (r,c) | r /= c = error "colmumns /= rows"
+              | otherwise = SymMMatrix r `liftM` GM.new ((r*(r+1)) `shiftR` 1)
+
+    replicate (r,c) x
+        | r /= c = error "colmumns /= rows"
+        | otherwise = SymMMatrix r `liftM` GM.replicate ((r*(r+1)) `shiftR` 1) x
+
+-- row major upper triangular indexing
+idx :: Int -> Int -> Int -> Int
+idx n i j | i <= j = (i * (2 * n - i - 1)) `shiftR` 1 + j
+          | otherwise = (j * (2 * n - j - 1)) `shiftR` 1 + i
+{-# INLINE idx #-}
diff --git a/src/Data/Matrix/Symmetric/Mutable.hs b/src/Data/Matrix/Symmetric/Mutable.hs
deleted file mode 100644
--- a/src/Data/Matrix/Symmetric/Mutable.hs
+++ /dev/null
@@ -1,46 +0,0 @@
-{-# LANGUAGE FlexibleInstances     #-}
-{-# LANGUAGE MultiParamTypeClasses #-}
-module Data.Matrix.Symmetric.Mutable
-   ( -- * Mutable Matrix
-     SymMMatrix(..)
-   , C.dim
-   , C.write
-   , C.unsafeWrite
-   , C.read
-   , C.unsafeRead
-   , C.new
-   , C.replicate
-   ) where
-
-import           Control.Monad               (liftM)
-import           Data.Bits                   (shiftR)
-import qualified Data.Vector.Generic.Mutable as GM
-import           Prelude                     hiding (read, replicate)
-
-import qualified Data.Matrix.Generic.Mutable as C
-
--- | mutable matrix
-data SymMMatrix v s a = SymMMatrix !Int !(v s a)
-
-instance GM.MVector v a => C.MMatrix SymMMatrix v a where
-    dim (SymMMatrix n _) = (n,n)
-    {-# INLINE dim #-}
-
-    unsafeRead (SymMMatrix n v) (i,j) = GM.unsafeRead v (idx n i j)
-    {-# INLINE unsafeRead #-}
-
-    unsafeWrite (SymMMatrix n v) (i,j) = GM.unsafeWrite v (idx n i j)
-    {-# INLINE unsafeWrite #-}
-
-    new (r,c) | r /= c = error "colmumns /= rows"
-              | otherwise = SymMMatrix r `liftM` GM.new ((r*(r+1)) `shiftR` 1)
-
-    replicate (r,c) x
-        | r /= c = error "colmumns /= rows"
-        | otherwise = SymMMatrix r `liftM` GM.replicate ((r*(r+1)) `shiftR` 1) x
-
--- row major upper triangular indexing
-idx :: Int -> Int -> Int -> Int
-idx n i j | i <= j = (i * (2 * n - i - 1)) `shiftR` 1 + j
-          | otherwise = (j * (2 * n - j - 1)) `shiftR` 1 + i
-{-# INLINE idx #-}
diff --git a/src/Data/Matrix/Unboxed.hs b/src/Data/Matrix/Unboxed.hs
--- a/src/Data/Matrix/Unboxed.hs
+++ b/src/Data/Matrix/Unboxed.hs
@@ -1,10 +1,384 @@
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE ConstraintKinds #-}
+{-# LANGUAGE KindSignatures #-}
+{-# LANGUAGE FlexibleContexts #-}
+
 module Data.Matrix.Unboxed
     ( Matrix
-    , module Data.Matrix.Dense.Generic
+
+    -- * Accessors
+    -- ** length information
+    , dim
+    , rows
+    , cols
+
+    -- ** Indexing
+    , unsafeIndex
+    , (!)
+    , takeRow
+    , takeColumn
+    , takeDiag
+
+    -- * Construction
+    , unsafeFromVector
+    , fromVector
+    , matrix
+    , fromList
+    , fromLists
+    , fromRows
+    , fromColumns
+    , empty
+
+    -- * Conversions
+    , flatten
+    , toRows
+    , toColumns
+    , toList
+    , toLists
+
+    , tr
+    , subMatrix
+    , ident
+    , diag
+    , diagRect
+    , fromBlocks
+    , isSymmetric
+    , force
+
+    , foldl
+
+    -- * Mapping
+    , map
+    , imap
+
+    -- * Monadic mapping
+    , mapM
+    , imapM
+    , mapM_
+    , imapM_
+    , forM
+    , forM_
+
+    -- * Zipping
+    , zipWith
+    , zipWith3
+    , zipWith4
+    , zipWith5
+    , zipWith6
+    , izipWith
+    , izipWith3
+    , izipWith4
+    , izipWith5
+    , izipWith6
+    , zip
+    , zip3
+    , zip4
+    , zip5
+    , zip6
+
+    -- * Monadic Zipping
+    , zipWithM
+    , zipWithM_
+
+    -- * Unzipping
+    , unzip
+    , unzip3
+    , unzip4
+    , unzip5
+    , unzip6
+
+    , generate
+
+    -- * Mutable matrix
+    , thaw
+    , unsafeThaw
+    , freeze
+    , unsafeFreeze
+    , create
     ) where
 
-import           Data.Matrix.Dense.Generic hiding (Matrix)
-import qualified Data.Matrix.Dense.Generic as MG
-import qualified Data.Vector.Unboxed       as V
+import GHC.Exts (Constraint)
+import Prelude hiding (sequence, sequence_, mapM_, zip, zip, zip3, zipWith, zipWith3, foldl, unzip, map, mapM, unzip3)
+import           Control.Monad.Primitive     (PrimMonad, PrimState)
+import           Control.Monad.ST            (ST)
+import Data.Foldable (Foldable)
+import Data.Vector.Unboxed (Vector, Unbox)
 
-type Matrix a = MG.Matrix V.Vector a
+import qualified Data.Matrix.Generic as MG
+import Data.Matrix.Unboxed.Mutable (MMatrix)
+
+type Matrix = MG.Matrix Vector
+type Context x = (Unbox x :: Constraint)
+
+dim :: Context a => Matrix a -> (Int, Int)
+dim = MG.dim
+
+rows :: Context a => Matrix a -> Int
+rows = MG.rows
+
+cols :: Context a => Matrix a -> Int
+cols = MG.cols
+
+unsafeIndex :: Context a => Matrix a -> (Int, Int) -> a
+unsafeIndex = MG.unsafeIndex
+
+(!) :: Context a => Matrix a -> (Int, Int) -> a
+(!) = (MG.!)
+
+takeRow :: Context a => Matrix a -> Int -> Vector a
+takeRow = MG.takeRow
+
+takeColumn :: Context a => Matrix a -> Int -> Vector a
+takeColumn = MG.takeColumn
+
+takeDiag :: Context a => Matrix a -> Vector a
+takeDiag = MG.takeDiag
+
+unsafeFromVector :: Context a => (Int, Int) -> Vector a -> Matrix a
+unsafeFromVector = MG.unsafeFromVector
+
+fromVector :: Context a => (Int, Int) -> Vector a -> Matrix a
+fromVector = MG.fromVector
+
+-- | O(m*n) Matrix construction
+matrix :: Context a => Int -> [a] -> Matrix a
+matrix = MG.matrix
+
+fromList :: Context a => (Int, Int) -> [a] -> Matrix a
+fromList = MG.fromList
+
+-- | O(m*n) Create matrix from list of lists, it doesn't check if the list of
+-- list is a valid matrix
+fromLists :: Context a => [[a]] -> Matrix a
+fromLists = MG.fromLists
+
+-- | O(m*n) Create matrix from rows
+fromRows :: Context a => [Vector a] -> Matrix a
+fromRows = MG.fromRows
+
+-- | O(m*n) Create matrix from columns
+fromColumns :: Context a => [Vector a] -> Matrix a
+fromColumns = MG.fromColumns
+
+empty :: Context a => Matrix a
+empty = MG.empty
+
+flatten :: Context a => Matrix a -> Vector a
+flatten = MG.flatten
+
+-- | O(m) Return the rows
+toRows :: Context a => Matrix a -> [Vector a]
+toRows = MG.toRows
+
+toColumns :: Context a => Matrix a -> [Vector a]
+toColumns = MG.toColumns
+
+-- | O(m*n) Create a list by concatenating rows
+toList :: Context a => Matrix a -> [a]
+toList = MG.toList
+
+-- | O(m*n) List of lists
+toLists :: Context a => Matrix a -> [[a]]
+toLists = MG.toLists
+
+-- | O(m*n) Matrix transpose
+tr :: Context a => Matrix a -> Matrix a
+tr = MG.tr
+
+-- | O(1) Extract sub matrix
+subMatrix :: Context a
+          => (Int, Int)  -- ^ upper left corner of the submatrix
+          -> (Int, Int)  -- ^ bottom right corner of the submatrix
+          -> Matrix a -> Matrix a
+subMatrix = MG.subMatrix
+
+-- | O(m*n) Create an identity matrix
+ident :: (Context a, Num a) => Int -> Matrix a
+ident = MG.ident
+
+-- | O(m*n) Create a square matrix with given diagonal, other entries default to 0
+diag :: (Context a, Num a, Foldable t)
+     => t a  -- ^ diagonal
+     -> Matrix a
+diag = MG.diag
+
+-- | O(m*n) Create a rectangular matrix with default values and given diagonal
+diagRect :: (Context a, Foldable t)
+         => a         -- ^ default value
+         -> (Int, Int)
+         -> t a       -- ^ diagonal
+         -> Matrix a
+diagRect = MG.diagRect
+
+fromBlocks :: Context a
+           => a    -- ^ default value
+           -> [[Matrix a]]
+           -> Matrix a
+fromBlocks = MG.fromBlocks
+
+isSymmetric :: (Context a, Eq a) => Matrix a -> Bool
+isSymmetric = MG.isSymmetric
+
+force :: Context a => Matrix a -> Matrix a
+force = MG.force
+
+foldl :: Context b => (a -> b -> a) -> a -> Matrix b -> a
+foldl = MG.foldl
+
+map :: (Context a, Context b) => (a -> b) -> Matrix a -> Matrix b
+map = MG.map
+
+imap :: (Context a, Context b) => ((Int, Int) -> a -> b) -> Matrix a -> Matrix b
+imap = MG.imap
+
+mapM :: (Context a, Context b, Monad m) => (a -> m b) -> Matrix a -> m (Matrix b)
+mapM = MG.mapM
+
+-- | O(m*n) Apply the monadic action to every element and its index,
+-- yielding a matrix of results.
+imapM :: (Context a, Context b, Monad m) => ((Int, Int) -> a -> m b) -> Matrix a -> m (Matrix b)
+imapM = MG.imapM
+
+mapM_ :: (Context a, Monad m) => (a -> m b) -> Matrix a -> m ()
+mapM_ = MG.mapM_
+
+-- | O(m*n) Apply the monadic action to every element and its index,
+-- ignoring the results.
+imapM_ :: (Context a, Monad m) => ((Int, Int) -> a -> m b) -> Matrix a -> m ()
+imapM_ = MG.imapM_
+
+forM :: (Context a, Context b, Monad m) => Matrix a -> (a -> m b) -> m (Matrix b)
+forM = MG.forM
+
+forM_ :: (Context a, Monad m) => Matrix a -> (a -> m b) -> m ()
+forM_ = MG.forM_
+
+zipWith :: ( Context a, Context b, Context c)
+        => (a -> b -> c) -> Matrix a -> Matrix b -> Matrix c
+zipWith = MG.zipWith
+
+zipWith3 :: ( Context a, Context b, Context c, Context d)
+         => (a -> b -> c -> d) -> Matrix a -> Matrix b -> Matrix c
+         -> Matrix d
+zipWith3 = MG.zipWith3
+
+zipWith4 :: ( Context a, Context b, Context c, Context d, Context e)
+         => (a -> b -> c -> d -> e) -> Matrix a -> Matrix b -> Matrix c
+         -> Matrix d -> Matrix e
+zipWith4 = MG.zipWith4
+
+zipWith5 :: ( Context a, Context b, Context c, Context d, Context e, Context f)
+         => (a -> b -> c -> d -> e -> f) -> Matrix a -> Matrix b
+         -> Matrix c -> Matrix d -> Matrix e -> Matrix f
+zipWith5 = MG.zipWith5
+
+zipWith6 :: ( Context a, Context b, Context c, Context d, Context e, Context f
+            , Context g )
+         => (a -> b -> c -> d -> e -> f -> g) -> Matrix a -> Matrix b
+         -> Matrix c -> Matrix d -> Matrix e -> Matrix f -> Matrix g
+zipWith6 = MG.zipWith6
+
+izipWith :: ( Context a, Context b, Context c)
+         => ((Int, Int) -> a -> b -> c) -> Matrix a -> Matrix b -> Matrix c
+izipWith = MG.izipWith
+
+izipWith3 :: ( Context a, Context b, Context c, Context d)
+          => ((Int, Int) -> a -> b -> c -> d) -> Matrix a -> Matrix b
+          -> Matrix c -> Matrix d
+izipWith3 = MG.izipWith3
+
+izipWith4 :: ( Context a, Context b, Context c, Context d, Context e)
+          => ((Int, Int) -> a -> b -> c -> d -> e) -> Matrix a -> Matrix b
+          -> Matrix c -> Matrix d -> Matrix e
+izipWith4 = MG.izipWith4
+
+izipWith5 :: ( Context a, Context b, Context c, Context d, Context e, Context f)
+          => ((Int, Int) -> a -> b -> c -> d -> e -> f) -> Matrix a
+          -> Matrix b -> Matrix c -> Matrix d -> Matrix e -> Matrix f
+izipWith5 = MG.izipWith5
+
+izipWith6 :: ( Context a, Context b, Context c, Context d, Context e, Context f
+             , Context g )
+          => ((Int, Int) -> a -> b -> c -> d -> e -> f -> g) -> Matrix a
+          -> Matrix b -> Matrix c -> Matrix d -> Matrix e -> Matrix f
+          -> Matrix g
+izipWith6 = MG.izipWith6
+
+zip :: ( Context a, Context b
+       , Context (a,b) )
+    => Matrix a -> Matrix b -> Matrix (a,b)
+zip = MG.zip
+
+zip3 :: ( Context a, Context b, Context c
+        , Context (a,b,c) )
+     => Matrix a -> Matrix b -> Matrix c -> Matrix (a,b,c)
+zip3 = MG.zip3
+
+zip4 :: ( Context a, Context b, Context c, Context d
+        , Context (a,b,c,d) )
+     => Matrix a -> Matrix b -> Matrix c -> Matrix d -> Matrix (a,b,c,d)
+zip4 = MG.zip4
+
+zip5 :: ( Context a, Context b, Context c, Context d, Context e
+        , Context (a,b,c,d,e) )
+     => Matrix a -> Matrix b -> Matrix c -> Matrix d -> Matrix e
+     -> Matrix (a,b,c,d,e)
+zip5 = MG.zip5
+
+zip6 :: ( Context a, Context b, Context c, Context d, Context e, Context f
+        , Context (a,b,c,d,e,f) )
+     => Matrix a -> Matrix b -> Matrix c -> Matrix d -> Matrix e
+     -> Matrix f -> Matrix (a,b,c,d,e,f)
+zip6 = MG.zip6
+
+zipWithM :: (Context a, Context b, Context c, Monad m)
+         => (a -> b -> m c) -> Matrix a -> Matrix b -> m (Matrix c)
+zipWithM = MG.zipWithM
+
+zipWithM_ :: (Context a, Context b, Monad m)
+          => (a -> b -> m c) -> Matrix a -> Matrix b -> m ()
+zipWithM_ = MG.zipWithM_
+
+unzip :: (Context a, Context b, Context (a,b))
+      => Matrix (a,b) -> (Matrix a, Matrix b )
+unzip = MG.unzip
+
+unzip3 :: ( Context a, Context b, Context c
+          , Context (a,b,c) )
+       => Matrix (a,b,c) -> (Matrix a, Matrix b, Matrix c)
+unzip3 = MG.unzip3
+
+unzip4 :: ( Context a, Context b, Context c, Context d
+          , Context (a,b,c,d) )
+       => Matrix (a,b,c,d) -> (Matrix a, Matrix b, Matrix c, Matrix d)
+unzip4 = MG.unzip4
+
+unzip5 :: ( Context a, Context b, Context c, Context d, Context e
+          , Context (a,b,c,d,e) )
+       => Matrix (a,b,c,d,e)
+       -> (Matrix a, Matrix b, Matrix c, Matrix d, Matrix e)
+unzip5 = MG.unzip5
+
+unzip6 :: ( Context a, Context b, Context c, Context d, Context e, Context f
+          , Context (a,b,c,d,e,f) )
+       => Matrix (a,b,c,d,e,f)
+       -> (Matrix a, Matrix b, Matrix c, Matrix d, Matrix e, Matrix f)
+unzip6 = MG.unzip6
+
+generate :: Context a => (Int, Int) -> ((Int, Int) -> a) -> Matrix a
+generate = MG.generate
+
+thaw :: (Context a, PrimMonad s) => Matrix a -> s (MMatrix (PrimState s) a)
+thaw = MG.thaw
+
+unsafeThaw :: (Context a, PrimMonad s) => Matrix a -> s (MMatrix (PrimState s) a)
+unsafeThaw = MG.unsafeThaw
+
+freeze :: (Context a, PrimMonad s) => MMatrix (PrimState s) a -> s (Matrix a)
+freeze = MG.freeze
+
+unsafeFreeze :: (Context a, PrimMonad s) => MMatrix (PrimState s) a -> s (Matrix a)
+unsafeFreeze = MG.unsafeFreeze
+
+create :: Context a => (forall s . ST s (MMatrix s a)) -> Matrix a
+create = MG.create
diff --git a/src/Data/Matrix/Unboxed/Mutable.hs b/src/Data/Matrix/Unboxed/Mutable.hs
--- a/src/Data/Matrix/Unboxed/Mutable.hs
+++ b/src/Data/Matrix/Unboxed/Mutable.hs
@@ -1,10 +1,49 @@
+{-# LANGUAGE ConstraintKinds #-}
+{-# LANGUAGE KindSignatures #-}
 module Data.Matrix.Unboxed.Mutable
-    ( MMatrix
-    , module Data.Matrix.Dense.Generic.Mutable
+    ( -- * Mutable Matrix
+      MMatrix
+    , dim
+    , takeRow
+    , write
+    , unsafeWrite
+    , read
+    , unsafeRead
+    , new
+    , replicate
     ) where
 
-import           Data.Matrix.Dense.Generic.Mutable hiding (MMatrix)
-import qualified Data.Matrix.Dense.Generic.Mutable as MG
-import qualified Data.Vector.Unboxed.Mutable       as VM
+import GHC.Exts (Constraint)
+import Prelude hiding (read, replicate)
+import           Control.Monad.Primitive     (PrimMonad, PrimState)
+import Data.Vector.Unboxed.Mutable (MVector, Unbox)
 
-type MMatrix a = MG.MMatrix VM.MVector a
+import qualified Data.Matrix.Generic.Mutable as MG
+
+type MMatrix a = MG.MMatrix MVector a
+type Context x = (Unbox x :: Constraint)
+
+dim :: Context a => MMatrix s a -> (Int, Int)
+dim = MG.dim
+
+takeRow :: Context a => MMatrix s a -> Int -> MVector s a
+takeRow = MG.takeRow
+
+write :: Context a => PrimMonad s => MMatrix (PrimState s) a -> (Int, Int) -> a -> s ()
+write = MG.write
+
+unsafeWrite :: (Context a, PrimMonad s) => MMatrix (PrimState s) a -> (Int, Int) -> a -> s ()
+unsafeWrite = MG.unsafeWrite
+
+read :: (Context a, PrimMonad s) => MMatrix (PrimState s) a -> (Int, Int) -> s a
+read = MG.read
+
+unsafeRead :: (Context a, PrimMonad s) => MMatrix (PrimState s) a -> (Int, Int) -> s a
+unsafeRead = MG.unsafeRead
+
+-- | Create a mutable matrix without initialization
+new :: (Context a, PrimMonad s) => (Int, Int) -> s (MMatrix (PrimState s) a)
+new = MG.new
+
+replicate :: (Context a, PrimMonad s) => (Int, Int) -> a -> s (MMatrix (PrimState s) a)
+replicate = MG.replicate
diff --git a/tests/test.hs b/tests/test.hs
--- a/tests/test.hs
+++ b/tests/test.hs
@@ -1,27 +1,54 @@
 import Test.Tasty
-import qualified Data.Matrix.Generic as MG
-import qualified Data.Matrix.Dense.Generic as MD
+import qualified Data.Matrix.Unboxed as MU
+import qualified Data.Matrix.Class as C
 import qualified Data.Matrix.Sparse.Generic as MS
 import qualified Data.Vector.Unboxed as U
 import Test.Tasty.HUnit
 
 main :: IO ()
 main = defaultMain $ testGroup "Main"
-    [ testCase "xx" testEqual ]
+    [ testCase "xx" testEqual
+    , subMatrixTest
+    , subMatrixEqual
+    ]
 
+subMatrixEqual :: TestTree
+subMatrixEqual = testCase "submatrix equal" $ assertEqual "x" mat submat
+  where
+    mat = MU.fromLists [[1]] :: MU.Matrix Int
+    submat = MU.subMatrix (0,0) (0,0) (MU.fromLists [[1,2],[3,4]]) :: MU.Matrix Int
 
 testEqual :: Assertion
 testEqual = do
     let xs = [0,0,0,0,1,2,3,0,0,0,0,0,4,5,67,0,0,2,40,0,2,0,0,20,0,0,0]
-        m1 = MG.fromList (3,9) xs :: MD.Matrix U.Vector Int
-        al = filter ((/=0) . snd) $ MD.toList $ MD.imap (\i v -> (i,v)) m1
-        m2 = MG.fromList (3,9) xs :: MS.CSR U.Vector Int
+        m1 = MU.fromList (3,9) xs
+        al = filter ((/=0) . snd) $ MU.toList $ MU.imap (\i v -> (i,v)) m1
+        m2 = C.fromList (3,9) xs :: MS.CSR U.Vector Int
         m3 = MS.fromAscAL (3,9) (length al) al :: MS.CSR U.Vector Int
 
-        row1 = MG.toRows m1
-        row2 = MG.toRows m2
-    
+        row1 = C.toRows m1
+        row2 = C.toRows m2
+
 --    assertEqual "x" (MG.flatten m1) (MG.flatten m2)
-    assertEqual "x" (MG.flatten m2) (MG.flatten m3)
+    assertEqual "x" (C.flatten m2) (C.flatten m3)
     assertEqual "x" row1 row2
 
+subMatrixTest :: TestTree
+subMatrixTest = testGroup "subMatrix"
+    [ testCase "case 1" $ [[5,6], [8,9]] @=? MU.toLists sub1
+    , testCase "case 2" $ [[5,6]] @=? MU.toLists sub2
+    , testCase "case 3" $ [[5], [8]] @=? MU.toLists sub3
+    , testCase "case 4" $ [[6], [9]] @=? MU.toLists sub4
+    , testCase "case 5" $ [[9]] @=? MU.toLists sub5
+    , testCase "case 6" $ [[5]] @=? MU.toLists sub6
+    , testCase "case 7" $ [[8]] @=? MU.toLists sub7
+    ]
+  where
+    ori = MU.fromLists [[1, 2, 3], [4, 5, 6], [7, 8, 9]] :: MU.Matrix Int
+    sub1 = MU.subMatrix (1,1) (2,2) ori
+    sub2 = MU.subMatrix (0,0) (0,1) sub1
+    sub3 = MU.subMatrix (0,0) (1,0) sub1
+    sub4 = MU.subMatrix (0,1) (1,1) sub1
+    sub5 = MU.subMatrix (1,1) (1,1) sub1
+    sub6 = MU.subMatrix (0,0) (0,0) sub1
+    sub7 = MU.subMatrix (1,0) (1,0) sub1