diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -37,7 +37,7 @@
 
 f :: (SingI n, (2 <= n - 2) ~ 'True)
   => Matrix n n Double -> Matrix n n (Complex Double)
-f m = let (d, v) = eigs (sing :: Sing 2) m
+f m = let (d, v) = eigS (sing :: Sing 2) m
       in v @@ S.diag d @@ G.transpose v
 
 main :: IO ()
diff --git a/cbits/eigen-runtime.h b/cbits/eigen-runtime.h
--- a/cbits/eigen-runtime.h
+++ b/cbits/eigen-runtime.h
@@ -39,6 +39,18 @@
 	GUARD_END\
 }
 
+#define API2(name,args,call) \
+extern "C" RET eigen_##name args {\
+	GUARD_START\
+	switch (code) {\
+		case 0: return name<T0, T0>call;\
+		case 1: return name<T1, T1>call;\
+		case 2: return name<T2, T0>call;\
+		case 3: return name<T3, T1>call;\
+	}\
+	GUARD_END\
+}
+
 typedef float T0;
 typedef double T1;
 typedef std::complex<float> T2;
diff --git a/cbits/eigen-solver.cpp b/cbits/eigen-solver.cpp
--- a/cbits/eigen-solver.cpp
+++ b/cbits/eigen-solver.cpp
@@ -1,4 +1,5 @@
 #include <Spectra/GenEigsSolver.h>
+#include <Spectra/SymEigsSolver.h>
 #include "eigen-runtime.h"
 #include <Eigen/Sparse>
 #include <Spectra/MatOp/SparseGenMatProd.h>
@@ -6,6 +7,20 @@
 using namespace Spectra;
 using namespace Eigen;
 
+extern "C" RET eigen_eig( 
+    void* d, void* v, 
+    const void* p, int n)
+{
+    typedef Map< Matrix<T1,Dynamic,Dynamic> > MapMatrix;
+    typedef Map< Matrix<T3,Dynamic,Dynamic> > MapComplexMatrix;
+    MapMatrix M((T1*)p, n, n);
+    MapComplexMatrix D((T3*)d, n, 1);
+    MapComplexMatrix V((T3*)v, n, n);
+    EigenSolver<MatrixXd> es(M);
+    D = es.eigenvalues();
+    V = es.eigenvectors();
+}
+
 extern "C" RET spectral_eigs( 
     int k,
     void* d, void* v, 
@@ -29,6 +44,28 @@
     return 0;
 }
 
+extern "C" RET spectral_eigsh( 
+    int k,
+    void* d, void* v, 
+    const void* p, int n)
+{
+    typedef Map< Matrix<T1,Dynamic,Dynamic> > MapMatrix;
+    MapMatrix M((T1*)p, n, n);
+    MapMatrix D((T1*)d, k, 1);
+    MapMatrix V((T1*)v, n, k);
+
+    DenseGenMatProd<double> op(M);
+    int ncv = 2 * k;
+    ncv = (ncv <= n) ? ncv : n;
+    SymEigsSolver< double, LARGEST_MAGN, DenseGenMatProd<double> > eigsh(&op, k, ncv);
+    eigsh.init();
+    int nconv = eigsh.compute();
+    if(eigsh.info() == 0)
+        D = eigsh.eigenvalues();
+        V = eigsh.eigenvectors();
+    return 0;
+}
+
 extern "C" RET spectral_seigs( 
     int k,
     void* d, void* v,
@@ -55,6 +92,33 @@
     return 0;
 }
 
+extern "C" RET spectral_seigsh( 
+    int k,
+    void* d, void* v,
+    const void* values,
+    const void* outerIndexPtr,
+    const void* innerIndices,
+    int n, int s)
+{
+    typedef Map< Matrix<T1,Dynamic,Dynamic> > MapMatrix;
+    typedef Map<const SparseMatrix<T1> > MapSparseMatrix;
+    MapSparseMatrix M(n, n, s, (int*)outerIndexPtr, (int*)innerIndices, (T1*)values);
+    MapMatrix D((T1*)d, k, 1);
+    MapMatrix V((T1*)v, n, k);
+
+    SparseGenMatProd<double> op(M);
+    int ncv = 2 * k;
+    ncv = (ncv <= n) ? ncv : n;
+    SymEigsSolver< double, LARGEST_MAGN, SparseGenMatProd<double> > eigsh(&op, k, ncv);
+    eigsh.init();
+    int nconv = eigsh.compute();
+    if(eigsh.info() == 0)
+        D = eigsh.eigenvalues();
+        V = eigsh.eigenvectors();
+    return 0;
+}
+
+
 template <class T>
 RET cholesky(void* px, const void* pa, int n)
 {
@@ -68,16 +132,23 @@
     void* px, const void* pa, int n), (px,pa,n));
 
 
-/*
-template <class T>
-RET bdcsvd(void* px, int r, int c
+template <class T, class TT>
+RET bdcsvd(
+    void* pu, void* ps, void* pv, 
+    const void* px, int r, int c)
 {
+    int m = r < c ? r : c;
     typedef Map< Matrix<T,Dynamic,Dynamic> > MapMatrix;
-    MapMatrix x((T*)px, n, n);
-    MapMatrix A((T*)pa, n, n);
-    x = A.llt().matrixL();
+    typedef Map< Matrix<TT,Dynamic,Dynamic> > MapMatrix2;
+    MapMatrix A((T*)px, r, c);
+    MapMatrix U((T*)pu, r, m);
+    MapMatrix2 s((TT*)ps, m, 1);
+    MapMatrix V((T*)pv, c, m);
+    BDCSVD< Matrix<T,Dynamic,Dynamic>> svd(A, ComputeThinU|ComputeThinV);
+    U = svd.matrixU();
+    V = svd.matrixV();
+    s = svd.singularValues();
     return 0;
 }
-API(cholesky, (int code,
-    void* px, const void* pa, int n), (px,pa,n));
-*/
+API2(bdcsvd, (int code,
+    void* pu, void* ps, void* pv, const void* px, int r, int c), (pu,ps,pv,px,r,c));
diff --git a/matrix-sized.cabal b/matrix-sized.cabal
--- a/matrix-sized.cabal
+++ b/matrix-sized.cabal
@@ -4,10 +4,10 @@
 --
 -- see: https://github.com/sol/hpack
 --
--- hash: ef7b97a49f99defe3cd58f2c5a3411fe8c958f0ad40f4b3f0bd763720b643db8
+-- hash: a94165a01f28a1a01026b95b46341900f92789d9a5b17aec1fe9ab414b560a26
 
 name:           matrix-sized
-version:        0.0.3
+version:        0.0.4
 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
@@ -392,6 +392,7 @@
       Data.Matrix.Static.Generic.Mutable
   other-modules:
       Data.Matrix.Static.Internal
+      Data.Matrix.Static.LinearAlgebra.Internal
   hs-source-dirs:
       src
   ghc-options: -Wall
@@ -408,6 +409,7 @@
       stdc++
   build-depends:
       base >=4.10 && <5
+    , conduit
     , primitive >=0.6.4.0
     , singletons
     , vector >=0.11
@@ -420,11 +422,18 @@
 test-suite test
   type: exitcode-stdio-1.0
   main-is: Main.hs
+  other-modules:
+      Test.Base
+      Test.LinearAlgebra
+      Test.Utils
   hs-source-dirs:
       tests
   ghc-options: -Wall -threaded
   build-depends:
       base
+    , conduit
+    , data-ordlist
+    , ieee754
     , matrix-sized
     , primitive >=0.6.4.0
     , singletons
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
@@ -42,6 +42,7 @@
     , C.fromColumns
     , C.unsafeFromVector
     , replicate
+    , diag
     , diagRect
 
     -- * Conversions
@@ -88,11 +89,17 @@
     , C.freeze
     , C.unsafeFreeze
     , C.create
+
+    , sum
+    , all
+    , any
     ) where
 
 import           Control.Monad                     (liftM)
 import qualified Data.Vector.Generic               as G
-import Prelude hiding (replicate, mapM, mapM_, zipWith, map, sequence, sequence_, zip, unzip, zipWith3, zip3, unzip3)
+import Prelude hiding ( replicate, mapM, mapM_, zipWith, map
+                      , sequence, sequence_, zip, unzip, zipWith3
+                      , zip3, unzip3, sum, all, any )
 import GHC.TypeLits (type (<=))
 import Data.Singletons
 import Data.Tuple (swap)
@@ -203,6 +210,17 @@
     c = fromIntegral $ fromSing (sing :: Sing c)
 {-# INLINE replicate #-}
 
+-- | O(m*n) Create a square matrix with default values and given diagonal
+diag :: (G.Vector v a, SingI n)
+     => a                    -- ^ default value
+     -> Matrix n 1 v a       -- ^ diagonal
+     -> Matrix n n v a
+diag z0 d = C.create $ do
+    mat <- DM.replicate z0
+    C.imapM_ (DM.unsafeWrite mat) d
+    return mat
+{-# INLINE diag #-}
+
 -- | O(m*n) Create a rectangular matrix with default values and given diagonal
 diagRect :: (G.Vector v a, SingI r, SingI c, n <= r, n <= c)
          => a                    -- ^ default value
@@ -352,8 +370,21 @@
 convert (Matrix vec) = Matrix $ G.convert vec
 {-# INLINE convert #-}
 
+sum :: (Num a, G.Vector v a) => Matrix r c v a -> a
+sum (Matrix vec) = G.sum vec
+{-# INLINE sum #-}
 
+all :: G.Vector v a => (a -> Bool) -> Matrix r c v a -> Bool
+all f (Matrix vec) = G.all f vec
+{-# INLINE all #-}
+
+any :: G.Vector v a => (a -> Bool) -> Matrix r c v a -> Bool
+any f (Matrix vec) = G.any f vec
+{-# INLINE any #-}
+
 -- Helper
-toIndex :: Int -> Int -> (Int, Int)
+toIndex :: Int   -- ^ Number of rows
+        -> Int   -- ^ 1-d index
+        -> (Int, Int)   -- ^ 2-d index
 toIndex r i = swap $ i `divMod` r
 {-# INLINE toIndex #-}
diff --git a/src/Data/Matrix/Static/Internal.hs b/src/Data/Matrix/Static/Internal.hs
--- a/src/Data/Matrix/Static/Internal.hs
+++ b/src/Data/Matrix/Static/Internal.hs
@@ -9,8 +9,12 @@
     , c_ss_plus
     , c_inverse
     , c_cholesky
+    , c_eig
     , c_eigs
+    , c_eigsh
     , c_seigs
+    , c_seigsh
+    , c_bdcsvd
     ) where
 
 import Data.Complex (Complex)
@@ -80,11 +84,27 @@
     c_cholesky :: CInt
                -> Ptr a -> Ptr a ->  CInt -> IO CString
 
+foreign import ccall "eigen_eig"
+    c_eig :: Ptr (Complex Double) -> Ptr (Complex Double)
+          -> Ptr Double -> CInt -> IO CString
+
 foreign import ccall "spectral_eigs"
     c_eigs :: CInt -> Ptr (Complex Double)
            -> Ptr (Complex Double) -> Ptr Double -> CInt -> IO CString
 
+foreign import ccall "spectral_eigsh"
+    c_eigsh :: CInt -> Ptr Double -> Ptr Double -> Ptr Double -> CInt -> IO CString
+
 foreign import ccall "spectral_seigs"
     c_seigs :: CInt -> Ptr (Complex Double) -> Ptr (Complex Double)
             -> Ptr Double -> Ptr CInt -> Ptr CInt
             -> CInt -> CInt -> IO CString
+
+foreign import ccall "spectral_seigsh"
+    c_seigsh :: CInt -> Ptr Double -> Ptr Double
+             -> Ptr Double -> Ptr CInt -> Ptr CInt
+             -> CInt -> CInt -> IO CString
+
+foreign import ccall "eigen_bdcsvd"
+    c_bdcsvd :: CInt -> Ptr a -> Ptr b -> Ptr a
+             -> Ptr a -> CInt -> CInt -> IO CString
diff --git a/src/Data/Matrix/Static/LinearAlgebra.hs b/src/Data/Matrix/Static/LinearAlgebra.hs
--- a/src/Data/Matrix/Static/LinearAlgebra.hs
+++ b/src/Data/Matrix/Static/LinearAlgebra.hs
@@ -6,15 +6,26 @@
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE FunctionalDependencies #-}
 {-# LANGUAGE GADTs #-}
+{-# LANGUAGE FlexibleContexts #-}
 {-# LANGUAGE ExplicitNamespaces #-}
+{-# LANGUAGE TypeFamilies #-}
 {-# LANGUAGE TypeOperators #-}
+{-# LANGUAGE InstanceSigs #-}
 module Data.Matrix.Static.LinearAlgebra
-    ( Arithmetic(..)
+    ( module Data.Matrix.Static.LinearAlgebra.Types
+    , Arithmetic(..)
     , Factorization(..)
-    , module Data.Matrix.Static.LinearAlgebra.Types
+    , LinearAlgebra(..)
+
+      -- * Dense matrix operation
+    , inverse
+    , eig
+    , svd
+    , cond
     ) where
 
 import qualified Data.Vector.Storable as VS
+import Data.Vector.Storable (Vector)
 import System.IO.Unsafe (unsafePerformIO)
 import Data.Complex (Complex)
 import Data.Singletons.Prelude hiding ((@@), type (==))
@@ -27,35 +38,39 @@
 import qualified Data.Matrix.Static.Generic as C
 import qualified Data.Matrix.Static.Internal as Internal
 import Data.Matrix.Static.LinearAlgebra.Types
+import Data.Matrix.Static.LinearAlgebra.Internal
 
 class Arithmetic (mat1 :: C.MatrixKind) (mat2 :: C.MatrixKind) where
-    -- | Matrix multiplication
+    -- | Matrix multiplication between different types of matrices.
     (@@) :: ( Numeric a, SingI n, SingI m
             , If (mat1 == mat2) mat1 D.Matrix ~ mat3 )
-         => mat1 n p VS.Vector a
-         -> mat2 p m VS.Vector a
-         -> mat3 n m VS.Vector a
+         => mat1 n p Vector a
+         -> mat2 p m Vector a
+         -> mat3 n m Vector a
     infixr 8 @@
 
+    -- | Element-wise addition between different types of matrices.
     (%+%) :: ( Numeric a, SingI n, SingI m
              , If (mat1 == mat2) mat1 D.Matrix ~ mat3 )
-          => mat1 n m VS.Vector a
-          -> mat2 n m VS.Vector a
-          -> mat3 n m VS.Vector a
+          => mat1 n m Vector a
+          -> mat2 n m Vector a
+          -> mat3 n m Vector a
     infixr 8 %+%
 
+    -- | Element-wise substraction between different types of matrices.
     (%-%) :: ( Numeric a, SingI n, SingI m
              , If (mat1 == mat2) mat1 D.Matrix ~ mat3 )
-          => mat1 n m VS.Vector a
-          -> mat2 n m VS.Vector a
-          -> mat3 n m VS.Vector a
+          => mat1 n m Vector a
+          -> mat2 n m Vector a
+          -> mat3 n m Vector a
     infixr 8 %-%
 
+    -- | Element-wise multiplication between different types of matrices.
     (%*%) :: ( Numeric a, SingI n, SingI m
              , If (mat1 == mat2) mat1 S.SparseMatrix ~ mat3 )
-          => mat1 n m VS.Vector a
-          -> mat2 n m VS.Vector a
-          -> mat3 n m VS.Vector a
+          => mat1 n m Vector a
+          -> mat2 n m Vector a
+          -> mat3 n m Vector a
     infixr 8 %*%
 
 instance Arithmetic D.Matrix D.Matrix where
@@ -82,45 +97,69 @@
     (%-%) a b = a %+% C.map negate b
     (%*%) = withSS Internal.c_ss_cmul
 
+class LinearAlgebra (mat :: C.MatrixKind) where
+    ident :: (Numeric a, SingI n) => mat n n Vector a
 
-class Factorization mat where
-    -- | Matrix inverse
-    inverse :: (SingI n, Numeric a) => mat n n VS.Vector a -> mat n n VS.Vector a
+instance LinearAlgebra D.Matrix where
+    ident = D.diag 0 $ D.replicate 1
 
-    -- | Eigenvalues (not ordered) and
+instance LinearAlgebra S.SparseMatrix where
+    ident = S.diag $ D.replicate 1
+
+class Factorization mat where
+    -- | Eigenvalues (from largest to smallest) and
     -- eigenvectors (as columns) of a general square matrix.
-    eigs :: (SingI k, SingI n, (k <= n - 2) ~ 'True)
+    eigS :: (SingI k, SingI n, (k <= n - 2) ~ 'True)
          => Sing k
-         -> mat n n VS.Vector Double
+         -> mat n n Vector Double
          -> (Matrix k 1 (Complex Double), Matrix n k (Complex Double))
 
+    -- | Eigenvalues (from largest to smallest) and
+    -- eigenvectors (as columns) of a symmetric square matrix.
+    eigSH :: (SingI k, SingI n, (k <= n - 1) ~ 'True)
+          => Sing k
+          -> mat n n Vector Double
+          -> (Matrix k 1 Double, Matrix n k Double)
+
     -- | Cholesky decomposition
-    cholesky :: (Numeric a, SingI n) => mat n n VS.Vector a -> mat n n VS.Vector a
+    cholesky :: (Numeric a, SingI n) => mat n n Vector a -> mat n n Vector a
 
-instance Factorization D.Matrix where
 
-    inverse = withFun1 Internal.c_inverse
+instance Factorization D.Matrix where
+    eigS s mat
+        | D.all (==0) mat = ( D.replicate 0, D.replicate 1)
+        | otherwise = unsafePerformIO $ do
+            m1 <- CM.new
+            m2 <- CM.new
+            _ <- unsafeWith' m1 $ \v1 _ _ -> unsafeWith' m2 $ \v2 _ _ ->
+                unsafeWith mat $ \v n _ -> Internal.c_eigs k v1 v2 v n
+            m1' <- C.unsafeFreeze m1
+            m2' <- C.unsafeFreeze m2
+            return (m1', m2')
+      where
+        k = fromIntegral $ fromSing s
+    {-# INLINE eigS #-}
 
-    eigs s mat = unsafePerformIO $ do
-        m1 <- CM.new
-        m2 <- CM.new
-        _ <- unsafeWith' m1 $ \v1 _ _ -> unsafeWith' m2 $ \v2 _ _ -> do
-            unsafeWith mat $ \v n _ -> Internal.c_eigs k v1 v2 v n
-        m1' <- C.unsafeFreeze m1
-        m2' <- C.unsafeFreeze m2
-        return (m1', m2')
+    eigSH s mat
+        | D.all (==0) mat = (D.replicate 0, D.replicate 1)
+        | otherwise = unsafePerformIO $ do
+            m1 <- CM.new
+            m2 <- CM.new
+            _ <- unsafeWith' m1 $ \v1 _ _ -> unsafeWith' m2 $ \v2 _ _ ->
+                unsafeWith mat $ \v n _ -> Internal.c_eigsh k v1 v2 v n
+            m1' <- C.unsafeFreeze m1
+            m2' <- C.unsafeFreeze m2
+            return (m1', m2')
       where
         k = fromIntegral $ fromSing s
-    {-# INLINE eigs #-}
+    {-# INLINE eigSH #-}
 
     cholesky mat = flip withFun1 mat $
         \code p1 c1 _ p2 _ _ -> Internal.c_cholesky code p1 p2 c1
     {-# INLINE cholesky #-}
 
 instance Factorization S.SparseMatrix where
-    inverse = undefined
-
-    eigs s mat = unsafePerformIO $ do
+    eigS s mat = unsafePerformIO $ do
         m1 <- CM.new
         m2 <- CM.new
         _ <- unsafeWith' m1 $ \v1 _ _ -> unsafeWith' m2 $ \v2 _ _ ->
@@ -131,6 +170,72 @@
         return (m1', m2')
       where
         k = fromIntegral $ fromSing s
-    {-# INLINE eigs #-}
+    {-# INLINE eigS #-}
 
+    eigSH s mat = unsafePerformIO $ do
+        m1 <- CM.new
+        m2 <- CM.new
+        _ <- unsafeWith' m1 $ \v1 _ _ -> unsafeWith' m2 $ \v2 _ _ ->
+            unsafeWithS mat $ \pv pin po n _ size ->
+                Internal.c_seigsh k v1 v2 pv po pin n size
+        m1' <- C.unsafeFreeze m1
+        m2' <- C.unsafeFreeze m2
+        return (m1', m2')
+      where
+        k = fromIntegral $ fromSing s
+    {-# INLINE eigSH #-}
+
     cholesky = undefined
+
+type family R a where
+    R Float = Float
+    R Double = Double
+    R (Complex Double) = Double
+    R (Complex Float) = Float
+    
+-- | The inverse of a dense matrix.
+inverse :: (SingI n, Numeric a) => Matrix n n a -> Matrix n n a
+inverse = withFun1 Internal.c_inverse
+{-# INLINE inverse #-}
+
+-- | Compute the full eigendecomposition for dense matrix.
+eig :: forall n . SingI n
+    => Matrix n n Double
+    -> (Matrix n 1 (Complex Double), Matrix n n (Complex Double))
+eig mat = unsafePerformIO $ do
+    m1 <- CM.new
+    m2 <- CM.new
+    _ <- unsafeWith' m1 $ \v1 _ _ -> unsafeWith' m2 $ \v2 _ _ ->
+        unsafeWith mat $ \v n _ -> Internal.c_eig v1 v2 v n
+    m1' <- C.unsafeFreeze m1
+    m2' <- C.unsafeFreeze m2
+    return (m1', m2')
+{-# INLINE eig #-}
+
+-- | Compute the full singular value decomposition for dense matrix.
+svd :: forall n p a m. (Numeric (R a), Numeric a, SingI n, SingI p, SingI m, m ~ Min n p)
+    => Matrix n p a
+    -> (Matrix n m a, Matrix m 1 (R a), Matrix p m a)
+svd mat = unsafePerformIO $ do
+    mu <- CM.new
+    ms <- CM.new
+    mv <- CM.new
+    checkResult $ unsafeWith' mu $ \pu _ _ -> unsafeWith' ms $ \ps _ _ ->
+        unsafeWith' mv $ \pv _ _ -> unsafeWith mat $ \px r c ->
+            Internal.c_bdcsvd (foreignType (undefined :: a))
+                pu ps pv px r c
+    u <- C.unsafeFreeze mu
+    s <- C.unsafeFreeze ms
+    v <- C.unsafeFreeze mv
+    return (u, s, v)
+{-# INLINE svd #-}
+
+-- | Condition number.
+cond :: ( Numeric a, Numeric (R a), Ord (R a), Fractional (R a)
+        , SingI n, SingI m, SingI (Min n m))
+     => Matrix n m a -> R a
+cond mat = VS.maximum val / VS.minimum val
+  where
+    val = VS.filter (/=0) $ D.flatten s
+    (_,s,_) = svd mat
+{-# INLINE cond #-}
diff --git a/src/Data/Matrix/Static/LinearAlgebra/Internal.hs b/src/Data/Matrix/Static/LinearAlgebra/Internal.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Matrix/Static/LinearAlgebra/Internal.hs
@@ -0,0 +1,189 @@
+{-# LANGUAGE FlexibleInstances     #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE RankNTypes #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE DataKinds #-}
+{-# OPTIONS_GHC -Wno-simplifiable-class-constraints #-}
+module Data.Matrix.Static.LinearAlgebra.Internal
+    ( withFun1
+    , withFun2
+    , withDS
+    , withSD
+    , withSS
+    , checkResult
+    , unsafeWith
+    , unsafeWith'
+    , unsafeWithS
+    ) where
+
+import Data.Vector.Storable (Storable)
+import qualified Data.Vector.Storable as VS
+import qualified Data.Vector.Storable.Mutable as VSM
+import System.IO.Unsafe (unsafePerformIO)
+import Control.Monad (when)
+import Control.Monad.ST (RealWorld)
+import Data.Singletons
+import Foreign
+import Foreign.C.Types
+import Foreign.C.String
+
+import qualified Data.Matrix.Static.Dense as D
+import qualified Data.Matrix.Static.Dense.Mutable as DM
+import qualified Data.Matrix.Static.Sparse as S
+import qualified Data.Matrix.Static.Generic.Mutable as CM
+import qualified Data.Matrix.Static.Generic as C
+import Data.Matrix.Static.LinearAlgebra.Types
+
+withFun1 :: forall r1 c1 r2 c2 a. (SingI r2, SingI c2, Numeric a)
+         => (CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString)
+         -> Matrix r1 c1 a -> Matrix r2 c2 a
+withFun1 f m1 = unsafePerformIO $ do
+    m0 <- CM.new
+    checkResult $ unsafeWith' m0 $ \vals0 rows0 cols0 ->
+        unsafeWith m1 $ \vals1 rows1 cols1 -> f (foreignType (undefined :: a))
+            vals0 rows0 cols0
+            vals1 rows1 cols1
+    C.unsafeFreeze m0
+{-# INLINE withFun1 #-}
+
+withFun2 :: forall r1 c1 r2 c2 r3 c3 a.
+            (SingI r3, SingI c3, Numeric a)
+         => ( CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt
+           -> Ptr a -> CInt -> CInt -> IO CString )
+         -> Matrix r1 c1 a
+         -> Matrix r2 c2 a
+         -> Matrix r3 c3 a
+withFun2 f m1 m2 = unsafePerformIO $ do
+    m0 <- CM.new
+    checkResult $ unsafeWith' m0 $ \vals0 rows0 cols0 ->
+        unsafeWith m1 $ \vals1 rows1 cols1 ->
+            unsafeWith m2 $ \vals2 rows2 cols2 ->
+                f (foreignType (undefined :: a))
+                    vals0 rows0 cols0
+                    vals1 rows1 cols1
+                    vals2 rows2 cols2
+    C.unsafeFreeze m0
+{-# INLINE withFun2 #-}
+
+withDS :: forall r1 c1 r2 c2 r3 c3 a.
+            (SingI r3, SingI c3, Numeric a)
+       => ( CInt
+         -> Ptr a -> CInt -> CInt
+         -> Ptr a -> CInt -> CInt
+         -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt
+         -> IO CString )
+       -> Matrix r1 c1 a
+       -> SparseMatrix r2 c2 a
+       -> Matrix r3 c3 a
+withDS f m1 m2 = unsafePerformIO $ do
+    m0 <- CM.new
+    checkResult $ unsafeWith' m0 $ \v0 r0 c0 ->
+        unsafeWith m1 $ \v1 r1 c1 ->
+            unsafeWithS m2 $ \v2 inner outer r2 c2 s ->
+                f (foreignType (undefined :: a))
+                    v0 r0 c0
+                    v1 r1 c1
+                    v2 outer inner r2 c2 s
+    C.unsafeFreeze m0
+{-# INLINE withDS #-}
+
+withSD :: forall r1 c1 r2 c2 r3 c3 a.
+            (SingI r3, SingI c3, Numeric a)
+       => ( CInt
+         -> Ptr a -> CInt -> CInt
+         -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt
+         -> Ptr a -> CInt -> CInt
+         -> IO CString )
+       -> SparseMatrix r2 c2 a
+       -> Matrix r1 c1 a
+       -> Matrix r3 c3 a
+withSD f m2 m1 = unsafePerformIO $ do
+    m0 <- CM.new
+    checkResult $ unsafeWith' m0 $ \v0 r0 c0 ->
+        unsafeWith m1 $ \v1 r1 c1 ->
+            unsafeWithS m2 $ \v2 inner outer r2 c2 s ->
+                f (foreignType (undefined :: a))
+                    v0 r0 c0
+                    v2 outer inner r2 c2 s
+                    v1 r1 c1
+    C.unsafeFreeze m0
+{-# INLINE withSD #-}
+
+mkSparseMatrix :: forall r c a. (Storable a, SingI r, SingI c)
+    => (Ptr (Ptr a) -> Ptr (Ptr CInt) -> Ptr CInt -> IO Int)
+    -> IO (SparseMatrix r c a)
+mkSparseMatrix f = do
+    outer' <- VSM.new $ c + 1
+    (n, pv, pinner) <- VSM.unsafeWith outer' $ \pouter -> alloca $ \ppv -> alloca $ \ppi -> do
+        n <- f ppv ppi pouter
+        pv <- peek ppv >>= newForeignPtr finalizerFree
+        pinner <- peek ppi >>= newForeignPtr finalizerFree
+        return (n, pv, pinner)
+    outer <- VS.unsafeFreeze outer'
+    return $ S.SparseMatrix (VS.unsafeFromForeignPtr0 pv n)
+        (VS.unsafeFromForeignPtr0 pinner n)
+        outer
+  where
+    c = fromIntegral $ fromSing (sing :: Sing c)
+{-# INLINE mkSparseMatrix #-}
+
+withSS :: forall r1 c1 r2 c2 r3 c3 a.
+            (SingI r3, SingI c3, Numeric a)
+       => ( CInt
+         -> Ptr (Ptr a) -> Ptr CInt -> Ptr (Ptr CInt) -> CInt -> CInt -> Ptr CInt
+         -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt
+         -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt
+         -> IO CString )
+       -> SparseMatrix r1 c1 a
+       -> SparseMatrix r2 c2 a
+       -> SparseMatrix r3 c3 a
+withSS f m1 m2 = unsafePerformIO $ mkSparseMatrix $ \v0 inner0 outer0 ->
+    alloca $ \pn -> unsafeWithS m1 $ \v1 inner1 outer1 r1 c1 s1 ->
+        unsafeWithS m2 $ \v2 inner2 outer2 r2 c2 s2 -> do
+            checkResult $ f (foreignType (undefined :: a))
+                v0 outer0 inner0 r c pn
+                v1 outer1 inner1 r1 c1 s1
+                v2 outer2 inner2 r2 c2 s2
+            fromIntegral <$> peek pn
+  where
+    r = fromIntegral $ fromSing (sing :: Sing r3)
+    c = fromIntegral $ fromSing (sing :: Sing c3)
+{-# INLINE withSS #-}
+
+checkResult :: IO CString -> IO ()
+checkResult func = func >>= \c_str -> when (c_str /= nullPtr) $
+    peekCString c_str >>= \str -> error str
+{-# INLINE checkResult #-}
+
+-------------------------------------------------------------------------------
+-- Raw pointers
+-------------------------------------------------------------------------------
+
+-- | Pass a pointer to the matrix's data to the IO action.
+-- The data may not be modified through the pointer.
+unsafeWith :: Storable a => Matrix n m a -> (Ptr a -> CInt -> CInt -> IO b) -> IO b
+unsafeWith mat@(D.Matrix vec) f = VS.unsafeWith vec $ \p ->
+    f p (fromIntegral r) $ fromIntegral c 
+  where
+    (r,c) = C.dim mat
+{-# INLINE unsafeWith #-}
+
+unsafeWith' :: Storable a => MMatrix n m RealWorld a -> (Ptr a -> CInt -> CInt -> IO b) -> IO b
+unsafeWith' mat@(DM.MMatrix vec) f = VSM.unsafeWith vec $ \p ->
+    f p (fromIntegral r) $ fromIntegral c
+  where
+    (r,c) = CM.dim mat
+{-# INLINE unsafeWith' #-}
+
+-- | Pass a pointer to the matrix's data to the IO action.
+-- The data may not be modified through the pointer.
+unsafeWithS :: (Storable a, S.Zero a)
+            => SparseMatrix n m a
+            -> (Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt -> IO b)
+            -> IO b
+unsafeWithS mat@(S.SparseMatrix val inner outer) f = VS.unsafeWith val $ \pval ->
+    VS.unsafeWith inner $ \pinner -> VS.unsafeWith outer $ \pouter ->
+        f pval pinner pouter (fromIntegral r) (fromIntegral c) (fromIntegral $ VS.length val)
+  where
+    (r,c) = C.dim mat
+{-# INLINE unsafeWithS #-}
diff --git a/src/Data/Matrix/Static/LinearAlgebra/Types.hs b/src/Data/Matrix/Static/LinearAlgebra/Types.hs
--- a/src/Data/Matrix/Static/LinearAlgebra/Types.hs
+++ b/src/Data/Matrix/Static/LinearAlgebra/Types.hs
@@ -9,39 +9,22 @@
     , Matrix
     , MMatrix
     , SparseMatrix
-    , withFun1
-    , withFun2
-    , withDS
-    , withSD
-    , withSS
-    , unsafeWith
-    , unsafeWith'
-    , unsafeWithS
     ) where
 
 import Data.Vector.Storable (Vector, Storable)
-import qualified Data.Vector.Storable as VS
-import qualified Data.Vector.Storable.Mutable as VSM
 import Data.Vector.Storable.Mutable (MVector)
-import System.IO.Unsafe (unsafePerformIO)
-import Control.Monad (when)
 import Data.Complex (Complex)
-import Control.Monad.ST (RealWorld)
-import Data.Singletons
-import Foreign
 import Foreign.C.Types
-import Foreign.C.String
 
 import qualified Data.Matrix.Static.Dense as D
 import qualified Data.Matrix.Static.Dense.Mutable as DM
 import qualified Data.Matrix.Static.Sparse as S
-import qualified Data.Matrix.Static.Generic.Mutable as CM
-import qualified Data.Matrix.Static.Generic as C
 
 class (S.Zero a, Storable a, Num a) => Numeric a where
     foreignType :: a -> CInt
 
 instance Numeric Float where foreignType _ = 0
+instance Numeric CFloat where foreignType _ = 0
 instance Numeric Double where foreignType _ =1
 instance Numeric (Complex Float) where foreignType _ = 2
 instance Numeric (Complex Double) where foreignType _ = 3
@@ -50,157 +33,3 @@
 type MMatrix r c s a = DM.MMatrix r c MVector s a
 
 type SparseMatrix r c a = S.SparseMatrix r c Vector a
-
-withFun1 :: forall r1 c1 r2 c2 a. (SingI r2, SingI c2, Numeric a)
-         => (CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt -> IO CString)
-         -> Matrix r1 c1 a -> Matrix r2 c2 a
-withFun1 f m1 = unsafePerformIO $ do
-    m0 <- CM.new
-    checkResult $ unsafeWith' m0 $ \vals0 rows0 cols0 ->
-        unsafeWith m1 $ \vals1 rows1 cols1 -> f (foreignType (undefined :: a))
-            vals0 rows0 cols0
-            vals1 rows1 cols1
-    C.unsafeFreeze m0
-{-# INLINE withFun1 #-}
-
-withFun2 :: forall r1 c1 r2 c2 r3 c3 a.
-            (SingI r3, SingI c3, Numeric a)
-         => ( CInt -> Ptr a -> CInt -> CInt -> Ptr a -> CInt -> CInt
-           -> Ptr a -> CInt -> CInt -> IO CString )
-         -> Matrix r1 c1 a
-         -> Matrix r2 c2 a
-         -> Matrix r3 c3 a
-withFun2 f m1 m2 = unsafePerformIO $ do
-    m0 <- CM.new
-    checkResult $ unsafeWith' m0 $ \vals0 rows0 cols0 ->
-        unsafeWith m1 $ \vals1 rows1 cols1 ->
-            unsafeWith m2 $ \vals2 rows2 cols2 ->
-                f (foreignType (undefined :: a))
-                    vals0 rows0 cols0
-                    vals1 rows1 cols1
-                    vals2 rows2 cols2
-    C.unsafeFreeze m0
-{-# INLINE withFun2 #-}
-
-withDS :: forall r1 c1 r2 c2 r3 c3 a.
-            (SingI r3, SingI c3, Numeric a)
-       => ( CInt
-         -> Ptr a -> CInt -> CInt
-         -> Ptr a -> CInt -> CInt
-         -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt
-         -> IO CString )
-       -> Matrix r1 c1 a
-       -> SparseMatrix r2 c2 a
-       -> Matrix r3 c3 a
-withDS f m1 m2 = unsafePerformIO $ do
-    m0 <- CM.new
-    checkResult $ unsafeWith' m0 $ \v0 r0 c0 ->
-        unsafeWith m1 $ \v1 r1 c1 ->
-            unsafeWithS m2 $ \v2 inner outer r2 c2 s ->
-                f (foreignType (undefined :: a))
-                    v0 r0 c0
-                    v1 r1 c1
-                    v2 outer inner r2 c2 s
-    C.unsafeFreeze m0
-{-# INLINE withDS #-}
-
-withSD :: forall r1 c1 r2 c2 r3 c3 a.
-            (SingI r3, SingI c3, Numeric a)
-       => ( CInt
-         -> Ptr a -> CInt -> CInt
-         -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt
-         -> Ptr a -> CInt -> CInt
-         -> IO CString )
-       -> SparseMatrix r2 c2 a
-       -> Matrix r1 c1 a
-       -> Matrix r3 c3 a
-withSD f m2 m1 = unsafePerformIO $ do
-    m0 <- CM.new
-    checkResult $ unsafeWith' m0 $ \v0 r0 c0 ->
-        unsafeWith m1 $ \v1 r1 c1 ->
-            unsafeWithS m2 $ \v2 inner outer r2 c2 s ->
-                f (foreignType (undefined :: a))
-                    v0 r0 c0
-                    v2 outer inner r2 c2 s
-                    v1 r1 c1
-    C.unsafeFreeze m0
-{-# INLINE withSD #-}
-
-mkSparseMatrix :: forall r c a. (Storable a, SingI r, SingI c)
-    => (Ptr (Ptr a) -> Ptr (Ptr CInt) -> Ptr CInt -> IO Int)
-    -> IO (SparseMatrix r c a)
-mkSparseMatrix f = do
-    outer' <- VSM.new $ c + 1
-    (n, pv, pinner) <- VSM.unsafeWith outer' $ \pouter -> alloca $ \ppv -> alloca $ \ppi -> do
-        n <- f ppv ppi pouter
-        pv <- peek ppv >>= newForeignPtr finalizerFree
-        pinner <- peek ppi >>= newForeignPtr finalizerFree
-        return (n, pv, pinner)
-    outer <- VS.unsafeFreeze outer'
-    return $ S.SparseMatrix (VS.unsafeFromForeignPtr0 pv n)
-        (VS.unsafeFromForeignPtr0 pinner n)
-        outer
-  where
-    c = fromIntegral $ fromSing (sing :: Sing c)
-{-# INLINE mkSparseMatrix #-}
-
-withSS :: forall r1 c1 r2 c2 r3 c3 a.
-            (SingI r3, SingI c3, Numeric a)
-       => ( CInt
-         -> Ptr (Ptr a) -> Ptr CInt -> Ptr (Ptr CInt) -> CInt -> CInt -> Ptr CInt
-         -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt
-         -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt
-         -> IO CString )
-       -> SparseMatrix r1 c1 a
-       -> SparseMatrix r2 c2 a
-       -> SparseMatrix r3 c3 a
-withSS f m1 m2 = unsafePerformIO $ mkSparseMatrix $ \v0 inner0 outer0 ->
-    alloca $ \pn -> unsafeWithS m1 $ \v1 inner1 outer1 r1 c1 s1 ->
-        unsafeWithS m2 $ \v2 inner2 outer2 r2 c2 s2 -> do
-            checkResult $ f (foreignType (undefined :: a))
-                v0 outer0 inner0 r c pn
-                v1 outer1 inner1 r1 c1 s1
-                v2 outer2 inner2 r2 c2 s2
-            fromIntegral <$> peek pn
-  where
-    r = fromIntegral $ fromSing (sing :: Sing r3)
-    c = fromIntegral $ fromSing (sing :: Sing c3)
-{-# INLINE withSS #-}
-
-checkResult :: IO CString -> IO ()
-checkResult func = func >>= \c_str -> when (c_str /= nullPtr) $
-    peekCString c_str >>= \str -> error str
-{-# INLINE checkResult #-}
-
--------------------------------------------------------------------------------
--- Raw pointers
--------------------------------------------------------------------------------
-
--- | Pass a pointer to the matrix's data to the IO action.
--- The data may not be modified through the pointer.
-unsafeWith :: Storable a => Matrix n m a -> (Ptr a -> CInt -> CInt -> IO b) -> IO b
-unsafeWith mat@(D.Matrix vec) f = VS.unsafeWith vec $ \p ->
-    f p (fromIntegral r) $ fromIntegral c 
-  where
-    (r,c) = C.dim mat
-{-# INLINE unsafeWith #-}
-
-unsafeWith' :: Storable a => MMatrix n m RealWorld a -> (Ptr a -> CInt -> CInt -> IO b) -> IO b
-unsafeWith' mat@(DM.MMatrix vec) f = VSM.unsafeWith vec $ \p ->
-    f p (fromIntegral r) $ fromIntegral c
-  where
-    (r,c) = CM.dim mat
-{-# INLINE unsafeWith' #-}
-
--- | Pass a pointer to the matrix's data to the IO action.
--- The data may not be modified through the pointer.
-unsafeWithS :: (Storable a, S.Zero a)
-            => SparseMatrix n m a
-            -> (Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt -> IO b)
-            -> IO b
-unsafeWithS mat@(S.SparseMatrix val inner outer) f = VS.unsafeWith val $ \pval ->
-    VS.unsafeWith inner $ \pinner -> VS.unsafeWith outer $ \pouter ->
-        f pval pinner pouter (fromIntegral r) (fromIntegral c) (fromIntegral $ VS.length val)
-  where
-    (r,c) = C.dim mat
-{-# INLINE unsafeWithS #-}
diff --git a/src/Data/Matrix/Static/Sparse.hs b/src/Data/Matrix/Static/Sparse.hs
--- a/src/Data/Matrix/Static/Sparse.hs
+++ b/src/Data/Matrix/Static/Sparse.hs
@@ -35,6 +35,8 @@
     -- * Construction
     , C.empty
     , fromTriplet
+    , fromTripletC
+    , toTriplet
     , C.fromVector
     , C.fromList
     , C.unsafeFromVector
@@ -54,10 +56,12 @@
 import qualified Data.Vector.Storable as S
 import qualified Data.Vector.Storable.Mutable as SM
 import Data.Singletons
-import Control.Monad
 import Control.Monad.ST (runST)
 import           Data.Bits                         (shiftR)
 import Text.Printf (printf)
+import Conduit
+import Data.Conduit.Internal (zipSinks)
+import Data.Tuple (swap)
 import GHC.TypeLits (type (<=))
 import Foreign.C.Types
 import Data.Complex
@@ -77,6 +81,9 @@
 instance Zero Float where
     zero = 0.0
 
+instance Zero CFloat where
+    zero = 0.0
+
 instance Zero Double where
     zero = 0.0
 
@@ -101,28 +108,16 @@
                                      -- non-zero in the previous two arrays.
                  -> SparseMatrix r c v a
 
+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'
+
 instance (G.Vector v a, Zero a, Show a) => Show (SparseMatrix r c v a) where
     show mat = printf "(%d x %d)\n%s" r c vals
       where
         (r,c) = C.dim mat
         vals = unlines $ map (unwords . map show . G.toList) $ C.toRows mat
 
-instance (SingI r, SingI c, G.Vector v a, Zero a, Num a) =>
-    Num (SparseMatrix r c v a) where
-        m1 + m2 = undefined
-        m1 - m2 = undefined
-        m1 * m2 = undefined
-        negate = C.map negate
-        abs = C.map abs
-        signum = undefined
-        fromInteger = undefined
-
-instance (SingI r, SingI c, G.Vector v a, Zero a, Fractional a) =>
-    Fractional (SparseMatrix r c v a) where
-        m1 / m2 = undefined
-        recip = C.map recip
-        fromRational = undefined
-
 instance (G.Vector v a, Zero a) => C.Matrix SparseMatrix v a where
     -- | O(1) Return the size of matrix.
     dim :: forall r c. SparseMatrix r c v a -> (Int, Int)
@@ -145,27 +140,18 @@
     -- | O(1) Create matrix from vector containing columns.
     unsafeFromVector :: forall r c. (G.Vector v a, SingI r, SingI c)
            => v a -> SparseMatrix r c v a
-    unsafeFromVector vec = SparseMatrix
-        (G.generate n (G.unsafeIndex vec . S.unsafeIndex nz))
-        inner outer
+    unsafeFromVector vec = fromTriplet vec'
       where
-        inner = S.map fromIntegral $ S.map (`mod` c) nz
-        outer = S.create $ do
-            v <- SM.replicate (c+1) 0
-            S.forM_ nz $ \x -> do
-                let i = x `div` r
-                SM.unsafeModify v succ (i+1)
-            forM_ [1..c] $ \i -> do
-                x <- SM.unsafeRead v (i-1)
-                SM.unsafeModify v (+x) i
-            return v
-        nz = S.filter (\i -> vec `G.unsafeIndex` i /= zero) $ S.enumFromN 0 (r*c)
-        n = S.length nz
+        vec' = map (\((a,b),c) -> (a,b,c)) $ filter ((/=zero) . snd) $
+            zipWith (\i x -> (toIndex i, x)) [0..] $ G.toList vec
+        toIndex i = swap $ i `divMod` r
         r = fromIntegral $ fromSing (sing :: Sing r)
-        c = fromIntegral $ fromSing (sing :: Sing c)
     {-# INLINE unsafeFromVector #-}
 
-    transpose (SparseMatrix val inner outer) = undefined 
+    transpose mat = runIdentity $ fromTripletC source
+      where
+        source = toTriplet mat .| mapC (\(i,j,x) -> (j,i,x))
+    {-# INLINE transpose #-}
 
     thaw = undefined
     {-# INLINE thaw #-}
@@ -180,9 +166,15 @@
     {-# INLINE unsafeFreeze #-}
 
     map f (SparseMatrix vec inner outer) = SparseMatrix (G.map f vec) inner outer
-    imap = undefined
     {-# INLINE map #-}
 
+    imap f mat@(SparseMatrix _ inner outer) = SparseMatrix vec' inner outer
+      where
+        vec' = runST $ runConduit $ toTriplet mat .| mapC g .| sinkVector
+        g (i,j,x) = f (i,j) x
+    {-# INLINE imap #-}
+    
+
 -- | O(n) Create matrix from triplet. row and column indices *are not* assumed to be ordered
 -- duplicate entries are carried over to the CSR represention
 fromTriplet :: forall t r c v a. (Traversable t, G.Vector v a, SingI r, SingI c)
@@ -192,28 +184,71 @@
     outer = S.scanl (+) 0 $ S.create $ do
         vec <- SM.replicate c 0
         _ <- flip mapM triplets $ \(_, j, _) -> 
-            SM.modify vec (+1) j
+            SM.unsafeModify vec (+1) j
         return vec
     (val, inner) = runST $ do
         outer' <- S.thaw outer
         val' <- GM.new nnz
         inner' <- SM.new nnz
         _ <- flip mapM triplets $ \(i, j, v) -> do
-            idx <- fromIntegral <$> SM.read outer' j
-            GM.write val' idx v
-            SM.write inner' idx $ fromIntegral i
-            SM.modify outer' (+1) j
+            idx <- fromIntegral <$> SM.unsafeRead outer' j
+            GM.unsafeWrite val' idx v
+            SM.unsafeWrite inner' idx $ fromIntegral i
+            SM.unsafeModify outer' (+1) j
         (,) <$> G.unsafeFreeze val' <*> S.unsafeFreeze inner'
     nnz = length triplets
     c = fromIntegral $ fromSing (sing :: Sing c)
 {-# INLINE fromTriplet #-}
 
-{-
-toTriplet :: (G.Vector v1 (Int, Int, a), G.Vector v2 a, SingI r, SingI c)
-          => SparseMatrix r c v2 a -> v1 (Int, Int, a)
-toTriplet mat = 
--}
+-- | O(n) Create matrix from triplet. row and column indices *are not* assumed to be ordered
+-- duplicate entries are carried over to the CSR represention
+fromTripletC :: forall m r c v a. (Monad m, G.Vector v a, SingI r, SingI c)
+             => ConduitT () (Int, Int, a) m ()
+             -> m (SparseMatrix r c v a)
+fromTripletC triplets = do
+    (nnz, outer) <- runConduit $ triplets .| zipSinks lengthC sinkOuter
+    (val, inner, _) <- runConduit $ triplets .| sinkValInner nnz (clone outer)
+    return $ SparseMatrix val inner outer
+  where
+    sinkOuter = S.scanl (+) 0 <$> foldlC f (S.replicate c 0)
+      where
+        f vec (_, j, _) = S.modify (\v -> SM.unsafeModify v (+1) j) vec
+    sinkValInner nnz outer0 = foldlC f (val0, inner0, outer0)
+      where
+        val0 = G.create $ GM.new nnz
+        inner0 = S.create $ SM.new nnz
+        f (val, inner, outer) (i, j, v) = (val', inner', outer')
+          where
+            idx = fromIntegral $ outer `S.unsafeIndex` j
+            val' = G.create $ do
+                vec <- G.unsafeThaw val
+                GM.unsafeWrite vec idx v
+                return vec
+            inner' = S.create $ do
+                vec <- S.unsafeThaw inner
+                SM.unsafeWrite vec idx $ fromIntegral i
+                return vec
+            outer' = S.create $ do
+                vec <- S.unsafeThaw outer
+                SM.unsafeModify vec (+1) j
+                return vec
+    c = fromIntegral $ fromSing (sing :: Sing c)
+    clone x = S.create $ S.thaw x
+{-# INLINE fromTripletC #-}
 
+toTriplet :: (Monad m, G.Vector v a, SingI r, SingI c)
+          => SparseMatrix r c v a -> ConduitT i (Int, Int, a) m ()
+toTriplet (SparseMatrix val inner outer) =
+    G.ifoldM_ go (fromIntegral $ G.head outer) outer
+  where
+    go start curC end = do
+        enumFromToC start (end'-1) .| mapC f
+        return end'
+      where
+        end' = fromIntegral end
+        f i = (fromIntegral $ inner `G.unsafeIndex` i, fromIntegral curC - 1, val `G.unsafeIndex` i)
+{-# INLINE toTriplet #-}
+
 -- | O(m*n) Create a rectangular matrix with default values and given diagonal
 diag :: (G.Vector v a, Zero a, SingI n)
      => D.Matrix n 1 v a       -- ^ diagonal
@@ -242,3 +277,12 @@
         k = (u+l) `shiftR` 1
         x' = vec `S.unsafeIndex` k
 {-# INLINE binarySearchByBounds #-}
+
+
+-------------------------------------------------------------------------------
+-- Helper
+-------------------------------------------------------------------------------
+
+--getIndex :: Int -> (Int, Int)
+--getIndex = 
+--{-# INLINE getIndex #-}
diff --git a/src/Data/Matrix/Static/Sparse/Mutable.hs b/src/Data/Matrix/Static/Sparse/Mutable.hs
--- a/src/Data/Matrix/Static/Sparse/Mutable.hs
+++ b/src/Data/Matrix/Static/Sparse/Mutable.hs
@@ -16,6 +16,7 @@
 import qualified Data.Vector.Generic.Mutable as GM
 import qualified Data.Vector.Storable as S
 import           Prelude                     hiding (read, replicate)
+import Foreign.C.Types
 import Data.Singletons
 
 import qualified Data.Matrix.Static.Generic.Mutable as C
@@ -25,16 +26,18 @@
     MSparseMatrix :: (SingI r, SingI c)
                   => (v s a)         -- ^ Values: stores the coefficient values
                                       -- of the non-zeros.
-                  -> (S.Vector Int)  -- ^ InnerIndices: stores the row
+                  -> (S.Vector CInt)  -- ^ InnerIndices: stores the row
                                       -- (resp. column) indices of the non-zeros.
-                  -> (S.Vector Int)  -- ^ OuterStarts: stores for each column
+                  -> (S.Vector CInt)  -- ^ OuterStarts: stores for each column
                                       -- (resp. row) the index of the first
                                       -- non-zero in the previous two arrays.
+                  -> (S.Vector CInt)  -- InnerNNZs: stores the number of non-zeros
+                                      -- of each column (resp. row). 
                   -> MSparseMatrix r c v s a
 
 instance GM.MVector v a => C.MMatrix MSparseMatrix v a where
     dim :: forall r c s. MSparseMatrix r c v s a -> (Int, Int)
-    dim (MSparseMatrix _ _ _) = (r,c)
+    dim (MSparseMatrix _ _ _ _) = (r,c)
       where
         r = fromIntegral $ fromSing (sing :: Sing r)
         c = fromIntegral $ fromSing (sing :: Sing c)
diff --git a/tests/Main.hs b/tests/Main.hs
--- a/tests/Main.hs
+++ b/tests/Main.hs
@@ -1,36 +1,13 @@
-{-# LANGUAGE DataKinds #-}
-{-# LANGUAGE ScopedTypeVariables #-}
-{-# LANGUAGE TypeSynonymInstances #-}
-{-# LANGUAGE FlexibleInstances #-}
-
 module Main where
 
 import Test.Tasty
-import Data.Matrix.Static.LinearAlgebra
-import qualified Data.Matrix.Static.Dense as D
-import qualified Data.Matrix.Static.Sparse as S
-import Data.Singletons
-import Data.Vector.Storable (Storable)
-import GHC.TypeNats (KnownNat)
-import Control.Monad.IO.Class (liftIO)
 
-import Test.Tasty.QuickCheck
-
-instance (Arbitrary a, KnownNat m, KnownNat n) => Arbitrary (Matrix m n a) where
-    arbitrary = D.fromList <$> vector (m*n)
-      where
-        m = fromIntegral $ fromSing (sing :: Sing m)
-        n = fromIntegral $ fromSing (sing :: Sing n)
-    shrink _v = []
-
-{-
-propTranspose :: Matrix 50 100 Double -> Bool
-propTranspose m = D.transpose (D.transpose m) == m && 
-    D.convertAny (S.transpose $ S.transpose (D.convertAny m)) == m
--}
+import Test.Base
+import Test.LinearAlgebra
 
 main :: IO ()
 main = defaultMain $ testGroup "Main"
-    [ testProperty "" square
+    [ base
+    , linearAlgebra
     ]
 
diff --git a/tests/Test/Base.hs b/tests/Test/Base.hs
new file mode 100644
--- /dev/null
+++ b/tests/Test/Base.hs
@@ -0,0 +1,61 @@
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE IncoherentInstances #-}
+
+module Test.Base (base) where
+
+import Test.Tasty
+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 Data.Singletons hiding ((@@))
+import Data.Singletons.Prelude (Min)
+import Data.Vector (Vector)
+import Control.Monad.IO.Class (liftIO)
+import Test.Tasty.QuickCheck
+import Conduit
+
+import Test.Utils
+
+base :: TestTree
+base = testGroup "Base"
+    [ pConversion
+    , pTranspose
+    ]
+
+pConversion = testGroup "Conversion"
+    [ testProperty "Dense: fromVector . flatten" t1
+    , testProperty "Sparse: fromVector . flatten" t2
+    , testProperty "Sparse -- fromTriplet . toTriplet" tTri
+    , testProperty "Sparse -- fromTripletC . toTriplet" tTriC
+    , testProperty "Sparse -- dense" t3 ]
+  where
+    t1 :: D.Matrix 80 60 Vector Int -> Bool
+    t1 mat = (D.fromVector $ D.flatten mat) == mat
+    t2 :: S.SparseMatrix 80 60 Vector Int -> Bool
+    t2 mat = mat == S.fromVector (S.flatten mat)
+    t3 :: D.Matrix 80 60 Vector Int -> Bool
+    t3 mat = (D.fromVector $ S.flatten mat') == mat
+      where
+        mat' = S.fromVector $ D.flatten mat :: S.SparseMatrix 80 60 Vector Int
+    tTri :: S.SparseMatrix 80 60 Vector Int -> Bool
+    tTri mat = S.fromTriplet xs == mat
+      where
+        xs = runIdentity $ runConduit $ S.toTriplet mat .| sinkList
+    tTriC :: S.SparseMatrix 80 60 Vector Int -> Bool
+    tTriC mat = mat == runIdentity (S.fromTripletC (S.toTriplet mat))
+
+pTranspose = testGroup "Transpose"
+    [ testProperty "Dense" t1
+    , testProperty "Sparse" tSp
+    ]
+  where
+    t1 :: D.Matrix 80 40 Vector Int -> Bool
+    t1 mat = D.transpose (D.transpose mat) == mat
+    tSp :: S.SparseMatrix 80 40 Vector Int -> Bool
+    tSp mat = G.transpose (G.transpose mat) == mat
diff --git a/tests/Test/LinearAlgebra.hs b/tests/Test/LinearAlgebra.hs
new file mode 100644
--- /dev/null
+++ b/tests/Test/LinearAlgebra.hs
@@ -0,0 +1,83 @@
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE IncoherentInstances #-}
+
+module Test.LinearAlgebra (linearAlgebra) where
+
+import Test.Tasty
+import Data.Matrix.Static.LinearAlgebra
+import qualified Data.Matrix.Static.Dense as D
+import qualified Data.Matrix.Static.Generic as G
+import qualified Data.Matrix.Static.Sparse as S
+import Data.Singletons hiding ((@@))
+import Data.Singletons.Prelude (Min)
+import Data.Complex
+import Data.Vector.Storable (Storable)
+import GHC.TypeNats (KnownNat)
+import Control.Monad.IO.Class (liftIO)
+import Test.Tasty.QuickCheck
+
+import Test.Utils
+
+linearAlgebra :: TestTree
+linearAlgebra = testGroup "Linear algebra"
+    [ svdTest
+    , eigenTest
+    ]
+
+svdTest = testGroup "SVD"
+    [ testProperty "SVD (Float)" svd1
+    , testProperty "SVD (Double)" svd2
+    ]
+  where
+    svd1 :: Matrix 50 30 Float -> Bool
+    svd1 m = m' ~= m
+      where
+        m' = u @@ S.diag d @@ D.transpose v
+        (u,d,v) = svd m
+    svd2 :: Matrix 50 30 Double -> Bool
+    svd2 m = m' ~= m
+      where
+        m' = u @@ S.diag d @@ D.transpose v
+        (u,d,v) = svd m
+
+eigenTest = testGroup "Eigendecomposition"
+    [ testProperty "Full" eigen1
+    , testProperty "Partial dense" eigen2
+    , testProperty "Partial symmetric dense" eigen3
+    , testProperty "Partial symmetric sparse" eigen4
+    ]
+  where
+    eigen1 :: Matrix 100 100 Double -> Bool
+    eigen1 m = (m' @@ v) ~= (v @@ S.diag d)
+      where
+        m' = D.map (\x -> mkPolar x 0) m
+        (d, v)= eig m
+    eigen2 :: Matrix 10 10 Double -> Bool
+    eigen2 m = m' @@ v ~= v @@ S.diag d
+      where
+        m' = D.map (\x -> mkPolar x 0) m
+        (d, v)= eigS (sing :: Sing 8) m
+    eigen3 :: Matrix 100 100 Double -> Bool
+    eigen3 raw = m @@ v ~= v @@ S.diag d
+      where
+        m = raw %+% D.transpose raw
+        (d, v)= eigSH (sing :: Sing 99) m
+    eigen4 :: SparseMatrix 100 100 Double -> Bool
+    eigen4 raw = m @@ v ~= v @@ S.diag d
+      where
+        m = raw %+% D.transpose raw
+        (d, v)= eigSH (sing :: Sing 99) m
+
+
+
+{-
+propTranspose :: Matrix 50 100 Double -> Bool
+propTranspose m = D.transpose (D.transpose m) == m && 
+    D.convertAny (S.transpose $ S.transpose (D.convertAny m)) == m
+-}
diff --git a/tests/Test/Utils.hs b/tests/Test/Utils.hs
new file mode 100644
--- /dev/null
+++ b/tests/Test/Utils.hs
@@ -0,0 +1,65 @@
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE TypeFamilies #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE TypeSynonymInstances #-}
+{-# LANGUAGE FlexibleContexts #-}
+{-# LANGUAGE FlexibleInstances #-}
+{-# LANGUAGE IncoherentInstances #-}
+
+module Test.Utils where
+
+import Test.Tasty.QuickCheck
+import Data.Complex
+import Data.Matrix.Static.LinearAlgebra
+import Data.Singletons
+import qualified Data.Matrix.Static.Dense as D
+import qualified Data.Matrix.Static.Sparse as S
+import Data.Vector.Storable (Storable)
+import Data.List.Ordered
+import Data.Ord
+import qualified Data.Vector.Generic as G
+import Data.AEq
+import Control.Monad
+
+class Approx a where
+    (~=) :: a -> a -> Bool
+    infixl 4 ~=
+
+instance Approx Double where
+    a ~= b = abs (a - b) < 1e-10
+
+instance Approx Float where
+    a ~= b = abs (a - b) < 1e-3
+
+instance Approx (Complex Double) where
+    a ~= b = r1 ~= r2 && i1 ~= i2
+      where
+        (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 (G.Vector v a, Arbitrary a, SingI m, SingI n)
+    => Arbitrary (D.Matrix m n v a) where
+        arbitrary = D.fromList <$> vector (m*n)
+          where
+            m = fromIntegral $ fromSing (sing :: Sing m)
+            n = fromIntegral $ fromSing (sing :: Sing n)
+        shrink _v = []
+
+instance (G.Vector v a, Arbitrary a, SingI m, SingI n, S.Zero a)
+    => Arbitrary (S.SparseMatrix m n v a) where
+        arbitrary = do
+            vals <- filter (/=S.zero) <$> vector p
+            xs <- fmap (nubSortBy (comparing fst)) $ forM vals $ \v -> do
+                i <- choose (0, m-1)
+                j <- choose (0, n-1)
+                return ((i,j),v)
+            return $ S.fromTriplet $ map (\((a,b),c) -> (a,b,c)) xs
+          where
+            p = (m * n) `div` 10
+            m = fromIntegral $ fromSing (sing :: Sing m)
+            n = fromIntegral $ fromSing (sing :: Sing n)
+        shrink _v = []