diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -8,16 +8,18 @@
 
 Contents :
 
-* Iterative linear solvers
+* Iterative linear solvers (`linSolve`)
 
-    * BiConjugate Gradient (`bcg`)
+    * Generalized Minimal Residual (GMRES) (non-Hermitian systems) 
 
-    * Conjugate Gradient Squared (`cgs`)
+    * BiConjugate Gradient (BCG)
 
-    * BiConjugate Gradient Stabilized (`bicgstab`) (non-Hermitian systems)
+    * Conjugate Gradient Squared (CGS)
 
-    * Transpose-Free Quasi-Minimal Residual (`tfqmr`)
+    * BiConjugate Gradient Stabilized (BiCGSTAB) (non-Hermitian systems)
 
+    * Transpose-Free Quasi-Minimal Residual (TFQMR)
+
 * Direct linear solvers
 
     * LU-based (`luSolve`)
@@ -64,7 +66,11 @@
     fromListSV :: Int -> [(Int, a)] -> SpVector a
 
 can be used to create sparse vectors.
+Alternatively, the user can copy the contents of a list to a (dense) SpVector using
 
+    fromListDenseSV :: Int -> [a] -> SpVector a
+
+
 ### Displaying sparse data
 
 Both sparse vectors and matrices can be pretty-printed using `prd`:
@@ -80,7 +86,7 @@
 
 ### Matrix operations
 
-Matrix factorizations are available as `lu` and `qr` respectively, and are straightforward to verify by using the matrix product `##`  :
+There are a few common matrix factorizations available; in the following example we compute the LU factorization of a matrix and verify it with the matrix-matrix product `##`  :
 
     > (l, u) = lu amat
     > prd $ l ## u
@@ -100,12 +106,33 @@
     [4.0,3.0,2.0]
     [0.0,0.0,5.0]
 
+A matrix is transposed using `transposeSM`.
 
+Sometimes we need to compute matrix-matrix transpose products, which is why the library offers the infix operators `#^#` (M^T N) and `##^` (M N^T):
+
+    > amat' = amat #^# amat
+    > prd amat'
+    ( 3 rows, 3 columns ) , 9 NZ ( sparsity 1.0 )
+
+    [20.0,12.0,8.0]
+    [12.0,9.0,6.0]
+    [8.0,6.0,29.0]
+
+    > l = chol amat'
+    > prd $ l ##^ l
+    ( 3 rows, 3 columns ) , 9 NZ ( sparsity 1.0 )
+
+    [20.000000000000004,12.0,8.0]
+    [12.0,9.0,10.8]
+    [8.0,10.8,29.0]
+
+In the above example we have also shown the Cholesky decomposition (M = L L^T where L is a lower-triangular matrix), which is only possible for symmetric positive-definite matrices.
+
 ### Linear systems
 
-Large sparse linear systems are best solved with iterative methods. `sparse-linear-algebra` provides a selection of these via the `linSolve` function, or alternatively `<\>` (which uses BiCGSTAB as default) :
+Large sparse linear systems are best solved with iterative methods. `sparse-linear-algebra` provides a selection of these via the `linSolve` function, or alternatively `<\>` (which uses GMRES as default solver method) :
 
-    > b = fromListSV 3 [(0,3),(1,2),(2,5)]
+    > b = fromListDenseSV 3 [3,2,5]
     > x = amat <\> b
     > prd x
     ( 3 elements ) ,  3 NZ ( sparsity 1.0 )
@@ -140,7 +167,7 @@
 
 * define appropriate data structures, and how they relate to those properties (sparse vectors and matrices, defined internally via `Data.IntMap`, are made instances of the VectorSpace and Additive classes respectively). This allows to decouple the algorithms from the actual implementation of the backend,
 
-* implement the algorithms, following 1:1 the textbook [1] 
+* implement the algorithms, following 1:1 the textbook [1, 2] 
 
 
 ## License
@@ -157,3 +184,5 @@
 ## References
 
 [1] : Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., 2000
+
+[2] : L. N. Trefethen, D. Bau, Numerical Linear Algebra, SIAM, 1997
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,17 @@
 name:                sparse-linear-algebra
-version:             0.2.1.1
+version:             0.2.2.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
+description:
+  /Overview/
+  .
+  The @sparse-linear-algebra@ library provides iterative linear solvers, matrix decompositions, eigenvalue computations and related utilities. The user interface is provided by the top-level module "Numeric.LinearAlgebra.Sparse":
+  .
+  @
+  import           Numeric.LinearAlgebra.Sparse
+  @
+  .
+  Please refer to the README file for usage examples.
+  .
 homepage:            https://github.com/ocramz/sparse-linear-algebra
 license:             GPL-3
 license-file:        LICENSE
@@ -20,10 +30,11 @@
   hs-source-dirs:      src
   exposed-modules:     Numeric.LinearAlgebra.Sparse
                        Numeric.LinearAlgebra.Class
-                       Data.Sparse.IntMap2.IntMap2
                        Data.Sparse.SpVector
                        Data.Sparse.SpMatrix
                        Data.Sparse.Common
+  other-modules:       Data.Sparse.Internal.IntMap2
+                       Data.Sparse.Internal.CSR
                        Data.Sparse.Utils
                        Data.Sparse.Types
                        Numeric.Eps
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
@@ -24,7 +24,7 @@
 
 import Data.Sparse.Utils as X
 import Data.Sparse.Types as X
-import Data.Sparse.IntMap2.IntMap2 -- as X
+import Data.Sparse.Internal.IntMap2 -- as X
 import Data.Sparse.SpMatrix as X
 import Data.Sparse.SpVector as X
 
diff --git a/src/Data/Sparse/IntMap2/IntMap2.hs b/src/Data/Sparse/IntMap2/IntMap2.hs
deleted file mode 100644
--- a/src/Data/Sparse/IntMap2/IntMap2.hs
+++ /dev/null
@@ -1,194 +0,0 @@
-module Data.Sparse.IntMap2.IntMap2 where
-
-import Numeric.LinearAlgebra.Class
-
-import qualified Data.IntMap.Strict as IM
-import Data.Sparse.Types
-
-import Data.Maybe
-
-
-
-instance Set IM.IntMap where
-  liftU2 = IM.unionWith
-  {-# INLINE liftU2 #-}
-  liftI2 = IM.intersectionWith
-  {-# INLINE liftI2 #-}
-
-instance Additive IM.IntMap where
-  zero = IM.empty
-  {-# INLINE zero #-}
-  (^+^) = liftU2 (+)
-  {-# INLINE (^+^) #-}
-
-
-instance VectorSpace IM.IntMap where
-  n .* im = IM.map (* n) im
-  
-instance Hilbert IM.IntMap where
-   a `dot` b = sum $ liftI2 (*) a b
-              
-
-instance Normed IM.IntMap where
-  norm p v | p==1 = norm1 v
-           | p==2 = norm2 v
-           | otherwise = normP p v
-
-
-
--- set-like brackets
-
-
-
-
-
-
-
--- * 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
-
--- |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 #-}  
-
--- | Lookup with default 0
-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)
-
-
-
-
--- |Populate 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)
-fromListIM2 iix sm = foldl ins sm iix where
-  ins t (i,j,x) = insertIM2 i j x t
-
-
--- * Folding
-
--- |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) ->  
-  IM.IntMap a
-ifoldlIM2 f m         = IM.foldlWithKey' accRow IM.empty m where
-  accRow    acc i row = IM.foldlWithKey' (accElem i) acc row
-  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
-  accElem acc x = f x acc
-{-# inline foldlIM2 #-}
-
-
--- | 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 #-}
-
--- specialized folds
-
--- -- extract diagonal elements
--- extractDiagonalIM2 :: IM.IntMap (IM.IntMap a) -> [a]
--- extractDiagonalIM2 = ifoldlIM2' (\i j x xs -> if i==j then x : xs else xs) []
-
-
-
-
--- * Filtering
-
--- |Map over outer IM and filter all inner IM's
-ifilterIM2 ::
-  (IM.Key -> IM.Key -> a -> Bool) ->
-  IM.IntMap (IM.IntMap a) ->
-  IM.IntMap (IM.IntMap a)
-ifilterIM2 f  =
-  IM.mapWithKey (\irow row -> IM.filterWithKey (f irow) row) 
-{-# inline ifilterIM2 #-}
-
--- |Specialized filtering : keep only sub-diagonal elements
-filterSubdiag :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)
-filterSubdiag = ifilterIM2 (\i j _ -> i>j)
-
-countSubdiagonalNZ :: IM.IntMap (IM.IntMap a) -> Int
-countSubdiagonalNZ im =
-  IM.size $ IM.filter (not . IM.null) (filterSubdiag im)
-
--- |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
-
-rpairs :: (a, [b]) -> [(a, b)]
-rpairs (i, jj@(_:_)) = zip (replicate (length jj) i) jj
-rpairs (_, []) = []
-
--- -- list of (row, col) indices of elements that satisfy a criterion
--- indicesThatIM2 ::
---   (IM.Key -> IM.IntMap a -> Bool) -> IM.IntMap (IM.IntMap a) -> [(IM.Key, IM.Key)]
--- indicesThatIM2 f im = concatMap rpairs $ IM.toList (IM.map IM.keys im') where
---   im' = IM.filterWithKey f im
-
-  
-
-
--- * 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
-imapIM2 ::
-  (IM.Key -> IM.Key -> a -> b) ->
-  IM.IntMap (IM.IntMap a) ->
-  IM.IntMap (IM.IntMap b)
-imapIM2 f im = IM.mapWithKey ff im where
-  ff j x = IM.mapWithKey (`f` j) x
-
-
-
--- |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
-  adjRows = IM.mapKeys fi im
-  adjCols = IM.mapKeys fj 
-
-
-
-
--- map over a single `column`
-
-mapColumnIM2 :: (b -> b) -> IM.IntMap (IM.IntMap b) -> Int -> IM.IntMap (IM.IntMap b)
-mapColumnIM2 f im jj = imapIM2 (\i j x -> if j == jj then f x else x) im
-
-
-
-
diff --git a/src/Data/Sparse/Internal/CSR.hs b/src/Data/Sparse/Internal/CSR.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Sparse/Internal/CSR.hs
@@ -0,0 +1,34 @@
+module Data.Sparse.Internal.CSR where
+
+import qualified Data.Vector as V
+import qualified Data.Vector.Unboxed as VU
+import qualified Data.Vector.Unboxed.Mutable as VM
+
+{-| Compressed Row Storage specification :
+
+   http://netlib.org/utk/people/JackDongarra/etemplates/node373.html
+
+   The compressed row storage (CRS) format puts the subsequent nonzeros of the matrix
+   rows in contiguous memory locations. Assuming we have a nonsymmetric sparse matrix
+   $A$, we create three vectors: one for floating point numbers (val) and the other
+   two for integers (col_ind, row_ptr).
+
+   The val vector stores the values of the nonzero elements of the matrix $A$ as
+   they are traversed in a row-wise fashion.
+ 
+   The col_ind vector stores the column indexes of the elements in the val vector,
+   that is, if val(k)=a_{i,j}, then  col_ind(k)=j$.
+
+   The row_ptr vector stores the locations in the val vector that start a row;
+   that is, if  val(k)=a_{i,j}, then row_ptr(i) <= k < row_ptr(i+1)
+
+-}
+
+data CsrMatrix a =
+  CsrMatrix { csrVal :: VU.Vector a,
+              csrColInd :: VU.Vector Int,
+              csrRowPtr :: VU.Vector Int,
+              csrNnz :: {-# UNPACK #-} !Int,
+              csrNrows :: {-# UNPACK #-} !Int,
+              csrNcols :: {-# UNPACK #-} !Int } deriving Eq
+
diff --git a/src/Data/Sparse/Internal/IntMap2.hs b/src/Data/Sparse/Internal/IntMap2.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Sparse/Internal/IntMap2.hs
@@ -0,0 +1,194 @@
+module Data.Sparse.Internal.IntMap2 where
+
+import Numeric.LinearAlgebra.Class
+
+import qualified Data.IntMap.Strict as IM
+import Data.Sparse.Types
+
+import Data.Maybe
+
+
+
+instance Set IM.IntMap where
+  liftU2 = IM.unionWith
+  {-# INLINE liftU2 #-}
+  liftI2 = IM.intersectionWith
+  {-# INLINE liftI2 #-}
+
+instance Additive IM.IntMap where
+  zero = IM.empty
+  {-# INLINE zero #-}
+  (^+^) = liftU2 (+)
+  {-# INLINE (^+^) #-}
+
+
+instance VectorSpace IM.IntMap where
+  n .* im = IM.map (* n) im
+  
+instance Hilbert IM.IntMap where
+   a `dot` b = sum $ liftI2 (*) a b
+              
+
+instance Normed IM.IntMap where
+  norm p v | p==1 = norm1 v
+           | p==2 = norm2 v
+           | otherwise = normP p v
+
+
+
+-- set-like brackets
+
+
+
+
+
+
+
+-- * 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
+
+-- |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 #-}  
+
+-- | Lookup with default 0
+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)
+
+
+
+
+-- |Populate 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)
+fromListIM2 iix sm = foldl ins sm iix where
+  ins t (i,j,x) = insertIM2 i j x t
+
+
+-- * Folding
+
+-- |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) ->  
+  IM.IntMap a
+ifoldlIM2 f m         = IM.foldlWithKey' accRow IM.empty m where
+  accRow    acc i row = IM.foldlWithKey' (accElem i) acc row
+  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
+  accElem acc x = f x acc
+{-# inline foldlIM2 #-}
+
+
+-- | 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 #-}
+
+-- specialized folds
+
+-- -- extract diagonal elements
+-- extractDiagonalIM2 :: IM.IntMap (IM.IntMap a) -> [a]
+-- extractDiagonalIM2 = ifoldlIM2' (\i j x xs -> if i==j then x : xs else xs) []
+
+
+
+
+-- * Filtering
+
+-- |Map over outer IM and filter all inner IM's
+ifilterIM2 ::
+  (IM.Key -> IM.Key -> a -> Bool) ->
+  IM.IntMap (IM.IntMap a) ->
+  IM.IntMap (IM.IntMap a)
+ifilterIM2 f  =
+  IM.mapWithKey (\irow row -> IM.filterWithKey (f irow) row) 
+{-# inline ifilterIM2 #-}
+
+-- |Specialized filtering : keep only sub-diagonal elements
+filterSubdiag :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)
+filterSubdiag = ifilterIM2 (\i j _ -> i>j)
+
+countSubdiagonalNZ :: IM.IntMap (IM.IntMap a) -> Int
+countSubdiagonalNZ im =
+  IM.size $ IM.filter (not . IM.null) (filterSubdiag im)
+
+-- |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
+
+rpairs :: (a, [b]) -> [(a, b)]
+rpairs (i, jj@(_:_)) = zip (replicate (length jj) i) jj
+rpairs (_, []) = []
+
+-- -- list of (row, col) indices of elements that satisfy a criterion
+-- indicesThatIM2 ::
+--   (IM.Key -> IM.IntMap a -> Bool) -> IM.IntMap (IM.IntMap a) -> [(IM.Key, IM.Key)]
+-- indicesThatIM2 f im = concatMap rpairs $ IM.toList (IM.map IM.keys im') where
+--   im' = IM.filterWithKey f im
+
+  
+
+
+-- * 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
+imapIM2 ::
+  (IM.Key -> IM.Key -> a -> b) ->
+  IM.IntMap (IM.IntMap a) ->
+  IM.IntMap (IM.IntMap b)
+imapIM2 f im = IM.mapWithKey ff im where
+  ff j x = IM.mapWithKey (`f` j) x
+
+
+
+-- |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
+  adjRows = IM.mapKeys fi im
+  adjCols = IM.mapKeys fj 
+
+
+
+
+-- map over a single `column`
+
+mapColumnIM2 :: (b -> b) -> IM.IntMap (IM.IntMap b) -> Int -> IM.IntMap (IM.IntMap b)
+mapColumnIM2 f im jj = imapIM2 (\i j x -> if j == jj then f x else x) im
+
+
+
+
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
@@ -16,7 +16,7 @@
 import Numeric.Eps
 import Numeric.LinearAlgebra.Class
 
-import Data.Sparse.IntMap2.IntMap2
+import Data.Sparse.Internal.IntMap2
 
 import qualified Data.IntMap as IM
 
@@ -58,12 +58,13 @@
 instance Set SpMatrix where
   liftU2 f2 (SM n1 x1) (SM n2 x2) = SM (maxTup n1 n2) ((liftU2.liftU2) f2 x1 x2)
   liftI2 f2 (SM n1 x1) (SM n2 x2) = SM (minTup n1 n2) ((liftI2.liftI2) f2 x1 x2)
-  
+
+-- | 'SpMatrix'es form a ring, in that they can be added and possess a zero element  
 instance Additive SpMatrix where
   zero = SM (0,0) IM.empty
   (^+^) = liftU2 (+)
 
-
+-- | 'SpMatrix'es are maps between finite-dimensional spaces
 instance FiniteDim SpMatrix where
   type FDSize SpMatrix = (Rows, Cols)
   dim = smDim
@@ -75,6 +76,7 @@
 instance Sparse SpMatrix a where
   spy = spySM
 
+-- | 'SpMatrix'es are sparse containers too, i.e. any specific component may be missing (so it is assumed to be 0)
 instance Num a => SpContainer SpMatrix a where
   type ScIx SpMatrix = (Rows, Cols)
   scInsert (i,j) = insertSpMatrix i j
@@ -288,6 +290,12 @@
 extractSubmatrix :: SpMatrix a -> (IxRow, IxRow) -> (IxCol, IxCol) -> SpMatrix a
 extractSubmatrix = extractSubmatrixSM id id
 
+
+takeRows :: IxRow -> SpMatrix a -> SpMatrix a
+takeRows n mm = extractSubmatrix mm (0, n-1) (0, ncols mm - 1)
+
+takeCols :: IxCol -> SpMatrix a -> SpMatrix a
+takeCols n mm = extractSubmatrix mm (0, nrows mm - 1) (0, n - 1)
 
 
 
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
@@ -12,7 +12,7 @@
 
 import Data.Sparse.Utils
 import Data.Sparse.Types
-import Data.Sparse.IntMap2.IntMap2
+import Data.Sparse.Internal.IntMap2
 
 import Numeric.Eps
 import Numeric.LinearAlgebra.Class
@@ -58,11 +58,11 @@
   (^+^) = liftU2 (+)
 
 
-                      
+-- | 'SpVector's form a vector space because they can be multiplied by a scalar
 instance VectorSpace SpVector where
   n .* v = scale n v
 
-
+-- | 'SpVector's are finite-dimensional vectors
 instance FiniteDim SpVector where
   type FDSize SpVector = Int
   dim = svDim  
@@ -75,7 +75,7 @@
   spy = spySV
 
 
-
+-- | 'SpVector's are sparse containers too, i.e. any specific component may be missing (so it is assumed to be 0)
 instance Num a => SpContainer SpVector a where
   type ScIx SpVector = Int
   scInsert = insertSpVector
@@ -84,14 +84,14 @@
 
 
 
-
+-- | 'SpVector's form a Hilbert space, in that we can define an inner product over them
 instance Hilbert SpVector where
   a `dot` b | dim a == dim b = dot (dat a) (dat b)
             | otherwise =
                      error $ "dot : sizes must coincide, instead we got " ++
                            show (dim a, dim b)
 
-
+-- | Since 'SpVector's form a Hilbert space, we can define a norm for them 
 instance Normed SpVector where
   norm p (SV _ v) = norm p v
 
@@ -104,17 +104,22 @@
 
 -- ** Creation
 
--- | empty sparse vector (length n, no entries)
+-- | Empty sparse vector (length n, no entries)
 zeroSV :: Int -> SpVector a
 zeroSV n = SV n IM.empty
 
 
--- | singleton sparse vector (length 1)
+-- | Singleton sparse vector (length 1)
 singletonSV :: a -> SpVector a
 singletonSV x = SV 1 (IM.singleton 0 x)
 
 
+-- | Canonical basis vector in R^n
+ei :: Num a => Int -> IM.Key -> SpVector a
+ei n i = SV n (IM.insert (i - 1) 1 IM.empty)
 
+
+
 -- | create a sparse vector from an association list while discarding all zero entries
 mkSpVector :: Epsilon a => Int -> IM.IntMap a -> SpVector a
 mkSpVector d im = SV d $ IM.filterWithKey (\k v -> isNz v && inBounds0 d k) im
@@ -203,6 +208,10 @@
 fromListSV :: Int -> [(Int, a)] -> SpVector a
 fromListSV d iix = SV d (IM.fromList (filter (inBounds0 d . fst) iix ))
 
+-- fromListSV' d iix = SV d (F.foldl' insf IM.empty iix') where
+--   insf im (i, x) = IM.insert i x im
+--   iix' = filter (inBounds0 d . fst) iix   -- filtering forces list as only instance
+ 
 -- ** toList
 toListSV :: SpVector a -> [(IM.Key, a)]
 toListSV sv = IM.toList (dat sv)
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
@@ -3,7 +3,8 @@
 module Numeric.LinearAlgebra.Sparse
        (
          -- * Matrix factorizations
-         qr, lu,
+         qr,
+         lu,
          chol,
          -- * Condition number
          conditionNumberSM,
@@ -12,18 +13,15 @@
          -- * Givens' rotation
          givens,
          -- * Arnoldi iteration
-         arnoldi,
+         arnoldi, 
          -- * Eigensolvers
          eigsQR, eigRayleigh,
          -- * Linear solvers
-         linSolve, LinSolveMethod, (<\>),
+         -- ** Iterative methods
+         linSolve, LinSolveMethod(..), (<\>),
+         pinv,
          -- ** Direct methods
          luSolve,
-         -- ** Iterative methods
-         cgne, tfqmr, bicgstab, cgs, bcg,
-         _xCgne, _xTfq, _xBicgstab, _x, _xBcg,
-         cgsStep, bicgstabStep,
-         CGNE, TFQMR, BICGSTAB, CGS, BCG,
          -- * Preconditioners
          ilu0, mSsor,
          -- * Matrix partitioning
@@ -49,7 +47,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)
@@ -85,7 +83,7 @@
 
 -- |uses the R matrix from the QR factorization
 conditionNumberSM :: (Epsilon a, RealFloat a) => SpMatrix a -> a
-conditionNumberSM m | isInfinite kappa = error "Infinite condition number : rank-deficient system"
+conditionNumberSM m | nearZero lmin = error "Infinite condition number : rank-deficient system"
                     | otherwise = kappa where
   kappa = lmax / lmin
   (_, r) = qr m
@@ -162,8 +160,8 @@
 {-# inline givens #-}
 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))
+  | isValidIxSM mm (i,j) && nrows mm >= ncols mm =
+       fromListSM' [(i,i,c),(j,j,c),(j,i,-s),(i,j,s)] (eye (nrows mm))
   | otherwise = error "givens : indices out of bounds"      
   where
     (c, s, _) = givensCoef a b
@@ -188,40 +186,33 @@
 
 
 
--- * QR decomposition
 
 
--- | 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 :: (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
- gm _ [] = []
+-- * QR decomposition
 
 
+-- | Given a matrix A, returns a pair of matrices (Q, R) such that Q R = A, where Q is orthogonal and R is upper triangular. Applies Givens rotation iteratively to zero out sub-diagonal elements.
+qr :: (Epsilon a, Ord a, Floating a) => SpMatrix a -> (SpMatrix a, SpMatrix a)
+qr mm = (transposeSM qt, r) where
+  (qt, r, _) = execState (modifyUntil qf stepf) gminit
+  qf (_, _, iis) = null iis
+  stepf (qmatt, m, iis) = (qmatt', m', tail iis) where
+    (i, j) = head iis
+    g = givens m i j
+    qmatt' = g #~# qmatt  -- update Q'
+    m' = g #~# m          -- update R
+  gminit = (eye (nrows mm), mm, subdiagIndicesSM mm)
 
+    
+  
 
 
--- -- | QR algorithm, state transformer version
--- gmatST0 (m, (i,j):is) = (m', is) where    -- WRONG, possible access to []
---   g = givens m i j                        
---   m' = g #~# m
--- gmatST0 (m, []) = (eye (nrows m), [])
 
--- gmatST m = gmatST0 (m, subdiagIndicesSM m)
 
 
 
 
 
-
 -- * Eigenvalue algorithms
 
 -- ** QR algorithm
@@ -353,7 +344,7 @@
 -- ** 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 -}
 
--- | Given a matrix A, returns a pair of matrices (L, U) such that L U = A
+-- | Given a matrix A, returns a pair of matrices (L, U) where L is lower triangular and U is upper triangular 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
@@ -429,61 +420,48 @@
 
 
 
--- * Incomplete LU
 
--- | 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
-  uh = sparsifyLU u aa
-  sparsifyLU m m2 = ifilterSM f m where
-    f i j _ = isJust (lookupSM m2 i j)
 
 
-
-
-
-
 -- * Arnoldi iteration
 
--- | Given a matrix A and a positive integer `n`, this procedure finds the basis of an order `n` Krylov subspace (as the columns of matrix Q), along with an upper Hessenberg matrix H, such that A = Q^T H Q.
+-- | Given a matrix A, a vector b and a positive integer `n`, this procedure finds the basis of an order `n` Krylov subspace (as the columns of matrix Q), along with an upper Hessenberg matrix H, such that A = Q^T H Q.
 -- At the i`th iteration, it finds (i + 1) coefficients (the i`th column of the Hessenberg matrix H) and the (i + 1)`th Krylov vector.
+
 arnoldi ::
-  (Floating a, Eq a) => SpMatrix a -> Int -> (SpMatrix a, SpMatrix a)
-arnoldi aa kn = (fromCols qvfin, hhfin) where
-  (qvfin, hhfin, _) = execState (modifyUntil tf arnoldiStep) arnInit 
-  tf (_, _, ii) = ii == kn -- termination criterion
+  (Epsilon a, Floating a, Eq a) => SpMatrix a -> SpVector a -> Int -> (SpMatrix a, SpMatrix a)
+arnoldi aa b kn = (fromCols qvfin, fromListSM (nmax + 1, nmax) hhfin)
+  where
+  (qvfin, hhfin, nmax, _) = execState (modifyUntil tf arnoldiStep) arnInit 
+  tf (_, _, ii, fbreak) = ii == kn || fbreak -- termination criterion
   (m, n) = dim aa
-  arnInit = (qv1, hh1, 1) where      
-      q0 = normalize 2 $ onesSV n -- starting basis vector
-      aq0 = aa #> q0
-      h11 = q0 `dot` aq0
+  arnInit = (qv1, hh1, 1, False) where
+      q0 = normalize 2 b   -- starting basis vector
+      aq0 = aa #> q0       -- A q0
+      h11 = q0 `dot` aq0          
       q1nn = (aq0 ^-^ (h11 .* q0))
-      hh1 = fromListSM (m + 1, n) [(0, 0, h11), (1, 0, h21)] where        
+      hh1 = V.fromList [(0, 0, h11), (1, 0, h21)] where        
         h21 = norm 2 q1nn
-      q1 = normalize 2 q1nn
+      q1 = normalize 2 q1nn       -- q1 `dot` q0 ~ 0
       qv1 = V.fromList [q0, q1]
-  arnoldiStep (qv, hh, i) = (qv', hh', i + 1)
-   where
+  arnoldiStep (qv, hh, i, _) = (qv', hh', i + 1, fb') where
     qi = V.last qv
     aqi = aa #> qi
     hhcoli = fmap (`dot` aqi) qv -- H_{1, i}, H_{2, i}, .. , H_{m + 1, i}
     zv = zeroSV m
     qipnn =
       aqi ^-^ (V.foldl' (^+^) zv (V.zipWith (.*) hhcoli qv)) -- unnormalized q_{i+1}
-    qipnorm = singletonSV $ norm 2 qipnn      -- normalization factor H_{i+1, i}
+    qipnorm = norm 2 qipnn      -- normalization factor H_{i+1, i}
     qip = normalize 2 qipnn              -- q_{i + 1}
-    hh' = insertCol hh (concatSV (fromVector hhcoli) qipnorm) i -- update H
+    hh' = (V.++) hh (indexed2 $ V.snoc hhcoli qipnorm) where -- update H
+      indexed2 v = V.zip3 ii jj v
+      ii = V.fromList [0 .. n]    -- nth col of upper Hessenberg has `n+1` nz
+      jj = V.replicate (n + 1) i  -- `n+1` replicas of `i`
     qv' = V.snoc qv qip        -- append q_{i+1} to Krylov basis Q_i
-
-
-  
-
-
-
+    fb' | nearZero qipnorm = True  -- breakdown condition
+        | otherwise = False
 
 
 
@@ -501,6 +479,25 @@
   f = extractSuperDiag aa
 
 
+-- -- ** Jacobi preconditioner
+
+-- -- | Returns the reciprocal of the diagonal 
+-- jacobiPreconditioner :: SpMatrix Double -> SpMatrix Double
+-- jacobiPreconditioner = reciprocal . extractDiag
+
+
+-- ** Incomplete LU
+
+-- | 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
+  uh = sparsifyLU u aa
+  sparsifyLU m m2 = ifilterSM f m where
+    f i j _ = isJust (lookupSM m2 i j)
+
+
 -- ** SSOR
 
 -- | `mSsor aa omega` : if `omega = 1` it returns the symmetric Gauss-Seidel preconditioner
@@ -524,18 +521,19 @@
 luSolve ::
   (Fractional a, Eq a, Epsilon a) => SpMatrix a -> SpMatrix a -> SpVector a -> SpVector a
 luSolve ll uu b
-  | isLowerTriSM ll && isUpperTriSM uu = lubwSolve uu (lufwSolve ll b)
+  | isLowerTriSM ll && isUpperTriSM uu = triUpperSolve uu (triLowerSolve ll b)
   | otherwise = error "luSolve : factors must be triangular matrices" 
 
-lufwSolve ll b = sparsifySV v where
+triLowerSolve :: (Epsilon a, Fractional a) => SpMatrix a -> SpVector a -> SpVector a
+triLowerSolve ll b = sparsifySV v where
   (v, _) = execState (modifyUntil q lStep) lInit where
   q (_, i) = i == dim b
-  lStep (ww, i) = (wwi, i + 1) where
+  lStep (ww, i) = (ww', i + 1) where
     lii = ll @@ (i, i)
     bi = b @@ i
     wi = (bi - r)/lii where
       r = extractSubRow ll i (0, i-1) `dot` takeSV i ww
-    wwi = insertSpVector i wi ww
+    ww' = insertSpVector i wi ww
   lInit = (ww0, 1) where
     l00 = ll @@ (0, 0)
     b0 = b @@ 0
@@ -543,15 +541,16 @@
     ww0 = insertSpVector 0 w0 $ zeroSV (dim b)  
 
 -- | NB in the computation of `xi` we must rebalance the subrow indices because `dropSV` does that too, in order to take the inner product with consistent index pairs
-lubwSolve uu w = sparsifySV x where
+triUpperSolve :: (Epsilon a, Fractional a) => SpMatrix a -> SpVector a -> SpVector a
+triUpperSolve uu w = sparsifySV x where
   (x, _) = execState (modifyUntil q uStep) uInit
   q (_, i) = i == (- 1)
-  uStep (xx, i) = (xxi, i - 1) where
+  uStep (xx, i) = (xx', i - 1) where
     uii = uu @@ (i, i)
     wi = w @@ i
     xi = (wi - r) / uii where
         r = extractSubRow_RK uu i (i + 1, dim w - 1) `dot` dropSV (i + 1) xx
-    xxi = insertSpVector i xi xx
+    xx' = insertSpVector i xi xx
   uInit = (xx0, i - 1) where
     i = dim w - 1
     u00 = uu @@ (i, i)
@@ -569,12 +568,42 @@
 
 -- ** GMRES
 
+-- | Given a linear system `A x = b` where `A` is an (m x m) real-valued matrix, the GMRES method finds an approximate solution `xhat` such that the Euclidean norm of the residual `A xhat - b` is minimized. `xhat` is spanned by the order-`n` Krylov subspace of (A, b).
+-- In this implementation:
+-- 1) the Arnoldi iteration is carried out until numerical breakdown (therefore yielding _at_most_ `m` Krylov basis vectors)
+-- 2) the resulting Hessenberg matrix is factorized in QR form
+-- 3) the Krylov-subspace solution `yhat` is found by backsubstitution
+-- 4) the approximate solution in the original space `xhat` is computed using the Krylov basis, `xhat = Q_n yhat`
+--
+-- Many optimizations are possible, for example interleaving the QR factorization with the Arnoldi process (and employing an updating QR factorization which only requires one Givens' rotation at every update). 
+
+gmres :: (Epsilon a, Ord a, Floating a) => SpMatrix a -> SpVector a -> SpVector a
+gmres aa b = qa' #> yhat where
+  m = ncols aa
+  (qa, ha) = arnoldi aa b m   -- at most m steps of Arnoldi (aa, b)
+  -- b' = (transposeSe qa) #> b
+  b' = norm2 b .* (ei mp1 1)  -- b rotated back to canonical basis by Q^T
+     where mp1 = nrows ha     -- 1 + # Arnoldi iterations
+  (qh, rh) = qr ha            -- QR fact.of H
+  yhat = triUpperSolve rh' rhs' where
+    rhs' = takeSV (dim b' - 1) (transposeSM qh #> b')
+    rh' = takeRows (nrows rh - 1) rh
+  qa' = takeCols (ncols qa - 1) qa
+
+
+
+
+
+
+
 -- *** Left-preconditioning
 
 
 
 
 
+
+
 -- ** CGNE
 
 cgneStep :: SpMatrix Double -> CGNE -> CGNE
@@ -587,6 +616,10 @@
 
 data CGNE =
   CGNE {_xCgne , _rCgne, _pCgne :: SpVector Double} deriving Eq
+instance Show CGNE where
+    show (CGNE x r p) = "x = " ++ show x ++ "\n" ++
+                       "r = " ++ show r ++ "\n" ++
+                       "p = " ++ show p ++ "\n"
 
 cgne :: SpMatrix Double -> SpVector Double -> SpVector Double -> CGNE
 cgne aa b x0 = execState (untilConverged _xCgne (cgneStep aa)) cgneInit where
@@ -645,7 +678,8 @@
           _mTfq :: Int,
           _tauTfq, _thetaTfq, _etaTfq, _rhoTfq, _alphaTfq :: Double}
   deriving Eq
-
+instance Show TFQMR where
+    show (TFQMR x _ _ _ _ _ _ _ _ _ _) = "x = " ++ show x ++ "\n"
 
 
 -- ** BCG
@@ -782,7 +816,7 @@
 
 -- * Linear solver interface
 
-data LinSolveMethod = CGNE_ | TFQMR_ | BCG_ | CGS_ | BICGSTAB_ deriving (Eq, Show) 
+data LinSolveMethod = GMRES_ | CGNE_ | TFQMR_ | BCG_ | CGS_ | BICGSTAB_ deriving (Eq, Show) 
 
 -- -- | Linear solve with _random_ starting vector
 -- linSolveM ::
@@ -805,18 +839,19 @@
       solve aa' b' | isDiagonalSM aa' = reciprocal aa' #> b' -- diagonal solve is easy
                    | otherwise = solveWith aa' b' 
       solveWith aa' b' = case method of
+        GMRES_ -> gmres aa' b'
         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)
+        BICGSTAB_ ->  _xBicgstab (bicgstab aa' b' x0 x0)
+        CGS_ -> _x (cgs aa' b' x0 x0)
       x0 = mkSpVectorD n $ replicate n 0.1 
       (m, n) = dim aa
       nb     = dim b
 
 -- | <\\> : linSolve using the BiCGSTAB method as default
 (<\>) :: SpMatrix Double -> SpVector Double -> SpVector Double      
-(<\>) = linSolve BICGSTAB_ 
+(<\>) = linSolve GMRES_ 
   
 
 
diff --git a/test/LibSpec.hs b/test/LibSpec.hs
--- a/test/LibSpec.hs
+++ b/test/LibSpec.hs
@@ -13,14 +13,13 @@
 import Numeric.LinearAlgebra.Sparse
 -- import Numeric.LinearAlgebra.Class
 
-import Control.Monad (liftM, liftM2)
+import Control.Applicative (liftA2)
+-- import Control.Monad (liftM, liftM2, replicateM)
 import Control.Monad.Primitive
 import Data.Foldable (foldrM)
 
 import Data.Sparse.Common
 
-
-import Control.Monad (replicateM)
 import Control.Monad.State.Strict (execState)
 
 import qualified System.Random.MWC as MWC
@@ -55,9 +54,9 @@
     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)])
+      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)])
+      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" $
@@ -77,61 +76,76 @@
   describe "Numeric.LinearAlgebra.Sparse : Iterative linear solvers" $ do
     -- it "TFQMR (2 x 2 dense)" $
     --   normSq (_xTfq (tfqmr aa0 b0 x0) ^-^ x0true) <= eps `shouldBe` True
+    it "GMRES (2 x 2 dense)" $
+      nearZero (normSq (linSolve GMRES_ aa0 b0 ^-^ x0true)) `shouldBe` True
+    it "GMRES (3 x 3 sparse, s.p.d.)" $
+      nearZero (normSq (linSolve GMRES_ aa2 b2 ^-^ x2)) `shouldBe` True
+    it "GMRES (4 x 4 sparse)" $
+      nearZero (normSq (linSolve GMRES_ aa1 b1 ^-^ x1)) `shouldBe` True      
     it "BCG (2 x 2 dense)" $
-      nearZero (normSq (_xBcg (bcg aa0 b0 x0) ^-^ x0true)) `shouldBe` True
-    it "BiCGSTAB (2 x 2 dense)" $ 
-      nearZero (normSq (aa0 <\> b0 ^-^ x0true)) `shouldBe` True
+      nearZero (normSq (linSolve BCG_ aa0 b0 ^-^ x0true)) `shouldBe` True
+    it "BCG (3 x 3 sparse, s.p.d.)" $
+      nearZero (normSq (linSolve BCG_ aa2 b2 ^-^ x2)) `shouldBe` True      
+    -- it "BiCGSTAB (2 x 2 dense)" $ 
+    --   nearZero (normSq (linSolve BICGSTAB_ aa0 b0 ^-^ x0true)) `shouldBe` True
+    it "BiCGSTAB (3 x 3 sparse, s.p.d.)" $ 
+      nearZero (normSq (linSolve BICGSTAB_ aa2 b2 ^-^ x2)) `shouldBe` True      
     it "CGS (2 x 2 dense)" $ 
-      nearZero (normSq (_x (cgs aa0 b0 x0 x0) ^-^ x0true)) `shouldBe` True
-  describe "Numeric.LinearAlgebra.Sparse : Direct linear solvers" $ do
-    it "LU (unoptimized) (4 x 4 sparse)" $ 
+      nearZero (normSq (linSolve CGS_ aa0 b0 ^-^ x0true)) `shouldBe` True
+    it "CGS (3 x 3 sparse, s.p.d.)" $ 
+      nearZero (normSq (linSolve CGS_ aa2 b2 ^-^ x2)) `shouldBe` True      
+  describe "Numeric.LinearAlgebra.Sparse : Direct linear solvers" $ 
+    it "luSolve (4 x 4 sparse)" $ 
       checkLuSolve aa1 b1 `shouldBe` True         
   describe "Numeric.LinearAlgebra.Sparse : QR decomposition" $ do    
-    it "QR (4 x 4 sparse)" $
+    it "qr (4 x 4 sparse)" $
       checkQr tm4 `shouldBe` True
-    it "QR (3 x 3 dense)" $ 
+    it "qr (3 x 3 dense)" $ 
       checkQr tm2 `shouldBe` True
   describe "Numeric.LinearAlgebra.Sparse : LU decomposition" $ do
-    it "LU (4 x 4 dense)" $
+    it "lu (4 x 4 dense)" $
       checkLu tm6 `shouldBe` True
-    it "LU (10 x 10 sparse)" $
+    it "lu (10 x 10 sparse)" $
       checkLu tm7 `shouldBe` True
-  describe "Numeric.LinearAlgebra.Sparse : Cholesky decomposition (PSD matrices only)" $ 
+  describe "Numeric.LinearAlgebra.Sparse : Cholesky decomposition (PSD matrices)" $ 
     it "chol (5 x 5 sparse)" $
       checkChol tm7 `shouldBe` True
-  describe "Numeric.LinearAlgebra.Sparse : Arnoldi iteration" $ do
-    it "Arnoldi iteration (3 x 3 dense)" $
-      checkArnoldi aa2 3 `shouldBe` True
-    it "Arnoldi iteration (5 x 5 sparse)" $
+  describe "Numeric.LinearAlgebra.Sparse : Arnoldi iteration, early breakdown detection" $ do      
+    it "arnoldi (4 x 4 dense)" $
+      checkArnoldi tm6 4 `shouldBe` True
+    it "arnoldi (5 x 5 sparse)" $
       checkArnoldi tm7 5 `shouldBe` True    
 
 
-
 {- QR-}
 
 
 checkQr :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool
-checkQr a = c1 && c2 where
+checkQr a = c1 && c2 && c3 where
   (q, r) = qr a
   c1 = nearZero $ normFrobenius ((q #~# r) ^-^ a)
   c2 = isOrthogonalSM q
+  c3 = isUpperTriSM r
 
 
 
 {- LU -}
 
 checkLu :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool
-checkLu a = lup == a where
+checkLu a = c1 && c2 where
   (l, u) = lu a
-  lup = l #~# u
+  c1 = nearZero (normFrobenius ((l #~# u) ^-^ a))
+  c2 = isUpperTriSM u && isLowerTriSM l
 
 
 
 {- Cholesky -}
 
 checkChol :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool
-checkChol a = nearZero $ normFrobenius ((l ##^ l) ^-^ a) where
+checkChol a = c1 && c2 where
   l = chol a
+  c1 = nearZero $ normFrobenius ((l ##^ l) ^-^ a)
+  c2 = isLowerTriSM l
 
 
 {- direct linear solver -}
@@ -145,12 +159,13 @@
   
 {- Arnoldi iteration -}
 checkArnoldi :: (Epsilon a, Floating a, Eq a) => SpMatrix a -> Int -> Bool
-checkArnoldi aa kn = nearZero $ normFrobenius $ (aa #~# qvprev) ^-^ (qv #~# hh) where
-  (qv, hh) = arnoldi aa kn
-  (m, n) = dim qv
-  qvprev = extractSubmatrix qv (0, m - 1) (0, n - 2)
-
-
+checkArnoldi aa kn = nearZero (normFrobenius $ lhs ^-^ rhs) where
+  b = onesSV (nrows aa)
+  (q, h) = arnoldi aa b kn
+  (m, n) = dim q
+  q' = extractSubmatrix q (0, m - 1) (0, n - 2) -- q' = all but one column of q
+  rhs = q #~# h
+  lhs = aa #~# q'
 
 
 {-
@@ -258,19 +273,20 @@
 --   let
 --     diags = [-bw .. bw - 1]
 
-randDiagMat :: PrimMonad m =>
-     Rows -> Double -> Double -> Int -> m (SpMatrix Double)
-randDiagMat n mu sig i = do
-  x <- randArray n mu sig
-  return $ mkSubDiagonal n i x
+-- randDiagMat :: PrimMonad m =>
+--      Rows -> Double -> Double -> Int -> m (SpMatrix Double)
+-- randDiagMat n mu sig i = do
+--   x <- randArray n mu sig
+--   return $ mkSubDiagonal n i x
 
 
-go (m:ms) mat =
-  m ^+^ go ms mat
-go [] mat = mat
+-- go (m:ms) mat =
+--   m ^+^ go ms mat
+-- go [] mat = mat
 
-  
-plusM x y = return $ x ^+^ y
+-- plusM ::
+--   (Additive f1, Applicative f, Num a) => f (f1 a) -> f (f1 a) -> f (f1 a)  
+-- plusM = liftA2 (^+^)
 
 
 
@@ -397,6 +413,11 @@
   b = mkSubDiagonal n 0 $ replicate n 2
   c = mkSubDiagonal n (-1) $ replicate n (-1)
 
+
+
+
+
+
 -- -- run N iterations 
 
 -- -- runNBiC :: Int -> SpMatrix Double -> SpVector Double -> BICGSTAB
@@ -434,3 +455,8 @@
 
 tm8' :: SpMatrix Double
 tm8' = fromListSM (2,2) [(0,0,1), (1,0,1), (1,1,1)]
+
+
+
+tm9 :: SpMatrix Double
+tm9 = fromListSM (4, 3) [(0,0,pi), (1,1, 3), (2,2,4), (3,2, 1), (3,1, 5)]