diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -18,11 +18,13 @@
 
     * Transpose-Free Quasi-Minimal Residual (TFQMR)
 
-* Matrix decompositions
+* Matrix factorization algorithms
 
-    * QR factorization
+    * QR
 
-    * LU factorization
+    * LU
+
+    * Cholesky
 
 * Eigenvalue algorithms
 
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.9
+version:             0.2.1.0
 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/Data/Sparse/Common.hs b/src/Data/Sparse/Common.hs
--- a/src/Data/Sparse/Common.hs
+++ b/src/Data/Sparse/Common.hs
@@ -1,6 +1,16 @@
+-----------------------------------------------------------------------------
+-- |
+-- Copyright   :  (C) 2016 Marco Zocca
+-- License     :  GPL-3 (see LICENSE)
+-- Maintainer  :  zocca.marco gmail
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-----------------------------------------------------------------------------
 module Data.Sparse.Common
        ( module X,
          insertRowWith, insertRow, insertColWith, insertCol,
+         diagonalSM,
          outerProdSV, (><), toSV, svToSM, 
          extractCol, extractRow,
          extractVectorDenseWith, extractRowDense, extractColDense,
@@ -58,6 +68,7 @@
 
 -- * Outer vector product
 
+-- | Outer product (all-with-all matrix)
 outerProdSV, (><) :: Num a => SpVector a -> SpVector a -> SpMatrix a
 outerProdSV v1 v2 = fromListSM (m, n) ixy where
   m = dim v1
@@ -65,6 +76,16 @@
   ixy = [(i,j, x * y) | (i,x) <- toListSV v1 , (j, y) <- toListSV v2]
 
 (><) = outerProdSV
+
+
+
+-- * Diagonal matrix
+
+-- | Fill the diagonal of a SpMatrix with the components of a SpVector
+diagonalSM :: SpVector a -> SpMatrix a
+diagonalSM sv = ifoldSV iins (zeroSM n n) sv where
+  n = dim sv
+  iins i = insertSpMatrix i i
 
 
 
diff --git a/src/Data/Sparse/IntMap2/IntMap2.hs b/src/Data/Sparse/IntMap2/IntMap2.hs
--- a/src/Data/Sparse/IntMap2/IntMap2.hs
+++ b/src/Data/Sparse/IntMap2/IntMap2.hs
@@ -35,10 +35,9 @@
 
 
 
--- * set-like brackets
+-- set-like brackets
 
 
--- unionWithKeyIM2 f im1 im2 = undefined where
 
 
 
@@ -110,7 +109,7 @@
 
 
 
--- * filtering
+-- * Filtering
 
 -- |Map over outer IM and filter all inner IM's
 ifilterIM2 ::
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
@@ -1,4 +1,13 @@
 {-# language FlexibleInstances, MultiParamTypeClasses, TypeFamilies #-}
+-----------------------------------------------------------------------------
+-- |
+-- Copyright   :  (C) 2016 Marco Zocca
+-- License     :  GPL-3 (see LICENSE)
+-- Maintainer  :  zocca.marco gmail
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-----------------------------------------------------------------------------
 module Data.Sparse.SpMatrix where
 
 import Data.Sparse.Utils
@@ -15,6 +24,20 @@
 
 
 
+-- *
+
+-- instance IxContainer SpMatrix a where
+--   type Ix SpMatrix = (Int, Int)
+--   type IxSz SpMatrix = (Int, Int)
+--   ixcLookup m (i,j) = lookupSM m i j
+--   ixcIfilter g = ifilterSM g' where g' i j = g (i, j)
+--   ixcInsert (i, j) = insertSpMatrix i j
+--   -- ixcFromList = fromListSM
+
+-- instance SMatrix SpMatrix a where
+
+
+
 -- * Sparse Matrix
 
 data SpMatrix a = SM {smDim :: (Rows, Cols),
@@ -88,7 +111,7 @@
 permutationSM n iis = permutPairsSM n (zip [0 .. n-1] iis)
 
 -- | Permutation matrix from a (possibly incomplete) list of row pair swaps
--- e.g. `permutPairs 5 [(2,4)]` swaps rows (2, 4) :
+-- e.g. `permutPairs 5 [(2,4)]` swaps rows 2 and 4 :
 --
 -- [1,0,0,0,0]
 -- [0,1,0,0,0]
@@ -159,7 +182,7 @@
 
 -- ** toList
 
--- |Populate list with SpMatrix contents and populate missing entries with 0
+-- | Populate list with SpMatrix contents and populate missing entries with 0
 toDenseListSM :: Num t => SpMatrix t -> [(IxRow, IxCol, t)]
 toDenseListSM m =
   [(i, j, m @@ (i, j)) | i <- [0 .. nrows m - 1], j <- [0 .. ncols m- 1]]
@@ -263,6 +286,7 @@
 
 
 -- *** Extract i'th row
+-- | Extract whole row
 extractRowSM :: SpMatrix a -> IxRow -> SpMatrix a
 extractRowSM sm i = extractSubmatrix sm (i, i) (0, ncols sm - 1)
 
@@ -279,7 +303,7 @@
 
 
 -- *** Extract j'th column
--- | Extract all column
+-- | Extract whole column
 extractColSM :: SpMatrix a -> IxCol -> SpMatrix a
 extractColSM sm j = extractSubmatrix sm (0, nrows sm - 1) (j, j)
 
@@ -323,8 +347,8 @@
   ff irow row = IM.size row == 1 &&
                 IM.size (IM.filterWithKey (\j _ -> j == irow) row) == 1
 
--- |is the matrix orthogonal? i.e. Q^t ## Q == I
-isOrthogonalSM :: SpMatrix Double -> Bool
+-- |Is the matrix orthogonal? i.e. Q^t ## Q == I
+isOrthogonalSM :: (Eq a, Epsilon a) => SpMatrix a -> Bool
 isOrthogonalSM sm@(SM (_,n) _) = rsm == eye n where
   rsm = roundZeroOneSM $ transposeSM sm ## sm
 
@@ -497,6 +521,11 @@
 
 -- ** Misc. SpMatrix operations
 
+ifilterSM :: (IM.Key -> IM.Key -> a -> Bool) -> SpMatrix a -> SpMatrix a
+ifilterSM f (SM d im) = SM d $ ifilterIM2 f im
+
+ 
+
 -- | Left fold over SpMatrix
 foldlSM :: (a -> b -> b) -> b -> SpMatrix a -> b
 foldlSM f n (SM _ m)= foldlIM2 f n m
@@ -536,11 +565,12 @@
 
 
 -- ** Sparsify : remove almost-0 elements (|x| < eps)
-sparsifyIM2 :: IM.IntMap (IM.IntMap Double) -> IM.IntMap (IM.IntMap Double)
-sparsifyIM2 = ifilterIM2 (\_ _ x -> abs x >= eps)
+sparsifyIM2 ::
+  Epsilon a => IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)
+sparsifyIM2 = ifilterIM2 (\_ _ x -> isNz x)
 
 -- | Sparsify an SpMatrix
-sparsifySM :: SpMatrix Double -> SpMatrix Double
+sparsifySM :: Epsilon a => SpMatrix a -> SpMatrix a
 sparsifySM (SM d im) = SM d $ sparsifyIM2 im
 
 
@@ -548,7 +578,7 @@
 
 -- ** Value rounding
 -- | Round almost-0 and almost-1 to 0 and 1 respectively
-roundZeroOneSM :: SpMatrix Double -> SpMatrix Double
+roundZeroOneSM :: Epsilon a => SpMatrix a -> SpMatrix a
 roundZeroOneSM (SM d im) = sparsifySM $ SM d $ mapIM2 roundZeroOne im  
 
 
@@ -560,14 +590,14 @@
 
 
 -- ** Matrix row swap
--- | swap two rows of a SpMatrix (bounds not checked)
+-- | Swap two rows of a SpMatrix (bounds not checked)
 swapRows :: IxRow -> IxRow -> SpMatrix a -> SpMatrix a
 swapRows i1 i2 (SM d im) = SM d $ IM.insert i1 ro2 im' where
   ro1 = im IM.! i1
   ro2 = im IM.! i2
   im' = IM.insert i2 ro1 im
 
--- | swap two rows of a SpMatrix (bounds checked)  
+-- | Swap two rows of a SpMatrix (bounds checked)  
 swapRowsSafe :: IxRow -> IxRow -> SpMatrix a -> SpMatrix a
 swapRowsSafe i1 i2 m
   | inBounds02 (nro, nro) (i1, i2) = swapRows i1 i2 m
@@ -581,11 +611,10 @@
 
 
 -- ** Matrix transpose
--- | transposeSM, (#^) : Matrix transpose
-transposeSM, (#^) :: SpMatrix a -> SpMatrix a
+-- | transposeSM : Matrix transpose
+transposeSM :: SpMatrix a -> SpMatrix a
 transposeSM (SM (m, n) im) = SM (n, m) (transposeIM2 im)
 
-(#^) = transposeSM
 
 
 
@@ -603,7 +632,7 @@
 
 
 -- ** Frobenius norm
-normFrobenius :: SpMatrix Double -> Double
+normFrobenius :: Floating a => SpMatrix a -> a
 normFrobenius m = sqrt $ foldlSM (+) 0 m' where
   m' | nrows m > ncols m = transposeSM m ## m
      | otherwise = m ## transposeSM m 
@@ -655,7 +684,7 @@
 
 -- ** Matrix-matrix product, sparsified
 -- | Removes all elements `x` for which `| x | <= eps`)
-matMatSparsified, (#~#)  :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double
+matMatSparsified, (#~#) :: Epsilon a => SpMatrix a -> SpMatrix a -> SpMatrix a
 matMatSparsified m1 m2 = sparsifySM $ matMat m1 m2
 
 (#~#) = matMatSparsified
@@ -666,12 +695,12 @@
 -- *** Sparsified matrix products of two matrices
 
 -- | A^T B
-(#^#) :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double
+(#^#) :: Epsilon a => SpMatrix a -> SpMatrix a -> SpMatrix a
 a #^# b = transposeSM a #~# b
 
 
 -- | A B^T
-(##^) :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double
+(##^) :: Epsilon a => SpMatrix a -> SpMatrix a -> SpMatrix a
 a ##^ b = a #~# transposeSM b
 
 
@@ -686,8 +715,8 @@
 
 
 
--- *** Matrix contraction
--- | Contract two matrices A and B up to an index `n`, i.e. summing over repeated indices: 
+-- ** 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 :: Num a => SpMatrix a -> SpMatrix a -> IxRow -> IxCol -> Int -> a
 contractSub a b i j n
diff --git a/src/Data/Sparse/SpVector.hs b/src/Data/Sparse/SpVector.hs
--- a/src/Data/Sparse/SpVector.hs
+++ b/src/Data/Sparse/SpVector.hs
@@ -1,12 +1,23 @@
 {-# language TypeFamilies, MultiParamTypeClasses, FlexibleInstances #-}
+-----------------------------------------------------------------------------
+-- |
+-- Copyright   :  (C) 2016 Marco Zocca
+-- License     :  GPL-3 (see LICENSE)
+-- Maintainer  :  zocca.marco gmail
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-----------------------------------------------------------------------------
 module Data.Sparse.SpVector where
 
 import Data.Sparse.Utils
 import Data.Sparse.Types
+import Data.Sparse.IntMap2.IntMap2
 
+import Numeric.Eps
 import Numeric.LinearAlgebra.Class
-import Data.Sparse.IntMap2.IntMap2
 
+
 import Data.Maybe
 
 import qualified Data.IntMap as IM
@@ -94,11 +105,11 @@
 
 
 -- | create a sparse vector from an association list while discarding all zero entries
-mkSpVector :: (Num a, Eq a) => Int -> IM.IntMap a -> SpVector a
-mkSpVector d im = SV d $ IM.filterWithKey (\k v -> v /= 0 && inBounds0 d k) im
+mkSpVector :: Epsilon a => Int -> IM.IntMap a -> SpVector a
+mkSpVector d im = SV d $ IM.filterWithKey (\k v -> isNz v && inBounds0 d k) im
 
 -- | ", from logically dense array (consecutive indices)
-mkSpVectorD :: (Num a, Eq a) => Int -> [a] -> SpVector a
+mkSpVectorD :: Epsilon a => Int -> [a] -> SpVector a
 mkSpVectorD d ll = mkSpVector d (IM.fromList $ denseIxArray (take d ll))
 
 -- ", don't filter zero elements
@@ -110,6 +121,18 @@
 fromListDenseSV d ll = SV d (IM.fromList $ denseIxArray (take d ll))
 
 
+-- | Map a function over a range of indices and filter the result (indices and values) to fit in a `n`-long SpVector
+spVectorDenseIx :: Epsilon a => (Int -> a) -> UB -> [Int] -> SpVector a
+spVectorDenseIx f n ix =
+  fromListSV n $ filter q $ zip ix $ map f ix where
+    q (i, v) = inBounds0 n i && isNz v
+
+-- | ", using just the integer bounds of the interval
+spVectorDenseLoHi :: Epsilon a => (Int -> a) -> UB -> Int -> Int -> SpVector a
+spVectorDenseLoHi f n lo hi = spVectorDenseIx f n [lo .. hi]    
+
+
+
 -- | one-hot encoding : `oneHotSV n k` produces a SpVector of length n having 1 at the k-th position
 oneHotSVU :: Num a => Int -> IxRow -> SpVector a
 oneHotSVU n k = SV n (IM.singleton k 1)
@@ -152,6 +175,11 @@
 toDenseListSV (SV d im) = fmap (\i -> IM.findWithDefault 0 i im) [0 .. d-1]
 
 
+
+
+-- | Indexed fold over SpVector
+ifoldSV :: (IM.Key -> a -> b -> b) -> b -> SpVector a -> b
+ifoldSV f e (SV d im) = IM.foldWithKey f e im
 
 
 
diff --git a/src/Numeric/Eps.hs b/src/Numeric/Eps.hs
--- a/src/Numeric/Eps.hs
+++ b/src/Numeric/Eps.hs
@@ -1,36 +1,75 @@
-module Numeric.Eps where
+-----------------------------------------------------------------------------
+-- |
+-- Copyright   :  (C) 2016 Marco Zocca, 2012-2015 Edward Kmett
+-- License     :  GPL-3 (see LICENSE)
+-- Maintainer  :  zocca.marco gmail
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-- Testing for values "near" zero
+-----------------------------------------------------------------------------
+module Numeric.Eps
+  ( Epsilon(..), nearZero, isNz, roundZero, roundOne, roundZeroOne
+  ) where
+import Foreign.C.Types (CFloat, CDouble)
 
--- * Numerical tolerance for "near-0" tests
--- | eps = 1e-8 
-eps :: Double
-eps = 1e-8
+-- | Provides a test to see if a quantity is near zero.
+--
+-- >>> nearZero (1e-11 :: Double)
+-- False
+--
+-- >>> nearZero (1e-17 :: Double)
+-- True
+--
+-- >>> nearZero (1e-5 :: Float)
+-- False
+--
+-- >>> nearZero (1e-7 :: Float)
+-- True
+class Num a => Epsilon a where
+  -- | Determine if a quantity is near zero.
+  nearZero :: a -> Bool
 
+-- | @'abs' a '<=' 1e-6@
+instance Epsilon Float where
+  nearZero a = abs a <= 1e-6
 
+-- | @'abs' a '<=' 1e-12@
+instance Epsilon Double where
+  nearZero a = abs a <= 1e-12
 
+-- | @'abs' a '<=' 1e-6@
+instance Epsilon CFloat where
+  nearZero a = abs a <= 1e-6
 
+-- | @'abs' a '<=' 1e-12@
+instance Epsilon CDouble where
+  nearZero a = abs a <= 1e-12
+
+
+
+
 -- * Rounding operations
 
--- | Rounding rule
-almostZero, almostOne :: Double -> Bool
-almostZero x = abs x <= eps
-almostOne x = x >= (1-eps) && x < (1+eps)
 
-isNz :: Double -> Bool
-isNz = not . almostZero
+-- | Rounding rule
+almostZero, almostOne, isNz :: Epsilon a => a -> Bool
+almostZero = nearZero
+almostOne x = nearZero (1 - x)
+isNz x = not (almostZero x)
 
 withDefault :: (t -> Bool) -> t -> t -> t
 withDefault q d x | q x = d
                   | otherwise = x
 
-roundZero, roundOne :: Double -> Double
-roundZero = withDefault almostZero 0
-roundOne = withDefault almostOne 1
+roundZero, roundOne, roundZeroOne :: Epsilon a => a -> a
+roundZero = withDefault almostZero (fromIntegral 0)
+roundOne = withDefault almostOne (fromIntegral 1)
 
 with2Defaults :: (t -> Bool) -> (t -> Bool) -> t -> t -> t -> t
 with2Defaults q1 q2 d1 d2 x | q1 x = d1
                             | q2 x = d2
                             | otherwise = x
 
--- | Round to respectively 0 or 1 within some predefined numerical precision eps
-roundZeroOne :: Double -> Double
-roundZeroOne = with2Defaults almostZero almostOne 0 1
+-- | Round to respectively 0 or 1
+roundZeroOne = with2Defaults almostZero almostOne (fromIntegral 0) (fromIntegral 1)
diff --git a/src/Numeric/LinearAlgebra/Class.hs b/src/Numeric/LinearAlgebra/Class.hs
--- a/src/Numeric/LinearAlgebra/Class.hs
+++ b/src/Numeric/LinearAlgebra/Class.hs
@@ -10,8 +10,12 @@
   (^+^) :: Num a => f a -> f a -> f a
 
 
+  one :: Num a => f a
 
+  (^*^) :: Num a => f a -> f a -> f a
 
+
+
 -- | negate the values in a functor
 negated :: (Num a, Functor f) => f a -> f a
 negated = fmap negate
@@ -185,15 +189,14 @@
 
 -- * IxContainer : indexed container types
 
-class IxContainer (c :: * -> *) a where
-  type Ix c :: *
-  ixcLookup :: Ix c -> c a -> Maybe a
-  ixcLookupDefault :: a -> Ix c -> c a -> a
-  ixcFilter :: (a -> Bool) -> c a -> c a
-  ixcIfilter :: (Ix c -> a -> Bool) -> c a -> c a
-  ixcInsert :: Ix c -> a -> c a -> c a
-  ixcFromList :: [(Ix c, a)] -> c a
-  ixcToList :: c a -> [(Ix c, a)]
+-- class IxContainer (c :: * -> *) a where
+--   type Ix c :: *
+--   type IxSz c :: *
+--   ixcLookup :: c a -> Ix c -> Maybe a
+--   ixcIfilter :: (Ix c -> a -> Bool) -> c a -> c a
+--   ixcInsert :: Ix c -> a -> c a -> c a
+--   ixcFromList :: Foldable t => IxSz c -> t (Ix c, a) -> c a
+--   ixcToList :: c a -> [(Ix c, a)]
 
 -- newtype IM_ a = IM (IM.IntMap a)
 
@@ -208,3 +211,23 @@
 -- instance IxContainer IM2 a where
 --   type Ix IM2 = (Int, Int)
 --   ixcIfilter f im2 = IM2 $ ifilterIM2 (curry f) (unIM2 im2)
+
+
+
+-- class Rank2 (c :: * -> *) a where
+--   type R2IxRow c :: *
+--   type R2IxCol c :: *
+--   type R2V c :: * -> *
+--   rank2Act :: c a -> R2V c a -> R2V c a
+--   extractRow :: c a -> R2IxRow c -> Maybe (R2V c a)
+--   extractCol :: c a -> R2IxCol c -> Maybe (R2V c a)
+--   -- extractRows :: c a -> [R2IxRow c] -> [R2V c]
+--   -- extractCols :: c a -> [R2IxCol c] -> [R2V c]
+  
+
+
+
+
+-- * SMatrix : sparse matrix types
+
+-- class (IxContainer c a, Sparse c a, Additive c) => SMatrix c a where
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
@@ -4,6 +4,7 @@
        (
          -- * Matrix factorizations
          qr, lu,
+         chol,
          -- * Incomplete LU
          ilu0,
          -- * Condition number
@@ -21,6 +22,8 @@
          _xCgne, _xTfq, _xBicgstab, _x, _xBcg,
          cgsStep, bicgstabStep,
          CGNE, TFQMR, BICGSTAB, CGS, BCG,
+         -- * Preconditioners
+         ilu0, mSsor,
          -- * Matrix partitioning
          diagPartitions,
          -- * Random arrays
@@ -69,8 +72,8 @@
   
 -- * Sparsify : remove almost-0 elements (|x| < eps)
 -- | Sparsify an SpVector
-sparsifySV :: SpVector Double -> SpVector Double
-sparsifySV (SV d im) = SV d $ IM.filter (\x -> abs x >= eps) im
+sparsifySV :: Epsilon a => SpVector a -> SpVector a
+sparsifySV (SV d im) = SV d $ IM.filter isNz im
 
 
 
@@ -79,7 +82,7 @@
 -- * Matrix condition number
 
 -- |uses the R matrix from the QR factorization
-conditionNumberSM :: SpMatrix Double -> Double
+conditionNumberSM :: (Epsilon a, RealFloat a) => SpMatrix a -> a
 conditionNumberSM m | isInfinite kappa = error "Infinite condition number : rank-deficient system"
                     | otherwise = kappa where
   kappa = lmax / lmin
@@ -104,8 +107,8 @@
 {-| a vector `x` uniquely defines an orthogonal plane; the Householder operator reflects any point `v` with respect to this plane:
  v' = (I - 2 x >< x) v
 -}
-hhRefl :: SpVector Double -> SpMatrix Double
-hhRefl = hhMat 2.0
+hhRefl :: Num a => SpVector a -> SpMatrix a
+hhRefl = hhMat (fromInteger 2)
 
 
 
@@ -151,7 +154,7 @@
 non-zero but A has zeros in row k for all columns less than j.
 -}
 
-givens :: SpMatrix Double -> IxRow -> IxCol -> SpMatrix Double
+givens :: (Floating a, Epsilon a, Ord a) => SpMatrix a -> IxRow -> IxCol -> SpMatrix a
 givens mm i j 
   | isValidIxSM mm (i,j) && isSquareSM mm =
        sparsifySM $ fromListSM' [(i,i,c),(j,j,c),(j,i,-s),(i,j,s)] (eye (nrows mm))
@@ -182,15 +185,15 @@
 -- * QR decomposition
 
 
--- | Applies Givens rotation iteratively to zero out sub-diagonal elements
-qr :: SpMatrix Double -> (SpMatrix Double, SpMatrix Double)
+-- | Given a matrix A, returns a pair of matrices (Q, R) such that Q R = A, Q is orthogonal and R is upper triangular. Applies Givens rotation iteratively to zero out sub-diagonal elements
+qr :: (Epsilon a, Floating a, Real a) => SpMatrix a -> (SpMatrix a, SpMatrix a)
 qr mm = (transposeSM qmatt, rmat)  where
   qmatt = F.foldl' (#~#) ee $ gmats mm -- Q^T = (G_n * G_n-1 ... * G_1)
   rmat = qmatt #~# mm                  -- R = Q^T A
   ee = eye (nrows mm)
       
 -- | Givens matrices in order [G1, G2, .. , G_N ]
-gmats :: SpMatrix Double -> [SpMatrix Double]
+gmats :: (Epsilon a, Real a, Floating a) => SpMatrix a -> [SpMatrix a]
 gmats mm = gm mm (subdiagIndicesSM mm) where
  gm m ((i,j):is) = let g = givens m i j
                    in g : gm (g #~# m) is
@@ -218,13 +221,13 @@
 -- ** QR algorithm
 
 -- | `eigsQR n mm` performs `n` iterations of the QR algorithm on matrix `mm`, and returns a SpVector containing all eigenvalues
-eigsQR :: Int -> SpMatrix Double -> SpVector Double
+eigsQR :: (Epsilon a, Real a, Floating a) => Int -> SpMatrix a -> SpVector a
 eigsQR nitermax m = extractDiagDense $ execState (convergtest eigsStep) m where
   eigsStep m = r #~# q where (q, r) = qr m
   convergtest g = modifyInspectN nitermax f g where
     f [m1, m2] = let dm1 = extractDiagDense m1
                      dm2 = extractDiagDense m2
-                 in norm2 (dm1 ^-^ dm2) <= eps
+                 in nearZero $ norm2 (dm1 ^-^ dm2)
 
 
 
@@ -234,13 +237,13 @@
 -- ** Rayleigh iteration
 
 -- | `eigsRayleigh n mm` performs `n` iterations of the Rayleigh algorithm on matrix `mm` and returns the eigenpair closest to the initialization. It displays cubic-order convergence, but it also requires an educated guess on the initial eigenpair
-eigRayleigh :: Int                -- max # iterations
-     -> SpMatrix Double           -- matrix
-     -> (SpVector Double, Double) -- initial guess of (eigenvector, eigenvalue)
-     -> (SpVector Double, Double) -- final estimate of (eigenvector, eigenvalue)
+-- eigRayleigh :: Int                -- max # iterations
+--      -> SpMatrix Double           -- matrix
+--      -> (SpVector Double, Double) -- initial guess of (eigenvector, eigenvalue)
+--      -> (SpVector Double, Double) -- final estimate of (eigenvector, eigenvalue)
 eigRayleigh nitermax m = execState (convergtest (rayleighStep m)) where
   convergtest g = modifyInspectN nitermax f g where
-    f [(b1, _), (b2, _)] = norm2 (b2 ^-^ b1) <= eps 
+    f [(b1, _), (b2, _)] = nearZero $ norm2 (b2 ^-^ b1)
   rayleighStep aa (b, mu) = (b', mu') where
       ii = eye (nrows aa)
       nom = (aa ^-^ (mu `matScale` ii)) <\> b
@@ -253,13 +256,13 @@
 -- * Householder vector 
 
 -- (Golub & Van Loan, Alg. 5.1.1, function `house`)
-hhV :: SpVector Double -> (SpVector Double, Double)
+hhV :: (Epsilon a, Real a, Floating a) => SpVector a -> (SpVector a, a)
 hhV x = (v, beta) where
   n = dim x
   tx = tailSV x
   sigma = tx `dot` tx
   vtemp = singletonSV 1 `concatSV` tx
-  (v, beta) | sigma <= eps = (vtemp, 0)
+  (v, beta) | nearZero sigma = (vtemp, 0)
             | otherwise = let mu = sqrt (headSV x**2 + sigma)
                               xh = headSV x
                               vh | xh <= 1 = xh - mu
@@ -298,109 +301,92 @@
 
 
 
+-- * Cholesky factorization
 
+-- ** Cholesky–Banachiewicz algorithm
 
+-- | Given a positive semidefinite matrix A, returns a lower-triangular matrix L such that L L^T = A
+chol :: (Epsilon a, Real a, Floating a) => SpMatrix a -> SpMatrix a
+chol aa = lfin where
+  (_, lfin) = execState (modifyUntil q cholUpd) cholInit
+  q (i, _) = i == nrows aa              -- stopping criterion
+  cholInit = cholUpd (0, zeroSM n n)    -- initialization
+  n = nrows aa
+  cholUpd (i, ll) = (i + 1, ll') where
+    ll' = cholDiagUpd (cholSDRowUpd ll) -- first upd subdiagonal entries in the row
+    cholSDRowUpd ll_ = insertRow ll_ lrs i where
+       lrs = fromListSV (i + 1) $ onRangeSparse (cholSubDiag ll i) [0 .. i-1]
+    cholDiagUpd ll_ = insertSpMatrix i i (cholDiag ll_ i) ll_ 
+  cholSubDiag ll i j = 1/ljj*(aij - inn) where
+    ljj = ll@@(j, j)
+    aij = aa@@(i, j)
+    inn = contractSub ll ll i j (j - 1)
+  cholDiag ll i | i == 0 = sqrt aai
+                | i > 0 = sqrt $ aai - sum (fmap (**2) lrow)
+                | otherwise = error "cholDiag : index must be nonnegative" where
+        aai = aa@@(i,i)
+        lrow = ifilterSV (\j _ -> j < i) (extractRow ll i) -- sub-diagonal elems of L
 
 
 
 
 
--- * LU factorization
--- ** Doolittle algorithm
-{- Doolittle algorithm for factoring A' = P A, where P is a permutation matrix such that A' has a nonzero as its (0, 0) entry -}
 
--- | LU factors
-lu :: SpMatrix Double -> (SpMatrix Double, SpMatrix Double)
-lu aa = (lfin, ufin) where
-  (ixf,lf,uf) = execState (modifyUntil q (luUpd aa)) (luInit aa)
-  lfin = lf
-  ufin = uUpd aa (ixf, lf, uf)
-  q (i, _, _) = i == (nrows aa - 1)
 
--- | First iteration of LU
-luInit ::
-  (Num t, Fractional a) => SpMatrix a -> (t, SpMatrix a, SpMatrix a)
-luInit aa = (1, l0, u0) where
-  n = nrows aa
-  l0 = insertCol (eye n) ((1/u00) .* extractSubCol aa 0 (1,n - 1)) 0  -- initial L
-  u0 = insertRow (zeroSM n n) (extractRow aa 0) 0               -- initial U
-  u00 = u0 @@ (0,0)  -- make sure this is non-zero by applying permutation
 
--- | LU update step
-luUpd :: SpMatrix Double
-     -> (Int, SpMatrix Double, SpMatrix Double)
-     -> (Int, SpMatrix Double, SpMatrix Double)
-luUpd aa (i, l, u) = (i', l', u') where
-  n = nrows aa  
-  u' = uUpdSparse aa (i, l, u)  -- update U
-  l' = lUpdSparse aa (i, l, u') -- update L
-  i' = i + 1     -- increment i
 
 
-uUpd' ::
-  Num a =>
-  ([(Int, a)] -> [(Int, a)]) ->
-  SpMatrix a ->
-  (Rows, SpMatrix a, SpMatrix a) ->
-  SpMatrix a
-uUpd' ff amat (ix, lmat, umat) = insertRow umat uv ix where
-  n = nrows amat
-  colsix = [ix .. n - 1]
-  us = ff $ zip colsix $ map (solveForUij amat lmat umat ix) colsix
-  uv = fromListSV n us
 
-uUpd :: Num a => SpMatrix a -> (Rows, SpMatrix a, SpMatrix a) -> SpMatrix a
-uUpd = uUpd' id
 
--- update U while sparsifying
-uUpdSparse ::
-  SpMatrix Double -> (Rows, SpMatrix Double, SpMatrix Double) -> SpMatrix Double
-uUpdSparse = uUpd' (filter (isNz . snd))
 
 
 
 
--- solve for element Uij
-solveForUij ::
-  Num a => SpMatrix a -> SpMatrix a -> SpMatrix a -> IxRow -> IxCol -> a
-solveForUij amat lmat umat i j = a - p where
-  a = amat @@! (i, j)
-  p = contractSub lmat umat i j (i - 1)
 
+-- * LU factorization
+-- ** Doolittle algorithm
+{- Doolittle algorithm for factoring A' = P A, where P is a permutation matrix such that A' has a nonzero as its (0, 0) entry -}
 
--- solve for element Lij
-solveForLij ::
-  SpMatrix Double -> SpMatrix Double -> SpMatrix Double -> IxRow -> IxCol -> Double
-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)
-   ujj = umat @@! (j , j)   -- NB this must be /= 0
-   p = contractSub lmat umat i j (i - 1)
+-- | Given a matrix A, returns a pair of matrices (L, U) such that L U = A
+lu :: (Epsilon a, Fractional a, Real a) => SpMatrix a -> (SpMatrix a, SpMatrix a)
+lu aa = (lf, ufin) where
+  (ixf, lf, uf) = execState (modifyUntil q luUpd) luInit
+  ufin = uUpdSparse (ixf, lf, uf) -- final U update
+  q (i, _, _) = i == (nrows aa - 1)
+  n = nrows aa
+  luInit = (1, l0, u0) where
+    l0 = insertCol (eye n) ((1/u00) .* extractSubCol aa 0 (1,n - 1)) 0  -- initial L
+    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) = (i + 1, l', u') where
+    u' = uUpdSparse (i, l, u)  -- update U
+    l' = lUpdSparse (i, l, u') -- update L
+  uUpdSparse (ix, lmat, umat) = insertRow umat (fromListSV n us) ix where
+    us = onRangeSparse (solveForUij ix) [ix .. n - 1]
+    solveForUij i j = a - p where
+      a = aa @@! (i, j)
+      p = contractSub lmat umat i j (i - 1)
+  lUpdSparse (ix, lmat, umat) = insertCol lmat (fromListSV n ls) ix where
+    ls = onRangeSparse (`solveForLij` ix) [ix + 1 .. n - 1]
+    solveForLij 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 = aa @@! (i, j)
+       ujj = umat @@! (j , j)   -- NB this must be /= 0
+       p = contractSub lmat umat i j (i - 1)
 
 
 
+-- | Apply a function over a range of integer indices, zip the result with it and filter out the almost-zero entries
+onRangeSparse :: (Epsilon b, Real b) => (Int -> b) -> [Int] -> [(Int, b)]
+onRangeSparse f ixs = filter (isNz . snd) $ zip ixs $ map f ixs
 
-lUpd' :: ([(Rows, Double)] -> [(Int, Double)])
-     -> SpMatrix Double
-     -> (Rows, SpMatrix Double, SpMatrix Double)
-     -> SpMatrix Double
-lUpd' ff amat (ix, lmat, umat) = insertCol lmat lv ix where
-  n = nrows amat
-  rowsix = [ix + 1 .. n - 1]
-  ls = ff $ zip rowsix $ map (\i -> solveForLij amat lmat umat i ix) rowsix
-  lv = fromListSV n ls
 
-lUpd :: SpMatrix Double -> (Rows, SpMatrix Double, SpMatrix Double) -> SpMatrix Double
-lUpd = lUpd' id
 
-lUpdSparse ::
-  SpMatrix Double -> (Rows, SpMatrix Double, SpMatrix Double) -> SpMatrix Double
-lUpdSparse = lUpd' (filter (isNz . snd))
 
 
 
@@ -442,7 +428,7 @@
 -- | used for Incomplete LU : remove entries in `m` corresponding to zero entries in `m2`
 
 
-
+ilu0 :: (Epsilon a, Real a, Fractional a) => SpMatrix a -> (SpMatrix a, SpMatrix a)
 ilu0 aa = (lh, uh) where
   (l, u) = lu aa
   lh = sparsifyLU l aa
@@ -471,13 +457,13 @@
 
 -- ** SSOR
 
--- | `mSsor aa omega` : if `omega = 1` it returns the diagonal of `aa`, 
-mSsor :: Fractional a => SpMatrix a -> a -> SpMatrix a
-mSsor aa omega = l ## r where
+-- | `mSsor aa omega` : if `omega = 1` it returns the symmetric Gauss-Seidel preconditioner
+mSsor :: Fractional a => SpMatrix a -> a -> (SpMatrix a, SpMatrix a)
+mSsor aa omega = (l, r) where
   (e, d, f) = diagPartitions aa
   n = nrows e
-  l = d ^-^ scale omega e
-  r = eye n ^-^ scale omega (reciprocal d ## f)
+  l = (eye n ^-^ scale omega e) ## reciprocal d
+  r = d ^-^ scale omega f 
 
 
 
@@ -486,11 +472,27 @@
 
 
 
+-- Linear solver, LU-based
+
+
+
+
+
+
+
+
+
 -- * Iterative linear solvers
 
 
+-- ** GMRES
 
+-- *** Left-preconditioning
 
+
+
+
+
 -- ** CGNE
 
 cgneStep :: SpMatrix Double -> CGNE -> CGNE
@@ -840,7 +842,7 @@
 -- | convergence check (FIXME)
 normDiffConverged :: (Foldable t, Functor t) =>
      (a -> SpVector Double) -> t a -> Bool
-normDiffConverged fp xx = normSq (foldrMap fp (^-^) (zeroSV 0) xx) <= eps
+normDiffConverged fp xx = nearZero $ normSq (foldrMap fp (^-^) (zeroSV 0) xx)
 
 
   
diff --git a/test/LibSpec.hs b/test/LibSpec.hs
--- a/test/LibSpec.hs
+++ b/test/LibSpec.hs
@@ -1,4 +1,13 @@
 {-# language ScopedTypeVariables #-}
+-----------------------------------------------------------------------------
+-- |
+-- Copyright   :  (C) 2016 Marco Zocca
+-- License     :  GPL-3 (see LICENSE)
+-- Maintainer  :  zocca.marco gmail
+-- Stability   :  provisional
+-- Portability :  portable
+--
+-----------------------------------------------------------------------------
 module LibSpec where
 
 import Numeric.LinearAlgebra.Sparse
@@ -40,9 +49,9 @@
     it "transposeSM : sparse matrix transpose" $
       transposeSM m1 `shouldBe` m1t
     it "matVec : matrix-vector product" $
-      normSq ((aa0 #> x0true) ^-^ b0 ) <= eps `shouldBe` True
+      nearZero ( normSq ((aa0 #> x0true) ^-^ b0 )) `shouldBe` True
     it "vecMat : vector-matrix product" $
-      normSq ((x0true <# aa0) ^-^ aa0tx0 ) <= eps `shouldBe` True  
+      nearZero ( normSq ((x0true <# aa0) ^-^ aa0tx0 ))`shouldBe` True  
     it "matMat : matrix-matrix product" $
       (m1 `matMat` m2) `shouldBe` m1m2
     it "eye : identity matrix" $
@@ -58,22 +67,22 @@
     it "countSubdiagonalNZ : # of nonzero elements below the diagonal" $
       countSubdiagonalNZSM m3 `shouldBe` 1
     it "permutPairsSM : permutation matrices are orthogonal" $ do
-      let pm0 = permutPairsSM 3 [(0,2), (1,2)]
+      let pm0 = permutPairsSM 3 [(0,2), (1,2)] :: SpMatrix Double
       pm0 ##^ pm0 `shouldBe` eye 3
       pm0 #^# pm0 `shouldBe` eye 3         
     it "modifyInspectN : early termination by iteration count" $
-      execState (modifyInspectN 2 ((< eps) . diffSqL) (/2)) 1 `shouldBe` 1/8
+      execState (modifyInspectN 2 (nearZero . diffSqL) (/2)) (1 :: Double) `shouldBe` 1/8
     it "modifyInspectN : termination by value convergence" $
-      execState (modifyInspectN (2^16) ((< eps) . head) (/2)) 1 < eps `shouldBe` True 
+      nearZero (execState (modifyInspectN (2^16) (nearZero . head) (/2)) (1 :: Double)) `shouldBe` True 
   describe "Numeric.LinearAlgebra.Sparse : Linear solvers" $ do
     -- it "TFQMR (2 x 2 dense)" $
     --   normSq (_xTfq (tfqmr aa0 b0 x0) ^-^ x0true) <= eps `shouldBe` True
     it "BCG (2 x 2 dense)" $
-      normSq (_xBcg (bcg aa0 b0 x0) ^-^ x0true) <= eps `shouldBe` True
+      nearZero (normSq (_xBcg (bcg aa0 b0 x0) ^-^ x0true)) `shouldBe` True
     it "BiCGSTAB (2 x 2 dense)" $ 
-      normSq (aa0 <\> b0 ^-^ x0true) <= eps `shouldBe` True
+      nearZero (normSq (aa0 <\> b0 ^-^ x0true)) `shouldBe` True
     it "CGS (2 x 2 dense)" $ 
-      normSq (_x (cgs aa0 b0 x0 x0) ^-^ x0true) <= eps `shouldBe` True
+      nearZero (normSq (_x (cgs aa0 b0 x0 x0) ^-^ x0true)) `shouldBe` True
   describe "Numeric.LinearAlgebra.Sparse : QR decomposition" $ do    
     it "QR (4 x 4 sparse)" $
       checkQr tm4 `shouldBe` True
@@ -84,6 +93,9 @@
       checkLu tm6 `shouldBe` True
     it "LU (10 x 10 sparse)" $
       checkLu tm7 `shouldBe` True
+  describe "Numeric.LinearAlgebra.Sparse : Cholesky decomposition (PSD matrices only)" $ do
+    it "chol (5 x 5 sparse)" $
+      checkChol tm7 `shouldBe` True
 
 
 {-
@@ -252,10 +264,11 @@
 
 {- QR-}
 
-checkQr :: SpMatrix Double -> Bool
+
+checkQr :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool
 checkQr a = c1 && c2 where
   (q, r) = qr a
-  c1 = normFrobenius ((q #~# r) ^-^ a) <= eps
+  c1 = nearZero $ normFrobenius ((q #~# r) ^-^ a)
   c2 = isOrthogonalSM q
 
 
@@ -266,13 +279,24 @@
 
 {- LU -}
 
-checkLu :: SpMatrix Double -> Bool
+checkLu :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool
 checkLu a = lup == a where
   (l, u) = lu a
   lup = l #~# u
 
 
 
+{- Cholesky -}
+
+checkChol :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool
+checkChol a = nearZero $ normFrobenius ((l ##^ l) ^-^ a) where
+  l = chol a
+
+
+
+
+
+
 {- eigenvalues -}
 
 
@@ -345,9 +369,9 @@
 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
+  a = mkSubDiagonal n 1 $ replicate n (-1)
+  b = mkSubDiagonal n 0 $ replicate n 2
+  c = mkSubDiagonal n (-1) $ replicate n (-1)
 
 -- -- run N iterations