diff --git a/app/Main.hs b/app/Main.hs
deleted file mode 100644
--- a/app/Main.hs
+++ /dev/null
@@ -1,8 +0,0 @@
-module Main where
-
-import Lib (ourAdd)
-
-import Text.Printf (printf)
-
-main :: IO ()
-main = printf "2 + 3 = %d\n" (ourAdd 2 3)
diff --git a/sparse-linear-algebra.cabal b/sparse-linear-algebra.cabal
--- a/sparse-linear-algebra.cabal
+++ b/sparse-linear-algebra.cabal
@@ -1,14 +1,14 @@
 name:                sparse-linear-algebra
-version:             0.1.0.0
-synopsis:            Sparse linear algebra datastructures and algorithms
+version:             0.1.0.1
+synopsis:            Sparse linear algebra datastructures and algorithms. Currently it provides iterative linear solvers, matrix decompositions, eigenvalue computations and related utilities.
 description:         Please see README.md
 homepage:            https://github.com/ocramz/sparse-linear-algebra
-license:             BSD3
+license:             GPL-3
 license-file:        LICENSE
 author:              Marco Zocca
 maintainer:          zocca.marco gmail
 copyright:           2016 Marco Zocca
-category:            Math
+category:            Numeric
 build-type:          Simple
 extra-source-files:  README.md
 cabal-version:       >=1.10
@@ -26,29 +26,28 @@
                      , containers
                      , hspec
                      , primitive >= 0.6.1.0
-                     , transformers >= 0.5.2.0
-                     -- , lens
                      , mtl >= 2.2.1
                      , mwc-random
-                     , monad-loops
 
-executable sparse-linear-algebra
-  default-language:    Haskell2010
-  ghc-options:         -threaded -rtsopts -with-rtsopts=-N
-  hs-source-dirs:      app
-  main-is:             Main.hs
-  build-depends:       base
-                     , mtl >= 2.2.1
-                     , mwc-random
-                     , primitive >= 0.6.1.0
-                     , sparse-linear-algebra
-                     , transformers >= 0.5.2.0
 
+-- executable sparse-linear-algebra
+--   default-language:    Haskell2010
+--   ghc-options:         -threaded -rtsopts -with-rtsopts=-N
+--   hs-source-dirs:      app
+--   main-is:             Main.hs
+--   build-depends:       base
+--                      , mtl >= 2.2.1
+--                      , mwc-random
+--                      , primitive >= 0.6.1.0
+--                      , sparse-linear-algebra
+--                      , transformers >= 0.5.2.0
+
 test-suite spec
   default-language:    Haskell2010
   ghc-options:         -Wall
   type:                exitcode-stdio-1.0
   hs-source-dirs:      test
+  other-modules:       LibSpec
   main-is:             Spec.hs
   build-depends:       base
                      , containers
@@ -57,7 +56,6 @@
                      , mwc-random
                      , primitive >= 0.6.1.0
                      , sparse-linear-algebra
-                     , transformers >= 0.5.2.0
                      , criterion
                      -- , QuickCheck
 
diff --git a/src/Math/Linear/Sparse.hs b/src/Math/Linear/Sparse.hs
--- a/src/Math/Linear/Sparse.hs
+++ b/src/Math/Linear/Sparse.hs
@@ -9,15 +9,13 @@
 import Control.Monad.Primitive
 
 import Control.Monad (mapM_, forM_, replicateM)
-import Control.Monad.Loops
 
-import Control.Monad.Cont
 import Control.Monad.State.Strict
-import Control.Monad.Writer
-import Control.Monad.Trans
+-- import Control.Monad.Writer
+-- import Control.Monad.Trans
 
-import Control.Monad.Trans.State (runStateT)
-import Control.Monad.Trans.Writer (runWriterT)
+-- import Control.Monad.Trans.State (runStateT)
+-- import Control.Monad.Trans.Writer (runWriterT)
 
 import qualified Data.IntMap.Strict as IM
 -- import Data.Utils.StrictFold (foldlStrict) -- hidden in `containers`
@@ -30,33 +28,35 @@
 import qualified Data.Traversable as T
 
 
-import qualified Data.List as L
+-- import qualified Data.List as L
 import Data.Maybe
 
 
 
--- | ========= CLASSES and common operations
+{-|  CLASSES and common operations -}
 
--- | Additive ring 
+-- * Additive ring 
 class Functor f => Additive f where
-  -- | zero element
+  -- | Ring zero element
   zero :: Num a => f a
   
-  -- | componentwise operations
+  -- | Ring +
   (^+^) :: Num a => f a -> f 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
 
-x `minus` y = x ^+^ negated y
+-- | subtract two Additive objects
+(^-^) :: (Additive f, Num a) => f a -> f a -> f a
+x ^-^ y = x ^+^ negated y
 
 
 
--- | Vector space
+-- * Vector space
 class Additive f => VectorSpace f where
   -- | multiplication by a scalar
   (.*) :: Num a => a -> f a -> f a
@@ -67,44 +67,45 @@
 lerp a u v = a .* u ^+^ ((1-a) .* v)
 
 
--- | Hilbert space (inner product)
+-- * Hilbert space (inner product)
 class VectorSpace f => Hilbert f where
   -- | inner product
   dot :: Num a => f a -> f a -> a
 
+-- * Normed vector space
 class Hilbert f => Normed f where
   norm :: (Floating a, Eq a) => a -> f a -> a
 
 
--- some norms and related results
+-- ** Norms and related results
 
--- squared norm 
+-- *** squared 2-norm
 normSq :: (Hilbert f, Num a) => f a -> a
 normSq v = v `dot` v
 
 
--- L1 norm
+-- *** L1 norm
 norm1 :: (Foldable t, Num a, Functor t) => t a -> a
 norm1 v = sum (fmap abs v)
 
--- Euclidean norm
+-- *** Euclidean norm
 norm2 :: (Hilbert f, Floating a) => f a -> a
 norm2 v = sqrt (normSq v)
 
--- Lp norm (p > 0)
+-- *** Lp norm (p > 0)
 normP :: (Foldable t, Functor t, Floating a) => a -> t a -> a
 normP p v = sum u**(1/p) where
   u = fmap (**p) v
 
--- infinity-norm
+-- *** infinity-norm
 normInfty :: (Foldable t, Ord a) => t a -> a
 normInfty = maximum
 
 
 
--- normalize
+-- *** normalize w.r.t. p-norm (p finite)
 normalize :: (Normed f, Floating a, Eq a) => a -> f a -> f a
-normalize n v = (1 / norm n v) .* v
+normalize p v = (1 / norm p v) .* v
 
 
 
@@ -112,19 +113,19 @@
 
 
 
--- -- Lp inner product (p > 0)
+-- *** Lp inner product (p > 0)
 dotLp :: (Set t, Foldable t, Floating a) => a -> t a -> t a ->  a
 dotLp p v1 v2 = sum u**(1/p) where
   f a b = (a*b)**p
   u = liftI2 f v1 v2
 
 
--- reciprocal
+-- *** reciprocal
 reciprocal :: (Functor f, Fractional b) => f b -> f b
 reciprocal = fmap recip
 
 
--- scale
+-- *** scale
 scale :: (Num b, Functor f) => b -> f b -> f b
 scale n = fmap (* n)
 
@@ -134,24 +135,20 @@
 
 
 
--- | FiniteDim : finite-dimensional objects
+-- ** FiniteDim : finite-dimensional objects
 
 class Additive f => FiniteDim f where
   type FDSize f :: *
   dim :: f a -> FDSize f
 
 
-instance FiniteDim SpVector where
-  type FDSize SpVector = Int
-  dim = svDim
 
 
-instance FiniteDim SpMatrix where
-  type FDSize SpMatrix = (Rows, Cols)
-  dim = smDim
 
 
--- unary dimension-checking bracket
+
+
+-- | unary dimension-checking bracket
 withDim :: (FiniteDim f, Show e) =>
      f a
      -> (FDSize f -> f a -> Bool)
@@ -162,7 +159,7 @@
 withDim x p f e ef | p (dim x) x = f x
                    | otherwise = error e' where e' = e ++ show (ef x)
 
--- binary dimension-checking bracket
+-- | binary dimension-checking bracket
 withDim2 :: (FiniteDim f, FiniteDim g, Show e) =>
      f a
      -> g b
@@ -179,40 +176,28 @@
 
 
 
--- | HasData : accessing inner data (do not export)
+-- ** HasData : accessing inner data (do not export)
 
 class Additive f => HasData f a where
   type HDData f a :: * 
   dat :: f a -> HDData f a
 
-instance HasData SpVector a where
-  type HDData SpVector a = IM.IntMap a
-  dat = svData
 
-instance HasData SpMatrix a where
-  type HDData SpMatrix a = IM.IntMap (IM.IntMap a)
-  dat = smData
 
 
-
-
-
-
--- | Sparse : sparse datastructures
+-- ** Sparse : sparse datastructures
 
 class (FiniteDim f, HasData f a) => Sparse f a where
   spy :: Fractional b => f a -> b
 
 
-instance Sparse SpVector a where
-  spy = spySV
 
-instance Sparse SpMatrix a where
-  spy = spySM
 
 
 
 
+-- ** Set : things that behave as sets (e.g. of which we can take the union and the intersection)
+
 class Functor f => Set f where
   -- |union binary lift
   liftU2 :: (a -> a -> a) -> f a -> f a -> f a
@@ -228,7 +213,7 @@
 
 
 
--- | =======================================================
+-- * IntMap implementation
 
 instance Set IM.IntMap where
   liftU2 = IM.unionWith
@@ -236,15 +221,13 @@
   liftI2 = IM.intersectionWith
   {-# INLINE liftI2 #-}
 
--- | IntMap implementation
 instance Additive IM.IntMap where
   zero = IM.empty
   {-# INLINE zero #-}
   (^+^) = liftU2 (+)
   {-# INLINE (^+^) #-}
-  x ^-^ y = x ^+^ negated y
-  {-# INLINE (^-^) #-}
 
+
 instance VectorSpace IM.IntMap where
   n .* im = IM.map (* n) im
   
@@ -259,27 +242,22 @@
 
 
 
--- | =======================================================
+-- * Sparse Vector
 
--- | Sparse Vector
 data SpVector a = SV { svDim :: Int ,
                        svData :: IM.IntMap a} deriving Eq
 
-dimSV :: SpVector a -> Int
-dimSV = svDim
-
+-- | SpVector sparsity
 spySV :: Fractional b => SpVector a -> b
-spySV s = fromIntegral (IM.size (dat s)) / fromIntegral (svDim s)
+spySV s = fromIntegral (IM.size (dat s)) / fromIntegral (dim s)
 
 
--- internal : projection functions, do not export
-imSV :: SpVector a -> IM.IntMap a
-imSV = svData
 
 
 
 
--- | instances for SpVector
+
+-- ** instances for SpVector
 instance Functor SpVector where
   fmap f (SV n x) = SV n (fmap f x)
 
@@ -293,12 +271,25 @@
 instance Additive SpVector where
   zero = SV 0 IM.empty
   (^+^) = liftU2 (+)
-  (^-^) = liftU2 (-)
 
+
                       
 instance VectorSpace SpVector where
   n .* v = scale n v
 
+
+instance FiniteDim SpVector where
+  type FDSize SpVector = Int
+  dim = svDim  
+
+instance HasData SpVector a where
+  type HDData SpVector a = IM.IntMap a
+  dat = svData
+
+instance Sparse SpVector a where
+  spy = spySV
+
+
 instance Hilbert SpVector where
   a `dot` b | dim a == dim b = dot (dat a) (dat b)
             | otherwise =
@@ -323,6 +314,7 @@
 singletonSV :: a -> SpVector a
 singletonSV x = SV 1 (IM.singleton 0 x)
 
+-- ** Create new sparse vector
 
 -- | create a sparse vector from an association list while discarding all zero entries
 mkSpVector :: (Num a, Eq a) => Int -> IM.IntMap a -> SpVector a
@@ -336,8 +328,9 @@
 mkSpVector1 :: Int -> IM.IntMap a -> SpVector a
 mkSpVector1 d ll = SV d $ IM.filterWithKey (\ k _ -> inBounds0 d k) ll
 
-mkSpVector1D :: Int -> [a] -> SpVector a
-mkSpVector1D d ll = mkSpVector1 d (IM.fromList $ denseIxArray (take d ll))
+-- | Create new sparse vector, assumin 0-based, contiguous indexing
+fromListDenseSV :: Int -> [a] -> SpVector a
+fromListDenseSV d ll = SV d (IM.fromList $ denseIxArray (take d ll))
 
 
 
@@ -352,7 +345,7 @@
 
 
 
--- insert
+-- |insert element `x` at index `i` in a preexisting SpVector
 insertSpVector :: Int -> a -> SpVector a -> SpVector a
 insertSpVector i x (SV d xim)
   | inBounds0 d i = SV d (IM.insert i x xim)
@@ -362,11 +355,11 @@
 fromListSV :: Int -> [(Int, a)] -> SpVector a
 fromListSV d iix = SV d (IM.fromList (filter (inBounds0 d . fst) iix ))
 
--- toList
+-- |toList
 toListSV :: SpVector a -> [(IM.Key, a)]
-toListSV sv = IM.toList (imSV sv)
+toListSV sv = IM.toList (dat sv)
 
--- to dense list (default = 0)
+-- |To dense list (default = 0)
 toDenseListSV :: Num b => SpVector b -> [b]
 toDenseListSV (SV d im) = fmap (\i -> IM.findWithDefault 0 i im) [0 .. d-1]
 
@@ -396,13 +389,14 @@
 
 -- | SV manipulation
 
+-- | Tail elements
 tailSV :: SpVector a -> SpVector a
 tailSV (SV n sv) = SV (n-1) ta where
   ta = IM.mapKeys (\i -> i - 1) $ IM.delete 0 sv
   
-
+-- | Head element
 headSV :: Num a => SpVector a -> a
-headSV sv = fromMaybe 0 (IM.lookup 0 (imSV sv))
+headSV sv = fromMaybe 0 (IM.lookup 0 (dat sv))
 
 
 
@@ -431,7 +425,7 @@
 
     
 
--- | outer vector product
+-- *** Outer vector product
 
 outerProdSV, (><) :: Num a => SpVector a -> SpVector a -> SpMatrix a
 outerProdSV v1 v2 = fromListSM (m, n) ixy where
@@ -449,17 +443,14 @@
 
 
 
--- | =======================================================
-
-
-
+-- * Sparse Matrix
 
 data SpMatrix a = SM {smDim :: (Rows, Cols),
                       smData :: IM.IntMap (IM.IntMap a)} deriving Eq
 
 
 
--- | instances for SpMatrix
+-- ** Instances for SpMatrix
 instance Show a => Show (SpMatrix a) where
   show sm@(SM _ x) = "SM " ++ sizeStr sm ++ " "++ show (IM.toList x)
 
@@ -473,9 +464,20 @@
 instance Additive SpMatrix where
   zero = SM (0,0) IM.empty
   (^+^) = liftU2 (+)
-  (^-^) = liftU2 (-)
 
 
+instance FiniteDim SpMatrix where
+  type FDSize SpMatrix = (Rows, Cols)
+  dim = smDim
+
+instance HasData SpMatrix a where
+  type HDData SpMatrix a = IM.IntMap (IM.IntMap a)
+  dat = smData
+
+instance Sparse SpMatrix a where
+  spy = spySM
+  
+
 -- | TODO : use semilattice properties instead
 maxTup, minTup :: Ord t => (t, t) -> (t, t) -> (t, t)
 maxTup (x1,y1) (x2,y2) = (max x1 x2, max y1 y2)
@@ -487,11 +489,11 @@
 
 
 
--- multiply matrix by a scalar
+-- *** multiply matrix by a scalar
 matScale :: Num a => a -> SpMatrix a -> SpMatrix a
 matScale a = fmap (*a)
 
--- Frobenius norm (sqrt of trace of M^T M)
+-- *** Frobenius norm (sqrt of trace of M^T M)
 normFrobenius :: SpMatrix Double -> Double
 normFrobenius m = sqrt $ foldlSM (+) 0 m' where
   m' | nrows m > ncols m = transposeSM m ## m
@@ -502,7 +504,7 @@
 
 
 
--- | ========= MATRIX METADATA
+-- ** MATRIX METADATA
 
 -- type synonyms
 type Rows = Int
@@ -511,16 +513,16 @@
 type IxRow = Int
 type IxCol = Int
 
--- -- predicates
--- are the supplied indices within matrix bounds?
+-- *** predicates
+-- |Are the supplied indices within matrix bounds?
 validIxSM :: SpMatrix a -> (Int, Int) -> Bool
 validIxSM mm = inBounds02 (dim mm)
 
--- is the matrix square?
+-- |Is the matrix square?
 isSquareSM :: SpMatrix a -> Bool
 isSquareSM m = nrows m == ncols m
 
--- is the matrix diagonal?
+-- |Is the matrix diagonal?
 isDiagonalSM :: SpMatrix a -> Bool
 isDiagonalSM m = IM.size d == nrows m where
   d = IM.filterWithKey ff (immSM m)
@@ -565,7 +567,7 @@
 spySM s = fromIntegral (nzSM s) / fromIntegral (nelSM s)
 
 
--- # NZ in row i
+-- *** Non-zero elements in a row
 
 nzRowU :: SpMatrix a -> IM.Key -> Int
 nzRowU s i = maybe 0 IM.size (IM.lookup i $ immSM s)
@@ -577,7 +579,7 @@
 
 
 
--- | bandwidth bounds (min, max)
+-- *** bandwidth bounds (min, max)
 
 bwMinSM :: SpMatrix a -> Int
 bwMinSM = fst . bwBoundsSM
@@ -603,7 +605,7 @@
 
 
 
--- | ========= SPARSE MATRIX BUILDERS
+-- ** SPARSE MATRIX BUILDERS
 
 zeroSM :: Int -> Int -> SpMatrix a
 zeroSM m n = SM (m,n) IM.empty 
@@ -617,15 +619,17 @@
       d = dim s
 
 
--- | from list (row, col, value)
+-- | Add to existing SpMatrix using data from list (row, col, value)
 fromListSM' :: Foldable t => t (IxRow, IxCol, a) -> SpMatrix a -> SpMatrix a
 fromListSM' iix sm = foldl ins sm iix where
   ins t (i,j,x) = insertSpMatrix i j x t
 
+-- | Create new SpMatrix using data from list (row, col, value)
 fromListSM :: Foldable t => (Int, Int) -> t (IxRow, IxCol, a) -> SpMatrix a
 fromListSM (m,n) iix = fromListSM' iix (zeroSM m n)
 
 
+-- | Create new SpMatrix assuming contiguous, 0-based indexing of elements
 fromListDenseSM :: Int -> [a] -> SpMatrix a
 fromListDenseSM m ll = fromListSM (m, n) $ denseIxArray2 m ll where
   n = length ll `div` m
@@ -634,7 +638,7 @@
 
 -- | to List
 
--- toDenseListSM : populate missing entries with 0
+-- |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]]
@@ -643,11 +647,11 @@
 
 
 
--- -- create diagonal and identity matrix
+-- ** Diagonal matrix
 mkDiagonal :: Int -> [a] -> SpMatrix a
 mkDiagonal n = mkSubDiagonal n 0
 
-
+-- *** Identity matrix
 eye :: Num a => Int -> SpMatrix a
 eye n = mkDiagonal n (ones n)
 
@@ -658,7 +662,7 @@
   
 
 
--- super- and sub- diagonal
+-- *** Create Super- or sub- diagonal matrix
 
 mkSubDiagonal :: Int -> Int -> [a] -> SpMatrix a
 mkSubDiagonal n o xx | abs o < n = if o >= 0
@@ -693,7 +697,7 @@
 
 
 
--- | ========= SUB-MATRICES
+-- ** SUB-MATRICES
 
 
 extractSubmatrixSM :: SpMatrix a -> (Int, Int) -> (Int, Int) -> SpMatrix a
@@ -714,16 +718,7 @@
       inBounds0 c j2 &&      
       i2 >= i1
 
--- extract row / column
-extractRowSM :: SpMatrix a -> Int -> SpMatrix a
-extractRowSM sm i = extractSubmatrixSM sm (i, i) (0, ncols sm - 1)
-
-extractColSM :: SpMatrix a -> Int -> SpMatrix a
-extractColSM sm j = extractSubmatrixSM sm (0, nrows sm - 1) (j, j)
-
-
-
--- demote (n x 1) or (1 x n) SpMatrix to SpVector
+-- |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
   d | m==1 && n==1 = 1
@@ -732,12 +727,22 @@
     | otherwise = error $ "toSV : incompatible dimensions " ++ show (m,n)
 
 
--- extract row or column and place into SpVector
+-- *** Extract jth column
+extractColSM :: SpMatrix a -> Int -> SpMatrix a
+extractColSM sm j = extractSubmatrixSM sm (0, nrows sm - 1) (j, j)
+
+-- |", and place into SpVector
 extractCol :: SpMatrix a -> Int -> SpVector a
-extractCol m i = toSV $ extractColSM m i
+extractCol m j = toSV $ extractColSM m j
 
+
+-- *** Extract ith row
+extractRowSM :: SpMatrix a -> Int -> SpMatrix a
+extractRowSM sm i = extractSubmatrixSM sm (i, i) (0, ncols sm - 1)
+
+-- |", and place into SpVector
 extractRow :: SpMatrix a -> Int -> SpVector a
-extractRow m j = toSV $ extractRowSM m j
+extractRow m i = toSV $ extractRowSM m i
 
 
 
@@ -745,8 +750,9 @@
 
 
 
--- | ========= MATRIX STACKING
+-- ** MATRIX STACKING
 
+-- | Vertical stacking
 vertStackSM, (-=-) :: SpMatrix a -> SpMatrix a -> SpMatrix a
 vertStackSM mm1 mm2 = SM (m, n) $ IM.union u1 u2 where
   nro1 = nrows mm1
@@ -757,7 +763,7 @@
 
 (-=-) = vertStackSM
 
-
+-- | Horizontal stacking
 horizStackSM, (-||-) :: SpMatrix a -> SpMatrix a -> SpMatrix a
 horizStackSM mm1 mm2 = t (t mm1 -=- t mm2) where
   t = transposeSM
@@ -772,12 +778,12 @@
 
 
 
--- | ========= LOOKUP
+-- ** MATRIX ELEMENT LOOKUP
 
 lookupSM :: SpMatrix a -> IM.Key -> IM.Key -> Maybe a
 lookupSM (SM _ im) i j = IM.lookup i im >>= IM.lookup j
 
--- | Looks up an element in the matrix (if not found, zero is returned)
+-- | 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 -> (IM.Key, IM.Key) -> a
 lookupWD_SM sm (i,j) =
@@ -786,6 +792,7 @@
 lookupWD_IM :: Num a => IM.IntMap (IM.IntMap a) -> (IM.Key, IM.Key) -> a
 lookupWD_IM im (i,j) = fromMaybe 0 (IM.lookup i im >>= IM.lookup j)
 
+-- | Zero-default lookup, infix form
 (@@) = lookupWD_SM
 
 
@@ -799,7 +806,7 @@
 
 
 
--- | ========= MISC SpMatrix OPERATIONS
+-- *** MISC SpMatrix OPERATIONS
 
 foldlSM :: (a -> b -> b) -> b -> SpMatrix a -> b
 foldlSM f n (SM _ m)= foldlIM2 f n m
@@ -830,9 +837,9 @@
 -- extractDiagonalSM :: (Num a, Eq a) => SpMatrix a -> SpVector a
 -- extractDiagonalSM (SM (m,n) im) = mkSpVectorD m $ extractDiagonalIM2 im
 
--- extract with default 0
+-- | extract with default 0
 extractDiagonalDSM :: Num a => SpMatrix a -> SpVector a
-extractDiagonalDSM mm = mkSpVector1D n $ foldr ins [] ll  where
+extractDiagonalDSM mm = fromListDenseSV n $ foldr ins [] ll  where
   ll = [0 .. n - 1]
   n = nrows mm
   ins i acc = mm@@(i,i) : acc
@@ -852,7 +859,7 @@
 
 
 
--- | sparsify : remove 0s (!!!)
+-- ** sparsify : remove 0s (!!!)
 
 sparsifyIM2 :: IM.IntMap (IM.IntMap Double) -> IM.IntMap (IM.IntMap Double)
 sparsifyIM2 = ifilterIM2 (\_ _ x -> abs x >= eps)
@@ -862,7 +869,7 @@
 
 
 
--- | ROUNDING operations (!!!)
+-- ** ROUNDING operations (!!!)
                               
 roundZeroOneSM :: SpMatrix Double -> SpMatrix Double
 roundZeroOneSM (SM d im) = sparsifySM $ SM d $ mapIM2 roundZeroOne im
@@ -877,10 +884,10 @@
 
 
 
--- | ========= ALGEBRAIC PRIMITIVE OPERATIONS
+-- * ALGEBRAIC PRIMITIVE OPERATIONS
 
 
--- | transpose
+-- ** Matrix transpose
 
 
 transposeSM, (#^) :: SpMatrix a -> SpMatrix a
@@ -904,7 +911,7 @@
 
 
 
--- | matrix action on a vector
+-- ** matrix action on a vector
 
 {- 
 FIXME : matVec is more general than SpVector's :
@@ -915,7 +922,7 @@
 
 
 
--- matrix on vector
+-- |matrix on vector
 matVec, (#>) :: Num a => SpMatrix a -> SpVector a -> SpVector a
 matVec (SM (nr, nc) mdata) (SV n sv)
   | nc == n = SV nr $ fmap (`dot` sv) mdata
@@ -923,7 +930,7 @@
 
 (#>) = matVec
 
--- vector on matrix (FIXME : transposes matrix: more costly than `matVec`)
+-- |vector on matrix (FIXME : transposes matrix: more costly than `matVec`, I think)
 vecMat, (<#) :: Num a => SpVector a -> SpMatrix a -> SpVector a  
 vecMat (SV n sv) (SM (nr, nc) mdata)
   | n == nr = SV nc $ fmap (`dot` sv) (transposeIM2 mdata)
@@ -940,17 +947,12 @@
 
 
 
--- matVec' mm vv =
---   withDim2 mm vv (\(nro, nco) nv _ _ -> nco == nv) matVecU "matVec : mismatching dimensions"
---    (\ m v -> unwords [show (dim m), show (dim v)])
 
 
--- asdfm ll = unwords (map (show . dim) ll)
 
 
 
-
--- | matrix-matrix product
+-- ** Matrix-matrix product
 
 -- unsafe matMat
 matMatU :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a
@@ -959,12 +961,6 @@
     im = fmap (\vm1 -> (`dot` vm1) <$> transposeIM2 (immSM m2)) (immSM m1)
 
 
--- matMat, (##) :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a
--- matMat (SM (nr1,nc1) m1) (SM (nr2,nc2) m2)
---   | nc1 == nr2 = SM (nr1, nc2) $
---       fmap (\vm1 -> fmap (`dot` vm1) (transposeIM2 m2)) m1
---   | otherwise = error "matMat : incompatible matrix sizes"
-
 matMat, (##) :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a
 matMat m1 m2
   | c1 == r2 = matMatU m1 m2
@@ -986,7 +982,7 @@
 
 
 
--- | sparsified matrix-matrix product (prunes all elements `x` for which `abs x <= eps`)
+-- ** Matrix-matrix product, sparsified (prunes all elements `x` for which `abs x <= eps`)
 matMatSparsified, (#~#)  :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double
 matMatSparsified m1 m2 = sparsifySM $ matMat m1 m2
 
@@ -997,9 +993,9 @@
 
 
 
--- | ========= predicates
+-- * Predicates
 
--- is the matrix orthogonal? i.e. Q^t ## Q == I
+-- |is the matrix orthogonal? i.e. Q^t ## Q == I
 isOrthogonalSM :: SpMatrix Double -> Bool
 isOrthogonalSM sm@(SM (_,n) _) = rsm == eye n where
   rsm = roundZeroOneSM $ transposeSM sm ## sm
@@ -1012,9 +1008,9 @@
 
 
 
--- | ========= condition number
+-- ** Condition number
 
--- uses the R matrix from the QR factorization
+-- |uses the R matrix from the QR factorization
 conditionNumberSM :: SpMatrix Double -> Double
 conditionNumberSM m | isInfinite kappa = error "Infinite condition number : rank-deficient system"
                     | otherwise = kappa where
@@ -1030,15 +1026,16 @@
 
 
 
--- | ========= Householder transformation
+-- ** Householder transformation
 
 hhMat :: Num a => a -> SpVector a -> SpMatrix a
 hhMat beta x = eye n ^-^ scale beta (x >< x) where
   n = dim x
 
 
--- 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 
+{-| 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
 
@@ -1052,7 +1049,7 @@
 
 
 
--- | ========= Givens rotation matrix
+-- ** Givens rotation matrix
 
 
 hypot :: Floating a => a -> a -> a
@@ -1064,6 +1061,7 @@
   | x == 0 = 0
   | otherwise = -1 
 
+-- | Givens coefficients (using stable algorithm shown in  Anderson, Edward (4 December 2000). "Discontinuous Plane Rotations and the Symmetric Eigenvalue Problem". LAPACK Working Note)
 givensCoef :: (Ord a, Floating a) => a -> a -> (a, a, a)
 givensCoef a b  -- returns (c, s, r) where r = norm (a, b)
   | b==0 = (sign a, 0, abs a)
@@ -1076,7 +1074,7 @@
                 in (t/u, - 1/u, b*u)
 
 
-{-
+{- |
 Givens method, row version: choose other row index i' s.t. i' is :
 * below the diagonal
 * corresponding element is nonzero
@@ -1096,12 +1094,12 @@
     a = mm @@ (i', j)
     b = mm @@ (i, j)   -- element to zero out
 
--- is the `k`th the first nonzero column in the row?
+-- |Is the `k`th the first nonzero column in the row?
 firstNonZeroColumn :: IM.IntMap a -> IM.Key -> Bool
 firstNonZeroColumn mm k = isJust (IM.lookup k mm) &&
                           isNothing (IM.lookupLT k mm)
 
--- returns a set of rows {k} that satisfy QR.C1
+-- |Returns a set of rows {k} that satisfy QR.C1
 candidateRows :: IM.IntMap (IM.IntMap a) -> IM.Key -> IM.Key -> Maybe [IM.Key]
 candidateRows mm i j | IM.null u = Nothing
                      | otherwise = Just (IM.keys u) where
@@ -1112,20 +1110,17 @@
 
 
 
--- | ========= QR algorithm
-
-{-
-applies Givens rotation iteratively to zero out sub-diagonal elements
--}
+-- ** QR decomposition
 
 
+-- | Applies Givens rotation iteratively to zero out sub-diagonal elements
 qr :: SpMatrix Double -> (SpMatrix Double, SpMatrix Double)
 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 ]
+-- | Givens matrices in order [G1, G2, .. , G_N ]
 gmats :: SpMatrix Double -> [SpMatrix Double]
 gmats mm = gm mm (subdiagIndicesSM mm) where
  gm m ((i,j):is) = let g = givens m i j
@@ -1148,9 +1143,12 @@
 
 
 
--- | ========= Eigenvalues, using QR
 
+-- ** Eigenvalue algorithms
 
+-- *** All eigenvalues using QR algorithm
+
+
 eigsQR :: Int -> SpMatrix Double -> SpVector Double
 eigsQR nitermax m = extractDiagonalDSM $ execState (convergtest eigsStep) m where
   eigsStep m = r #~# q where (q, r) = qr m
@@ -1164,11 +1162,12 @@
 
 
 
--- | ========= Eigenvalues, using Rayleigh iteration
+-- *** One eigenvalue and corresponding eigenvector, using Rayleigh iteration 
 
+-- | Cubic-order convergence, but it requires a mildly educated guess on the initial eigenpair
 rayleighStep ::
   SpMatrix Double ->
-  (SpVector Double, Double) ->
+  (SpVector Double, Double) -> 
   (SpVector Double, Double)    -- updated estimate of (eigenvector, eigenvalue)
 rayleighStep aa (b, mu) = (b', mu') where
   ii = eye (nrows aa)
@@ -1188,9 +1187,8 @@
 
 
 
--- | ========= Householder vector (G & VL Alg. 5.1.1, function `house`)
+-- ** Householder vector (G & VL Alg. 5.1.1, function `house`)
 
--- hhV :: (Ord a, Floating a) => SpVector a -> (SpVector a, a)
 hhV :: SpVector Double -> (SpVector Double, Double)
 hhV x = (v, beta) where
   n = dim x
@@ -1212,7 +1210,7 @@
 
 
 
--- | ========= SVD
+-- * SVD
 
 {- Golub & Van Loan, sec 8.6.2 (p 452 segg.)
 
@@ -1240,10 +1238,10 @@
 
 
 
--- | =======================================================
 
--- | LINEAR SOLVERS : solve A x = b
 
+-- * LINEAR SOLVERS : solve A x = b
+
 -- | numerical tolerance for e.g. solution convergence
 eps :: Double
 eps = 1e-8
@@ -1257,7 +1255,7 @@
 
 
 
--- | CGS
+-- ** CGS
 
 -- | one step of CGS
 cgsStep :: SpMatrix Double -> SpVector Double -> CGS -> CGS
@@ -1277,15 +1275,14 @@
                  _p :: SpVector Double,
                  _u :: SpVector Double } deriving Eq
 
+-- | iterate solver until convergence or until max # of iterations is reached
 cgs ::
   SpMatrix Double ->
   SpVector Double ->
   SpVector Double ->
   SpVector Double ->
-  -- Int ->
   CGS
 cgs aa b x0 rhat =
-  -- execState (replicateM n (modify (cgsStep aa rhat))) cgsInit where
   execState (untilConverged _x (cgsStep aa rhat)) cgsInit where
   r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution
   p0 = r0
@@ -1303,7 +1300,7 @@
 
   
 
--- | BiCSSTAB
+-- ** BiCSSTAB
 
 -- _aa :: SpMatrix Double,    -- matrix
 -- _b :: SpVector Double,     -- rhs
@@ -1327,21 +1324,18 @@
                            _rBicgstab :: SpVector Double,
                            _pBicgstab :: SpVector Double} deriving Eq
 
-
+-- | iterate solver until convergence or until max # of iterations is reached
 bicgstab
   :: SpMatrix Double
      -> SpVector Double
      -> SpVector Double
      -> SpVector Double
-     -- -> Int
      -> BICGSTAB
 bicgstab aa b x0 r0hat =
-  -- execState (replicateM n (modify (bicgstabStep aa r0hat))) bicgsInit where
   execState (untilConverged _xBicgstab (bicgstabStep aa r0hat)) bicgsInit where
    r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution
    p0 = r0
    bicgsInit = BICGSTAB x0 r0 p0
-   -- q (BICGSTAB xi _ _) = nor
 
 instance Show BICGSTAB where
   show (BICGSTAB x r p) = "x = " ++ show x ++ "\n" ++
@@ -1358,11 +1352,11 @@
 
 
 
--- | ========= LINEAR SOLVERS INTERFACE
+-- * LINEAR SOLVERS INTERFACE
 
 data LinSolveMethod = CGS_ | BICGSTAB_ deriving (Eq, Show) 
 
--- random starting vector
+-- | linear solve with _random_ starting vector
 linSolveM ::
   PrimMonad m =>
     LinSolveMethod -> SpMatrix Double -> SpVector Double -> m (SpVector Double)
@@ -1374,7 +1368,7 @@
     case method of CGS_ -> return $ _xBicgstab (bicgstab aa b x0 x0)
                    BICGSTAB_ -> return $ _x (cgs aa b x0 x0)
 
--- deterministic starting vector (every component at 0.1) 
+-- | linear solve with _deterministic_ starting vector (every component at 0.1) 
 linSolve ::
   LinSolveMethod -> SpMatrix Double -> SpVector Double -> SpVector Double
 linSolve method aa b
@@ -1422,9 +1416,8 @@
 
 
 
--- | ========= PRETTY PRINTING
+-- * PRETTY PRINTING
 
--- | Show details and contents of sparse matrix
 
 sizeStr :: SpMatrix a -> String
 sizeStr sm =
@@ -1484,6 +1477,7 @@
       rr_' | dim sv > nco = unwords [take (nco - 2) rr_ , " ... " , [last rr_]]
            | otherwise = rr_
 
+-- ** Pretty printer typeclass
 class PrintDense a where
   prd :: a -> IO ()
 
@@ -1550,9 +1544,13 @@
                            else go (i + 1) (take 2 $ y:ll) y
                 where y = f xx
 
--- modify state and append, until max # of iterations is reached
+-- | keep state `x` in a moving window of length 2 to assess convergence, stop when either a condition on that list is satisfied or when max # of iterations is reached
 modifyInspectN ::
-  MonadState s m => Int -> ([s] -> Bool) -> (s -> s) -> m s
+  MonadState s m =>
+    Int ->           -- iteration budget
+    ([s] -> Bool) -> -- convergence criterion
+    (s -> s) ->      -- state stepping function
+    m s
 modifyInspectN nitermax q f 
   | nitermax > 0 = go 0 []
   | otherwise = error "modifyInspectN : n must be > 0" where
@@ -1569,13 +1567,14 @@
                        go (i + 1) (take 2 $ y : ll)
 
 
+-- helper functions for estimating convergence
 meanl :: (Foldable t, Fractional a) => t a -> a
 meanl xx = 1/fromIntegral (length xx) * sum xx
 
 norm2l :: (Foldable t, Functor t, Floating a) => t a -> a
 norm2l xx = sqrt $ sum (fmap (**2) xx)
 
-
+diffSqL :: Floating a => [a] -> a
 diffSqL xx = (x1 - x2)**2 where [x1, x2] = [head xx, xx!!1]
 
 
@@ -1594,9 +1593,12 @@
 normDiffConverged :: (Foldable t, Functor t) =>
      (a -> SpVector Double) -> t a -> Bool
 normDiffConverged fp xx = normSq (foldrMap fp (^-^) (zeroSV 0) xx) <= eps
-              
 
 
+  
+
+
+
 -- run `niter` iterations and append the state `x` to a list `xs`, stop when either the `xs` satisfies a predicate `q` or when the counter reaches 0
 
 runAppendN :: ([t] -> Bool) -> (t -> t) -> Int -> t -> [t]
@@ -1616,8 +1618,6 @@
     if n <= 0 then xs
               else go f (n-1) x (x : xs)
 
--- runN :: Int -> (a -> a) -> a -> a
--- runN n stepf x0 = runAppendN' stepf n x0
   
 
 
@@ -1739,104 +1739,15 @@
 
 --
 
-tm0, tm1, tm2, tm3, tm4 :: SpMatrix Double
-tm0 = fromListSM (2,2) [(0,0,pi), (1,0,sqrt 2), (0,1, exp 1), (1,1,sqrt 5)]
 
-tv0, tv1 :: SpVector Double
-tv0 = mkSpVectorD 2 [5, 6]
 
 
-tv1 = SV 2 $ IM.singleton 0 1
 
--- wikipedia test matrix for Givens rotation
 
-tm1 = sparsifySM $ fromListDenseSM 3 [6,5,0,5,1,4,0,4,3]
 
-tm1g1 = givens tm1 1 0
-tm1a2 = tm1g1 ## tm1
 
-tm1g2 = givens tm1a2 2 1
-tm1a3 = tm1g2 ## tm1a2
 
-tm1q = transposeSM (tm1g2 ## tm1g1)
 
-
--- wp test matrix for QR decomposition via Givens rotation
-
-tm2 = fromListDenseSM 3 [12, 6, -4, -51, 167, 24, 4, -68, -41]
-
-
-
-
-tm3 = transposeSM $ fromListDenseSM 3 [1 .. 9]
-
-tm3g1 = fromListDenseSM 3 [1, 0,0, 0,c,-s, 0, s, c]
-  where c= 0.4961
-        s = 0.8682
-
-
---
-
-tm4 = sparsifySM $ fromListDenseSM 4 [1,0,0,0,2,5,0,10,3,6,8,11,4,7,9,12]
-
-
--- playground
-
--- | terminate after n iterations or when q becomes true, whichever comes first
-untilC :: (a -> Bool) -> Int ->  (a -> a) -> a -> a
-untilC p n f = go n
-  where
-    go m x | p x || m <= 0 = x
-           | otherwise     = True `seq` go (m-1) (f x)
-
-
-
-
-
--- testing State
-
-
--- data T0 = T0 {unT :: Int} deriving Eq
--- instance Show T0 where
---   show (T0 x) = show x
-
--- -- modifyT :: MonadState T0 m => (Int -> Int) -> m String
--- modifyT f = state (\(T0 i) -> (i, T0 (f i)))
-  
-
--- t00 = T0 0
-
--- testT n = execState $ replicateM n (modifyT (+1)) 
-
-
--- testT2 = execState $ when 
-  
-
--- replicateSwitch p m f = loop m where
---       loop n | n <= 0 || p = pure (#)
---              | otherwise = f *> loop (n-1)
-
-
-
--- testing Writer
-               
--- asdfw n = runWriter $ do
---   tell $ "potato " ++ show n
---   tell "jam"
---   return (n+1)
-
-
--- --
-
-
--- testing State and Writer
-
-
-
--- runMyApp runA k maxDepth =
---     let config = maxDepth
---         state =  0
---     in runStateT (runWriterT (runA k) config) state
 
 
   
diff --git a/test/LibSpec.hs b/test/LibSpec.hs
new file mode 100644
--- /dev/null
+++ b/test/LibSpec.hs
@@ -0,0 +1,320 @@
+{-# language ScopedTypeVariables #-}
+module LibSpec where
+
+import qualified Data.IntMap as IM
+
+import Control.Monad (replicateM)
+import Control.Monad.State.Strict (execState)
+
+import qualified System.Random.MWC as MWC
+import qualified System.Random.MWC.Distributions as MWC
+       
+import Test.Hspec
+-- import Test.Hspec.QuickCheck
+
+import Lib
+import Math.Linear.Sparse
+
+
+
+main :: IO ()
+main = hspec spec
+
+-- niter = 5
+
+spec :: Spec
+spec = do
+  describe "Math.Linear.Sparse : library" $ do
+    -- prop "subtraction is cancellative" $ \(x :: SpVector Double) ->
+    --   x ^-^ x `shouldBe` zero
+    it "dot : inner product" $
+      tv0 `dot` tv0 `shouldBe` 61
+    it "transposeSM : sparse matrix transpose" $
+      transposeSM m1 `shouldBe` m1t
+    it "matVec : matrix-vector product" $
+      normSq ((aa0 #> x0true) ^-^ b0 ) <= eps `shouldBe` True
+    it "vecMat : vector-matrix product" $
+      normSq ((x0true <# aa0) ^-^ aa0tx0 ) <= eps `shouldBe` True  
+    it "matMat : matrix-matrix product" $
+      (m1 `matMat` m2) `shouldBe` m1m2
+    it "eye : identity matrix" $
+      infoSM (eye 10) `shouldBe` SMInfo 10 0.1
+    it "countSubdiagonalNZ : # of nonzero elements below the diagonal" $
+      countSubdiagonalNZSM m3 `shouldBe` 1
+    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    
+    it "BiCGSTAB (2 x 2 dense)" $ 
+      -- normSq (_xBicgstab (bicgstab aa0 b0 x0 x0) ^-^ x0true) <= eps `shouldBe` True
+      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    
+    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
+
+
+
+-- -- 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
+
+[1 2] [2] = [8]
+[3 4] [3]   [18]
+
+
+[1 3] [2] = [11]
+[2 4] [3]   [16]
+
+
+-}
+
+aa0 :: SpMatrix Double
+aa0 = SM (2,2) im where
+  im = IM.fromList [(0, aa0r0), (1, aa0r1)]
+
+aa0r0, aa0r1 :: IM.IntMap Double
+aa0r0 = IM.fromList [(0,1),(1,2)]
+aa0r1 = IM.fromList [(0,3),(1,4)]
+
+
+-- b0, x0 : r.h.s and initial solution resp.
+b0, x0, x0true :: SpVector Double
+b0 = mkSpVectorD 2 [8,18]
+x0 = mkSpVectorD 2 [0.3,1.4]
+
+
+-- x0true : true solution
+x0true = mkSpVectorD 2 [2,3]
+
+
+
+aa0tx0 = mkSpVectorD 2 [11,16]
+
+
+
+
+
+
+
+{- 4x4 system -}
+
+aa1 :: SpMatrix Double
+aa1 = sparsifySM $ fromListDenseSM 4 [1,0,0,0,2,5,0,10,3,6,8,11,4,7,9,12]
+
+x1, b1 :: SpVector Double
+x1 = mkSpVectorD 4 [1,2,3,4]
+
+b1 = mkSpVectorD 4 [30,56,60,101]
+
+
+
+{- 3x3 system -}
+aa2 :: SpMatrix Double
+aa2 = sparsifySM $ fromListDenseSM 3 [2, -1, 0, -1, 2, -1, 0, -1, 2]
+x2, b2 :: SpVector Double
+x2 = mkSpVectorD 3 [3,2,3]
+
+b2 = mkSpVectorD 3 [4,-2,4]
+
+
+
+-- --
+
+{-
+example 1 : random linear system
+
+-}
+
+
+
+-- dense
+solveRandom n = do
+  aa0 <- randMat n
+  let aa = aa0 ^+^ eye n
+  xtrue <- randVec n
+  -- x0 <- randVec n
+  let b = aa #> xtrue
+      dx = aa <\> b ^-^ xtrue
+  return $ normSq dx
+  -- let xhatB = _xBicgstab (bicgstab aa b x0 x0)
+  --     xhatC = _x (cgs aa b x0 x0)
+  -- return (aa, x, x0, b, xhatB, xhatC)
+
+-- sparse
+solveSpRandom :: Int -> Int -> IO Double
+solveSpRandom n nsp = do
+  aa0 <- randSpMat n nsp
+  let aa = aa0 ^+^ eye n
+  xtrue <- randSpVec n nsp
+  let b = (aa ^+^ eye n) #> xtrue
+      dx = aa <\> b ^-^ xtrue
+  return $ normSq dx
+
+
+
+-- `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))
+
+--
+
+{-
+matMat
+
+[1, 2] [5, 6] = [19, 22]
+[3, 4] [7, 8]   [43, 50]
+-}
+
+m1 = fromListDenseSM 2 [1,3,2,4]
+m2 = fromListDenseSM 2 [5, 7, 6, 8]     
+m1m2 = fromListDenseSM 2 [19, 43, 22, 50]
+
+-- transposeSM
+
+m1t = fromListDenseSM 2 [1,2,3,4]
+
+
+--
+
+{-
+countSubdiagonalNZ
+-}
+
+m3 = fromListSM (3,3) [(0,2,3),(2,0,4),(1,1,3)] 
+
+
+
+
+{- mkSubDiagonal -}
+
+testLaplacian1 :: Int -> SpMatrix Double
+testLaplacian1 n = m where
+  m :: SpMatrix Double
+  m = mksd (-1) l1 ^+^
+       mksd 0 l2 ^+^
+       mksd 1 l3
+    where
+    mksd = mkSubDiagonal n
+    l1 = replicate n (-1)
+    l2 = replicate n 2
+    l3 = l1
+  -- x :: SpVector Double
+  -- x = mkSpVectorD n (replicate n 2)
+  -- b = m #> x
+
+-- t3 n = normSq $ (aa <\> b) ^-^ xhat where
+--   aa = testLaplacian1 n :: SpMatrix Double
+--   xhat = mkSpVectorD n (concat $ replicate 20 [1,2,3,4,5]) :: SpVector Double
+--   b = aa #> xhat
+
+{- QR-}
+
+checkQr :: SpMatrix Double -> Bool
+checkQr a = c1 && c2 where
+  (q, r) = qr a
+  c1 = normFrobenius ((q #~# r) ^-^ a) <= eps
+  c2 = isOrthogonalSM q
+
+
+aa22 = fromListDenseSM 2 [2,1,1,2] :: SpMatrix Double
+
+
+
+
+{- eigenvalues -}
+
+
+aa3 = fromListDenseSM 3 [1,1,3,2,2,2,3,1,1] :: SpMatrix Double
+
+b3 = mkSpVectorD 3 [1,1,1] :: SpVector Double
+
+
+
+
+
+
+-- test data
+
+tm0, tm1, tm2, tm3, tm4 :: SpMatrix Double
+tm0 = fromListSM (2,2) [(0,0,pi), (1,0,sqrt 2), (0,1, exp 1), (1,1,sqrt 5)]
+
+tv0, tv1 :: SpVector Double
+tv0 = mkSpVectorD 2 [5, 6]
+
+
+tv1 = SV 2 $ IM.singleton 0 1
+
+-- wikipedia test matrix for Givens rotation
+
+tm1 = sparsifySM $ fromListDenseSM 3 [6,5,0,5,1,4,0,4,3]
+
+tm1g1 = givens tm1 1 0
+tm1a2 = tm1g1 ## tm1
+
+tm1g2 = givens tm1a2 2 1
+tm1a3 = tm1g2 ## tm1a2
+
+tm1q = transposeSM (tm1g2 ## tm1g1)
+
+
+-- wp test matrix for QR decomposition via Givens rotation
+
+tm2 = fromListDenseSM 3 [12, 6, -4, -51, 167, 24, 4, -68, -41]
+
+
+
+
+tm3 = transposeSM $ fromListDenseSM 3 [1 .. 9]
+
+tm3g1 = fromListDenseSM 3 [1, 0,0, 0,c,-s, 0, s, c]
+  where c= 0.4961
+        s = 0.8682
+
+
+--
+
+tm4 = sparsifySM $ fromListDenseSM 4 [1,0,0,0,2,5,0,10,3,6,8,11,4,7,9,12]