diff --git a/README.md b/README.md
--- a/README.md
+++ b/README.md
@@ -77,7 +77,7 @@
     [0.0,0.0,5.0]
 
 Notice that the result is _dense_, i.e. certain entries are numerically zero but have been inserted into the result along with all the others (thus taking up memory!).
-To preserve sparsity, we can use a sparsifying matrix-matrix product `#~#`, which filters out all the elements x for which `|x| <= eps`, where `eps` (defined) in `Numeric.Eps`, is fixed at 10^-8.
+To preserve sparsity, we can use a sparsifying matrix-matrix product `#~#`, which filters out all the elements x for which `|x| <= eps`, where `eps` (defined in `Numeric.Eps`) is fixed at 10^-8.
 
     > prd $ l #~# u
     ( 3 rows, 3 columns ) , 5 NZ ( sparsity 0.5555555555555556 )
diff --git a/sparse-linear-algebra.cabal b/sparse-linear-algebra.cabal
--- a/sparse-linear-algebra.cabal
+++ b/sparse-linear-algebra.cabal
@@ -1,5 +1,5 @@
 name:                sparse-linear-algebra
-version:             0.2.0.8
+version:             0.2.0.9
 synopsis:            Numerical computation in native Haskell 
 description:         Currently the library provides iterative linear solvers, matrix decompositions, eigenvalue computations and related utilities. Please see README.md for details
 homepage:            https://github.com/ocramz/sparse-linear-algebra
@@ -20,7 +20,7 @@
   hs-source-dirs:      src
   exposed-modules:     Numeric.LinearAlgebra.Sparse
                        Numeric.LinearAlgebra.Class
-                       Numeric.LinearAlgebra.Sparse.IntMap
+                       Data.Sparse.IntMap2.IntMap2
                        Data.Sparse.SpVector
                        Data.Sparse.SpMatrix
                        Data.Sparse.Common
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
@@ -12,13 +12,12 @@
 
 import Data.Sparse.Utils as X
 import Data.Sparse.Types as X
+import Data.Sparse.IntMap2.IntMap2 as X
 import Data.Sparse.SpMatrix as X
 import Data.Sparse.SpVector as X
 
 import Numeric.Eps as X
 import Numeric.LinearAlgebra.Class as X
-
-import Numeric.LinearAlgebra.Sparse.IntMap as X
 
 import qualified Data.IntMap as IM
 
diff --git a/src/Data/Sparse/IntMap2/IntMap2.hs b/src/Data/Sparse/IntMap2/IntMap2.hs
new file mode 100644
--- /dev/null
+++ b/src/Data/Sparse/IntMap2/IntMap2.hs
@@ -0,0 +1,187 @@
+module Data.Sparse.IntMap2.IntMap2 where
+
+import Numeric.LinearAlgebra.Class
+
+import qualified Data.IntMap.Strict as IM
+
+
+
+
+
+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
+
+
+-- 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
+
+-- |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 #-}  
+
+-- |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)
+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
@@ -7,7 +7,7 @@
 import Numeric.Eps
 import Numeric.LinearAlgebra.Class
 
-import Numeric.LinearAlgebra.Sparse.IntMap
+import Data.Sparse.IntMap2.IntMap2
 
 import qualified Data.IntMap as IM
 
@@ -209,7 +209,20 @@
 
 -- ** Sub-matrices
 
+-- | Indexed filtering function
+filterSM :: (IM.Key -> IM.Key -> a -> Bool) -> SpMatrix a -> SpMatrix a
+filterSM f sm = SM (dim sm) $ ifilterIM2 f (dat sm)
 
+-- | Diagonal, subdiagonal, superdiagonal partitions of a SpMatrix (useful for writing preconditioners)
+extractDiag, extractSuperDiag, extractSubDiag :: SpMatrix a -> SpMatrix a
+extractSubDiag = filterSM (\i j _ -> i > j)
+
+extractSuperDiag = filterSM (\i j _ -> i < j)
+
+extractDiag = filterSM (\i j _ -> i == j)
+
+
+
 -- | Extract a submatrix given the specified index bounds, rebalancing keys with the two supplied functions
 extractSubmatrixSM ::
   (IM.Key -> IM.Key) ->   -- row index function
@@ -582,6 +595,12 @@
 -- ** Multiply matrix by a scalar
 matScale :: Num a => a -> SpMatrix a -> SpMatrix a
 matScale a = fmap (*a)
+
+
+
+
+
+
 
 -- ** Frobenius norm
 normFrobenius :: SpMatrix Double -> Double
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
@@ -5,7 +5,7 @@
 import Data.Sparse.Types
 
 import Numeric.LinearAlgebra.Class
-import Numeric.LinearAlgebra.Sparse.IntMap
+import Data.Sparse.IntMap2.IntMap2
 
 import Data.Maybe
 
diff --git a/src/Numeric/LinearAlgebra/Sparse.hs b/src/Numeric/LinearAlgebra/Sparse.hs
--- a/src/Numeric/LinearAlgebra/Sparse.hs
+++ b/src/Numeric/LinearAlgebra/Sparse.hs
@@ -4,6 +4,8 @@
        (
          -- * Matrix factorizations
          qr, lu,
+         -- * Incomplete LU
+         ilu0,
          -- * Condition number
          conditionNumberSM,
          -- * Householder reflection
@@ -19,6 +21,10 @@
          _xCgne, _xTfq, _xBicgstab, _x, _xBcg,
          cgsStep, bicgstabStep,
          CGNE, TFQMR, BICGSTAB, CGS, BCG,
+         -- * Matrix partitioning
+         diagPartitions,
+         -- * Random arrays
+         randArray,
          -- * Random matrices and vectors
          randMat, randVec, 
          -- ** Sparse "
@@ -434,12 +440,15 @@
 -- * 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)
 
 
 
+ilu0 aa = (lh, uh) where
+  (l, u) = lu aa
+  lh = sparsifyLU l aa
+  uh = sparsifyLU u aa
+  sparsifyLU m m2 = SM (dim m) $ ifilterIM2 f (dat m) where
+    f i j _ = isJust (lookupSM m2 i j)
 
 
 
@@ -450,7 +459,33 @@
 
 
 
+-- * Preconditioning
 
+-- | Partition a matrix into strictly subdiagonal, diagonal and strictly superdiagonal parts
+diagPartitions :: SpMatrix a -> (SpMatrix a, SpMatrix a, SpMatrix a)
+diagPartitions aa = (e,d,f) where
+  e = extractSubDiag aa
+  d = extractDiag aa
+  f = extractSuperDiag aa
+
+
+-- ** SSOR
+
+-- | `mSsor aa omega` : if `omega = 1` it returns the diagonal of `aa`, 
+mSsor :: Fractional a => SpMatrix a -> a -> SpMatrix a
+mSsor aa omega = l ## r where
+  (e, d, f) = diagPartitions aa
+  n = nrows e
+  l = d ^-^ scale omega e
+  r = eye n ^-^ scale omega (reciprocal d ## f)
+
+
+
+
+
+
+
+
 -- * Iterative linear solvers
 
 
@@ -843,8 +878,16 @@
 
 
 
+-- * Random arrays
 
+randArray :: PrimMonad m => Int -> Double -> Double -> m [Double]
+randArray n mu sig = do
+  g <- MWC.create
+  replicateM n (MWC.normal mu sig g)
+  
 
+
+
 -- * Random matrices and vectors
 
 -- |Dense SpMatrix
@@ -884,4 +927,6 @@
   aav <- replicateM nsp (MWC.normal 0 1 g)
   ii <- replicateM nsp (MWC.uniformR (0, n-1) g :: IO Int)
   return $ fromListSV n $ zip ii aav
+
+
 
diff --git a/src/Numeric/LinearAlgebra/Sparse/IntMap.hs b/src/Numeric/LinearAlgebra/Sparse/IntMap.hs
deleted file mode 100644
--- a/src/Numeric/LinearAlgebra/Sparse/IntMap.hs
+++ /dev/null
@@ -1,184 +0,0 @@
-module Numeric.LinearAlgebra.Sparse.IntMap where
-
-import Numeric.LinearAlgebra.Class
-
-import qualified Data.IntMap.Strict as IM
-
-
-
-
-
-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
-
-
--- 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
-
--- |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 #-}  
-
--- |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)
-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/test/LibSpec.hs b/test/LibSpec.hs
--- a/test/LibSpec.hs
+++ b/test/LibSpec.hs
@@ -4,6 +4,10 @@
 import Numeric.LinearAlgebra.Sparse
 -- import Numeric.LinearAlgebra.Class
 
+import Control.Monad (liftM, liftM2)
+import Control.Monad.Primitive
+import Data.Foldable (foldrM)
+
 import Data.Sparse.Common
 
 
@@ -155,33 +159,63 @@
 
 
 
--- 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)
+-- -- 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
+-- -- 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
 
 
 
 
+-- solveRandomBanded n bw mu sig = do
+--   let ndiags = 2*bw
+--   bands <- replicateM (ndiags + 1) (randArray n mu sig)
+--   xtrue <- randVec n
+--   b <- randVec n
+--   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
+
+
+go (m:ms) mat =
+  m ^+^ go ms mat
+go [] mat = mat
+
+  
+plusM x y = return $ x ^+^ y
+
+
+
+
+
+
+
+
+
+
 --
 
 {-
@@ -213,25 +247,8 @@
 
 {- 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-}