diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c) 2014, Kai Zhang
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Kai Zhang nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/matrices.cabal b/matrices.cabal
new file mode 100644
--- /dev/null
+++ b/matrices.cabal
@@ -0,0 +1,40 @@
+-- Initial matrices.cabal generated by cabal init.  For further 
+-- documentation, see http://haskell.org/cabal/users-guide/
+
+name:                matrices
+version:             0.1.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.
+license:             BSD3
+license-file:        LICENSE
+author:              Kai Zhang
+maintainer:          kai@kzhang.org
+copyright:           <c> 2014 Kai Zhang
+category:            Data
+build-type:          Simple
+-- extra-source-files:  
+cabal-version:       >=1.10
+
+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.Unboxed
+    , Data.Matrix.Unboxed.Mutable
+
+  -- other-modules:       
+
+  build-depends:
+      base >=4.0 && <4.8
+    , vector >=0.9
+    , primitive
+
+  hs-source-dirs:      src
+  default-language:    Haskell2010
diff --git a/src/Data/Matrix.hs b/src/Data/Matrix.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix.hs
@@ -0,0 +1,10 @@
+module Data.Matrix
+    ( Matrix
+    , module Data.Matrix.Generic.Base
+    ) where
+
+import qualified Data.Matrix.Generic.Types as MG
+import Data.Matrix.Generic.Base
+import qualified Data.Vector as V
+
+type Matrix a = MG.Matrix V.Vector a
diff --git a/src/Data/Matrix/Generic.hs b/src/Data/Matrix/Generic.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Generic.hs
@@ -0,0 +1,21 @@
+{-# LANGUAGE GADTs #-}
+--------------------------------------------------------------------------------
+-- |
+-- Module      :  $Header$
+-- Copyright   :  (c) 2014 Kai Zhang
+-- License     :  BSD3
+
+-- Maintainer  :  kai@kzhang.org
+-- Stability   :  experimental
+-- Portability :  portable
+
+--
+--------------------------------------------------------------------------------
+
+module Data.Matrix.Generic
+    ( module Data.Matrix.Generic.Base
+    , module Data.Matrix.Generic.Types
+    ) where
+
+import Data.Matrix.Generic.Base
+import Data.Matrix.Generic.Types
diff --git a/src/Data/Matrix/Generic/Base.hs b/src/Data/Matrix/Generic/Base.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Generic/Base.hs
@@ -0,0 +1,162 @@
+{-# LANGUAGE GADTs #-}
+{-# 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
+    , (!)
+    , matrix
+    , flatten
+    , fromVector
+    , toRows
+    , toList
+    , toLists
+    , tr
+    , subMatrix
+    , ident
+    , diag
+    , diagRect
+    , fromBlocks
+    , isSymmetric
+    ) 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
+
+rows :: G.Vector v a => Matrix v a -> Int
+rows (Matrix m _ _ _ _) = m
+
+cols :: G.Vector v a => Matrix v a -> Int
+cols (Matrix _ n _ _ _) = n
+
+(!) :: 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 (!) #-}
+
+matrix :: G.Vector v a => Int -> [a] -> 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 #-}
+
+flatten :: 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 #-}
+
+fromVector :: G.Vector v a => Int -> Int -> v a -> Matrix v a
+fromVector r c = Matrix r c c 0
+{-# INLINE fromVector #-}
+
+toList :: G.Vector v a => Matrix v a -> [a]
+toList = G.toList . flatten
+{-# INLINE toList #-}
+
+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 #-}
+
+toLists :: G.Vector v a => Matrix v a -> [[a]]
+toLists = map G.toList . toRows
+{-# INLINE toLists #-}
+
+subMatrix :: G.Vector v a => (Int, Int) -> (Int, Int) -> Matrix v a -> Matrix v a
+subMatrix (ri, rj) (ci, cj) (Matrix _ n tda offset vec) =
+    Matrix m' n' tda offset' vec
+  where
+    m' = rj - ri + 1
+    n' = cj - ci + 1
+    offset' = offset + ri * n + ci
+{-# INLINE subMatrix #-}
+
+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 #-}
+
+ident :: (Num a, G.Vector v a) => Int -> Matrix v a
+ident n = diagRect 0 n n $ replicate n 1
+{-# INLINE ident #-}
+
+-- | 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 #-}
+
+-- | create a rectangular matrix with default values and given diagonal
+diagRect :: (G.Vector v a, F.Foldable t)
+         => a         -- ^ default value
+         -> Int       -- ^ number of rows
+         -> Int       -- ^ number of columns
+         -> 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 #-}
+
+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 . map ((maximum *** sum) . unzip . 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 #-}
diff --git a/src/Data/Matrix/Generic/Mutable.hs b/src/Data/Matrix/Generic/Mutable.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Generic/Mutable.hs
@@ -0,0 +1,62 @@
+{-# LANGUAGE GADTs #-}
+module Data.Matrix.Generic.Mutable
+   ( thaw
+   , unsafeThaw
+   , freeze
+   , unsafeFreeze
+   , write
+   , unsafeWrite
+   , read
+   , unsafeRead
+   ) where
+
+import Prelude hiding (read)
+import Control.Monad
+import Data.Matrix.Generic.Types
+import qualified Data.Vector.Generic as G
+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
+
+thaw :: 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 #-}
+
+unsafeThaw :: 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 #-}
+
+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 #-}
+
+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 #-}
+
+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
+{-# 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
+{-# 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 #-}
diff --git a/src/Data/Matrix/Generic/Types.hs b/src/Data/Matrix/Generic/Types.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Generic/Types.hs
@@ -0,0 +1,63 @@
+{-# LANGUAGE GADTs #-}
+{-# 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 qualified Data.Vector.Generic as G
+import qualified Data.Vector.Generic.Mutable as GM
+
+-- | row-major matrix supporting efficient slice
+data Matrix v a where
+    Matrix :: G.Vector v a
+           => !Int    -- ^ number of rows
+           -> !Int    -- ^ number of cols
+           -> !Int    -- ^ physical row dimension
+           -> !Int    -- ^ offset
+           -> !(v a)  -- ^ flat matrix
+           -> Matrix v a
+
+-- | mutable matrix
+data MMatrix v m a where
+    MMatrix :: GM.MVector v a
+            => !Int
+            -> !Int
+            -> !Int
+            -> !Int
+            -> !(v m a)
+            -> MMatrix 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
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Mutable.hs
@@ -0,0 +1,10 @@
+module Data.Matrix.Mutable
+    ( MMatrix
+    , module Data.Matrix.Generic.Mutable
+    ) where
+
+import qualified Data.Matrix.Generic.Types as MG
+import Data.Matrix.Generic.Mutable
+import qualified Data.Vector.Mutable as VM
+
+type MMatrix a = MG.MMatrix VM.MVector a
diff --git a/src/Data/Matrix/Unboxed.hs b/src/Data/Matrix/Unboxed.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Unboxed.hs
@@ -0,0 +1,10 @@
+module Data.Matrix.Unboxed
+    ( Matrix
+    , module Data.Matrix.Generic.Base
+    ) where
+
+import qualified Data.Matrix.Generic.Types as MG
+import Data.Matrix.Generic.Base
+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
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Unboxed/Mutable.hs
@@ -0,0 +1,10 @@
+module Data.Matrix.Unboxed.Mutable
+    ( MMatrix
+    , module Data.Matrix.Generic.Mutable
+    ) where
+
+import qualified Data.Matrix.Generic.Types as MG
+import Data.Matrix.Generic.Mutable
+import qualified Data.Vector.Unboxed.Mutable as UM
+
+type MMatrix a = MG.MMatrix UM.MVector a