diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -39,8 +39,10 @@
 
 ## Examples
 
-The module `Numeric.LinearAlgebra.Sparse` contains the interface functions:
+The module `Numeric.LinearAlgebra.Sparse` contains the user interface.
 
+### Creation and pretty-printing
+
 To create a sparse matrix from an array of its entries we use `fromListSM` :
 
     fromListSM :: Foldable t => (Int, Int) -> t (IxRow, IxCol, a) -> SpMatrix a
@@ -60,7 +62,10 @@
     [4,3,2]
     [0,0,5]
 
+The zeros are just added at pretty printing time; sparse vectors and matrices should only contain non-zero entries.
 
+### Matrix operations
+
 Matrix factorizations are available as `lu` and `qr` respectively, and are straightforward to verify by using the matrix product `##`  :
 
     > (l, u) = lu amat
@@ -70,6 +75,19 @@
     [2.0,0.0,0.0]
     [4.0,3.0,2.0]
     [0.0,0.0,5.0]
+
+Notice that the result is _dense_, i.e. certain entries are numerically zero but have been inserted into the result along with all the others (thus taking up memory!).
+To preserve sparsity, we can use a sparsifying matrix-matrix product `#~#`, which filters out all the elements x for which `|x| <= eps`, where `eps` (defined) in `Numeric.Eps`, is fixed at 10^-8.
+
+    > prd $ l #~# u
+    ( 3 rows, 3 columns ) , 5 NZ ( sparsity 0.5555555555555556 )
+
+    [2.0,0.0,0.0]
+    [4.0,3.0,2.0]
+    [0.0,0.0,5.0]
+
+
+### Linear systems
 
 Linear systems can be solved with either `linSolve` (which also requires choosing a method) or with `<\>` (which uses BiCGSTAB as default) :
 
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.2.0.7
+version:             0.2.0.8
 synopsis:            Numerical computation in native Haskell 
 description:         Currently the library provides iterative linear solvers, matrix decompositions, eigenvalue computations and related utilities. Please see README.md for details
 homepage:            https://github.com/ocramz/sparse-linear-algebra
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
@@ -2,20 +2,32 @@
 -- {-# OPTIONS_GHC -O2 -rtsopts -with-rtsopts=-K32m -prof#-}
 module Numeric.LinearAlgebra.Sparse
        (
-         sparsifySV,
+         -- * Matrix factorizations
+         qr, lu,
+         -- * Condition number
          conditionNumberSM,
+         -- * Householder reflection
          hhMat, hhRefl,
+         -- * Givens' rotation
          givens,
+         -- * Eigensolvers
          eigsQR, eigRayleigh,
+         -- * Linear solvers
+         linSolve, LinSolveMethod, (<\>),
+         -- ** Methods
          cgne, tfqmr, bicgstab, cgs, bcg,
          _xCgne, _xTfq, _xBicgstab, _x, _xBcg,
          cgsStep, bicgstabStep,
          CGNE, TFQMR, BICGSTAB, CGS, BCG,
-         linSolve, LinSolveMethod, (<\>),
-         randMat, randVec, randSpMat, randSpVec,
+         -- * Random matrices and vectors
+         randMat, randVec, 
+         -- ** Sparse "
+         randSpMat, randSpVec,
+         -- * Sparsify data
+         sparsifySV,
+         -- * Iteration combinators
          modifyInspectN, runAppendN',
-         diffSqL,
-         qr, lu
+         diffSqL
        )
        where
 
@@ -353,11 +365,15 @@
 -- solve for element Lij
 solveForLij ::
   SpMatrix Double -> SpMatrix Double -> SpMatrix Double -> IxRow -> IxCol -> Double
-solveForLij amat lmat umat i j | isNz uii = (a - p)/uii
-                               | otherwise = error $ unwords ["solveForLij : U",show (i,i)," is close to 0. Permute rows in order to have a nonzero diagonal of U"]
+solveForLij amat lmat umat i j
+  | isNz ujj = (a - p)/ujj
+  | otherwise =
+     error $ unwords ["solveForLij : U",
+                      show (j ,j ),
+                      "is close to 0. Permute rows in order to have a nonzero diagonal of U"]
   where
    a = amat @@! (i, j)
-   uii = umat @@! (i-1, i-1) -- NB this must be /= 0
+   ujj = umat @@! (j , j)   -- NB this must be /= 0
    p = contractSub lmat umat i j (i - 1)
 
 
@@ -389,20 +405,6 @@
 
 
 
-
-
-
-
-
-
-
-
-
-
-
--- test data
-
-aa0 = fromListDenseSM 3 [2, -4, -4, -1, 6, -2, -2, 3, 8] :: SpMatrix Double
 
 
 
diff --git a/test/LibSpec.hs b/test/LibSpec.hs
--- a/test/LibSpec.hs
+++ b/test/LibSpec.hs
@@ -76,9 +76,10 @@
     it "QR (3 x 3 dense)" $ 
       checkQr tm2 `shouldBe` True
   describe "Numeric.LinearAlgebra.Sparse : LU decomposition" $ do
-    it "LU (3 x 3 dense)" $
-      checkLu tm2 `shouldBe` True
-
+    it "LU (4 x 4 dense)" $
+      checkLu tm6 `shouldBe` True
+    it "LU (10 x 10 sparse)" $
+      checkLu tm7 `shouldBe` True
 
 
 {-
@@ -248,9 +249,10 @@
 
 {- LU -}
 
+checkLu :: SpMatrix Double -> Bool
 checkLu a = lup == a where
   (l, u) = lu a
-  lup = l ## u
+  lup = l #~# u
 
 
 
@@ -318,10 +320,17 @@
 tm4 = sparsifySM $ fromListDenseSM 4 [1,0,0,0,2,5,0,10,3,6,8,11,4,7,9,12]
 
 
-
+tm5 = fromListDenseSM 3 [2, -4, -4, -1, 6, -2, -2, 3, 8] :: SpMatrix Double
 
 
+tm6 = fromListDenseSM 4 [1,3,4,2,2,5,2,10,3,6,8,11,4,7,9,12] :: SpMatrix Double
 
+tm7 :: SpMatrix Double
+tm7 = a ^+^ b ^+^ c where
+  n = 5
+  a = mkSubDiagonal n 1 $ replicate n 1
+  b = mkSubDiagonal n 0 $ replicate n (-2)
+  c = mkSubDiagonal n (-1) $ replicate n 1
 
 -- -- run N iterations