diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -18,18 +18,39 @@
 -------
 
 ```haskell
-let mat = D.matrix [ [1,0,3]
-                   , [0,5,6]
-                   , [0,0,0] ] :: Matrix 3 3 Double
-    mat' = D.convertAny mat :: SparseMatrix 3 3 Double
+{-# LANGUAGE DataKinds #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE Rank2Types #-}
+{-# LANGUAGE GADTs #-}
+{-# LANGUAGE TypeOperators #-}
 
-print mat
-print mat'
+import Data.Matrix.Static.LinearAlgebra
+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.Prelude hiding ((@@))
+import Data.Singletons.TypeLits
+import Data.Complex
+import System.Random
+import Control.Monad
+import Data.Type.Equality
 
-print $ eigs (sing :: Sing 1) mat == eigs (sing :: Sing 1) mat'
+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
+      in v @@ S.diag d @@ G.transpose v
 
-print $ cholesky mat
+main :: IO ()
+main = do
+    n <- randomRIO (2, 6)
+    vals <- replicateM (n*n) $ randomRIO (-100,100) :: IO [Double]
 
-print $ mat %*% mat %*% mat
-print $ mat' %*% mat' %*% mat
+    withSomeSing (fromIntegral n) $ \sn@(SNat :: Sing n) ->
+        let s0 = SNat :: Sing 2
+            sn2 = sn %- s0
+        in case testEquality (s0 %<= sn2) STrue of
+            Just Refl -> do
+                let mat = G.fromList vals :: Matrix n n Double
+                print $ f mat
+            Nothing -> error $ "Requiring Matrix size >= 4, but got: " <> show n
 ```
diff --git a/cbits/eigen-basic.cpp b/cbits/eigen-basic.cpp
--- a/cbits/eigen-basic.cpp
+++ b/cbits/eigen-basic.cpp
@@ -70,8 +70,9 @@
     const void* val, const void* outer, const void* inner, int r2, int c2, int s,
     const void* p1, int r1, int c1), (p,r,c,val,outer,inner,r2,c2,s,p1,r1,c1));
 
+
 template <class T>
-RET ss_mul( void** v, void** o, void** i, int r, int c, int* s,
+RET ss_mul( void** v, void* o, void** i, int r, int c, int* s,
     const void* v1, const void* o1, const void* i1, int r1, int c1, int s1,
     const void* v2, const void* o2, const void* i2, int r2, int c2, int s2)
 {
@@ -79,24 +80,96 @@
     MapSparseMatrix a(r1, c1, s1, (int*)o1, (int*)i1, (T*)v1);
     MapSparseMatrix b(r2, c2, s2, (int*)o2, (int*)i2, (T*)v2);
     SparseMatrix<T> M = (a * b).pruned();
+
+    memcpy(o, M.outerIndexPtr(), (c+1) * sizeof(int));
+
     *s = M.nonZeros();
-    T* p1 = (T*) malloc(*s * sizeof(T));
-    memcpy(p1, M.valuePtr(), *s * sizeof(T));
+    T* p1 = (T*) malloc((*s) * sizeof(T));
+    memcpy(p1, M.valuePtr(), (*s) * sizeof(T));
     *v = p1;
-    int* p2 = (int*) malloc(*s * sizeof(int));
-    memcpy(p2, M.innerIndexPtr(), *s * sizeof(T));
+    int* p2 = (int*) malloc((*s) * sizeof(int));
+    memcpy(p2, M.innerIndexPtr(), (*s) * sizeof(int));
     *i = p2;
-    int* p3 = (int*) malloc(c * sizeof(int));
-    memcpy(p3, M.outerIndexPtr(), c * sizeof(int));
-    *o = p3;
     return 0;
 }
 API(ss_mul, (int code,
-    void** v, void** o, void** i, int r, int c, int* s,
+    void** v, void* o, void** i, int r, int c, int* s,
     const void* v1, const void* o1, const void* i1, int r1, int c1, int s1,
     const void* v2, const void* o2, const void* i2, int r2, int c2, int s2),
     (v,o,i,r,c,s,v1,o1,i1,r1,c1,s1,v2,o2,i2,r2,c2,s2));
 
+
+template <class T>
+RET ss_plus( void** v, void* o, void** i, int r, int c, int* s,
+    const void* v1, const void* o1, const void* i1, int r1, int c1, int s1,
+    const void* v2, const void* o2, const void* i2, int r2, int c2, int s2)
+{
+    typedef Map<SparseMatrix<T> > MapSparseMatrix;
+    MapSparseMatrix a(r1, c1, s1, (int*)o1, (int*)i1, (T*)v1);
+    MapSparseMatrix b(r2, c2, s2, (int*)o2, (int*)i2, (T*)v2);
+    SparseMatrix<T> M = a + b;
+
+    memcpy(o, M.outerIndexPtr(), (c+1) * sizeof(int));
+
+    *s = M.nonZeros();
+    T* p1 = (T*) malloc((*s) * sizeof(T));
+    memcpy(p1, M.valuePtr(), (*s) * sizeof(T));
+    *v = p1;
+    int* p2 = (int*) malloc((*s) * sizeof(int));
+    memcpy(p2, M.innerIndexPtr(), (*s) * sizeof(int));
+    *i = p2;
+    return 0;
+}
+API(ss_plus, (int code,
+    void** v, void* o, void** i, int r, int c, int* s,
+    const void* v1, const void* o1, const void* i1, int r1, int c1, int s1,
+    const void* v2, const void* o2, const void* i2, int r2, int c2, int s2),
+    (v,o,i,r,c,s,v1,o1,i1,r1,c1,s1,v2,o2,i2,r2,c2,s2));
+
+template <class T>
+RET ss_cmul( void** v, void* o, void** i, int r, int c, int* s,
+    const void* v1, const void* o1, const void* i1, int r1, int c1, int s1,
+    const void* v2, const void* o2, const void* i2, int r2, int c2, int s2)
+{
+    typedef Map<SparseMatrix<T> > MapSparseMatrix;
+    MapSparseMatrix a(r1, c1, s1, (int*)o1, (int*)i1, (T*)v1);
+    MapSparseMatrix b(r2, c2, s2, (int*)o2, (int*)i2, (T*)v2);
+
+    SparseMatrix<T> M = a.cwiseProduct(b);
+
+    memcpy(o, M.outerIndexPtr(), (c+1) * sizeof(int));
+
+    *s = M.nonZeros();
+    T* p1 = (T*) malloc((*s) * sizeof(T));
+    memcpy(p1, M.valuePtr(), (*s) * sizeof(T));
+    *v = p1;
+    int* p2 = (int*) malloc((*s) * sizeof(int));
+    memcpy(p2, M.innerIndexPtr(), (*s) * sizeof(int));
+    *i = p2;
+    return 0;
+}
+API(ss_cmul, (int code,
+    void** v, void* o, void** i, int r, int c, int* s,
+    const void* v1, const void* o1, const void* i1, int r1, int c1, int s1,
+    const void* v2, const void* o2, const void* i2, int r2, int c2, int s2),
+    (v,o,i,r,c,s,v1,o1,i1,r1,c1,s1,v2,o2,i2,r2,c2,s2));
+
+template <class T>
+RET sd_plus( void* p, int r, int c,
+    const void* val, const void* outer, const void* inner,
+    int r2, int c2, int s,
+    const void* p1, int r1, int c1)
+{
+    typedef Map< Matrix<T,Dynamic,Dynamic> > MapMatrix;
+    MapMatrix x((T*)p, r, c);
+    x = matrix<T>(p1,r1,c1);
+    x += smatrix<T>(val, outer, inner, r2, c2, s);
+    return 0;
+}
+API(sd_plus, (int code,
+    void* p, int r, int c,
+    const void* val, const void* outer, const void* inner, int r2, int c2, int s,
+    const void* p1, int r1, int c1), (p,r,c,val,outer,inner,r2,c2,s,p1,r1,c1));
 
 #define UNOP(name) \
 extern "C" RET __attribute__((noinline)) eigen_##name(int code, void* p, int r, int c, const void* p1, int r1, int c1) {\
diff --git a/cbits/eigen-solver.cpp b/cbits/eigen-solver.cpp
--- a/cbits/eigen-solver.cpp
+++ b/cbits/eigen-solver.cpp
@@ -18,7 +18,9 @@
     MapComplexMatrix V((T3*)v, n, k);
 
     DenseGenMatProd<double> op(M);
-    GenEigsSolver< double, LARGEST_MAGN, DenseGenMatProd<double> > eigs(&op, k, 2*k+1);
+    int ncv = 2 * k + 1;
+    ncv = (ncv <= n) ? ncv : n;
+    GenEigsSolver< double, LARGEST_MAGN, DenseGenMatProd<double> > eigs(&op, k, ncv);
     eigs.init();
     int nconv = eigs.compute();
     if(eigs.info() == 0)
@@ -42,7 +44,9 @@
     MapComplexMatrix V((T3*)v, n, k);
 
     SparseGenMatProd<double> op(M);
-    GenEigsSolver< double, LARGEST_MAGN, SparseGenMatProd<double> > eigs(&op, k, 2*k+1);
+    int ncv = 2 * k + 1;
+    ncv = (ncv <= n) ? ncv : n;
+    GenEigsSolver< double, LARGEST_MAGN, SparseGenMatProd<double> > eigs(&op, k, ncv);
     eigs.init();
     int nconv = eigs.compute();
     if(eigs.info() == 0)
@@ -62,3 +66,18 @@
 }
 API(cholesky, (int code,
     void* px, const void* pa, int n), (px,pa,n));
+
+
+/*
+template <class T>
+RET bdcsvd(void* px, int r, int c
+{
+    typedef Map< Matrix<T,Dynamic,Dynamic> > MapMatrix;
+    MapMatrix x((T*)px, n, n);
+    MapMatrix A((T*)pa, n, n);
+    x = A.llt().matrixL();
+    return 0;
+}
+API(cholesky, (int code,
+    void* px, const void* pa, int n), (px,pa,n));
+*/
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: 6910822a984416806867bac316c8f9685ade7b6a7672d4437c67fd4770977af1
+-- hash: ef7b97a49f99defe3cd58f2c5a3411fe8c958f0ad40f4b3f0bd763720b643db8
 
 name:           matrix-sized
-version:        0.0.2
+version:        0.0.3
 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
@@ -375,6 +375,11 @@
   type: git
   location: https://github.com/kaizhang/matrix-sized
 
+flag parallel
+  description: Enable multithreading
+  manual: True
+  default: False
+
 library
   exposed-modules:
       Data.Matrix.Static.Dense
@@ -403,10 +408,27 @@
       stdc++
   build-depends:
       base >=4.10 && <5
-    , binary
-    , bytestring
-    , deepseq
     , primitive >=0.6.4.0
     , singletons
+    , vector >=0.11
+  if flag(parallel)
+    cc-options: -fopenmp
+    cxx-options: -fopenmp -std=c++11
+    ld-options: -fopenmp
+  default-language: Haskell2010
+
+test-suite test
+  type: exitcode-stdio-1.0
+  main-is: Main.hs
+  hs-source-dirs:
+      tests
+  ghc-options: -Wall -threaded
+  build-depends:
+      base
+    , matrix-sized
+    , primitive >=0.6.4.0
+    , singletons
+    , tasty
+    , tasty-quickcheck
     , vector >=0.11
   default-language: Haskell2010
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
@@ -41,8 +41,7 @@
     , C.fromRows
     , C.fromColumns
     , C.unsafeFromVector
-
-    , diag
+    , replicate
     , diagRect
 
     -- * Conversions
@@ -91,10 +90,9 @@
     , C.create
     ) where
 
-import           Control.DeepSeq                   hiding (force)
 import           Control.Monad                     (liftM)
 import qualified Data.Vector.Generic               as G
-import Prelude hiding (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)
 import GHC.TypeLits (type (<=))
 import Data.Singletons
 import Data.Tuple (swap)
@@ -136,9 +134,6 @@
         recip = C.map recip
         fromRational = undefined
 
-instance NFData (v a) => NFData (Matrix r c v a) where
-    rnf (Matrix vec) = rnf vec
-
 instance G.Vector v a => C.Matrix Matrix v a where
     -- | O(1) Return the size of matrix.
     dim :: forall r c. Matrix r c v a -> (Int, Int)
@@ -199,20 +194,14 @@
     sequence_ (Matrix vec) = G.sequence_ vec
     {-# INLINE sequence_ #-}
 
-{-
--- | O(m*n) Create an identity matrix
-ident :: (Num a, G.Vector v a)
-      => Matrix n n v a
-ident = diagRect 0 $ replicate 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, SingI n)
-     => Matrix n 1 v a       -- ^ diagonal
-     -> Matrix n n v a
-diag = diagRect 0
-{-# INLINE diag #-}
+-- | O(m*n) Create a constant matrix.
+replicate :: forall r c v a. (G.Vector v a, SingI r, SingI c)
+          => a -> Matrix r c v a
+replicate = C.unsafeFromVector . G.replicate (r*c)
+  where
+    r = fromIntegral $ fromSing (sing :: Sing r)
+    c = fromIntegral $ fromSing (sing :: Sing c)
+{-# INLINE replicate #-}
 
 -- | 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)
diff --git a/src/Data/Matrix/Static/Dense/Mutable.hs b/src/Data/Matrix/Static/Dense/Mutable.hs
--- a/src/Data/Matrix/Static/Dense/Mutable.hs
+++ b/src/Data/Matrix/Static/Dense/Mutable.hs
@@ -18,7 +18,6 @@
    , C.replicate
    ) where
 
-import           Control.DeepSeq
 import qualified Data.Vector.Generic.Mutable as GM
 import           Prelude                     hiding (read, replicate)
 import Data.Singletons
@@ -29,9 +28,6 @@
 -- | Column-major mutable matrix.
 data MMatrix :: C.MMatrixKind where
     MMatrix :: (SingI r, SingI c) => v s a -> MMatrix r c v s a
-
-instance (NFData (v s a)) => NFData (MMatrix r c v s a) where
-    rnf (MMatrix vec) = rnf vec
 
 instance GM.MVector v a => C.MMatrix MMatrix v a where
     dim :: forall r c s. MMatrix r c v s a -> (Int, Int)
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
@@ -41,6 +41,7 @@
 import Text.Printf (printf)
 import Prelude hiding (map, mapM, mapM_, sequence, sequence_)
 import qualified Data.List as L
+import Data.Tuple (swap)
 import Data.Kind (Type)
 import GHC.TypeLits (Nat, type (<=))
 import Data.Singletons (SingI, Sing, fromSing, sing, withSomeSing)
@@ -62,7 +63,7 @@
     -- | Convert matrix to vector in column order.
     -- Default algorithm is O((m*n) * O(unsafeIndex)).
     flatten :: mat r c v a -> v a
-    flatten mat = G.generate (r*c) $ \i -> unsafeIndex mat (i `divMod` r)
+    flatten mat = G.generate (r*c) $ \i -> unsafeIndex mat (swap $ i `divMod` r)
       where
         (r,c) = dim mat
     {-# INLINE flatten #-}
@@ -90,7 +91,7 @@
 
     transpose :: (SingI r, SingI c) => mat r c v a -> mat c r v a
     transpose mat = unsafeFromVector $ G.generate (r*c) $ \x ->
-        unsafeIndex mat $ x `divMod` r
+        unsafeIndex mat $ x `divMod` c
       where
        (r, c) = dim mat
     {-# INLINE transpose #-}
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
@@ -4,6 +4,9 @@
     , c_ds_mul
     , c_sd_mul
     , c_ss_mul
+    , c_ss_cmul
+    , c_sd_plus
+    , c_ss_plus
     , c_inverse
     , c_cholesky
     , c_eigs
@@ -39,13 +42,33 @@
           -> Ptr a -> CInt -> CInt
           -> IO CString
 
+foreign import ccall "eigen_sd_plus"
+    c_sd_plus :: CInt
+          -> Ptr a -> CInt -> CInt
+          -> Ptr a -> Ptr CInt -> Ptr CInt -> CInt -> CInt -> CInt
+          -> Ptr a -> CInt -> CInt
+          -> IO CString
+
 foreign import ccall "eigen_ss_mul"
     c_ss_mul :: CInt
-          -> Ptr (Ptr a) -> Ptr (Ptr CInt) -> Ptr (Ptr CInt) -> CInt -> CInt -> Ptr 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
 
+foreign import ccall "eigen_ss_cmul"
+    c_ss_cmul :: 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
+
+foreign import ccall "eigen_ss_plus"
+    c_ss_plus :: 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
              
 foreign import ccall "eigen_inverse"
     c_inverse :: CInt
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
@@ -11,15 +11,15 @@
 module Data.Matrix.Static.LinearAlgebra
     ( Arithmetic(..)
     , Factorization(..)
-    , inverse
     , module Data.Matrix.Static.LinearAlgebra.Types
     ) where
 
 import qualified Data.Vector.Storable as VS
 import System.IO.Unsafe (unsafePerformIO)
 import Data.Complex (Complex)
-import Data.Singletons
-import GHC.TypeLits (type (<=), type (-))
+import Data.Singletons.Prelude hiding ((@@), type (==))
+import Data.Type.Bool (If)
+import Data.Type.Equality (type (==))
 
 import qualified Data.Matrix.Static.Dense as D
 import qualified Data.Matrix.Static.Sparse as S
@@ -28,35 +28,68 @@
 import qualified Data.Matrix.Static.Internal as Internal
 import Data.Matrix.Static.LinearAlgebra.Types
 
-class Arithmetic (mat1 :: C.MatrixKind)
-                 (mat2 :: C.MatrixKind)
-                 (mat3 :: C.MatrixKind) |
-                 mat1 mat2 -> mat3 where
-    (%*%) :: (Numeric a, SingI n, SingI m)
-          => mat1 n p VS.Vector a
-          -> mat2 p m VS.Vector a
+class Arithmetic (mat1 :: C.MatrixKind) (mat2 :: C.MatrixKind) where
+    -- | Matrix multiplication
+    (@@) :: ( 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
+    infixr 8 @@
+
+    (%+%) :: ( 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
+    infixr 8 %+%
+
+    (%-%) :: ( 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
+    infixr 8 %-%
+
+    (%*%) :: ( 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
     infixr 8 %*%
 
-instance Arithmetic D.Matrix D.Matrix D.Matrix where
-    (%*%) = withFun2 Internal.c_dd_mul
+instance Arithmetic D.Matrix D.Matrix where
+    (@@) = withFun2 Internal.c_dd_mul
+    (%+%) = (+)
+    (%-%) = (-)
+    (%*%) = (*)
 
-instance Arithmetic D.Matrix S.SparseMatrix D.Matrix where
-    (%*%) = withDS Internal.c_ds_mul
+instance Arithmetic D.Matrix S.SparseMatrix where
+    (@@) = withDS Internal.c_ds_mul
+    (%+%) = flip (%+%)
+    (%-%) a b = a %+% C.map negate b
+    (%*%) = undefined
 
-instance Arithmetic S.SparseMatrix D.Matrix D.Matrix where
-    (%*%) = withSD Internal.c_sd_mul
+instance Arithmetic S.SparseMatrix D.Matrix where
+    (@@) = withSD Internal.c_sd_mul
+    (%+%) = withSD Internal.c_sd_plus
+    (%-%) a b = a %+% C.map negate b
+    (%*%) = undefined
 
-instance Arithmetic S.SparseMatrix S.SparseMatrix S.SparseMatrix where
-    (%*%) = withSS Internal.c_ss_mul
+instance Arithmetic S.SparseMatrix S.SparseMatrix where
+    (@@) = withSS Internal.c_ss_mul
+    (%+%) = withSS Internal.c_ss_plus
+    (%-%) a b = a %+% C.map negate b
+    (%*%) = withSS Internal.c_ss_cmul
 
-inverse :: (SingI n, Numeric a) => Matrix n n a -> Matrix n n a
-inverse = withFun1 Internal.c_inverse
 
 class Factorization mat where
+    -- | Matrix inverse
+    inverse :: (SingI n, Numeric a) => mat n n VS.Vector a -> mat n n VS.Vector a
+
     -- | Eigenvalues (not ordered) and
     -- eigenvectors (as columns) of a general square matrix.
-    eigs :: (SingI k, SingI n, k <= n - 2)
+    eigs :: (SingI k, SingI n, (k <= n - 2) ~ 'True)
          => Sing k
          -> mat n n VS.Vector Double
          -> (Matrix k 1 (Complex Double), Matrix n k (Complex Double))
@@ -65,6 +98,9 @@
     cholesky :: (Numeric a, SingI n) => mat n n VS.Vector a -> mat n n VS.Vector a
 
 instance Factorization D.Matrix where
+
+    inverse = withFun1 Internal.c_inverse
+
     eigs s mat = unsafePerformIO $ do
         m1 <- CM.new
         m2 <- CM.new
@@ -82,6 +118,8 @@
     {-# INLINE cholesky #-}
 
 instance Factorization S.SparseMatrix where
+    inverse = undefined
+
     eigs s mat = unsafePerformIO $ do
         m1 <- CM.new
         m2 <- CM.new
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
@@ -38,7 +38,7 @@
 import qualified Data.Matrix.Static.Generic.Mutable as CM
 import qualified Data.Matrix.Static.Generic as C
 
-class (S.Zero a, Storable a) => Numeric a where
+class (S.Zero a, Storable a, Num a) => Numeric a where
     foreignType :: a -> CInt
 
 instance Numeric Float where foreignType _ = 0
@@ -127,16 +127,19 @@
 {-# INLINE withSD #-}
 
 mkSparseMatrix :: forall r c a. (Storable a, SingI r, SingI c)
-    => (Ptr (Ptr a) -> Ptr (Ptr CInt) -> Ptr (Ptr CInt) -> IO Int)
+    => (Ptr (Ptr a) -> Ptr (Ptr CInt) -> Ptr CInt -> IO Int)
     -> IO (SparseMatrix r c a)
-mkSparseMatrix f = alloca $ \ppv -> alloca $ \ppi -> alloca $ \ppo -> do
-    n <- f ppv ppi ppo
-    pv <- peek ppv >>= newForeignPtr finalizerFree
-    pinner <- peek ppi >>= newForeignPtr finalizerFree
-    pouter <- peek ppo >>= newForeignPtr finalizerFree
+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)
-        (VS.unsafeFromForeignPtr0 pouter $ c + 1)
+        outer
   where
     c = fromIntegral $ fromSing (sing :: Sing c)
 {-# INLINE mkSparseMatrix #-}
@@ -144,7 +147,7 @@
 withSS :: forall r1 c1 r2 c2 r3 c3 a.
             (SingI r3, SingI c3, Numeric a)
        => ( CInt
-         -> Ptr (Ptr a) -> Ptr (Ptr CInt) -> Ptr (Ptr CInt) -> CInt -> CInt -> Ptr 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 )
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
@@ -34,10 +34,10 @@
 
     -- * Construction
     , C.empty
+    , fromTriplet
     , C.fromVector
     , C.fromList
     , C.unsafeFromVector
-
     , diag
     , diagRect
 
@@ -49,12 +49,13 @@
     , C.convertAny
    ) where
 
-import           Control.DeepSeq
 import qualified Data.Vector.Generic as G
+import qualified Data.Vector.Generic.Mutable as GM
 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 GHC.TypeLits (type (<=))
@@ -106,9 +107,22 @@
         (r,c) = C.dim mat
         vals = unlines $ map (unwords . map show . G.toList) $ C.toRows mat
 
-instance (NFData (v a)) => NFData (SparseMatrix r c v a) where
-    rnf (SparseMatrix vec inner outer) = rnf vec
+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)
@@ -151,6 +165,8 @@
         c = fromIntegral $ fromSing (sing :: Sing c)
     {-# INLINE unsafeFromVector #-}
 
+    transpose (SparseMatrix val inner outer) = undefined 
+
     thaw = undefined
     {-# INLINE thaw #-}
 
@@ -167,7 +183,38 @@
     imap = undefined
     {-# INLINE map #-}
 
--- | O(m*n) Create a square matrix with given diagonal.
+-- | 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)
+            => t (Int, Int, a) -> SparseMatrix r c v a
+fromTriplet triplets = SparseMatrix val inner outer
+  where
+    outer = S.scanl (+) 0 $ S.create $ do
+        vec <- SM.replicate c 0
+        _ <- flip mapM triplets $ \(_, j, _) -> 
+            SM.modify 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
+        (,) <$> 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(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
      -> SparseMatrix n n v a
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
@@ -13,12 +13,10 @@
      MSparseMatrix(..)
    ) where
 
-import           Control.DeepSeq
 import qualified Data.Vector.Generic.Mutable as GM
 import qualified Data.Vector.Storable as S
 import           Prelude                     hiding (read, replicate)
 import Data.Singletons
-import Control.Monad.Primitive     (PrimMonad, PrimState)
 
 import qualified Data.Matrix.Static.Generic.Mutable as C
 
@@ -33,9 +31,6 @@
                                       -- (resp. row) the index of the first
                                       -- non-zero in the previous two arrays.
                   -> MSparseMatrix r c v s a
-
-instance (NFData (v s a)) => NFData (MSparseMatrix r c v s a) where
-    rnf (MSparseMatrix vec inner outer) = rnf vec
 
 instance GM.MVector v a => C.MMatrix MSparseMatrix v a where
     dim :: forall r c s. MSparseMatrix r c v s a -> (Int, Int)
diff --git a/tests/Main.hs b/tests/Main.hs
new file mode 100644
--- /dev/null
+++ b/tests/Main.hs
@@ -0,0 +1,36 @@
+{-# 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
+-}
+
+main :: IO ()
+main = defaultMain $ testGroup "Main"
+    [ testProperty "" square
+    ]
+