diff --git a/matrices.cabal b/matrices.cabal
--- a/matrices.cabal
+++ b/matrices.cabal
@@ -2,17 +2,17 @@
 -- documentation, see http://haskell.org/cabal/users-guide/
 
 name:                matrices
-version:             0.3.2
+version:             0.4.0
 synopsis:            native matrix based on vector
 description:         This library provide the APIs for creating, indexing,
-                     modifying matrices (2d arrays). The underling data
-                     structure is vector. It's not intended to be a linear
-                     algebra library.
+                     modifying matrices (2d arrays), including dense and
+                     sparse representations. The underling data
+                     structure is vectors.
 license:             BSD3
 license-file:        LICENSE
 author:              Kai Zhang
 maintainer:          kai@kzhang.org
-copyright:           (c) 2014 Kai Zhang
+copyright:           (c) 2015 Kai Zhang
 category:            Data
 build-type:          Simple
 -- extra-source-files:
@@ -20,16 +20,19 @@
 
 library
   exposed-modules:
-      Data.Matrix
-    , Data.Matrix.Mutable
-    , Data.Matrix.Generic
-    , Data.Matrix.Generic.Mutable
-    , Data.Matrix.Generic.Base
-    , Data.Matrix.Generic.Types
-    , Data.Matrix.Storable
-    , Data.Matrix.Storable.Mutable
-    , Data.Matrix.Unboxed
-    , Data.Matrix.Unboxed.Mutable
+    Data.Matrix
+    Data.Matrix.Mutable
+    Data.Matrix.Storable
+    Data.Matrix.Storable.Mutable
+    Data.Matrix.Unboxed
+    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
 
   ghc-options: -Wall -fwarn-tabs -funbox-strict-fields
 
@@ -39,7 +42,6 @@
       base >=4.0 && <5
     , binary
     , vector >=0.9
-    , vector-binary-instances
     , primitive
 
   hs-source-dirs:      src
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,10 @@
 module Data.Matrix
     ( Matrix
-    , module Data.Matrix.Generic.Base
-    ) where
+    , module Data.Matrix.Dense.Generic
+    )where
 
-import qualified Data.Matrix.Generic.Types as MG
-import Data.Matrix.Generic.Base
+import qualified Data.Matrix.Dense.Generic as MG
+import Data.Matrix.Dense.Generic hiding (Matrix)
 import qualified Data.Vector as V
 
-type Matrix a = MG.Matrix V.Vector a
+type Matrix = MG.Matrix V.Vector
diff --git a/src/Data/Matrix/Dense/Generic.hs b/src/Data/Matrix/Dense/Generic.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Dense/Generic.hs
@@ -0,0 +1,280 @@
+{-# LANGUAGE BangPatterns #-}
+{-# 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
+    , imap
+    , Data.Matrix.Dense.Generic.map
+    -- * Monadic mapping
+    , mapM
+    , mapM_
+    , forM
+    , forM_
+
+    -- * 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 Prelude hiding (mapM_, mapM)
+import qualified Data.Vector.Generic as G
+import Control.Arrow ((***), (&&&))
+import Control.Monad (liftM, foldM, foldM_)
+import qualified Data.Foldable as F
+import qualified Data.Vector.Generic.Mutable as GM
+
+import qualified Data.Matrix.Generic as MG
+import Data.Matrix.Dense.Generic.Mutable (MMatrix(..))
+
+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)
+
+
+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.
+    takeRow (Matrix _ c tda offset vec) i = G.slice i' c vec
+      where
+        i' = offset + i * tda
+    {-# INLINE takeRow #-}
+
+    -- | 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) . unzip . Prelude.map ((maximum *** sum) .
+                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 #-}
+
+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 #-}
+
+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 #-}
+
+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 #-}
+
+mapM_ :: (G.Vector v a, Monad m) => (a -> m b) -> Matrix v a -> m ()
+mapM_ f = G.mapM_ f . MG.flatten
+{-# INLINE mapM_ #-}
+
+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_ #-}
+
+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
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Dense/Generic/Mutable.hs
@@ -0,0 +1,47 @@
+{-# LANGUAGE FlexibleInstances #-}
+{-# 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 Prelude hiding (read, replicate)
+import Control.Monad (liftM)
+import qualified Data.Vector.Generic.Mutable as GM
+
+import qualified Data.Matrix.Generic.Mutable as C
+
+-- | mutable matrix
+data MMatrix v s a = MMatrix !Int !Int !Int !Int !(v s a)
+
+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,20 +1,168 @@
---------------------------------------------------------------------------------
--- |
--- Module      :  $Header$
--- Copyright   :  (c) 2014 Kai Zhang
--- License     :  BSD3
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE Rank2Types #-}
+module Data.Matrix.Generic
+    ( Mutable
+    , Matrix(..)
 
--- Maintainer  :  kai@kzhang.org
--- Stability   :  experimental
--- Portability :  portable
+    -- * Derived mothods
+    , rows
+    , cols
+    , (!)
+    , fromVector
+    , empty
+    , toList
+    , fromLists
+    , matrix
+    , fromRows
+    , toRows
+    , toColumns
+    , toLists
+    , create
+    ) where
 
---
---------------------------------------------------------------------------------
+import Control.Monad.ST (runST, ST)
+import Control.Monad.Primitive (PrimMonad, PrimState)
+import qualified Data.Vector.Generic as G
 
-module Data.Matrix.Generic
-    ( module Data.Matrix.Generic.Base
-    , module Data.Matrix.Generic.Types
-    ) where
+import qualified Data.Matrix.Generic.Mutable as MM
 
-import Data.Matrix.Generic.Base
-import Data.Matrix.Generic.Types
+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 `div` c, i `mod` c)
+      where
+        (r,c) = dim mat
+    {-# INLINE flatten #-}
+
+    -- | Extract a row. Default algorithm is O(n * O(unsafeIndex)).
+    takeRow :: m v a -> Int -> v a
+    takeRow mat i = G.generate c $ \j -> unsafeIndex mat (i,j)
+      where
+        (_,c) = dim mat
+    {-# INLINE takeRow #-}
+
+    -- | Extract a column. Default algorithm is O(m * O(unsafeIndex)).
+    takeColumn :: m v a -> Int -> v a
+    takeColumn mat j = G.generate r $ \i -> unsafeIndex mat (i,j)
+      where
+        (r,_) = dim mat
+    {-# INLINE takeColumn #-}
+
+    -- | 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 >= r || 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 /= G.length vec = error "incorrect length"
+                     | otherwise = unsafeFromVector (r,c) vec
+{-# INLINE fromVector #-}
+
+-- | 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 #-}
+
+-- | O(m) Return the rows
+toRows :: Matrix m v a => m v a -> [v a]
+toRows mat = map (takeRow mat) [0..r-1]
+  where
+    (r,_) = dim mat
+{-# INLINE toRows #-}
+
+-- | O(m*n) Return the columns
+toColumns :: Matrix m v a => m v a -> [v a]
+toColumns mat = map (takeColumn 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/Generic/Base.hs b/src/Data/Matrix/Generic/Base.hs
deleted file mode 100644
--- a/src/Data/Matrix/Generic/Base.hs
+++ /dev/null
@@ -1,331 +0,0 @@
-{-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE FlexibleContexts #-}
---------------------------------------------------------------------------------
--- |
--- Module      :  $Header$
--- Copyright   :  (c) 2014 Kai Zhang
--- License     :  BSD3
-
--- Maintainer  :  kai@kzhang.org
--- Stability   :  experimental
--- Portability :  portable
-
---
---------------------------------------------------------------------------------
-module Data.Matrix.Generic.Base
-    ( rows
-    , cols
-    , dim
-    , (!)
-    , unsafeIndex
-    , empty
-
-    -- * Conversions
-    , matrix
-    , flatten
-    , fromVector
-    , toRows
-    , toColumns
-    , fromRows
-    , fromColumns
-    , toList
-    , toLists
-    , fromLists
-
-    -- * Different matrix types
-    , convert
-
-    , tr
-    , takeRow
-    , takeColumn
-    , subMatrix
-    , ident
-    , diag
-    , diagRect
-    , takeDiag
-    , fromBlocks
-    , isSymmetric
-    , force
-    , Data.Matrix.Generic.Base.foldl
-
-    -- * Mapping
-    , imap
-    , Data.Matrix.Generic.Base.map
-
-    -- * Monadic mapping
-    , Data.Matrix.Generic.Base.mapM
-    , Data.Matrix.Generic.Base.mapM_
-    , Data.Matrix.Generic.Base.forM
-    , Data.Matrix.Generic.Base.forM_
-
-    , generate
-
-    ) where
-
-import Control.Arrow ((***), (&&&))
-import Control.Monad
-import qualified Data.Foldable as F
-import qualified Data.Vector.Generic as G
-import qualified Data.Vector.Generic.Mutable as GM
-import Data.Matrix.Generic.Types
-
--- | O(1) Return the number of rows
-rows :: G.Vector v a => Matrix v a -> Int
-rows (Matrix m _ _ _ _) = m
-
--- | O(1) Return the number of columns
-cols :: G.Vector v a => Matrix v a -> Int
-cols (Matrix _ n _ _ _) = n
-
--- | O(1) Return the size of matrix
-dim :: G.Vector v a => Matrix v a -> (Int, Int)
-dim (Matrix r c _ _ _) = (r,c)
-{-# INLINE dim #-}
-
--- | O(1) Indexing
-(!) :: G.Vector v a => Matrix v a -> (Int, Int) -> a
-(!) (Matrix _ _ tda offset vec) (i, j) = vec G.! idx
-  where
-    idx = offset + i * tda + j
-{-# INLINE (!) #-}
-
--- | O(1) Unsafe indexing without bound check
-unsafeIndex :: G.Vector v a => Matrix v a -> (Int, Int) -> a
-unsafeIndex (Matrix _ _ tda offset vec) (i,j) = vec `G.unsafeIndex` idx
-  where
-    idx = offset + i * tda + j
-{-# INLINE unsafeIndex #-}
-
---reshape :: G.Vector v a => Matrix v a -> (Int, Int) -> Matrix v a
-
--- | O(1) Return an empty matrix
-empty :: G.Vector v a => Matrix v a
-empty = Matrix 0 0 0 0 G.empty
-{-# INLINE empty #-}
-
--- | O(m*n) Matrix construction
-matrix :: G.Vector v a
-       => Int  -- ^ number of columns
-       -> [a]  -- ^ row list
-       -> Matrix v a
-matrix ncol xs | n `mod` ncol /= 0 = error "incorrect length"
-               | otherwise = fromVector (nrow,ncol) vec
-  where
-    vec = G.fromList xs
-    nrow = n `div` ncol
-    n = G.length vec
-{-# INLINE matrix #-}
-
--- | Create a vector by concatenating rows
-flatten :: G.Vector v a => Matrix v a -> v a
-flatten (Matrix m n tda offset vec)
-    | n == tda = G.slice offset (m * n) vec
-    | otherwise = G.generate (m * n) f
-  where
-    f i = (G.!) vec $ offset + (i `div` n) * tda + i `mod` n
-{-# INLINE flatten #-}
-
--- | O(1) Create matrix from vector
-fromVector :: G.Vector v a => (Int, Int) -> v a -> Matrix v a
-fromVector (r,c) = Matrix r c c 0
-{-# INLINE fromVector #-}
-
--- | O(m*n) Create a list by concatenating rows
-toList :: G.Vector v a => Matrix v a -> [a]
-toList = G.toList . flatten
-{-# INLINE toList #-}
-
--- | O(m) Return the rows
-toRows :: G.Vector v a => Matrix v a -> [v a]
-toRows (Matrix m n tda offset vec) = loop 0
-  where
-    loop !i | i < m = G.slice (f i) n vec : loop (i+1)
-            | otherwise = []
-    f i = offset + i * tda
-{-# INLINE toRows #-}
-
--- | O(m*n) Return the columns
-toColumns :: G.Vector v a => Matrix v a -> [v a]
-toColumns m = Prelude.map (takeColumn m) [0 .. c-1]
-  where c = cols m
-{-# INLINE toColumns #-}
-
--- | O(m*n) Create matrix from rows
-fromRows :: G.Vector v a => [v a] -> Matrix v a
-fromRows xs | any (\x -> G.length x /= c) xs = error "inequal length"
-            | otherwise = fromVector (r,c) . G.concat $ xs
-  where
-    r = length xs
-    c = G.length . head $ xs
-{-# INLINE fromRows #-}
-
--- | O(m*n) Create matrix from columns
-fromColumns :: G.Vector v a => [v a] -> Matrix v a
-fromColumns = tr . fromRows
-{-# INLINE fromColumns #-}
-
--- | O(m*n) List of lists
-toLists :: G.Vector v a => Matrix v a -> [[a]]
-toLists = Prelude.map G.toList . toRows
-{-# INLINE toLists #-}
-
--- | O(m*n) Create matrix from list of lists, it doesn't check if the list of
--- list is a valid matrix
-fromLists :: G.Vector v a => [[a]] -> Matrix v a
-fromLists xs = fromVector (r,c) . G.fromList . concat $ xs
-  where
-    r = length xs
-    c = length .head $ xs
-{-# INLINE fromLists #-}
-
----- | 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 a row
-takeRow :: G.Vector v a => Matrix v a -> Int -> v a
-takeRow (Matrix _ c tda offset vec) i = G.slice i' c vec
-  where
-    i' = offset + i * tda
-{-# INLINE takeRow #-}
-
--- | O(m) Extract a column
-takeColumn :: G.Vector v a => Matrix v a -> Int -> v a
-takeColumn (Matrix r _ tda offset vec) j = G.create $ GM.new r >>= go idx vec r 0
-  where
-    go f vec' r' !i v | i >= r' = return v
-                      | otherwise = do GM.unsafeWrite v i $ vec' G.! f i
-                                       go f vec' r' (i+1) v
-    idx i = offset + i * tda + j
-{-# INLINE takeColumn #-}
-
--- | 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 = 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) = 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 = 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 #-}
-
--- | extracts the diagonal from a rectangular matrix
-takeDiag :: G.Vector v a => Matrix v a -> v a
-takeDiag mat@(Matrix r c _ _ _) = G.generate n $ \i -> unsafeIndex mat (i,i)
-  where
-    n = min r c
-{-# INLINE takeDiag #-}
-
-fromBlocks :: G.Vector v a
-           => a               -- ^ default value
-           -> [[Matrix v a]]
-           -> Matrix v a
-fromBlocks d ms = 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 c = cols x
-                    r = rows x
-                    vec = 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) . unzip . Prelude.map ((maximum *** sum) .
-                unzip . Prelude.map (rows &&& 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 ! (i,j) == m ! (j,i)
-{-# INLINE isSymmetric #-}
-
-force :: G.Vector v a => Matrix v a -> Matrix v a
-force m@(Matrix r c _ _ _) = fromVector (r,c) . G.force . flatten $ m
-{-# INLINE force #-}
-
-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 _ _ _) = fromVector (r,c) $ G.imap f' . flatten $ m
-  where
-    f' i = f (i `div` c, i `mod` c)
-{-# INLINE imap #-}
-
-map :: (G.Vector v a, G.Vector v b) => (a -> b) -> Matrix v a -> Matrix v b
-map f m@(Matrix r c _ _ _) = fromVector (r,c) $ G.map f . flatten $ m
-{-# INLINE map #-}
-
-foldl :: G.Vector v b => (a -> b -> a) -> a -> Matrix v b -> a
-foldl f acc m = G.foldl f acc . 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 (fromVector (r,c)) . G.mapM f . flatten $ m
-{-# INLINE mapM #-}
-
-mapM_ :: (G.Vector v a, Monad m) => (a -> m b) -> Matrix v a -> m ()
-mapM_ f = G.mapM_ f . flatten
-{-# INLINE mapM_ #-}
-
-forM :: (G.Vector v a, G.Vector v b, Monad m) => Matrix v a -> (a -> m b) -> m (Matrix v b)
-forM = flip Data.Matrix.Generic.Base.mapM
-{-# INLINE forM #-}
-
-forM_ :: (G.Vector v a, Monad m) => Matrix v a -> (a -> m b) -> m ()
-forM_ = flip Data.Matrix.Generic.Base.mapM_
-{-# INLINE forM_ #-}
-
-generate :: G.Vector v a => (Int, Int) -> ((Int, Int) -> a) -> Matrix v a
-generate (r,c) f = fromVector (r,c) . G.generate (r*c) $ \i -> f (i `div` c, i `mod` 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,112 +1,51 @@
-{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE TypeFamilies #-}
+
 module Data.Matrix.Generic.Mutable
-   ( fromMVector
-   , dim
-   , flatten
-   , takeRow
-   , thaw
-   , unsafeThaw
-   , freeze
-   , unsafeFreeze
-   , write
-   , unsafeWrite
-   , read
-   , unsafeRead
-   , replicate
-   , new
-   , create
-   ) where
+    ( MMatrix(..)
+    , write
+    , read
+    ) where
 
-import Prelude hiding (read, replicate)
-import Control.Monad
-import Control.Monad.ST
-import Data.Matrix.Generic.Types
-import qualified Data.Vector.Generic as G
+import Prelude hiding (read)
+import Control.Monad.Primitive (PrimMonad, PrimState)
 import qualified Data.Vector.Generic.Mutable as GM
-import Control.Monad.Primitive
 
--- to be removed in GHC-7.10
-(<$>) :: Monad m => (a -> b) -> m a -> m b
-(<$>) = liftM
-
-dim :: MMatrix v m a -> (Int, Int)
-dim (MMatrix r c _ _ _) = (r,c)
-{-# INLINE dim #-}
-
-fromMVector :: GM.MVector v a => (Int, Int) -> v m a -> MMatrix v m a
-fromMVector (r,c) = MMatrix r c c 0
-{-# INLINE fromMVector #-}
+class GM.MVector v a => MMatrix m v a where
+    dim ::  m v s a -> (Int, Int)
 
-flatten :: (GM.MVector v a, PrimMonad m)
-        => MMatrix v (PrimState m) a -> m (v (PrimState m) a)
-flatten (MMatrix m n tda offset vec)
-    | n == tda = return $ GM.slice offset (m * n) vec
-    | otherwise = do
-        vec' <- GM.new (m*n)
-        forM_ [0 .. m*n] $ \i ->
-            GM.unsafeRead vec (offset + (i `div` n) * tda + i `mod` n) >>=
-                GM.unsafeWrite vec' i
-        return vec'
-{-# INLINE flatten #-}
+    unsafeRead :: PrimMonad s => m v (PrimState s) a -> (Int, Int) -> s a
 
-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 #-}
+    unsafeWrite :: PrimMonad s => m v (PrimState s) a -> (Int, Int) -> a -> s ()
 
-thaw :: (G.Vector v a, PrimMonad m)
-     => Matrix v a -> m (MMatrix (G.Mutable v) (PrimState m) a)
-thaw (Matrix r c tda offset v) = MMatrix r c tda offset <$> G.thaw v
-{-# INLINE thaw #-}
+    -- | Create a mutable matrix without initialization
+    new :: PrimMonad s => (Int, Int) -> s (m v (PrimState s) a)
 
-unsafeThaw :: (G.Vector v a, PrimMonad m)
-           => Matrix v a -> m (MMatrix (G.Mutable v) (PrimState m) a)
-unsafeThaw (Matrix r c tda offset v) = MMatrix r c tda offset <$> G.unsafeThaw v
-{-# INLINE unsafeThaw #-}
+    replicate :: PrimMonad s => (Int, Int) -> a -> s (m v (PrimState s) a)
 
-freeze :: (PrimMonad m, G.Vector v a) => MMatrix (G.Mutable v) (PrimState m) a -> m (Matrix v a)
-freeze (MMatrix r c tda offset v) = Matrix r c tda offset <$> G.freeze v
-{-# INLINE freeze #-}
+    {-# MINIMAL dim, unsafeRead, unsafeWrite, new, replicate #-}
 
-unsafeFreeze :: (PrimMonad m, G.Vector v a) => MMatrix (G.Mutable v) (PrimState m) a -> m (Matrix v a)
-unsafeFreeze (MMatrix r c tda offset v) = Matrix r c tda offset <$> G.unsafeFreeze v
-{-# INLINE unsafeFreeze #-}
+-- | Derived methods
 
-write :: (PrimMonad m, GM.MVector v a)
-      => MMatrix v (PrimState m) a -> (Int, Int) -> a -> m ()
-write (MMatrix _ _ tda offset v) (i,j) = GM.write v idx
-  where idx = offset + i * tda + j
+write :: (PrimMonad s, MMatrix m v a)
+      => m v (PrimState s) a -> (Int, Int) -> a -> s ()
+write mat (i,j) | i >= r || j >= c = error "Index out of bounds"
+                | otherwise = unsafeWrite mat (i,j)
+  where
+    (r,c) = dim mat
 {-# INLINE write #-}
 
-unsafeWrite :: (PrimMonad m, GM.MVector v a)
-            => MMatrix v (PrimState m) a -> (Int, Int) -> a -> m ()
-unsafeWrite (MMatrix _ _ tda offset v) (i,j) = GM.unsafeWrite v idx
-  where idx = offset + i * tda + j
-{-# INLINE unsafeWrite #-}
-
-read :: (PrimMonad m, GM.MVector v a)
-     => MMatrix v (PrimState m) a -> (Int, Int) -> m a
-read (MMatrix _ _ tda offset v) (i,j) = GM.read v idx
-  where idx = offset + i * tda + j
+read :: (PrimMonad s, MMatrix m v a)
+     => m v (PrimState s) a -> (Int, Int) -> s a
+read mat (i,j) | i >= r || j >= c = error "Index out of bounds"
+               | otherwise = unsafeRead mat (i,j)
+  where
+    (r,c) = dim mat
 {-# INLINE read #-}
 
-unsafeRead :: (PrimMonad m, GM.MVector v a)
-           => MMatrix v (PrimState m) a -> (Int, Int) -> m a
-unsafeRead (MMatrix _ _ tda offset v) (i,j) = GM.unsafeRead v idx
-  where idx = offset + i * tda + j
-{-# INLINE unsafeRead #-}
-
-replicate :: (PrimMonad m, GM.MVector v a)
-          => (Int, Int) -> a -> m (MMatrix v (PrimState m) a)
-replicate (r,c) x = fromMVector (r,c) <$> GM.replicate (r*c) x
-{-# INLINE replicate #-}
-
-new :: (PrimMonad m, GM.MVector v a)
-    => (Int, Int) -> m (MMatrix v (PrimState m) a)
-new (r,c) = fromMVector (r,c) <$> GM.new (r*c)
-{-# INLINE new #-}
-
-create :: G.Vector v a => (forall s . ST s (MMatrix (G.Mutable v) s a)) -> Matrix v a
-create m = runST $ unsafeFreeze =<< m
-{-# INLINE create #-}
+{-
+replicate :: (PrimMonad s, MMatrix m v a)
+          => (Int, Int) -> a -> s (m v (PrimState s) a)
+replicate (r,c) x = do
+    mat <- new (r,c)
+ -}   
diff --git a/src/Data/Matrix/Generic/Types.hs b/src/Data/Matrix/Generic/Types.hs
deleted file mode 100644
--- a/src/Data/Matrix/Generic/Types.hs
+++ /dev/null
@@ -1,85 +0,0 @@
-{-# LANGUAGE BangPatterns #-}
-{-# LANGUAGE FlexibleInstances #-}
-{-# LANGUAGE FlexibleContexts #-}
-{-# LANGUAGE OverlappingInstances #-}
---------------------------------------------------------------------------------
--- |
--- Module      :  $Header$
--- Copyright   :  (c) 2014 Kai Zhang
--- License     :  BSD3
-
--- Maintainer  :  kai@kzhang.org
--- Stability   :  experimental
--- Portability :  portable
-
---
---------------------------------------------------------------------------------
-
-module Data.Matrix.Generic.Types
-    ( Matrix(..)
-    , MMatrix(..)
-    ) where
-
-import Data.Binary
-import qualified Data.Vector as V
-import qualified Data.Vector.Unboxed as U
-import qualified Data.Vector.Storable as S
-import qualified Data.Vector.Binary ()
-import qualified Data.Vector.Generic as G
-
--- | 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
-
-instance Binary a => Binary (Matrix V.Vector a) where
-    put = putGeneric
-    get = getGeneric
-
-instance (U.Unbox a, Binary a) => Binary (Matrix U.Vector a) where
-    put = putGeneric
-    get = getGeneric
-
-instance (S.Storable a, Binary a) => Binary (Matrix S.Vector a) where
-    put = putGeneric
-    get = getGeneric
-
-getGeneric :: (Binary (v a), G.Vector v a) => Get (Matrix v a)
-getGeneric = do
-    r <- get
-    c <- get
-    tda <- get
-    offset <- get
-    vec <- get
-    return $ Matrix r c tda offset vec
-
-putGeneric :: (Binary (v a), G.Vector v a) => Matrix v a -> Put
-putGeneric (Matrix r c tda offset vec) = do
-    put r
-    put c
-    put tda
-    put offset
-    put vec
-
--- | mutable matrix
-data MMatrix v m a = MMatrix !Int !Int !Int !Int !(v m a)
-
-instance (G.Vector v a, Show a) => Show (Matrix v a) where
-    show mat = unlines . map (unwords . map show) . toLists $ mat
-
-instance G.Vector v Bool => Show (Matrix v Bool) where
-    show mat = unlines . map (unwords . map showBool) . toLists $ mat
-      where
-        showBool x = if x then "1" else "0"
-
-toRows :: G.Vector v a => Matrix v a -> [v a]
-toRows (Matrix m n tda offset vec) = loop 0
-    where
-    loop !i | i < m = G.slice (f i) n vec : loop (i+1)
-            | otherwise = []
-    f i = offset + i * tda
-
-toLists :: G.Vector v a => Matrix v a -> [[a]]
-toLists = map G.toList . toRows
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,10 @@
 module Data.Matrix.Mutable
     ( MMatrix
-    , module Data.Matrix.Generic.Mutable
+    , module Data.Matrix.Dense.Generic.Mutable
     ) where
 
-import qualified Data.Matrix.Generic.Types as MG
-import Data.Matrix.Generic.Mutable
+import qualified Data.Matrix.Dense.Generic.Mutable as MG
+import Data.Matrix.Dense.Generic.Mutable hiding (MMatrix)
 import qualified Data.Vector.Mutable as VM
 
 type MMatrix a = MG.MMatrix VM.MVector a
diff --git a/src/Data/Matrix/Sparse/Generic.hs b/src/Data/Matrix/Sparse/Generic.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Sparse/Generic.hs
@@ -0,0 +1,125 @@
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE TypeFamilies #-}
+module Data.Matrix.Sparse.Generic
+    ( Zero(..)
+    , CSR(..)
+    , fromAscAL
+    , fromAscStream
+    , (!)
+    ) where
+
+import Control.Monad (when, forM_)
+import Control.Monad.ST (runST)
+import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Generic.Mutable as GM
+import qualified Data.Vector.Unboxed as U
+import qualified Data.Vector.Fusion.Stream as S
+import Data.Bits (shiftR)
+
+import Data.Matrix.Generic
+import Data.Matrix.Dense.Generic.Mutable (MMatrix)
+
+class Eq a => Zero a where
+    zero :: a
+
+instance Zero Int where
+    zero = 0
+
+instance Zero Double where
+    zero = 0.0
+
+instance Eq a => Zero ([] a) where
+    zero = []
+
+-- | mutable sparse matrix not implemented
+type instance Mutable CSR = MMatrix
+
+-- | Compressed Sparse Row (CSR) matrix
+data CSR v a = CSR !Int  -- rows
+                   !Int  -- cols
+                   !(v a)  -- values
+                   !(U.Vector Int)  -- column index of values
+                   !(U.Vector Int)  -- row pointer
+    deriving (Show)
+
+instance (Zero a, G.Vector v a) => Matrix CSR v a where
+    dim (CSR r c _ _ _) = (r,c)
+    {-# INLINE dim #-}
+
+    unsafeIndex (CSR _ _ vec ci rp) (i,j) =
+        case binarySearchByBounds ci j r0 r1 of
+            Nothing -> zero
+            Just k -> vec `G.unsafeIndex` k
+      where
+        r0 = rp `U.unsafeIndex` i
+        r1 = rp `U.unsafeIndex` (i+1) - 1
+    {-# INLINE unsafeIndex #-}
+
+    unsafeFromVector (r,c) vec =
+        CSR r c (G.generate n (G.unsafeIndex vec . U.unsafeIndex nz))
+                (U.map (`mod` c) nz)
+                (U.fromList . g . U.foldr f ((r-1,n-1), [n]) $ nz)
+      where
+        nz = U.filter (\i -> vec `G.unsafeIndex` i /= zero) . U.enumFromN 0 $ (r*c)
+        f i ((!prev,!acc), xs) | stride == 0 = ((prev, acc-1), xs)
+                               | otherwise = ((current, acc-1), replicate stride (acc+1) ++ xs)
+          where
+            stride = prev - current
+            current = i `div` c
+        g ((a, _), xs) | a == 0 = 0 : xs
+                       | otherwise = replicate (a+1) 0 ++ xs
+        n = U.length nz
+    {-# INLINE unsafeFromVector #-}
+
+    thaw = undefined
+    unsafeThaw = undefined
+    freeze = undefined
+    unsafeFreeze = undefined
+
+type AssocList a = [((Int, Int), a)]
+
+fromAscAL :: G.Vector v a => (Int, Int) -> Int -> AssocList a -> CSR v a
+fromAscAL (r,c) n al = fromAscStream (r,c) n . S.fromList $ al
+{-# INLINE fromAscAL #-}
+
+fromAscStream :: (GM.MVector (G.Mutable v) a, G.Vector v a) => (Int, Int) -> Int -> S.Stream ((Int,Int),a) -> CSR v a
+fromAscStream (r,c) n al = CSR r c values ci rp
+  where
+    (values, ci, rp) = runST $ do
+        v <- GM.new n
+        col <- GM.new n
+        row <- GM.new (r+1)
+
+        let f (i',acc) ((i,j),x) = do
+                GM.write v acc x
+                GM.write col acc j
+
+                let stride = i - i'
+                when (stride > 0) $ forM_ [0..stride-1] $ \s -> GM.write row (i-s) acc
+                
+                return (i,acc+1)
+
+        _ <- S.foldM f (0,0) al
+        GM.write row r n
+
+        v' <- G.unsafeFreeze v
+        col' <- G.unsafeFreeze col
+        row' <- G.unsafeFreeze row
+        return (v', col', row')
+{-# INLINE fromAscStream #-}
+    
+binarySearchByBounds :: U.Vector Int -> Int -> Int -> Int -> Maybe Int
+binarySearchByBounds vec x = loop
+  where
+    loop !l !u
+        | l > u = Nothing
+        | x == x' = Just k
+        | x < x' = loop l (k-1)
+        | otherwise = loop (k+1) u
+      where
+        k = (u+l) `shiftR` 1
+        x' = vec `U.unsafeIndex` k
+{-# INLINE binarySearchByBounds #-}
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,10 @@
 module Data.Matrix.Storable
     ( Matrix
-    , module Data.Matrix.Generic.Base
+    , module Data.Matrix.Dense.Generic
     ) where
 
-import qualified Data.Matrix.Generic.Types as MG
-import Data.Matrix.Generic.Base
-import qualified Data.Vector.Storable as S
+import qualified Data.Matrix.Dense.Generic as MG
+import Data.Matrix.Dense.Generic hiding (Matrix)
+import qualified Data.Vector.Storable as V
 
-type Matrix a = MG.Matrix S.Vector a
+type Matrix a = MG.Matrix V.Vector a
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,10 @@
 module Data.Matrix.Storable.Mutable
     ( MMatrix
-    , module Data.Matrix.Generic.Mutable
+    , module Data.Matrix.Dense.Generic.Mutable
     ) where
 
-import qualified Data.Matrix.Generic.Types as MG
-import Data.Matrix.Generic.Mutable
-import qualified Data.Vector.Storable.Mutable as SM
+import qualified Data.Matrix.Dense.Generic.Mutable as MG
+import Data.Matrix.Dense.Generic.Mutable hiding (MMatrix)
+import qualified Data.Vector.Storable.Mutable as VM
 
-type MMatrix a = MG.MMatrix SM.MVector a
+type MMatrix a = MG.MMatrix VM.MVector a
diff --git a/src/Data/Matrix/Symmetric.hs b/src/Data/Matrix/Symmetric.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Symmetric.hs
@@ -0,0 +1,96 @@
+{-# LANGUAGE MultiParamTypeClasses #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE TypeFamilies #-}
+module Data.Matrix.Symmetric
+    ( SymMatrix(..)
+    , dim
+    , rows
+    , cols
+    , unsafeIndex
+    , (!)
+    , flatten
+    , unsafeFromVector
+    , fromVector
+    , takeRow
+    , thaw
+    , unsafeThaw
+    , freeze
+    , unsafeFreeze
+    , zip
+    , zipWith
+    ) where
+
+import Prelude hiding (zip, zipWith)
+import Control.Monad (liftM)
+import Data.Bits (shiftR)
+import Data.Binary
+import qualified Data.Vector.Generic as G
+
+import Data.Matrix.Generic
+import Data.Matrix.Symmetric.Mutable (SymMMatrix(..))
+
+type instance Mutable SymMatrix = SymMMatrix
+
+-- | Symmetric square matrix
+data SymMatrix v a = SymMatrix !Int !(v a)
+    deriving (Show)
+
+instance (G.Vector v a, Binary a) => Binary (SymMatrix v a) where
+    put (SymMatrix n vec) = do put (52 :: Int)
+                               put n
+                               G.mapM_ put vec
+    get = do
+        magic <- get :: Get Int
+        if magic == 52
+           then do
+               n <- get 
+               vec <- G.replicateM ((n+1)*n `shiftR` 1) get
+               return $ SymMatrix n vec
+           else error "not a valid file"
+
+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 . 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 #-}
+
+-- 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 #-}
+
+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 #-}
diff --git a/src/Data/Matrix/Symmetric/Mutable.hs b/src/Data/Matrix/Symmetric/Mutable.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Symmetric/Mutable.hs
@@ -0,0 +1,46 @@
+{-# 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 Prelude hiding (read, replicate)
+import Control.Monad (liftM)
+import Data.Bits (shiftR)
+import qualified Data.Vector.Generic.Mutable as GM
+
+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,10 @@
 module Data.Matrix.Unboxed
     ( Matrix
-    , module Data.Matrix.Generic.Base
+    , module Data.Matrix.Dense.Generic
     ) where
 
-import qualified Data.Matrix.Generic.Types as MG
-import Data.Matrix.Generic.Base
+import qualified Data.Matrix.Dense.Generic as MG
+import Data.Matrix.Dense.Generic hiding (Matrix)
 import qualified Data.Vector.Unboxed as V
 
 type Matrix a = MG.Matrix V.Vector a
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,10 @@
 module Data.Matrix.Unboxed.Mutable
     ( MMatrix
-    , module Data.Matrix.Generic.Mutable
+    , module Data.Matrix.Dense.Generic.Mutable
     ) where
 
-import qualified Data.Matrix.Generic.Types as MG
-import Data.Matrix.Generic.Mutable
-import qualified Data.Vector.Unboxed.Mutable as UM
+import qualified Data.Matrix.Dense.Generic.Mutable as MG
+import Data.Matrix.Dense.Generic.Mutable hiding (MMatrix)
+import qualified Data.Vector.Unboxed.Mutable as VM
 
-type MMatrix a = MG.MMatrix UM.MVector a
+type MMatrix a = MG.MMatrix VM.MVector a