diff --git a/CHANGELOG.markdown b/CHANGELOG.markdown
--- a/CHANGELOG.markdown
+++ b/CHANGELOG.markdown
@@ -1,3 +1,10 @@
+	0.3.1
+	* Changed `SpMatrix` to use `foldlWithKey'` for efficiency (Joshua Moerman)
+	
+	* Bumped LTS to 11.3 (GHC 8.2.2)
+	
+	* Removed unneeded dependencies from stack.yaml
+
 	0.3
 	* Fixed a number of instances, un-commented tests (Joshua Moerman)
 	
diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md
--- a/CONTRIBUTORS.md
+++ b/CONTRIBUTORS.md
@@ -1,9 +1,9 @@
 # Original author 
 
-Marco Zocca (@ocramz)
+[Marco Zocca](https://github.com/ocramz)
 
 # Contributors
 
-Gregory Schwartz (@GregorySchwartz)
+[Gregory Schwartz](https://github.com/GregorySchwartz)
 
-Joshua Moerman (@Jaxan)
+[Joshua Moerman](https://github.com/Jaxan)
diff --git a/sparse-linear-algebra.cabal b/sparse-linear-algebra.cabal
--- a/sparse-linear-algebra.cabal
+++ b/sparse-linear-algebra.cabal
@@ -1,5 +1,5 @@
 name:                sparse-linear-algebra
-version:             0.3
+version:             0.3.1
 synopsis:            Numerical computing in native Haskell
 description:
   /Overview/
@@ -17,7 +17,7 @@
 license-file:        LICENSE
 author:              Marco Zocca
 maintainer:          zocca.marco gmail
-copyright:           2016-2017 Marco Zocca
+copyright:           2016-2018 Marco Zocca
 category:            Numeric
 build-type:          Simple
 extra-source-files:  README.md
@@ -72,18 +72,20 @@
                        Numeric.LinearAlgebra.LinearSolvers.Experimental
                        Numeric.LinearAlgebra.EigenSolvers.Experimental
   build-depends:       base >= 4.7 && < 5
-                     , primitive
                      -- , deepseq
                      , containers
+                     , exceptions
+                     , mtl >= 2.2.1                                          
+                     , primitive
+                     , transformers                                      
                      , vector
                      , vector-algorithms
                      -- , hybrid-vectors
-                     , exceptions
-                     , mtl >= 2.2.1                     
-                     , transformers
                      -- , writer-cps-transformers
                      -- , monad-log == 0.1.1.0
                      -- , monad-par
+
+                     -- -- DEBUG
                      , QuickCheck >= 2.8.2
                      , hspec
                      -- , primitive >= 0.6.1.0
diff --git a/src/Data/Sparse/SpMatrix.hs b/src/Data/Sparse/SpMatrix.hs
--- a/src/Data/Sparse/SpMatrix.hs
+++ b/src/Data/Sparse/SpMatrix.hs
@@ -843,13 +843,15 @@
 
 -- ** Partial inner product
 -- | Contract row `i` of A with column `j` of B up to an index `n`, i.e. summing over repeated indices: 
--- Aij Bjk , for j in [0 .. n] 
-contractSub :: Elt a => SpMatrix a -> SpMatrix a -> IxRow -> IxCol -> Int -> a
+-- Aik Bkj , for k in [0 .. n] 
+contractSub :: Num a => SpMatrix a -> SpMatrix a -> IxRow -> IxCol -> Int -> a
 contractSub a b i j n
   | ncols a == nrows b &&
     isValidIxSM a (i,j) &&
-    n <= ncols a = sum $ map (\i' -> (a@@!(i,i'))*b@@!(i',j)) [0 .. n]
-  | otherwise = error "contractSub : n must be <= i"
+    n <= ncols a = case I.lookup i (smData a) of
+                     Nothing -> 0
+                     Just r  -> I.foldlWithKey' (\acc k x -> if k > n then acc else acc + x * b @@! (k, j)) 0 r
+  | otherwise = error "contractSub: matrices have incompatible dimensions or i,j,n are too big"
 
 
 
diff --git a/src/Numeric/LinearAlgebra/Sparse.hs b/src/Numeric/LinearAlgebra/Sparse.hs
--- a/src/Numeric/LinearAlgebra/Sparse.hs
+++ b/src/Numeric/LinearAlgebra/Sparse.hs
@@ -478,7 +478,7 @@
       luInit | isNz u00 = return (1, l0, u0)
              | otherwise = oops (0 :: Int)
         where
-          l0 = insertCol (eye n) ((extractSubCol aa 0 (1, n-1)) ./ u00 ) 0
+          l0 = insertCol (eye n) (extractSubCol aa 0 (1, n-1) ./ u00 ) 0
           u0 = insertRow (zeroSM n n) (extractRow aa 0) 0   -- initial U
           u00 = u0 @@! (0,0)  -- make sure this is non-zero by applying permutation
       luUpd (i, l, u) = do -- (i + 1, l', u')