diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,8 +1,10 @@
 # Changelog
 
-## [0.1.0.0] - 2019-08-16
-Initial release.
 
+## [0.2.1] - 2019-08-22
+- fixed haddock build for ghc 8.6
+- fixed numerical instability in gaussian elimination algorithm
+
 ## [0.2] - 2019-08-21
 - removed singletons from dependencies
 - removed eigen from dependencies
@@ -17,3 +19,6 @@
   error when using haddock with ghc 8.6
 
 - ensured backwards compatibility down to ghc 8.0.2 (stackage lts-9.21)
+
+## [0.1.0.0] - 2019-08-16
+Initial release.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,31 @@
+{-
+Disable some errors and warnings during the haddock pass
+  (caused by compiler plugins and hs-boot)
+ -}
+{-# OPTIONS_GHC -Wall #-}
+module Main (main) where
+
+import Distribution.Simple
+import Distribution.Simple.Setup
+
+main :: IO ()
+main = defaultMainWithHooks simpleUserHooks
+  { confHook = \a -> confHook simpleUserHooks a . tweakFlags }
+
+tweakFlags :: ConfigFlags -> ConfigFlags
+tweakFlags flags = flags
+  { configProgramArgs = addHaddockArgs (configProgramArgs flags) }
+
+addHaddockArgs :: [(String, [String])] -> [(String, [String])]
+addHaddockArgs []
+  = [("haddock", newHaddockGhcArgs)]
+addHaddockArgs (("haddock", args):otherProgsArgs)
+  = ("haddock", args ++ newHaddockGhcArgs) : otherProgsArgs
+addHaddockArgs (progArgs:otherProgsArgs)
+  = progArgs : addHaddockArgs otherProgsArgs
+
+newHaddockGhcArgs :: [String]
+newHaddockGhcArgs =
+  [ "--optghc=-fdefer-type-errors"
+  , "--optghc=-fno-warn-deferred-type-errors"
+  ]
diff --git a/cabal.project b/cabal.project
deleted file mode 100644
--- a/cabal.project
+++ /dev/null
@@ -1,4 +0,0 @@
-packages: ./
-
-package sparse-tensor
-  haddock-options: --optghc=-fdefer-type-errors
diff --git a/sparse-tensor.cabal b/sparse-tensor.cabal
--- a/sparse-tensor.cabal
+++ b/sparse-tensor.cabal
@@ -1,5 +1,5 @@
 name:                sparse-tensor
-version:             0.2
+version:             0.2.1
 synopsis:            typesafe tensor algebra library
 description:
     .
@@ -17,13 +17,16 @@
 maintainer:          tobi.reinhart@fau.de, nils.alex@fau.de
 copyright:           2019 Tobias Reinhart and Nils Alex
 category:            Data, Math, Algebra
-build-type:          Simple
-cabal-version:       >=1.10
+build-type:          Custom
+cabal-version:       1.24
 extra-source-files:
     README.md
     CHANGELOG.md
-    stack.yaml
-    cabal.project
+
+custom-setup
+  setup-depends:
+    base >= 4.9 && < 5,
+    Cabal >= 1.24 && < 3.0
 
 source-repository head
   type:     git
diff --git a/src/Math/Tensor.hs b/src/Math/Tensor.hs
--- a/src/Math/Tensor.hs
+++ b/src/Math/Tensor.hs
@@ -387,24 +387,26 @@
 
 isZero :: forall n. KnownNat n => Proxy n -> IsZero n
 isZero _ = case sameNat (Proxy @n) (Proxy @0)
-                        of Nothing   -> unsafeCoerce (NonZero @1)
-                           Just Refl -> Zero
+             of Nothing   -> unsafeCoerce (NonZero @1)
+                Just Refl -> Zero
 
-fromList' :: forall n. KnownNat n => Proxy n -> forall (a :: *). [a] -> Maybe (IndList n a)
-fromList' (n :: Proxy n') xs = case isZero n
-                               of Zero    -> case xs
-                                               of [] -> Just Empty
-                                                  _  -> Nothing
-                                  NonZero -> case xs
-                                               of []    -> Nothing
-                                                  x:xs' -> case fromList' (undefined :: (Proxy (n' - 1))) xs'
-                                                             of Just v  -> Just (x `Append` v)
-                                                                Nothing -> Nothing
+isZero' :: forall n. KnownNat n => IsZero n
+isZero' = case sameNat (Proxy @n) (Proxy @0)
+            of Nothing   -> unsafeCoerce (NonZero @1)
+               Just Refl -> Zero
 
 -- | Construction of a length typed @'IndList'@ from an untyped list.
 
-fromList :: forall n.KnownNat n => forall (a :: *). [a] -> Maybe (IndList n a)
-fromList = fromList' (undefined :: Proxy n)
+fromList :: forall n. KnownNat n => forall (a :: *). [a] -> Maybe (IndList n a)
+fromList xs = case isZero' @n
+                of Zero    -> case xs
+                                of [] -> Just Empty
+                                   _  -> Nothing
+                   NonZero -> case xs
+                                of []    -> Nothing
+                                   x:xs' -> case fromList xs'
+                                              of Just v  -> Just (x `Append` v)
+                                                 Nothing -> Nothing
 
 -- | Construction of a length typed @'IndList'@ (partial function).
 
diff --git a/src/Math/Tensor/Internal/LinearAlgebra.hs b/src/Math/Tensor/Internal/LinearAlgebra.hs
--- a/src/Math/Tensor/Internal/LinearAlgebra.hs
+++ b/src/Math/Tensor/Internal/LinearAlgebra.hs
@@ -47,82 +47,81 @@
 pivotsU mat = go (0,0)
   where
     go (i,j)
-      = case findPivot mat (i,j) of
+      = case findPivot mat e (i,j) of
           Nothing       -> []
           Just (i', j') -> j' : go (i'+1, j'+1)
+    maxAbs = maximum $ map (maximum . map abs) $ toLists mat
+    e = eps * maxAbs
 
 
--- naive check for numerical zero
-
-eps :: Double -> Bool
-eps = (< 1e-12) . abs
+eps :: Double
+eps = 1e-12
 
 -- find next pivot in upper triangular matrix
 
-findPivot :: Matrix Double -> (Int, Int) -> Maybe (Int, Int)
-findPivot mat (i, j)
+findPivot :: Matrix Double -> Double -> (Int, Int) -> Maybe (Int, Int)
+findPivot mat e (i, j)
     | n == j = Nothing
     | m == i = Nothing
     | otherwise = case nonZeros of
                     []          -> if n == j+1
                                    then Nothing
-                                   else findPivot mat (i, j+1)
+                                   else findPivot mat e (i, j+1)
                     (pi, pj):_  -> Just (pi, pj+j)
     where
         m = rows mat
         n = cols mat
         col = mat ¿ [j]
-        nonZeros = filter (\(i', _) -> i' >= i) $ find (not . eps) col
+        nonZeros = filter (\(i', _) -> i' >= i) $ find (not . (< e) . abs) col
 
--- | Find pivot element below position (i, j) with greatest absolute value.
+-- | Find pivot element below position (i, j) with greatest absolute value in the ST monad.
 
-findPivotMax :: Matrix Double -> Int -> Int -> Maybe (Int, Int)
-findPivotMax mat i j
-    | n == j = Nothing
-    | m == i = Nothing
-    | otherwise = case nonZeros of
-                    [] -> if n == j+1
-                          then Nothing
-                          else findPivotMax mat i (j+1)
-                    _  -> Just (pi, pj+j)
-    where
-        m = rows mat
-        n = cols mat
-        col = mat ¿ [j]
-        nonZeros = filter (\(i', _) -> i' >= i) $ find (not . eps) col
-        (pi, pj) = maximumBy (\ix jx -> abs (col `atIndex` ix)
-                                        `compare`
-                                        abs (col `atIndex` jx))
-                             nonZeros
+findPivotMax :: Int -> Int -> Int -> Int -> STMatrix s Double -> ST s (Maybe (Int, Int))
+findPivotMax m n i j mat
+    | n == j = return Nothing
+    | m == i = return Nothing
+    | otherwise =
+        do
+          col      <- mapM (\i' -> do
+                                    x <- readMatrix mat i' j
+                                    return (i', abs x))
+                      [i..m-1]
+          let nonZeros = filter (not . (<eps) . abs . snd) col
+          let (pi, _) = maximumBy (\(_, x) (_, y) -> x `compare` y) nonZeros
+          case nonZeros of
+            [] -> if n == j+1
+                  then return Nothing
+                  else findPivotMax m n i (j+1) mat
+            _  -> return $ Just (pi, j)
 
 -- gaussian elimination of sub matrix below position (i, j)
 
-gaussian' :: Int -> Int -> STMatrix s Double -> ST s ()
-gaussian' i j mat = do
-    m       <- liftSTMatrix rows mat
-    n       <- liftSTMatrix cols mat
-    iPivot' <- liftSTMatrix (\x -> findPivotMax x i j) mat
+gaussian' :: Int -> Int -> Int -> Int -> STMatrix s Double -> ST s ()
+gaussian' m n i j mat = do
+    iPivot' <- findPivotMax m n i j mat
     case iPivot' of
         Nothing     -> return ()
         Just (r, p) -> do
                           rowOper (SWAP i r (FromCol j)) mat
-                          pv <- liftSTMatrix (`atIndex` (i, p)) mat
+                          pv <- readMatrix mat i p
                           mapM_ (reduce pv p) [i+1 .. m-1]
-                          gaussian' (i+1) (p+1) mat
+                          gaussian' m n (i+1) (p+1) mat
   where
     reduce pv p r = do
-                      rv <- liftSTMatrix (`atIndex` (r, p)) mat
-                      if eps rv
+                      rv <- readMatrix mat r p
+                      if abs rv < eps
                         then return ()
                         else
                          let frac = -rv / pv
                              op   = AXPY frac i r (FromCol p)
-                         in  rowOper op mat
+                         in do
+                             rowOper op mat
+                             mapM_ (\j' -> modifyMatrix mat r j' (\x -> if abs x < eps then 0 else x)) [p..n-1]
 
 -- | Gaussian elimination perfomed in-place in the @'ST'@ monad.
 
-gaussianST :: STMatrix s Double -> ST s ()
-gaussianST = gaussian' 0 0
+gaussianST :: Int -> Int -> STMatrix s Double -> ST s ()
+gaussianST m n = gaussian' m n 0 0
 
 
 -- | Gaussian elimination as pure function. Involves a copy of the input matrix.
@@ -144,9 +143,12 @@
 
 gaussian :: Matrix Double -> Matrix Double
 gaussian mat = runST $ do
-    m <- thawMatrix mat
-    gaussianST m
-    freezeMatrix m
+    matST <- thawMatrix mat
+    gaussianST m n matST
+    freezeMatrix matST
+  where
+    m = rows mat
+    n = cols mat
 
 -- | Returns the indices of a maximal linearly independent subset of the columns
 --   in the matrix.
diff --git a/stack.yaml b/stack.yaml
deleted file mode 100644
--- a/stack.yaml
+++ /dev/null
@@ -1,6 +0,0 @@
-resolver: lts-14.1
-
-build:
-  haddock-arguments:
-    haddock-args:
-    - --optghc=-fdefer-type-errors
diff --git a/test/LinearAlgebra.hs b/test/LinearAlgebra.hs
--- a/test/LinearAlgebra.hs
+++ b/test/LinearAlgebra.hs
@@ -13,7 +13,7 @@
 import Numeric.LinearAlgebra (rank)
 import qualified Numeric.LinearAlgebra.Data as Matrix
 
-import Math.Tensor.Internal.LinearAlgebra (independentColumnsMat)
+import Math.Tensor.Internal.LinearAlgebra (independentColumns, independentColumnsMat)
 
 data SmallInt = S0 | S1 deriving (Show, Ord, Eq, Enum, Bounded)
 
@@ -59,8 +59,16 @@
     mat  = (rows Matrix.>< cols) xs
     mat' = independentColumnsMat mat
 
+prop_consec :: Positive Int -> Int -> Bool
+prop_consec (Positive dim') start =
+    independentColumns mat == [0,1]
+  where
+    dim = dim' + 100
+    mat = (dim Matrix.>< dim) $ map fromIntegral [start..]
+
 testCase1 = testProperty "prop_smallValues" prop_smallValues
 testCase2 = testProperty "prop_ints" prop_ints
 testCase3 = testProperty "prop_doubles" prop_doubles
+testCase4 = testProperty "prop_consec" prop_consec
 
-linearAlgebraTest = testGroup "LinearAlgebraTest" [testCase1, testCase2, testCase3]
+linearAlgebraTest = testGroup "LinearAlgebraTest" [testCase1, testCase2, testCase3, testCase4]