diff --git a/matrix-sized.cabal b/matrix-sized.cabal
--- a/matrix-sized.cabal
+++ b/matrix-sized.cabal
@@ -1,13 +1,13 @@
 cabal-version: 2.2
 
--- This file has been generated from package.yaml by hpack version 0.31.2.
+-- This file has been generated from package.yaml by hpack version 0.33.0.
 --
 -- see: https://github.com/sol/hpack
 --
--- hash: 2b1169a59b262a2ee39daffd72cfe59b2b4f716e82ef7df496b186c155e1b003
+-- hash: 856dabd5492c97e6a433e61a9b3b39337897ded09cf66a9702ba3010a9ee0a83
 
 name:           matrix-sized
-version:        0.1.0
+version:        0.1.1
 synopsis:       Haskell matrix library with interface to C++ linear algebra libraries.
 description:    A Haskell implementation of matrices with statically known sizes. The library also comes with the bindings to high performance C++ linear algebra libraries such as Eigen and Spectra.
 category:       Math
@@ -391,6 +391,7 @@
       Data.Matrix.Static.Generic
       Data.Matrix.Static.Generic.Mutable
       Data.Matrix.Static.IO
+      Data.Matrix.Dynamic
   other-modules:
       Data.Matrix.Static.Internal
       Data.Matrix.Static.LinearAlgebra.Internal
@@ -413,7 +414,7 @@
     , bytestring
     , bytestring-lexing
     , conduit
-    , conduit-extra
+    , double-conversion
     , primitive >=0.6.4.0
     , singletons
     , store
diff --git a/src/Data/Matrix/Dynamic.hs b/src/Data/Matrix/Dynamic.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Dynamic.hs
@@ -0,0 +1,93 @@
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE Rank2Types            #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE KindSignatures #-}
+{-# LANGUAGE FlexibleContexts #-}
+module Data.Matrix.Dynamic
+    ( Dynamic(..)
+    , withDyn
+    , matrix
+    , fromList
+    , fromVector
+    , fromColumns
+    , fromRows
+    , fromTriplet
+    , decodeSparse
+    )where
+
+import Data.ByteString (ByteString)
+import qualified Data.Matrix.Static.Sparse as S
+import qualified Data.Matrix.Static.Generic as C
+import qualified Data.Vector.Generic         as G
+import Data.Kind (Type)
+import Data.Singletons
+import Data.Singletons.TypeLits
+import Data.Store (Store(..), decodeExWith)
+
+data Dynamic (m :: C.MatrixKind) (v :: Type -> Type) a where
+    Dynamic :: m r c v a -> Dynamic m v a
+
+withDyn :: Dynamic m v a -> (forall r c. m r c v a -> b) -> b
+withDyn (Dynamic x) f = f x
+{-# INLINE withDyn #-}
+
+matrix :: forall m v a. C.Matrix m v a => [[a]] -> Dynamic m v a
+matrix lists = withSomeSing (fromIntegral r) $ \(SNat :: Sing r) ->
+        withSomeSing (fromIntegral c) $ \(SNat :: Sing c) ->
+            Dynamic (C.matrix lists :: m r c v a)
+  where
+    r = length lists
+    c = length $ head lists
+{-# INLINE matrix #-}
+
+-- | Construct matrix from a list containg columns.
+fromList :: forall v a m. C.Matrix m v a
+         => (Int, Int) -> [a] -> Dynamic m v a
+fromList d = fromVector d . G.fromList
+{-# INLINE fromList #-}
+
+-- | Construct matrix from a vector containg columns.
+fromVector :: forall v a m. C.Matrix m v a
+           => (Int, Int) -> v a -> Dynamic m v a
+fromVector (r, c) vec = withSomeSing (fromIntegral r) $ \(SNat :: Sing r) ->
+        withSomeSing (fromIntegral c) $ \(SNat :: Sing c) ->
+            Dynamic (C.fromVector vec :: m r c v a)
+{-# INLINE fromVector #-}
+
+-- | O(m*n) Create matrix from rows
+fromRows :: forall m v a. C.Matrix m v a
+          => [v a] -> Dynamic m v a
+fromRows xs = withSomeSing r $ \(SNat :: Sing r) -> 
+    withSomeSing c $ \(SNat :: Sing c) -> Dynamic (C.fromRows xs :: m r c v a)
+  where
+    r = fromIntegral $ length xs
+    c = fromIntegral $ G.length $ head xs
+{-# INLINE fromRows #-}
+
+fromColumns :: forall m v a. C.Matrix m v a
+            => [v a] -> Dynamic m v a
+fromColumns xs = withSomeSing r $ \(SNat :: Sing r) -> 
+    withSomeSing c $ \(SNat :: Sing c) -> Dynamic (C.fromColumns xs :: m r c v a)
+  where
+    c = fromIntegral $ length xs
+    r = fromIntegral $ G.length $ head xs
+{-# INLINE fromColumns #-}
+
+fromTriplet :: forall u v a. (G.Vector u (Int, Int, a), G.Vector v a)
+            => (Int, Int) -> u (Int, Int, a) -> Dynamic S.SparseMatrix v a
+fromTriplet (r, c) triplets = withSomeSing (fromIntegral r) $ \(SNat :: Sing r) -> 
+    withSomeSing (fromIntegral c) $ \(SNat :: Sing c) ->
+        Dynamic (S.fromTriplet triplets :: S.SparseMatrix r c v a)
+{-# INLINE fromTriplet #-}
+
+decodeSparse :: forall v a. (Store (v a), G.Vector v a)
+             => ByteString -> Dynamic S.SparseMatrix v a
+decodeSparse bs = withSomeSing (fromIntegral (r :: Int)) $ \(SNat :: Sing r) ->
+    withSomeSing (fromIntegral (c :: Int)) $ \(SNat :: Sing c) ->
+        Dynamic (S.SparseMatrix nnz inner outer :: S.SparseMatrix r c v a)
+  where
+    (r,c,nnz,inner,outer) = decodeExWith 
+        ((,,,,) <$> peek <*> peek <*> peek <*> peek <*> peek)
+        bs
+{-# INLINE decodeSparse #-}
diff --git a/src/Data/Matrix/Static/Dense.hs b/src/Data/Matrix/Static/Dense.hs
--- a/src/Data/Matrix/Static/Dense.hs
+++ b/src/Data/Matrix/Static/Dense.hs
@@ -105,7 +105,7 @@
 import Data.Tuple (swap)
 import qualified Data.List as L
 import Text.Printf (printf)
-import Data.Store (Store(..), Size(..), decodeExWith)
+import Data.Store (Store(..), Size(..))
 import Foreign.Storable (sizeOf)
 
 import           Data.Matrix.Static.Dense.Mutable (MMatrix (..))
@@ -123,6 +123,7 @@
         size = VarSize $ \(Matrix vec) -> case size of
             VarSize f  ->
                 2 * sizeOf (0 :: Int) + f vec
+            _ -> undefined
 
         poke mat@(Matrix vec) = poke r >> poke c >> poke vec
           where
diff --git a/src/Data/Matrix/Static/Generic.hs b/src/Data/Matrix/Static/Generic.hs
--- a/src/Data/Matrix/Static/Generic.hs
+++ b/src/Data/Matrix/Static/Generic.hs
@@ -11,7 +11,7 @@
     , Matrix(..)
     , MatrixKind
 
-    -- * Derived mothods
+    -- * Matrix query
     , rows
     , cols
     , (!)
@@ -19,15 +19,15 @@
     , takeRow
     , toRows
     , toColumns
+
+    -- * Matrix Construction
     , empty
     , matrix
-    , withMatrix
     , fromRows
-    , withRows
     , fromColumns
-    , withColumns
     , fromVector
     , fromList
+
     , toList
     , create
     , convertAny
@@ -44,8 +44,7 @@
 import Data.Tuple (swap)
 import Data.Kind (Type)
 import GHC.TypeLits (Nat, type (<=))
-import Data.Singletons (SingI, Sing, fromSing, sing, withSomeSing)
-import Data.Singletons.TypeLits
+import Data.Singletons (SingI, Sing, fromSing, sing)
 
 import Data.Matrix.Static.Generic.Mutable (MMatrix, MMatrixKind)
 
@@ -158,15 +157,6 @@
 matrix = fromList . concat . L.transpose
 {-# INLINE matrix #-}
 
-withMatrix :: forall mat v a b. Matrix mat v a
-           => [[a]] -> (forall r c. mat r c v a -> b) -> b
-withMatrix xs f = withSomeSing n $ \(SNat :: Sing n) -> 
-    withSomeSing m $ \(SNat :: Sing m) -> f (matrix xs :: mat n m v a)
-  where
-    n = fromIntegral $ length xs
-    m = fromIntegral $ length $ head xs
-{-# INLINE withMatrix #-}
-
 -- | Construct matrix from a list containg columns.
 fromList :: (SingI r, SingI c, Matrix m v a)
          => [a] -> m r c v a
@@ -178,29 +168,11 @@
 fromRows = transpose . fromColumns
 {-# INLINE fromRows #-}
 
-withRows :: forall mat v a b. Matrix mat v a
-          => [v a] -> (forall r c. mat r c v a -> b) -> b
-withRows xs f = withSomeSing n $ \(SNat :: Sing n) -> 
-    withSomeSing m $ \(SNat :: Sing m) -> f (fromRows xs :: mat n m v a)
-  where
-    n = fromIntegral $ length xs
-    m = fromIntegral $ G.length $ head xs
-{-# INLINE withRows #-}
-
 -- | O(m*n) Create matrix from columns
 fromColumns :: (Matrix m v a, SingI r, SingI c)
             => [v a] -> m r c v a
 fromColumns = fromVector . G.concat
 {-# INLINE fromColumns #-}
-
-withColumns :: forall mat v a b. Matrix mat v a
-            => [v a] -> (forall r c. mat r c v a -> b) -> b
-withColumns xs f = withSomeSing n $ \(SNat :: Sing n) -> 
-    withSomeSing m $ \(SNat :: Sing m) -> f (fromRows xs :: mat n m v a)
-  where
-    m = fromIntegral $ length xs
-    n = fromIntegral $ G.length $ head xs
-{-# INLINE withColumns #-}
 
 -- | O(m*n) Create a list by concatenating columns
 toList :: Matrix m v a => m r c v a -> [a]
diff --git a/src/Data/Matrix/Static/IO.hs b/src/Data/Matrix/Static/IO.hs
--- a/src/Data/Matrix/Static/IO.hs
+++ b/src/Data/Matrix/Static/IO.hs
@@ -15,18 +15,24 @@
 
 module Data.Matrix.Static.IO
     ( fromMM
-    , withMM
+    , fromMM'
+    , toMM
+    , IOElement(..)
     ) where
 
 import qualified Data.ByteString.Char8 as B
 import Conduit
 import Control.Monad (when)
 import qualified Data.Vector.Generic as G
+import Data.Matrix.Dynamic (Dynamic(..))
+import qualified Data.Vector.Unboxed as U
 import           Data.ByteString.Lex.Fractional (readExponential, readSigned)
 import           Data.ByteString.Lex.Integral   (readDecimal, readDecimal_)
 import Data.Singletons
 import Data.Maybe
 import Data.Singletons.TypeLits
+import Data.Double.Conversion.ByteString (toShortest)
+import Text.Printf (printf)
 
 import qualified Data.Matrix.Static.Sparse as S
 
@@ -36,63 +42,84 @@
             | MMPattern
             deriving (Eq)
 
-class IOElement a where
+class U.Unbox a => IOElement a where
     decodeElem :: B.ByteString -> a
+    encodeElem :: a -> B.ByteString
     elemType :: Proxy a -> MMElem
 
 instance IOElement Int where
     decodeElem x = fst . fromMaybe errMsg . readSigned readDecimal $ x
       where
         errMsg = error $ "readInt: Fail to cast ByteString to Int:" ++ show x
+    encodeElem = B.pack . show
     elemType _ = MMInteger
 
 instance IOElement Double where
     decodeElem x = fst . fromMaybe errMsg . readSigned readExponential $ x
       where
         errMsg = error $ "readDouble: Fail to cast ByteString to Double:" ++ show x
+    encodeElem = toShortest
     elemType _ = MMReal
 
-withMM :: forall m v a b. (Monad m, G.Vector v a, IOElement a)
-       => ConduitT () B.ByteString m ()
-       -> (forall r c. S.SparseMatrix r c v a -> m b)
-       -> m b
-withMM conduit f = do
-    (ty, (r,c,nnz)) <- parseHeader conduit
+fromMM' :: forall o m v a. (PrimMonad m, G.Vector v a, IOElement a)
+        => ConduitT B.ByteString o m (Dynamic S.SparseMatrix v a)
+fromMM' = linesUnboundedAsciiC .| do
+    (ty, (r,c,nnz)) <- parseHeader
     when (elemType (Proxy :: Proxy a) /= ty) $ error "Element types do not match"
+    vec <- streamTriplet .| sinkVector 
+    when (U.length vec /= nnz) $ error $
+        "number of non-zeros do not match: " <> show nnz <> "/=" <> show (U.length vec)
     withSomeSing (fromIntegral (r :: Int)) $ \(SNat :: Sing r) ->
-        withSomeSing (fromIntegral (c :: Int)) $ \(SNat :: Sing c) -> do
-            mat@(S.SparseMatrix v _ _) <- S.fromTripletC triplets
-            let n = G.length v
-            when (n /= nnz) $ error $
-                "number of non-zeros do not match: " <> show nnz <> "/=" <> show n
-            f (mat :: S.SparseMatrix r c v a)
-  where
-    triplets = conduit .| linesUnboundedAsciiC .| filterC (not . (=='%') . B.head) .|
-        (dropC 1 >> streamTriplet)
+        withSomeSing (fromIntegral (c :: Int)) $ \(SNat :: Sing c) ->
+            return $ Dynamic (S.fromTriplet vec :: S.SparseMatrix r c v a)
 
-fromMM :: forall m r c v a. (Monad m, SingI r, SingI c, G.Vector v a, IOElement a)
-       => ConduitT () B.ByteString m () -> m (S.SparseMatrix r c v a)
-fromMM conduit = do
-    (ty, (r,c,nnz)) <- parseHeader conduit
+fromMM :: forall o m r c v a. (PrimMonad m, SingI r, SingI c, G.Vector v a, IOElement a)
+       => ConduitT B.ByteString o m (S.SparseMatrix r c v a)
+fromMM = linesUnboundedAsciiC .| do
+    (ty, (r,c,nnz)) <- parseHeader
     mat@(S.SparseMatrix v _ _) <- case () of
         _ | elemType (Proxy :: Proxy a) /= ty -> error "Element types do not match"
           | (r, c) /= (nrow, ncol) -> error $ "Dimensions do not match: " <>
                 show (r,c) <> "/=" <> show (nrow,ncol)
-          | otherwise -> S.fromTripletC triplets
+          | otherwise -> do
+              vec <- streamTriplet .| sinkVector 
+              return $ S.fromTriplet (vec :: U.Vector (Int, Int, a))
     let n = G.length v
     if n /= nnz
         then error $ "number of non-zeros do not match: " <> show nnz <> "/=" <> show n
         else return mat
   where
-    triplets = conduit .| linesUnboundedAsciiC .| filterC (not . (=='%') . B.head) .|
-        (dropC 1 >> streamTriplet)
     nrow = fromIntegral $ fromSing (sing :: Sing r) :: Int
     ncol = fromIntegral $ fromSing (sing :: Sing c) :: Int
 
-parseHeader :: Monad m => ConduitT () B.ByteString m () -> m (MMElem, (Int, Int, Int))
-parseHeader conduit = runConduit $
-    conduit .| linesUnboundedAsciiC .| headerParser
+toMM :: forall m r c v a i. (Monad m, S.Zero a, IOElement a, G.Vector v a)
+     => S.SparseMatrix r c v a -> ConduitT i B.ByteString m ()
+toMM mat@(S.SparseMatrix vec _ _) = ( do
+    yield header
+    yield "%"
+    yield $ B.pack $ printf "%d %d %d" r c n
+    S.toTriplet mat .| mapC f ) .| unlinesAsciiC
   where
+    f (i, j, x) = B.unwords [B.pack $ show (i+1), B.pack $ show (j+1), encodeElem x]
+    header = case elemType (Proxy :: Proxy a) of
+        MMReal -> "%%MatrixMarket matrix coordinate real general"
+        MMInteger -> "%%MatrixMarket matrix coordinate integer general"
+        _ -> undefined
+    (r, c) = S.dim mat
+    n = G.length vec
+
+parseHeader :: Monad m => ConduitT B.ByteString o m (MMElem, (Int, Int, Int))
+parseHeader = do
+    ty <- headC >>= \case
+        Nothing -> error "Empty file"
+        Just header -> return $ parse header
+    dropWhileC $ (=='%') . B.head 
+    headC >>= \case
+        Nothing -> error "Empty file"
+        Just x ->
+            let [r, c, nnz] = map decodeElem $ B.words x
+            in return (ty, (r, c, nnz))
+  where
     parse x
         | "%%MatrixMarket" `B.isPrefixOf` x = case B.words x of
             [_, _, format, ty, form] -> 
@@ -105,21 +132,11 @@
                 in ty'
             _ -> error $ "Cannot parse header: " <> show x
         | otherwise = error $ "Cannot parse header: " <> show x
-    headerParser = do
-        ty <- headC >>= \case
-            Nothing -> error "Empty file"
-            Just header -> return $ parse header
-        line <- filterC (not . (=='%') . B.head) .| headC
-        case line of
-            Nothing -> error "Empty file"
-            Just x ->
-                let [r, c, nnz] = map decodeElem $ B.words x
-                in return (ty, (r, c, nnz))
 {-# INLINE parseHeader #-}
 
 streamTriplet :: (Monad m, IOElement a) => ConduitT B.ByteString (Int, Int, a) m ()
 streamTriplet = mapC (f . B.words)
   where
-    f [i,j,x] = (readDecimal_ i, readDecimal_ j, decodeElem x)
+    f [i,j,x] = (readDecimal_ i - 1, readDecimal_ j - 1, decodeElem x)
     f x = error $ "Formatting error: " <> show x
 {-# INLINE streamTriplet #-}
diff --git a/src/Data/Matrix/Static/Sparse.hs b/src/Data/Matrix/Static/Sparse.hs
--- a/src/Data/Matrix/Static/Sparse.hs
+++ b/src/Data/Matrix/Static/Sparse.hs
@@ -38,7 +38,6 @@
     , fromTriplet
     , fromTripletC
     , toTriplet
-    , withDecodedMatrix
     , C.fromVector
     , C.fromList
     , C.unsafeFromVector
@@ -69,10 +68,8 @@
 import GHC.TypeLits (type (<=))
 import Foreign.C.Types
 import Data.Complex
-import Data.Store (Store(..), Size(..), decodeExWith)
+import Data.Store (Store(..), Size(..))
 import Foreign.Storable (sizeOf)
-import Data.ByteString (ByteString)
-import Data.Singletons.TypeLits
 
 import qualified Data.Matrix.Static.Dense as D
 import qualified Data.Matrix.Static.Dense.Mutable as DM
@@ -123,6 +120,7 @@
             case (size, size) of
                 (VarSize f, VarSize g) ->
                     2 * sizeOf (0 :: Int) + f nnz + g inner + g outer
+                _ -> undefined
 
         poke mat@(SparseMatrix nnz inner outer) = poke r >> poke c >> poke nnz >>
             poke inner >> poke outer
@@ -138,17 +136,6 @@
             r = fromIntegral $ fromSing (sing :: Sing r) :: Int
             c = fromIntegral $ fromSing (sing :: Sing c) :: Int
 
-withDecodedMatrix :: forall v a b. (G.Vector v a, Store (v a))
-                  => ByteString -> (forall r c. SparseMatrix r c v a -> b) -> b
-withDecodedMatrix bs f = withSomeSing (fromIntegral (r :: Int)) $ \(SNat :: Sing r) ->
-    withSomeSing (fromIntegral (c :: Int)) $ \(SNat :: Sing c) ->
-        f (SparseMatrix nnz inner outer :: SparseMatrix r c v a)
-  where
-    (r,c,nnz,inner,outer) = decodeExWith 
-        ((,,,,) <$> peek <*> peek <*> peek <*> peek <*> peek)
-        bs
-{-# INLINE withDecodedMatrix #-}
-
 instance (G.Vector v a, Eq (v a)) => Eq (SparseMatrix r c v a) where
     (==) (SparseMatrix a b c) (SparseMatrix a' b' c') =
         a == a' && b == b' && c == c'
@@ -261,8 +248,9 @@
   where
     outer = S.scanl (+) 0 $ S.create $ do
         vec <- SM.replicate c 0
-        G.forM_ triplets $ \(_, j, _) -> 
-            SM.unsafeModify vec (+1) j
+        G.forM_ triplets $ \(i, j, _) -> if i < r && j < c
+            then SM.unsafeModify vec (+1) j
+            else error $ printf "Index out of bound: (%d, %d) >= (%d, %d)" i j r c 
         return vec
     (val, inner) = runST $ do
         outer' <- S.thaw outer
@@ -275,6 +263,7 @@
             SM.unsafeModify outer' (+1) j
         (,) <$> G.unsafeFreeze val' <*> S.unsafeFreeze inner'
     nnz = G.length triplets
+    r = fromIntegral $ fromSing (sing :: Sing r)
     c = fromIntegral $ fromSing (sing :: Sing c)
 {-# INLINE fromTriplet #-}
 
diff --git a/tests/Test/Base.hs b/tests/Test/Base.hs
--- a/tests/Test/Base.hs
+++ b/tests/Test/Base.hs
@@ -13,11 +13,12 @@
 import qualified Data.Matrix.Static.Generic as G
 import qualified Data.Matrix.Static.Dense as D
 import qualified Data.Matrix.Static.Sparse as S
+import qualified Data.Matrix.Dynamic as Dyn
+import Control.Monad.ST (runST)
+import Data.Matrix.Static.IO
 import Data.Singletons hiding ((@@))
-import Data.Singletons.Prelude (Min)
 import Data.Vector (Vector)
 import qualified Data.Vector as V
-import Control.Monad.IO.Class (liftIO)
 import Test.Tasty.QuickCheck
 import Conduit
 import Data.Store
@@ -35,14 +36,18 @@
 pSerialization = testGroup "Serialization"
     [ testProperty "Dense: id == decode . encode" tStoreD
     , testProperty "Sparse: id == decode . encode" tStoreS
-    , testProperty "Sparse: id == decode . encode" tStoreS' ]
+    , testProperty "Sparse: id == decode . encode" tStoreS'
+    , testProperty "Sparse: id ~= fromMM . toMM" tMM ]
   where
     tStoreD :: D.Matrix 80 60 Vector Double -> Bool
     tStoreD mat = mat == decodeEx (encode mat)
     tStoreS :: S.SparseMatrix 80 60 Vector Double -> Bool
     tStoreS mat = mat == decodeEx (encode mat)
     tStoreS' :: S.SparseMatrix 80 60 Vector Double -> Bool
-    tStoreS' mat = G.flatten mat == S.withDecodedMatrix (encode mat) G.flatten
+    tStoreS' mat = G.flatten mat ==
+        Dyn.withDyn (Dyn.decodeSparse $ encode mat) G.flatten
+    tMM :: S.SparseMatrix 80 60 Vector Int -> Bool
+    tMM mat = runST (runConduit $ toMM mat .| fromMM) == mat
 
 pConversion :: TestTree
 pConversion = testGroup "Conversion"
diff --git a/tests/Test/Utils.hs b/tests/Test/Utils.hs
--- a/tests/Test/Utils.hs
+++ b/tests/Test/Utils.hs
@@ -39,8 +39,9 @@
         (r1, i1) = polar a
         (r2, i2) = polar b
 
-instance (SingI m, SingI n, Storable a, Approx a) => Approx (Matrix m n a) where
-    m1 ~= m2 = D.all id $ D.zipWith (~=) m1 m2
+instance (SingI m, SingI n, G.Vector v Bool, G.Vector v a, Approx a) =>
+    Approx (D.Matrix m n v a) where
+        m1 ~= m2 = D.all id $ D.zipWith (~=) m1 m2
 
 instance (G.Vector v a, Arbitrary a, SingI m, SingI n)
     => Arbitrary (D.Matrix m n v a) where