diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -16,10 +16,14 @@
 
     * BiConjugate Gradient Stabilized (BiCGSTAB) (non-Hermitian systems)
 
+    * Transpose-Free Quasi-Minimal Residual (TFQMR)
+
 * Matrix decompositions
 
     * QR factorization
 
+    * LU factorization
+
 * Eigenvalue algorithms
 
     * QR algorithm
@@ -29,6 +33,60 @@
 * Utilities : Vector and matrix norms, matrix condition number, Givens rotation, Householder reflection
 
 * Predicates : Matrix orthogonality test (A^T A ~= I)
+
+
+---------
+
+## Examples
+
+The module `Numeric.LinearAlgebra.Sparse` contains the interface functions:
+
+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
+
+e.g.
+
+    > amat = fromListSM (3,3) [(0,0,2),(1,0,4),(1,1,3),(1,2,2),(2,2,5)]
+
+And similarly for sparse vectors : `fromListSV :: Int -> [(Int, a)] -> SpVector a`.
+
+Both sparse vectors and matrices can be pretty-printed using `prd`:
+
+    > prd amat
+    ( 3 rows, 3 columns ) , 5 NZ ( sparsity 0.5555555555555556 )
+
+    [2,0,0]
+    [4,3,2]
+    [0,0,5]
+
+
+Matrix factorizations are available as `lu` and `qr` respectively, and are straightforward to verify by using the matrix product `##`  :
+
+    > (l, u) = lu amat
+    > prd $ l ## u
+    ( 3 rows, 3 columns ) , 9 NZ ( sparsity 1.0 )
+
+    [2.0,0.0,0.0]
+    [4.0,3.0,2.0]
+    [0.0,0.0,5.0]
+
+Linear systems can be solved with either `linSolve` (which also requires choosing a method) or with `<\>` (which uses BiCGSTAB as default) :
+
+    > b = fromListSV 3 [(0,3),(1,2),(2,5)]
+    > x = amat <\> b
+    > prd x
+    ( 3 elements ) ,  3 NZ ( sparsity 1.0 )
+
+    [1.4999999999999998,-1.9999999999999998,0.9999999999999998]
+
+The result can be verified by computing the matrix-vector action `amat #> x`, which should (ideally) be very close to the right-hand side `b` :
+
+    > prd $ amat #> x
+    ( 3 elements ) ,  3 NZ ( sparsity 1.0 )
+
+    [2.9999999999999996,1.9999999999999996,4.999999999999999]
+
 
 
 
diff --git a/sparse-linear-algebra.cabal b/sparse-linear-algebra.cabal
--- a/sparse-linear-algebra.cabal
+++ b/sparse-linear-algebra.cabal
@@ -1,7 +1,7 @@
 name:                sparse-linear-algebra
-version:             0.2.0.5
+version:             0.2.0.7
 synopsis:            Numerical computation in native Haskell 
-description:         Currently it provides iterative linear solvers, matrix decompositions, eigenvalue computations and related utilities. Please see README.md for details
+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
 license:             GPL-3
 license-file:        LICENSE
@@ -20,12 +20,12 @@
   hs-source-dirs:      src
   exposed-modules:     Numeric.LinearAlgebra.Sparse
                        Numeric.LinearAlgebra.Class
-                       Numeric.LinearAlgebra.Data
                        Numeric.LinearAlgebra.Sparse.IntMap
                        Data.Sparse.SpVector
                        Data.Sparse.SpMatrix
                        Data.Sparse.Common
                        Data.Sparse.Utils
+                       Data.Sparse.Types
                        Numeric.Eps
   build-depends:       QuickCheck
                      , base >= 4.7 && < 5
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,28 +1,64 @@
 module Data.Sparse.Common
        ( module X,
-         svToSM, outerProdSV, (><), toSV, extractCol, extractRow,
-         extractDiagonalDSM,
-         matVec, (#>), vecMat, (<#)) where
+         insertRowWith, insertRow, insertColWith, insertCol,
+         outerProdSV, (><), toSV, svToSM, 
+         extractCol, extractRow,
+         extractVectorDenseWith, extractRowDense, extractColDense,
+         extractDiagDense,
+         extractSubRow, extractSubCol,
+         extractSubRow_RK, extractSubCol_RK,
+         matVec, (#>), vecMat, (<#),
+         prd) where
 
 import Data.Sparse.Utils as X
+import Data.Sparse.Types as X
 import Data.Sparse.SpMatrix as X
 import Data.Sparse.SpVector as X
 
 import Numeric.Eps as X
-import Numeric.LinearAlgebra.Data as X
 import Numeric.LinearAlgebra.Class as X
 
 import Numeric.LinearAlgebra.Sparse.IntMap as X
 
 import qualified Data.IntMap as IM
 
--- | promote a SV to SM
-svToSM :: SpVector a -> SpMatrix a
-svToSM (SV n d) = SM (n, 1) $ IM.singleton 0 d
 
 
--- ** Outer vector product
+-- * Insert row/column vector in matrix
 
+-- | Insert row , using the provided row index transformation function
+insertRowWith :: (IxCol -> IxCol) -> SpMatrix a -> SpVector a -> IM.Key -> SpMatrix a
+insertRowWith fj (SM (m,n) im) (SV d sv) i
+  | not (inBounds0 m i) = error "insertRowSM : index out of bounds"
+  | n >= d = SM (m,n) $ IM.insert i (insertOrUnion i sv' im) im
+  | otherwise = error $ "insertRowSM : incompatible dimensions " ++ show (n, d)
+    where sv' = IM.mapKeys fj sv
+          insertOrUnion i sv im = maybe sv (IM.union sv) (IM.lookup i im)
+
+-- | Insert row          
+insertRow :: SpMatrix a -> SpVector a -> IM.Key -> SpMatrix a
+insertRow = insertRowWith id
+
+-- | Insert column, using the provided row index transformation function
+insertColWith :: (IxRow -> IxRow) -> SpMatrix a -> SpVector a -> IxCol -> SpMatrix a
+insertColWith fi smm sv j
+  | not (inBounds0 n j) = error "insertColSM : index out of bounds"
+  | m >= mv = insIM2 smm vl j
+  | otherwise = error $ "insertColSM : incompatible dimensions " ++ show (m,mv) where
+      (m, n) = dim smm
+      mv = dim sv
+      vl = toListSV sv
+      insIM2 im2 ((i,x):xs) j = insIM2 (insertSpMatrix (fi i) j x im2) xs j
+      insIM2 im2 [] _ = im2
+
+-- | Insert column
+insertCol :: SpMatrix a -> SpVector a -> IxCol -> SpMatrix a
+insertCol = insertColWith id
+
+
+
+-- * Outer vector product
+
 outerProdSV, (><) :: Num a => SpVector a -> SpVector a -> SpMatrix a
 outerProdSV v1 v2 = fromListSM (m, n) ixy where
   m = dim v1
@@ -33,39 +69,79 @@
 
 
 
+-- * Matrix-vector conversions
+
+-- | promote a SV to SM
+svToSM :: SpVector a -> SpMatrix a
+svToSM (SV n d) = SM (n, 1) $ IM.singleton 0 d
+
 -- |Demote (n x 1) or (1 x n) SpMatrix to SpVector
 toSV :: SpMatrix a -> SpVector a
-toSV (SM (m,n) im) = SV d $ snd . head $ IM.toList im where
+toSV (SM (m,n) im) = SV d (ff im) where
+  ff | m > n = IM.map g  -- column case
+     | otherwise = g
+  g = snd . head . IM.toList
   d | m==1 && n==1 = 1
     | m==1 && n>1 = n 
     | n==1 && m>1 = m
-    | otherwise = error $ "toSV : incompatible dimensions " ++ show (m,n)
+    | otherwise = error $ "toSV : incompatible matrix dimension " ++ show (m,n)
 
--- |Extract jth column, and place into SpVector
-extractCol :: SpMatrix a -> IxCol -> SpVector a
-extractCol m j = toSV $ extractColSM m j      
 
+-- * Extract a SpVector from an SpMatrix
+-- ** Sparse extract
 
--- |Extract ith row, and place into SpVector
+-- |Extract ith row
 extractRow :: SpMatrix a -> IxRow -> SpVector a
 extractRow m i = toSV $ extractRowSM m i
 
+-- |Extract jth column
+extractCol :: SpMatrix a -> IxCol -> SpVector a
+extractCol m j = toSV $ extractColSM m j      
 
 
 
+-- ** Dense extract (default == 0)
+-- | Generic extraction function
+extractVectorDenseWith ::
+  Num a => (Int -> (IxRow, IxCol)) -> SpMatrix a -> SpVector a
+extractVectorDenseWith f mm = fromListDenseSV n $ foldr ins [] ll  where
+  ll = [0 .. n - 1]
+  (m, n) = dim mm
+  ins i acc = mm @@ f i : acc
 
+-- | Extract ith row (dense)
+extractRowDense :: Num a => SpMatrix a -> IxRow -> SpVector a
+extractRowDense mm iref = extractVectorDenseWith (\j -> (iref, j)) mm
 
+-- | Extract jth column
+extractColDense :: Num a => SpMatrix a -> IxCol -> SpVector a
+extractColDense mm jref = extractVectorDenseWith (\i -> (i, jref)) mm
 
--- | Extract the diagonal as a SpVector (with default 0)
-extractDiagonalDSM :: Num a => SpMatrix a -> SpVector a
-extractDiagonalDSM mm = fromListDenseSV n $ foldr ins [] ll  where
-  ll = [0 .. n - 1]
-  n = nrows mm
-  ins i acc = mm@@(i,i) : acc
+-- | Extract the diagonal
+extractDiagDense :: Num a => SpMatrix a -> SpVector a
+extractDiagDense = extractVectorDenseWith (\i -> (i, i))
 
 
 
+-- | extract row interval
+extractSubRow :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpVector a
+extractSubRow m i (j1, j2)  = toSV $ extractSubRowSM m i (j1, j2)
 
+-- | extract row interval, rebalance keys by subtracting lowest one
+extractSubRow_RK :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpVector a
+extractSubRow_RK m i (j1, j2)  = toSV $ extractSubRowSM_RK m i (j1, j2)
+
+-- | extract column interval
+extractSubCol :: SpMatrix a -> IxCol -> (IxRow, IxRow) -> SpVector a
+extractSubCol m j (i1, i2)  = toSV $ extractSubColSM m j (i1, i2)
+
+-- | extract column interval, rebalance keys by subtracting lowest one
+extractSubCol_RK :: SpMatrix a -> IxCol -> (IxRow, IxRow) -> SpVector a
+extractSubCol_RK m j (i1, i2)  = toSV $ extractSubColSM_RK m j (i1, i2)
+
+
+
+
 -- ** Matrix action on a vector
 
 {- 
@@ -92,3 +168,99 @@
   | otherwise = error $ "vecMat : mismatching dimensions " ++ show (n, nr)
 
 (<#) = vecMat  
+
+
+
+
+
+
+-- * Pretty printing
+
+
+showNonZero :: (Show a, Num a, Eq a) => a -> String
+showNonZero x  = if x == 0 then " " else show x
+
+
+toDenseRow :: Num a => SpMatrix a -> IM.Key -> [a]
+toDenseRow (SM (_,ncol) im) irow =
+  fmap (\icol -> im `lookupWD_IM` (irow,icol)) [0..ncol-1]
+
+toDenseRowClip :: (Show a, Num a) => SpMatrix a -> IM.Key -> Int -> String
+toDenseRowClip sm irow ncomax
+  | ncols sm > ncomax = unwords (map show h) ++  " ... " ++ show t
+  | otherwise = show dr
+     where dr = toDenseRow sm irow
+           h = take (ncomax - 2) dr
+           t = last dr
+
+newline :: IO ()
+newline = putStrLn ""
+
+printDenseSM :: (Show t, Num t) => SpMatrix t -> IO ()
+printDenseSM sm = do
+  newline
+  putStrLn $ sizeStr sm
+  newline
+  printDenseSM' sm 5 5
+  newline
+  where    
+    printDenseSM' :: (Show t, Num t) => SpMatrix t -> Int -> Int -> IO ()
+    printDenseSM' sm'@(SM (nr,_) _) nromax ncomax = mapM_ putStrLn rr_' where
+      rr_ = map (\i -> toDenseRowClip sm' i ncomax) [0..nr - 1]
+      rr_' | nrows sm > nromax = take (nromax - 2) rr_ ++ [" ... "] ++[last rr_]
+           | otherwise = rr_
+
+
+toDenseListClip :: (Show a, Num a) => SpVector a -> Int -> String
+toDenseListClip sv ncomax
+  | dim sv > ncomax = unwords (map show h) ++  " ... " ++ show t
+  | otherwise = show dr
+     where dr = toDenseListSV sv
+           h = take (ncomax - 2) dr
+           t = last dr
+
+printDenseSV :: (Show t, Num t) => SpVector t -> IO ()
+printDenseSV sv = do
+  newline
+  putStrLn $ sizeStrSV sv
+  newline
+  printDenseSV' sv 5
+  newline where
+    printDenseSV' v nco = putStrLn rr_' where
+      rr_ = toDenseListClip v nco :: String
+      rr_' | dim sv > nco = unwords [take (nco - 2) rr_ , " ... " , [last rr_]]
+           | otherwise = rr_
+
+-- ** Pretty printer typeclass
+class PrintDense a where
+  prd :: a -> IO ()
+
+instance (Show a, Num a) => PrintDense (SpVector a) where
+  prd = printDenseSV
+
+instance (Show a, Num a) => PrintDense (SpMatrix a) where
+  prd = printDenseSM
+
+
+
+  
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+  
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
@@ -2,12 +2,10 @@
 module Data.Sparse.SpMatrix where
 
 import Data.Sparse.Utils
-
-import Data.Sparse.Utils
+import Data.Sparse.Types
 
 import Numeric.Eps
 import Numeric.LinearAlgebra.Class
-import Numeric.LinearAlgebra.Data
 
 import Numeric.LinearAlgebra.Sparse.IntMap
 
@@ -184,19 +182,23 @@
 
 -- | Looks up an element in the matrix with a default (if the element is not found, zero is returned)
 
-lookupWD_SM, (@@) :: Num a => SpMatrix a -> (IxRow, IxCol) -> a
+lookupWD_SM, (@@!), (@@) :: Num a => SpMatrix a -> (IxRow, IxCol) -> a
 lookupWD_SM sm (i,j) =
   fromMaybe 0 (lookupSM sm i j)
 
 lookupWD_IM :: Num a => IM.IntMap (IM.IntMap a) -> (IxRow, IxCol) -> a
 lookupWD_IM im (i,j) = fromMaybe 0 (IM.lookup i im >>= IM.lookup j)
 
--- | Zero-default lookup, infix form
-(@@) = lookupWD_SM
+-- | Zero-default lookup, infix form (no bound checking)
+(@@!) = lookupWD_SM
 
+-- | Zero-default lookup, infix form ("safe" : throws exception if lookup is outside matrix bounds)
+m @@ d | isValidIxSM m d = m @@! d
+       | otherwise = error $ "@@ : incompatible indices : matrix size is " ++ show (dim m) ++ ", but user looked up " ++ show d
 
 
 
+
 -- FIXME : to throw an exception or just ignore the out-of-bound access ?
 
 
@@ -207,12 +209,18 @@
 
 -- ** Sub-matrices
 
--- | Extract a submatrix given the specified index bounds
-extractSubmatrixSM :: SpMatrix a -> (IxRow, IxCol) -> (IxRow, IxCol) -> SpMatrix a
-extractSubmatrixSM (SM (r, c) im) (i1, i2) (j1, j2)
+
+-- | Extract a submatrix given the specified index bounds, rebalancing keys with the two supplied functions
+extractSubmatrixSM ::
+  (IM.Key -> IM.Key) ->   -- row index function
+  (IM.Key -> IM.Key) ->   -- column "  "
+  SpMatrix a ->
+  (IxRow, IxRow) -> (IxCol, IxCol) ->
+  SpMatrix a
+extractSubmatrixSM fi gi (SM (r, c) im) (i1, i2) (j1, j2)
   | q = SM (m', n') imm'
-  | otherwise = error $ "extractSubmatrixSM : invalid indexing " ++ show (i1, i2) ++ ", " ++ show (j1, j2) where
-  imm' = mapKeysIM2 (\i -> i - i1) (\j -> j - j1) $  -- rebalance keys
+  | otherwise = error $ "extractSubmatrixSM : invalid index " ++ show (i1, i2) ++ ", " ++ show (j1, j2) where
+  imm' = mapKeysIM2 fi gi $                          -- rebalance keys
           IM.filter (not . IM.null) $                -- remove all-null rows
           ifilterIM2 ff im                           -- keep `submatrix`
   ff i j _ = i1 <= i &&
@@ -226,19 +234,52 @@
       inBounds0 c j2 &&      
       i2 >= i1
 
+-- | Extract a submatrix given the specified index bounds
+-- NB : subtracts (i1, j1) from the indices
+extractSubmatrixRebalanceKeys ::
+  SpMatrix a -> (IxRow, IxRow) -> (IxCol, IxCol) -> SpMatrix a
+extractSubmatrixRebalanceKeys mm (i1,i2) (j1,j2) =
+  extractSubmatrixSM (\i -> i - i1) (\j -> j - j1) mm (i1,i2) (j1,j2)
 
+-- | Extract a submatrix given the specified index bounds
+-- NB : submatrix indices are _preserved_
+extractSubmatrix :: SpMatrix a -> (IxRow, IxRow) -> (IxCol, IxCol) -> SpMatrix a
+extractSubmatrix = extractSubmatrixSM id id
 
 
+
+
+-- *** Extract i'th row
+extractRowSM :: SpMatrix a -> IxRow -> SpMatrix a
+extractRowSM sm i = extractSubmatrix sm (i, i) (0, ncols sm - 1)
+
+
+-- | Extract column within a row range
+extractSubRowSM :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpMatrix a
+extractSubRowSM sm i (j1, j2) = extractSubmatrix sm (i, i) (j1, j2)
+
+-- | Extract column within a row range, rebalance keys
+extractSubRowSM_RK :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpMatrix a
+extractSubRowSM_RK sm i =
+  extractSubmatrixRebalanceKeys sm (i, i) 
+
+
+
 -- *** Extract j'th column
+-- | Extract all column
 extractColSM :: SpMatrix a -> IxCol -> SpMatrix a
-extractColSM sm j = extractSubmatrixSM sm (0, nrows sm - 1) (j, j)
+extractColSM sm j = extractSubmatrix sm (0, nrows sm - 1) (j, j)
 
+-- | Extract column within a row range
+extractSubColSM :: SpMatrix a -> IxCol -> (IxRow, IxRow) -> SpMatrix a
+extractSubColSM sm j (i1, i2) = extractSubmatrix sm (i1, i2) (j, j)
 
 
+-- | Extract column within a row range, rebalance keys
+extractSubColSM_RK :: SpMatrix a -> IxCol -> (IxRow, IxRow) -> SpMatrix a
+extractSubColSM_RK sm j (i1, i2) =
+  extractSubmatrixRebalanceKeys sm (i1, i2) (j, j) 
 
--- *** Extract i'th row
-extractRowSM :: SpMatrix a -> IxRow -> SpMatrix a
-extractRowSM sm i = extractSubmatrixSM sm (i, i) (0, ncols sm - 1)
 
 
 
@@ -255,8 +296,8 @@
 
 -- ** Predicates
 -- |Are the supplied indices within matrix bounds?
-validIxSM :: SpMatrix a -> (Int, Int) -> Bool
-validIxSM mm = inBounds02 (dim mm)
+isValidIxSM :: SpMatrix a -> (Int, Int) -> Bool
+isValidIxSM mm = inBounds02 (dim mm)
 
 -- |Is the matrix square?
 isSquareSM :: SpMatrix a -> Bool
@@ -626,4 +667,12 @@
 
 
 
-
+-- *** Matrix contraction
+-- | Contract two matrices A and 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
+  | 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"
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
@@ -2,9 +2,9 @@
 module Data.Sparse.SpVector where
 
 import Data.Sparse.Utils
+import Data.Sparse.Types
 
 import Numeric.LinearAlgebra.Class
-import Numeric.LinearAlgebra.Data
 import Numeric.LinearAlgebra.Sparse.IntMap
 
 import Data.Maybe
@@ -21,8 +21,13 @@
 spySV :: Fractional b => SpVector a -> b
 spySV s = fromIntegral (IM.size (dat s)) / fromIntegral (dim s)
 
+-- | Number of nonzeros
+nzSV :: SpVector a -> Int
+nzSV sv = IM.size (dat sv)
 
 
+sizeStrSV :: SpVector a -> String
+sizeStrSV sv = unwords ["(",show (dim sv),"elements ) , ",show (nzSV sv),"NZ ( sparsity", show (spy sv),")"]
 
 
 
diff --git a/src/Data/Sparse/Types.hs b/src/Data/Sparse/Types.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Sparse/Types.hs
@@ -0,0 +1,7 @@
+module Data.Sparse.Types where
+
+type Rows = Int
+type Cols = Int
+
+type IxRow = Int
+type IxCol = Int
diff --git a/src/Numeric/Eps.hs b/src/Numeric/Eps.hs
--- a/src/Numeric/Eps.hs
+++ b/src/Numeric/Eps.hs
@@ -15,6 +15,9 @@
 almostZero x = abs x <= eps
 almostOne x = x >= (1-eps) && x < (1+eps)
 
+isNz :: Double -> Bool
+isNz = not . almostZero
+
 withDefault :: (t -> Bool) -> t -> t -> t
 withDefault q d x | q x = d
                   | otherwise = x
diff --git a/src/Numeric/LinearAlgebra/Data.hs b/src/Numeric/LinearAlgebra/Data.hs
deleted file mode 100644
--- a/src/Numeric/LinearAlgebra/Data.hs
+++ /dev/null
@@ -1,7 +0,0 @@
-module Numeric.LinearAlgebra.Data where
-
-type Rows = Int
-type Cols = Int
-
-type IxRow = Int
-type IxCol = Int
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
@@ -1,6 +1,23 @@
 {-# LANGUAGE FlexibleContexts, TypeFamilies, MultiParamTypeClasses, FlexibleInstances  #-}
 -- {-# OPTIONS_GHC -O2 -rtsopts -with-rtsopts=-K32m -prof#-}
-module Numeric.LinearAlgebra.Sparse where
+module Numeric.LinearAlgebra.Sparse
+       (
+         sparsifySV,
+         conditionNumberSM,
+         hhMat, hhRefl,
+         givens,
+         eigsQR, eigRayleigh,
+         cgne, tfqmr, bicgstab, cgs, bcg,
+         _xCgne, _xTfq, _xBicgstab, _x, _xBcg,
+         cgsStep, bicgstabStep,
+         CGNE, TFQMR, BICGSTAB, CGS, BCG,
+         linSolve, LinSolveMethod, (<\>),
+         randMat, randVec, randSpMat, randSpVec,
+         modifyInspectN, runAppendN',
+         diffSqL,
+         qr, lu
+       )
+       where
 
 
 import Data.Sparse.Common
@@ -9,7 +26,7 @@
 import Control.Monad.Primitive
 import Control.Monad (mapM_, forM_, replicateM)
 import Control.Monad.State.Strict
--- import Control.Monad.Writer
+import Control.Monad.Writer
 -- import Control.Monad.Trans
 -- import Control.Monad.Trans.State (runStateT)
 -- import Control.Monad.Trans.Writer (runWriterT)
@@ -49,7 +66,7 @@
                     | otherwise = kappa where
   kappa = lmax / lmin
   (_, r) = qr m
-  u = extractDiagonalDSM r  -- FIXME : need to extract with default element 0 
+  u = extractDiagDense r  -- FIXME : need to extract with default element 0 
   lmax = abs (maximum u)
   lmin = abs (minimum u)
 
@@ -118,7 +135,7 @@
 
 givens :: SpMatrix Double -> IxRow -> IxCol -> SpMatrix Double
 givens mm i j 
-  | validIxSM mm (i,j) && isSquareSM mm =
+  | isValidIxSM mm (i,j) && isSquareSM mm =
        sparsifySM $ fromListSM' [(i,i,c),(j,j,c),(j,i,-s),(i,j,s)] (eye (nrows mm))
   | otherwise = error "givens : indices out of bounds"      
   where
@@ -180,15 +197,15 @@
 
 -- * Eigenvalue algorithms
 
--- ** All eigenvalues (QR algorithm)
+-- ** QR algorithm
 
--- | `eigsQR n mm` performs `n` iterations of the QR algorithm on matrix `mm` 
+-- | `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 nitermax m = extractDiagonalDSM $ execState (convergtest eigsStep) m where
+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 = extractDiagonalDSM m1
-                     dm2 = extractDiagonalDSM m2
+    f [m1, m2] = let dm1 = extractDiagDense m1
+                     dm2 = extractDiagDense m2
                  in norm2 (dm1 ^-^ dm2) <= eps
 
 
@@ -196,10 +213,9 @@
 
 
 
--- ** One eigenvalue and eigenvector (Rayleigh iteration)
-
+-- ** Rayleigh iteration
 
--- | `eigsRayleigh n mm` performs `n` iterations of the Rayleigh algorithm on matrix `mm`. Cubic-order convergence, but it requires a mildly educated guess on the initial eigenpair
+-- | `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)
@@ -272,22 +288,131 @@
 
 
 -- * 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)
 
--- lu aa | isSquareSM aa = undefined
---       | otherwise = error "LU factorization not currently defined for rectangular matrices" where
---           n = nrows aa
---           l0 = eye n
---           aa0 = 
+-- | 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
 
--- luStep aa i l u 
+-- | 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)
+
+
+-- 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"]
+  where
+   a = amat @@! (i, j)
+   uii = umat @@! (i-1, i-1) -- NB this must be /= 0
+   p = contractSub lmat umat i j (i - 1)
+
+
+
+
+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))
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+-- test data
+
+aa0 = fromListDenseSM 3 [2, -4, -4, -1, 6, -2, -2, 3, 8] :: SpMatrix Double
+
+
+
+
+
+
+
 -- Produces the permutation matrix necessary to have a nonzero in position (iref, jref). This is used in the LU factorization
-permutAA :: Num b => SpMatrix a -> IxRow -> IxCol -> Maybe (SpMatrix b)
-permutAA (SM (nro,_) mm) iref jref
+permutAA :: Num b => IxRow -> IxCol -> SpMatrix a -> Maybe (SpMatrix b)
+permutAA iref jref (SM (nro,_) mm) 
   | isJust (lookupIM2 iref jref mm) = Nothing -- eye nro
   | otherwise = Just $ permutationSM nro [head u] where
       u = IM.keys (ifilterIM2 ff mm)
@@ -302,19 +427,106 @@
 
 
 
+
+
+-- * Incomplete LU
+
+-- | used for Incomplete LU : remove entries in `m` corresponding to zero entries in `m2`
+ripHoles :: SpMatrix t -> SpMatrix a -> SpMatrix a
+ripHoles (SM d m) m2 = SM d $ ifilterIM2 f (dat m2) where
+  f i j _ = isJust (lookupSM m2 i j)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 -- * Iterative linear solvers
 
 
 
--- -- | residual of candidate solution x0 of a linear system
--- residual :: Num a => SpMatrix a -> SpVector a -> SpVector a -> SpVector a
--- residual aa b x0 = b ^-^ (aa #> x0)
 
--- converged :: SpMatrix Double -> SpVector Double -> SpVector Double -> Bool
--- converged aa b x0 = normSq (residual aa b x0) <= eps
+-- ** CGNE
 
+cgneStep :: SpMatrix Double -> CGNE -> CGNE
+cgneStep aa (CGNE x r p) = CGNE x1 r1 p1 where
+  alphai = r `dot` r / (p `dot` p)
+  x1 = x ^+^ (alphai .* p)
+  r1 = r ^-^ (alphai .* (aa #> p))
+  beta = r1 `dot` r1 / (r `dot` r)
+  p1 = transposeSM aa #> r ^+^ (beta .* p)
 
+data CGNE =
+  CGNE {_xCgne , _rCgne, _pCgne :: SpVector Double} deriving Eq
 
+cgne :: SpMatrix Double -> SpVector Double -> SpVector Double -> CGNE
+cgne aa b x0 = execState (untilConverged _xCgne (cgneStep aa)) cgneInit where
+  r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution
+  p0 = transposeSM aa #> r0
+  cgneInit = CGNE x0 r0 p0
+
+
+-- ** TFQMR
+
+-- | one step of TFQMR
+tfqmrStep :: SpMatrix Double -> SpVector Double -> TFQMR -> TFQMR
+tfqmrStep aa r0hat (TFQMR x w u v d m tau theta eta rho alpha) =
+  TFQMR x1 w1 u1 v1 d1 (m+1) tau1 theta1 eta1 rho1 alpha1
+  where
+  -- alpham = alpha
+  w1 = w ^-^ (alpha .* (aa #> u))
+  d1 = u ^+^ ((theta**2/alpha*eta) .* d)
+  theta1 = norm2 w1 / tau
+  c = recip $ sqrt (1 + theta1**2)
+  tau1 = tau * theta1 * c
+  eta1 = c**2 * alpha
+  x1 = x^+^ (eta1 .* d1)
+  (alpha1, u1, rho1, v1)
+    | even m = let
+                   alpha' = rho / (v `dot` r0hat)
+                   u' = u ^-^ (alpha' .* v)
+               in
+                   (alpha', u', rho, v)
+    | otherwise = let
+                   rho' = w1 `dot` r0hat
+                   beta = rho'/rho
+                   u' = w1 ^+^ (beta .* u)
+                   v' = (aa #> u') ^+^ (beta .* (aa #> u ^+^ (beta .* v)) )
+                  in (alpha, u', rho', v')
+
+tfqmr :: SpMatrix Double -> SpVector Double -> SpVector Double -> TFQMR  
+tfqmr aa b x0 = execState (untilConverged _xTfq (tfqmrStep aa r0)) tfqmrInit where
+  n = dim b
+  r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution
+  w0 = r0
+  u0 = r0
+  v0 = aa #> u0
+  d0 = zeroSV n
+  r0hat = r0
+  rho0 = r0hat `dot` r0
+  alpha0 = rho0 / (v0 `dot` r0hat)
+  m = 0
+  tau0 = norm2 r0
+  theta0 = 0
+  eta0 = 0
+  tfqmrInit = TFQMR x0 w0 u0 v0 d0 m tau0 theta0 eta0 rho0 alpha0
+
+data TFQMR =
+  TFQMR { _xTfq, _wTfq, _uTfq, _vTfq, _dTfq :: SpVector Double,
+          _mTfq :: Int,
+          _tauTfq, _thetaTfq, _etaTfq, _rhoTfq, _alphaTfq :: Double}
+  deriving Eq
+
+
+
 -- ** BCG
 
 -- | one step of BCG
@@ -449,7 +661,7 @@
 
 -- * Linear solver interface
 
-data LinSolveMethod = BCG_ | CGS_ | BICGSTAB_ deriving (Eq, Show) 
+data LinSolveMethod = CGNE_ | TFQMR_ | BCG_ | CGS_ | BICGSTAB_ deriving (Eq, Show) 
 
 -- -- | Linear solve with _random_ starting vector
 -- linSolveM ::
@@ -472,9 +684,11 @@
       solve aa' b' | isDiagonalSM aa' = reciprocal aa' #> b' -- diagonal solve is easy
                    | otherwise = solveWith aa' b' 
       solveWith aa' b' = case method of
-                                BCG_ -> _xBcg (bcg aa' b' x0)
-                                CGS_ ->  _xBicgstab (bicgstab aa' b' x0 x0)
-                                BICGSTAB_ -> _x (cgs aa' b' x0 x0)
+        CGNE_ -> _xCgne (cgne aa' b' x0)
+        TFQMR_ -> _xTfq (tfqmr aa' b' x0)
+        BCG_ -> _xBcg (bcg aa' b' x0)
+        CGS_ ->  _xBicgstab (bicgstab aa' b' x0 x0)
+        BICGSTAB_ -> _x (cgs aa' b' x0 x0)
       x0 = mkSpVectorD n $ replicate n 0.1 
       (m, n) = dim aa
       nb     = dim b
@@ -584,7 +798,7 @@
 
 -- | iterate until convergence is verified or we run out of a fixed iteration budget
 untilConverged :: MonadState a m => (a -> SpVector Double) -> (a -> a) -> m a
-untilConverged fproj = modifyInspectN 100 (normDiffConverged fproj)
+untilConverged fproj = modifyInspectN 200 (normDiffConverged fproj)
 
 -- | convergence check (FIXME)
 normDiffConverged :: (Foldable t, Functor t) =>
@@ -669,103 +883,3 @@
   ii <- replicateM nsp (MWC.uniformR (0, n-1) g :: IO Int)
   return $ fromListSV n $ zip ii aav
 
-
-
-
-
--- * Pretty printing
-
-
-
-
-
-showNonZero :: (Show a, Num a, Eq a) => a -> String
-showNonZero x  = if x == 0 then " " else show x
-
-
-toDenseRow :: Num a => SpMatrix a -> IM.Key -> [a]
-toDenseRow (SM (_,ncol) im) irow =
-  fmap (\icol -> im `lookupWD_IM` (irow,icol)) [0..ncol-1]
-
-toDenseRowClip :: (Show a, Num a) => SpMatrix a -> IM.Key -> Int -> String
-toDenseRowClip sm irow ncomax
-  | ncols sm > ncomax = unwords (map show h) ++  " ... " ++ show t
-  | otherwise = show dr
-     where dr = toDenseRow sm irow
-           h = take (ncomax - 2) dr
-           t = last dr
-
-newline :: IO ()
-newline = putStrLn ""
-
-printDenseSM :: (Show t, Num t) => SpMatrix t -> IO ()
-printDenseSM sm = do
-  newline
-  putStrLn $ sizeStr sm
-  newline
-  printDenseSM' sm 5 5
-  newline
-  where    
-    printDenseSM' :: (Show t, Num t) => SpMatrix t -> Int -> Int -> IO ()
-    printDenseSM' sm'@(SM (nr,_) _) nromax ncomax = mapM_ putStrLn rr_' where
-      rr_ = map (\i -> toDenseRowClip sm' i ncomax) [0..nr - 1]
-      rr_' | nrows sm > nromax = take (nromax - 2) rr_ ++ [" ... "] ++[last rr_]
-           | otherwise = rr_
-
-
-toDenseListClip :: (Show a, Num a) => SpVector a -> Int -> String
-toDenseListClip sv ncomax
-  | dim sv > ncomax = unwords (map show h) ++  " ... " ++ show t
-  | otherwise = show dr
-     where dr = toDenseListSV sv
-           h = take (ncomax - 2) dr
-           t = last dr
-
-printDenseSV :: (Show t, Num t) => SpVector t -> IO ()
-printDenseSV sv = do
-  newline
-  printDenseSV' sv 5
-  newline where
-    printDenseSV' v nco = putStrLn rr_' where
-      rr_ = toDenseListClip v nco :: String
-      rr_' | dim sv > nco = unwords [take (nco - 2) rr_ , " ... " , [last rr_]]
-           | otherwise = rr_
-
--- ** Pretty printer typeclass
-class PrintDense a where
-  prd :: a -> IO ()
-
-instance (Show a, Num a) => PrintDense (SpVector a) where
-  prd = printDenseSV
-
-instance (Show a, Num a) => PrintDense (SpMatrix a) where
-  prd = printDenseSM
-
-
-
-  
-
-
-
-
-
-
-
-
-
-
--- * Type synonyms for SpMatrix
-
-
-
-
-
-
-
-
-
-
-
-
-
-  
diff --git a/src/Numeric/LinearAlgebra/Sparse/IntMap.hs b/src/Numeric/LinearAlgebra/Sparse/IntMap.hs
--- a/src/Numeric/LinearAlgebra/Sparse/IntMap.hs
+++ b/src/Numeric/LinearAlgebra/Sparse/IntMap.hs
@@ -35,23 +35,33 @@
 
 
 
--- | ========= IntMap-of-IntMap (IM2) stuff
+-- * set-like brackets
 
 
--- insert an element
+-- unionWithKeyIM2 f im1 im2 = undefined where
+
+
+
+
+
+-- * Insertion
+
+-- | Insert an element
 insertIM2 ::
   IM.Key -> IM.Key -> a -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)
 insertIM2 i j x imm = IM.insert i ro imm where
   ro = maybe (IM.singleton j x) (IM.insert j x) (IM.lookup i imm)
 {-# inline insertIM2 #-}  
 
--- lookup a key
+-- * Lookup
+
+-- |Lookup a key
 lookupIM2 ::
   IM.Key -> IM.Key -> IM.IntMap (IM.IntMap a) -> Maybe a
 lookupIM2 i j imm = IM.lookup i imm >>= IM.lookup j
 {-# inline lookupIM2 #-}  
 
--- populate an IM2 from a list of (row index, column index, value)  
+-- |Ppopulate an IM2 from a list of (row index, column index, value)  
 fromListIM2 ::
   Foldable t =>
      t (IM.Key, IM.Key, a) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)
@@ -59,15 +69,16 @@
   ins t (i,j,x) = insertIM2 i j x t
 
 
--- | folding
+-- * folding
 
--- indexed fold over an IM2
+-- |Indexed left fold over an IM2, with general accumulator
 ifoldlIM2' :: (IM.Key -> IM.Key -> a -> b -> b) -> b -> IM.IntMap (IM.IntMap a) -> b
 ifoldlIM2' f empty mm = IM.foldlWithKey' accRow empty mm where
   accRow acc i r = IM.foldlWithKey' (accElem i) acc r
   accElem i acc j x = f i j x acc
 {-# inline ifoldlIM2' #-}
 
+-- |Indexed left fold over an IM2
 ifoldlIM2 ::
   (IM.Key -> IM.Key -> t -> IM.IntMap a -> IM.IntMap a) ->
   IM.IntMap (IM.IntMap t) ->  
@@ -77,6 +88,7 @@
   accElem i acc j x   = f i j x acc
 {-# inline ifoldlIM2 #-}  
 
+-- |Left fold over an IM2, with general accumulator
 foldlIM2 :: (a -> b -> b) -> b -> IM.IntMap (IM.IntMap a) -> b
 foldlIM2 f empty mm = IM.foldl accRow empty mm where
   accRow acc r = IM.foldl accElem acc r
@@ -84,7 +96,7 @@
 {-# inline foldlIM2 #-}
 
 
--- transposeIM2 : inner indices become outer ones and vice versa. No loss of information because both inner and outer IntMaps are nubbed.
+-- | Inner indices become outer ones and vice versa. No loss of information because both inner and outer IntMaps are nubbed.
 transposeIM2 :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)
 transposeIM2 = ifoldlIM2 (flip insertIM2)
 {-# inline transposeIM2 #-}
@@ -98,9 +110,9 @@
 
 
 
--- | filtering
+-- * filtering
 
--- map over outer IM and filter all inner IM's
+-- |Map over outer IM and filter all inner IM's
 ifilterIM2 ::
   (IM.Key -> IM.Key -> a -> Bool) ->
   IM.IntMap (IM.IntMap a) ->
@@ -109,9 +121,7 @@
   IM.mapWithKey (\irow row -> IM.filterWithKey (f irow) row) 
 {-# inline ifilterIM2 #-}
 
--- specialized filtering function
-
--- keep only sub-diagonal elements
+-- |Specialized filtering : keep only sub-diagonal elements
 filterSubdiag :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)
 filterSubdiag = ifilterIM2 (\i j _ -> i>j)
 
@@ -119,7 +129,7 @@
 countSubdiagonalNZ im =
   IM.size $ IM.filter (not . IM.null) (filterSubdiag im)
 
--- list of (row, col) indices of (nonzero) subdiagonal elements
+-- |List of (row, col) indices of (nonzero) subdiagonal elements
 subdiagIndices :: IM.IntMap (IM.IntMap a) -> [(IM.Key, IM.Key)]
 subdiagIndices im = concatMap rpairs $ IM.toList (IM.map IM.keys im') where
   im' = filterSubdiag im
@@ -137,15 +147,14 @@
   
 
 
--- | mapping
-
--- map over IM2
+-- * mapping
 
+-- |Map over IM2
 mapIM2 :: (a -> b) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap b)
 mapIM2 = IM.map . IM.map   -- imapIM2 (\_ _ x -> f x)
 
 
--- indexed map over IM2
+-- |Indexed map over IM2
 imapIM2 ::
   (IM.Key -> IM.Key -> a -> b) ->
   IM.IntMap (IM.IntMap a) ->
@@ -155,8 +164,7 @@
 
 
 
--- mapping keys
-
+-- |Mapping keys
 mapKeysIM2 ::
   (IM.Key -> IM.Key) -> (IM.Key -> IM.Key) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)
 mapKeysIM2 fi fj im = IM.map adjCols adjRows where
@@ -173,7 +181,4 @@
 
 
 
-
-
--- sparsification :
 
diff --git a/test/LibSpec.hs b/test/LibSpec.hs
--- a/test/LibSpec.hs
+++ b/test/LibSpec.hs
@@ -28,7 +28,7 @@
 
 spec :: Spec
 spec = do
-  describe "Math.Linear.Sparse : library" $ do
+  describe "Numeric.LinearAlgebra.Sparse : library" $ do
     -- prop "subtraction is cancellative" $ \(x :: SpVector Double) ->
     --   x ^-^ x `shouldBe` zero
     it "dot : inner product" $
@@ -43,58 +43,44 @@
       (m1 `matMat` m2) `shouldBe` m1m2
     it "eye : identity matrix" $
       infoSM (eye 10) `shouldBe` SMInfo 10 0.1
+    it "insertCol : insert a column in a SpMatrix" $
+      insertCol (eye 3) (fromListDenseSV 3 [2,2,2]) 0 `shouldBe` (fromListSM (3,3) [(0,0,2),(1,0,2),(1,1,1),(2,0,2),(2,2,1)])
+    it "insertRow : insert a row in a SpMatrix" $
+      insertRow (eye 3) (fromListDenseSV 3 [2,2,2]) 1 `shouldBe` (fromListSM (3,3) [(0,0,1), (1,0,2), (1,1,2), (1,2,2), (2,2,1)])
+    it "extractCol -> insertCol : identity" $
+      insertCol (eye 3) (extractCol (eye 3) 1) 1 `shouldBe` eye 3
+    it "extractRow -> insertRow : identity" $
+      insertRow (eye 3) (extractRow (eye 3) 1) 1 `shouldBe` eye 3      
     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)]
+      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
     it "modifyInspectN : termination by value convergence" $
       execState (modifyInspectN (2^16) ((< eps) . head) (/2)) 1 < eps `shouldBe` True 
-  describe "Math.Linear.Sparse : Linear solvers" $ do
+  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
     it "BiCGSTAB (2 x 2 dense)" $ 
       normSq (aa0 <\> b0 ^-^ x0true) <= eps `shouldBe` True
     it "CGS (2 x 2 dense)" $ 
       normSq (_x (cgs aa0 b0 x0 x0) ^-^ x0true) <= eps `shouldBe` True
-  describe "Math.Linear.Sparse : QR decomposition" $ do    
+  describe "Numeric.LinearAlgebra.Sparse : QR decomposition" $ do    
     it "QR (4 x 4 sparse)" $
       checkQr tm4 `shouldBe` True
     it "QR (3 x 3 dense)" $ 
       checkQr tm2 `shouldBe` True
-    
-  -- let n = 10
-  --     nsp = 3
-  -- describe ("random sparse linear system of size " ++ show n ++ " and sparsity " ++ show (fromIntegral nsp/fromIntegral n)) $ it "<\\>" $ do
-  --   aa <- randSpMat n nsp
-  --   xtrue <- randSpVec n nsp
-  --   b <- randSpVec n nsp    
-  --   let b = aa #> xtrue
-  --   printDenseSM aa
-  --   normSq (aa <\> b ^-^ xtrue) <= eps `shouldBe` True
-  -- --     normSq (_xBicgstab (bicgstab aa b x0 x0) ^-^ x) <= eps `shouldBe` True
+  describe "Numeric.LinearAlgebra.Sparse : LU decomposition" $ do
+    it "LU (3 x 3 dense)" $
+      checkLu tm2 `shouldBe` True
 
 
 
--- -- run N iterations 
-
--- runNBiC :: Int -> SpMatrix Double -> SpVector Double -> BICGSTAB
-runNBiC n aa b = map _xBicgstab $ runAppendN' (bicgstabStep aa x0) n bicgsInit where
-   x0 = mkSpVectorD nd $ replicate nd 0.9
-   nd = dim r0
-   r0 = b ^-^ (aa #> x0)    
-   p0 = r0
-   bicgsInit = BICGSTAB x0 r0 p0
-
--- runNCGS :: Int -> SpMatrix Double -> SpVector Double -> CGS
-runNCGS n aa b = map _x $ runAppendN' (cgsStep aa x0) n cgsInit where
-  x0 = mkSpVectorD nd $ replicate nd 0.1
-  nd = dim r0
-  r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution
-  p0 = r0
-  u0 = r0
-  cgsInit = CGS x0 r0 p0 u0  
-
-
 {-
 
 example 0 : 2x2 linear system
@@ -193,17 +179,7 @@
 
 
 
--- `ndim` iterations
 
-solveRandomN ndim nsp niter = do
-  aa0 <- randSpMat ndim (nsp ^ 2)
-  let aa = aa0 ^+^ eye ndim
-  xtrue <- randSpVec ndim nsp
-  let b = aa #> xtrue
-      xhatB = head $ runNBiC niter aa b
-      xhatC = head $ runNCGS niter aa b
-  printDenseSM aa    
-  return (normSq (xhatB ^-^ xtrue), normSq (xhatC ^-^ xtrue))
 
 --
 
@@ -270,6 +246,14 @@
 
 
 
+{- LU -}
+
+checkLu a = lup == a where
+  (l, u) = lu a
+  lup = l ## u
+
+
+
 {- eigenvalues -}
 
 
@@ -332,3 +316,40 @@
 --
 
 tm4 = sparsifySM $ fromListDenseSM 4 [1,0,0,0,2,5,0,10,3,6,8,11,4,7,9,12]
+
+
+
+
+
+
+
+-- -- run N iterations 
+
+-- -- runNBiC :: Int -> SpMatrix Double -> SpVector Double -> BICGSTAB
+-- runNBiC n aa b = map _xBicgstab $ runAppendN' (bicgstabStep aa x0) n bicgsInit where
+--    x0 = mkSpVectorD nd $ replicate nd 0.9
+--    nd = dim r0
+--    r0 = b ^-^ (aa #> x0)    
+--    p0 = r0
+--    bicgsInit = BICGSTAB x0 r0 p0
+
+-- -- runNCGS :: Int -> SpMatrix Double -> SpVector Double -> CGS
+-- runNCGS n aa b = map _x $ runAppendN' (cgsStep aa x0) n cgsInit where
+--   x0 = mkSpVectorD nd $ replicate nd 0.1
+--   nd = dim r0
+--   r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution
+--   p0 = r0
+--   u0 = r0
+--   cgsInit = CGS x0 r0 p0 u0
+
+
+
+-- solveRandomN ndim nsp niter = do
+--   aa0 <- randSpMat ndim (nsp ^ 2)
+--   let aa = aa0 ^+^ eye ndim
+--   xtrue <- randSpVec ndim nsp
+--   let b = aa #> xtrue
+--       xhatB = head $ runNBiC niter aa b
+--       xhatC = head $ runNCGS niter aa b
+--   -- printDenseSM aa    
+--   return (normSq (xhatB ^-^ xtrue), normSq (xhatC ^-^ xtrue))