sparse-linear-algebra 0.2.0.5 → 0.2.0.7

raw patch · 11 files changed

+673/−246 lines, 11 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• README.md +58/−0
• sparse-linear-algebra.cabal +3/−3
• src/Data/Sparse/Common.hs +192/−20
• src/Data/Sparse/SpMatrix.hs +67/−18
• src/Data/Sparse/SpVector.hs +6/−1
• src/Data/Sparse/Types.hs +7/−0
• src/Numeric/Eps.hs +3/−0
• src/Numeric/LinearAlgebra/Data.hs +0/−7
• src/Numeric/LinearAlgebra/Sparse.hs +245/−131
• src/Numeric/LinearAlgebra/Sparse/IntMap.hs +27/−22
• test/LibSpec.hs +65/−44

@@ -16,10 +16,14 @@      * BiConjugate Gradient Stabilized (BiCGSTAB) (non-Hermitian systems) +    * Transpose-Free Quasi-Minimal Residual (TFQMR)+ * Matrix decompositions      * QR factorization +    * LU factorization+ * Eigenvalue algorithms      * QR algorithm@@ -29,6 +33,60 @@ * Utilities : Vector and matrix norms, matrix condition number, Givens rotation, Householder reflection  * Predicates : Matrix orthogonality test (A^T A ~= I)+++---------++## Examples++The module `Numeric.LinearAlgebra.Sparse` contains the interface functions:++To create a sparse matrix from an array of its entries we use `fromListSM` :++    fromListSM :: Foldable t => (Int, Int) -> t (IxRow, IxCol, a) -> SpMatrix a++e.g.++    > amat = fromListSM (3,3) [(0,0,2),(1,0,4),(1,1,3),(1,2,2),(2,2,5)]++And similarly for sparse vectors : `fromListSV :: Int -> [(Int, a)] -> SpVector a`.++Both sparse vectors and matrices can be pretty-printed using `prd`:++    > prd amat+    ( 3 rows, 3 columns ) , 5 NZ ( sparsity 0.5555555555555556 )++    [2,0,0]+    [4,3,2]+    [0,0,5]+++Matrix factorizations are available as `lu` and `qr` respectively, and are straightforward to verify by using the matrix product `##`  :++    > (l, u) = lu amat+    > prd $ l ## u+    ( 3 rows, 3 columns ) , 9 NZ ( sparsity 1.0 )++    [2.0,0.0,0.0]+    [4.0,3.0,2.0]+    [0.0,0.0,5.0]++Linear systems can be solved with either `linSolve` (which also requires choosing a method) or with `<\>` (which uses BiCGSTAB as default) :++    > b = fromListSV 3 [(0,3),(1,2),(2,5)]+    > x = amat <\> b+    > prd x+    ( 3 elements ) ,  3 NZ ( sparsity 1.0 )++    [1.4999999999999998,-1.9999999999999998,0.9999999999999998]++The result can be verified by computing the matrix-vector action `amat #> x`, which should (ideally) be very close to the right-hand side `b` :++    > prd $ amat #> x+    ( 3 elements ) ,  3 NZ ( sparsity 1.0 )++    [2.9999999999999996,1.9999999999999996,4.999999999999999]+   

@@ -1,7 +1,7 @@ name:                sparse-linear-algebra-version:             0.2.0.5+version:             0.2.0.7 synopsis:            Numerical computation in native Haskell -description:         Currently it provides iterative linear solvers, matrix decompositions, eigenvalue computations and related utilities. Please see README.md for details+description:         Currently the library provides iterative linear solvers, matrix decompositions, eigenvalue computations and related utilities. Please see README.md for details homepage:            https://github.com/ocramz/sparse-linear-algebra license:             GPL-3 license-file:        LICENSE@@ -20,12 +20,12 @@   hs-source-dirs:      src   exposed-modules:     Numeric.LinearAlgebra.Sparse                        Numeric.LinearAlgebra.Class-                       Numeric.LinearAlgebra.Data                        Numeric.LinearAlgebra.Sparse.IntMap                        Data.Sparse.SpVector                        Data.Sparse.SpMatrix                        Data.Sparse.Common                        Data.Sparse.Utils+                       Data.Sparse.Types                        Numeric.Eps   build-depends:       QuickCheck                      , base >= 4.7 && < 5

@@ -1,28 +1,64 @@ module Data.Sparse.Common        ( module X,-         svToSM, outerProdSV, (><), toSV, extractCol, extractRow,-         extractDiagonalDSM,-         matVec, (#>), vecMat, (<#)) where+         insertRowWith, insertRow, insertColWith, insertCol,+         outerProdSV, (><), toSV, svToSM, +         extractCol, extractRow,+         extractVectorDenseWith, extractRowDense, extractColDense,+         extractDiagDense,+         extractSubRow, extractSubCol,+         extractSubRow_RK, extractSubCol_RK,+         matVec, (#>), vecMat, (<#),+         prd) where  import Data.Sparse.Utils as X+import Data.Sparse.Types as X import Data.Sparse.SpMatrix as X import Data.Sparse.SpVector as X  import Numeric.Eps as X-import Numeric.LinearAlgebra.Data as X import Numeric.LinearAlgebra.Class as X  import Numeric.LinearAlgebra.Sparse.IntMap as X  import qualified Data.IntMap as IM --- | promote a SV to SM-svToSM :: SpVector a -> SpMatrix a-svToSM (SV n d) = SM (n, 1) $ IM.singleton 0 d  --- ** Outer vector product+-- * Insert row/column vector in matrix +-- | Insert row , using the provided row index transformation function+insertRowWith :: (IxCol -> IxCol) -> SpMatrix a -> SpVector a -> IM.Key -> SpMatrix a+insertRowWith fj (SM (m,n) im) (SV d sv) i+  | not (inBounds0 m i) = error "insertRowSM : index out of bounds"+  | n >= d = SM (m,n) $ IM.insert i (insertOrUnion i sv' im) im+  | otherwise = error $ "insertRowSM : incompatible dimensions " ++ show (n, d)+    where sv' = IM.mapKeys fj sv+          insertOrUnion i sv im = maybe sv (IM.union sv) (IM.lookup i im)++-- | Insert row          +insertRow :: SpMatrix a -> SpVector a -> IM.Key -> SpMatrix a+insertRow = insertRowWith id++-- | Insert column, using the provided row index transformation function+insertColWith :: (IxRow -> IxRow) -> SpMatrix a -> SpVector a -> IxCol -> SpMatrix a+insertColWith fi smm sv j+  | not (inBounds0 n j) = error "insertColSM : index out of bounds"+  | m >= mv = insIM2 smm vl j+  | otherwise = error $ "insertColSM : incompatible dimensions " ++ show (m,mv) where+      (m, n) = dim smm+      mv = dim sv+      vl = toListSV sv+      insIM2 im2 ((i,x):xs) j = insIM2 (insertSpMatrix (fi i) j x im2) xs j+      insIM2 im2 [] _ = im2++-- | Insert column+insertCol :: SpMatrix a -> SpVector a -> IxCol -> SpMatrix a+insertCol = insertColWith id++++-- * Outer vector product+ outerProdSV, (><) :: Num a => SpVector a -> SpVector a -> SpMatrix a outerProdSV v1 v2 = fromListSM (m, n) ixy where   m = dim v1@@ -33,39 +69,79 @@   +-- * Matrix-vector conversions++-- | promote a SV to SM+svToSM :: SpVector a -> SpMatrix a+svToSM (SV n d) = SM (n, 1) $ IM.singleton 0 d+ -- |Demote (n x 1) or (1 x n) SpMatrix to SpVector toSV :: SpMatrix a -> SpVector a-toSV (SM (m,n) im) = SV d $ snd . head $ IM.toList im where+toSV (SM (m,n) im) = SV d (ff im) where+  ff | m > n = IM.map g  -- column case+     | otherwise = g+  g = snd . head . IM.toList   d | m==1 && n==1 = 1     | m==1 && n>1 = n      | n==1 && m>1 = m-    | otherwise = error $ "toSV : incompatible dimensions " ++ show (m,n)+    | otherwise = error $ "toSV : incompatible matrix dimension " ++ show (m,n) --- |Extract jth column, and place into SpVector-extractCol :: SpMatrix a -> IxCol -> SpVector a-extractCol m j = toSV $ extractColSM m j       +-- * Extract a SpVector from an SpMatrix+-- ** Sparse extract --- |Extract ith row, and place into SpVector+-- |Extract ith row extractRow :: SpMatrix a -> IxRow -> SpVector a extractRow m i = toSV $ extractRowSM m i +-- |Extract jth column+extractCol :: SpMatrix a -> IxCol -> SpVector a+extractCol m j = toSV $ extractColSM m j         +-- ** Dense extract (default == 0)+-- | Generic extraction function+extractVectorDenseWith ::+  Num a => (Int -> (IxRow, IxCol)) -> SpMatrix a -> SpVector a+extractVectorDenseWith f mm = fromListDenseSV n $ foldr ins [] ll  where+  ll = [0 .. n - 1]+  (m, n) = dim mm+  ins i acc = mm @@ f i : acc +-- | Extract ith row (dense)+extractRowDense :: Num a => SpMatrix a -> IxRow -> SpVector a+extractRowDense mm iref = extractVectorDenseWith (\j -> (iref, j)) mm +-- | Extract jth column+extractColDense :: Num a => SpMatrix a -> IxCol -> SpVector a+extractColDense mm jref = extractVectorDenseWith (\i -> (i, jref)) mm --- | Extract the diagonal as a SpVector (with default 0)-extractDiagonalDSM :: Num a => SpMatrix a -> SpVector a-extractDiagonalDSM mm = fromListDenseSV n $ foldr ins [] ll  where-  ll = [0 .. n - 1]-  n = nrows mm-  ins i acc = mm@@(i,i) : acc+-- | Extract the diagonal+extractDiagDense :: Num a => SpMatrix a -> SpVector a+extractDiagDense = extractVectorDenseWith (\i -> (i, i))   +-- | extract row interval+extractSubRow :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpVector a+extractSubRow m i (j1, j2)  = toSV $ extractSubRowSM m i (j1, j2) +-- | extract row interval, rebalance keys by subtracting lowest one+extractSubRow_RK :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpVector a+extractSubRow_RK m i (j1, j2)  = toSV $ extractSubRowSM_RK m i (j1, j2)++-- | extract column interval+extractSubCol :: SpMatrix a -> IxCol -> (IxRow, IxRow) -> SpVector a+extractSubCol m j (i1, i2)  = toSV $ extractSubColSM m j (i1, i2)++-- | extract column interval, rebalance keys by subtracting lowest one+extractSubCol_RK :: SpMatrix a -> IxCol -> (IxRow, IxRow) -> SpVector a+extractSubCol_RK m j (i1, i2)  = toSV $ extractSubColSM_RK m j (i1, i2)++++ -- ** Matrix action on a vector  {- @@ -92,3 +168,99 @@   | otherwise = error $ "vecMat : mismatching dimensions " ++ show (n, nr)  (<#) = vecMat  +++++++-- * Pretty printing+++showNonZero :: (Show a, Num a, Eq a) => a -> String+showNonZero x  = if x == 0 then " " else show x+++toDenseRow :: Num a => SpMatrix a -> IM.Key -> [a]+toDenseRow (SM (_,ncol) im) irow =+  fmap (\icol -> im `lookupWD_IM` (irow,icol)) [0..ncol-1]++toDenseRowClip :: (Show a, Num a) => SpMatrix a -> IM.Key -> Int -> String+toDenseRowClip sm irow ncomax+  | ncols sm > ncomax = unwords (map show h) ++  " ... " ++ show t+  | otherwise = show dr+     where dr = toDenseRow sm irow+           h = take (ncomax - 2) dr+           t = last dr++newline :: IO ()+newline = putStrLn ""++printDenseSM :: (Show t, Num t) => SpMatrix t -> IO ()+printDenseSM sm = do+  newline+  putStrLn $ sizeStr sm+  newline+  printDenseSM' sm 5 5+  newline+  where    +    printDenseSM' :: (Show t, Num t) => SpMatrix t -> Int -> Int -> IO ()+    printDenseSM' sm'@(SM (nr,_) _) nromax ncomax = mapM_ putStrLn rr_' where+      rr_ = map (\i -> toDenseRowClip sm' i ncomax) [0..nr - 1]+      rr_' | nrows sm > nromax = take (nromax - 2) rr_ ++ [" ... "] ++[last rr_]+           | otherwise = rr_+++toDenseListClip :: (Show a, Num a) => SpVector a -> Int -> String+toDenseListClip sv ncomax+  | dim sv > ncomax = unwords (map show h) ++  " ... " ++ show t+  | otherwise = show dr+     where dr = toDenseListSV sv+           h = take (ncomax - 2) dr+           t = last dr++printDenseSV :: (Show t, Num t) => SpVector t -> IO ()+printDenseSV sv = do+  newline+  putStrLn $ sizeStrSV sv+  newline+  printDenseSV' sv 5+  newline where+    printDenseSV' v nco = putStrLn rr_' where+      rr_ = toDenseListClip v nco :: String+      rr_' | dim sv > nco = unwords [take (nco - 2) rr_ , " ... " , [last rr_]]+           | otherwise = rr_++-- ** Pretty printer typeclass+class PrintDense a where+  prd :: a -> IO ()++instance (Show a, Num a) => PrintDense (SpVector a) where+  prd = printDenseSV++instance (Show a, Num a) => PrintDense (SpMatrix a) where+  prd = printDenseSM++++  ++++++++++++++++++++  

@@ -2,12 +2,10 @@ module Data.Sparse.SpMatrix where  import Data.Sparse.Utils--import Data.Sparse.Utils+import Data.Sparse.Types  import Numeric.Eps import Numeric.LinearAlgebra.Class-import Numeric.LinearAlgebra.Data  import Numeric.LinearAlgebra.Sparse.IntMap @@ -184,19 +182,23 @@  -- | Looks up an element in the matrix with a default (if the element is not found, zero is returned) -lookupWD_SM, (@@) :: Num a => SpMatrix a -> (IxRow, IxCol) -> a+lookupWD_SM, (@@!), (@@) :: Num a => SpMatrix a -> (IxRow, IxCol) -> a lookupWD_SM sm (i,j) =   fromMaybe 0 (lookupSM sm i j)  lookupWD_IM :: Num a => IM.IntMap (IM.IntMap a) -> (IxRow, IxCol) -> a lookupWD_IM im (i,j) = fromMaybe 0 (IM.lookup i im >>= IM.lookup j) --- | Zero-default lookup, infix form-(@@) = lookupWD_SM+-- | Zero-default lookup, infix form (no bound checking)+(@@!) = lookupWD_SM +-- | Zero-default lookup, infix form ("safe" : throws exception if lookup is outside matrix bounds)+m @@ d | isValidIxSM m d = m @@! d+       | otherwise = error $ "@@ : incompatible indices : matrix size is " ++ show (dim m) ++ ", but user looked up " ++ show d   + -- FIXME : to throw an exception or just ignore the out-of-bound access ?  @@ -207,12 +209,18 @@  -- ** Sub-matrices --- | Extract a submatrix given the specified index bounds-extractSubmatrixSM :: SpMatrix a -> (IxRow, IxCol) -> (IxRow, IxCol) -> SpMatrix a-extractSubmatrixSM (SM (r, c) im) (i1, i2) (j1, j2)++-- | Extract a submatrix given the specified index bounds, rebalancing keys with the two supplied functions+extractSubmatrixSM ::+  (IM.Key -> IM.Key) ->   -- row index function+  (IM.Key -> IM.Key) ->   -- column "  "+  SpMatrix a ->+  (IxRow, IxRow) -> (IxCol, IxCol) ->+  SpMatrix a+extractSubmatrixSM fi gi (SM (r, c) im) (i1, i2) (j1, j2)   | q = SM (m', n') imm'-  | otherwise = error $ "extractSubmatrixSM : invalid indexing " ++ show (i1, i2) ++ ", " ++ show (j1, j2) where-  imm' = mapKeysIM2 (\i -> i - i1) (\j -> j - j1) $  -- rebalance keys+  | otherwise = error $ "extractSubmatrixSM : invalid index " ++ show (i1, i2) ++ ", " ++ show (j1, j2) where+  imm' = mapKeysIM2 fi gi $                          -- rebalance keys           IM.filter (not . IM.null) $                -- remove all-null rows           ifilterIM2 ff im                           -- keep `submatrix`   ff i j _ = i1 <= i &&@@ -226,19 +234,52 @@       inBounds0 c j2 &&             i2 >= i1 +-- | Extract a submatrix given the specified index bounds+-- NB : subtracts (i1, j1) from the indices+extractSubmatrixRebalanceKeys ::+  SpMatrix a -> (IxRow, IxRow) -> (IxCol, IxCol) -> SpMatrix a+extractSubmatrixRebalanceKeys mm (i1,i2) (j1,j2) =+  extractSubmatrixSM (\i -> i - i1) (\j -> j - j1) mm (i1,i2) (j1,j2) +-- | Extract a submatrix given the specified index bounds+-- NB : submatrix indices are _preserved_+extractSubmatrix :: SpMatrix a -> (IxRow, IxRow) -> (IxCol, IxCol) -> SpMatrix a+extractSubmatrix = extractSubmatrixSM id id  +++-- *** Extract i'th row+extractRowSM :: SpMatrix a -> IxRow -> SpMatrix a+extractRowSM sm i = extractSubmatrix sm (i, i) (0, ncols sm - 1)+++-- | Extract column within a row range+extractSubRowSM :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpMatrix a+extractSubRowSM sm i (j1, j2) = extractSubmatrix sm (i, i) (j1, j2)++-- | Extract column within a row range, rebalance keys+extractSubRowSM_RK :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpMatrix a+extractSubRowSM_RK sm i =+  extractSubmatrixRebalanceKeys sm (i, i) +++ -- *** Extract j'th column+-- | Extract all column extractColSM :: SpMatrix a -> IxCol -> SpMatrix a-extractColSM sm j = extractSubmatrixSM sm (0, nrows sm - 1) (j, j)+extractColSM sm j = extractSubmatrix sm (0, nrows sm - 1) (j, j) +-- | Extract column within a row range+extractSubColSM :: SpMatrix a -> IxCol -> (IxRow, IxRow) -> SpMatrix a+extractSubColSM sm j (i1, i2) = extractSubmatrix sm (i1, i2) (j, j)  +-- | Extract column within a row range, rebalance keys+extractSubColSM_RK :: SpMatrix a -> IxCol -> (IxRow, IxRow) -> SpMatrix a+extractSubColSM_RK sm j (i1, i2) =+  extractSubmatrixRebalanceKeys sm (i1, i2) (j, j)  --- *** Extract i'th row-extractRowSM :: SpMatrix a -> IxRow -> SpMatrix a-extractRowSM sm i = extractSubmatrixSM sm (i, i) (0, ncols sm - 1)   @@ -255,8 +296,8 @@  -- ** Predicates -- |Are the supplied indices within matrix bounds?-validIxSM :: SpMatrix a -> (Int, Int) -> Bool-validIxSM mm = inBounds02 (dim mm)+isValidIxSM :: SpMatrix a -> (Int, Int) -> Bool+isValidIxSM mm = inBounds02 (dim mm)  -- |Is the matrix square? isSquareSM :: SpMatrix a -> Bool@@ -626,4 +667,12 @@   -+-- *** Matrix contraction+-- | Contract two matrices A and B up to an index `n`, i.e. summing over repeated indices: +-- Aij Bjk , for j in [0 .. n] +contractSub :: Num a => SpMatrix a -> SpMatrix a -> IxRow -> IxCol -> Int -> a+contractSub a b i j n+  | ncols a == nrows b &&+    isValidIxSM a (i,j) &&+    n <= ncols a = sum $ map (\i' -> a@@!(i,i')*b@@!(i',j)) [0 .. n]+  | otherwise = error "contractSub : n must be <= i"

@@ -2,9 +2,9 @@ module Data.Sparse.SpVector where  import Data.Sparse.Utils+import Data.Sparse.Types  import Numeric.LinearAlgebra.Class-import Numeric.LinearAlgebra.Data import Numeric.LinearAlgebra.Sparse.IntMap  import Data.Maybe@@ -21,8 +21,13 @@ spySV :: Fractional b => SpVector a -> b spySV s = fromIntegral (IM.size (dat s)) / fromIntegral (dim s) +-- | Number of nonzeros+nzSV :: SpVector a -> Int+nzSV sv = IM.size (dat sv)  +sizeStrSV :: SpVector a -> String+sizeStrSV sv = unwords ["(",show (dim sv),"elements ) , ",show (nzSV sv),"NZ ( sparsity", show (spy sv),")"]   

@@ -0,0 +1,7 @@+module Data.Sparse.Types where++type Rows = Int+type Cols = Int++type IxRow = Int+type IxCol = Int

@@ -15,6 +15,9 @@ almostZero x = abs x <= eps almostOne x = x >= (1-eps) && x < (1+eps) +isNz :: Double -> Bool+isNz = not . almostZero+ withDefault :: (t -> Bool) -> t -> t -> t withDefault q d x | q x = d                   | otherwise = x

@@ -1,7 +0,0 @@-module Numeric.LinearAlgebra.Data where--type Rows = Int-type Cols = Int--type IxRow = Int-type IxCol = Int

@@ -1,6 +1,23 @@ {-# LANGUAGE FlexibleContexts, TypeFamilies, MultiParamTypeClasses, FlexibleInstances  #-} -- {-# OPTIONS_GHC -O2 -rtsopts -with-rtsopts=-K32m -prof#-}-module Numeric.LinearAlgebra.Sparse where+module Numeric.LinearAlgebra.Sparse+       (+         sparsifySV,+         conditionNumberSM,+         hhMat, hhRefl,+         givens,+         eigsQR, eigRayleigh,+         cgne, tfqmr, bicgstab, cgs, bcg,+         _xCgne, _xTfq, _xBicgstab, _x, _xBcg,+         cgsStep, bicgstabStep,+         CGNE, TFQMR, BICGSTAB, CGS, BCG,+         linSolve, LinSolveMethod, (<\>),+         randMat, randVec, randSpMat, randSpVec,+         modifyInspectN, runAppendN',+         diffSqL,+         qr, lu+       )+       where   import Data.Sparse.Common@@ -9,7 +26,7 @@ import Control.Monad.Primitive import Control.Monad (mapM_, forM_, replicateM) import Control.Monad.State.Strict--- import Control.Monad.Writer+import Control.Monad.Writer -- import Control.Monad.Trans -- import Control.Monad.Trans.State (runStateT) -- import Control.Monad.Trans.Writer (runWriterT)@@ -49,7 +66,7 @@                     | otherwise = kappa where   kappa = lmax / lmin   (_, r) = qr m-  u = extractDiagonalDSM r  -- FIXME : need to extract with default element 0 +  u = extractDiagDense r  -- FIXME : need to extract with default element 0    lmax = abs (maximum u)   lmin = abs (minimum u) @@ -118,7 +135,7 @@  givens :: SpMatrix Double -> IxRow -> IxCol -> SpMatrix Double givens mm i j -  | validIxSM mm (i,j) && isSquareSM mm =+  | isValidIxSM mm (i,j) && isSquareSM mm =        sparsifySM $ fromListSM' [(i,i,c),(j,j,c),(j,i,-s),(i,j,s)] (eye (nrows mm))   | otherwise = error "givens : indices out of bounds"         where@@ -180,15 +197,15 @@  -- * Eigenvalue algorithms --- ** All eigenvalues (QR algorithm)+-- ** QR algorithm --- | `eigsQR n mm` performs `n` iterations of the QR algorithm on matrix `mm` +-- | `eigsQR n mm` performs `n` iterations of the QR algorithm on matrix `mm`, and returns a SpVector containing all eigenvalues eigsQR :: Int -> SpMatrix Double -> SpVector Double-eigsQR nitermax m = extractDiagonalDSM $ execState (convergtest eigsStep) m where+eigsQR nitermax m = extractDiagDense $ execState (convergtest eigsStep) m where   eigsStep m = r #~# q where (q, r) = qr m   convergtest g = modifyInspectN nitermax f g where-    f [m1, m2] = let dm1 = extractDiagonalDSM m1-                     dm2 = extractDiagonalDSM m2+    f [m1, m2] = let dm1 = extractDiagDense m1+                     dm2 = extractDiagDense m2                  in norm2 (dm1 ^-^ dm2) <= eps  @@ -196,10 +213,9 @@   --- ** One eigenvalue and eigenvector (Rayleigh iteration)-+-- ** Rayleigh iteration --- | `eigsRayleigh n mm` performs `n` iterations of the Rayleigh algorithm on matrix `mm`. Cubic-order convergence, but it requires a mildly educated guess on the initial eigenpair+-- | `eigsRayleigh n mm` performs `n` iterations of the Rayleigh algorithm on matrix `mm` and returns the eigenpair closest to the initialization. It displays cubic-order convergence, but it also requires an educated guess on the initial eigenpair eigRayleigh :: Int                -- max # iterations      -> SpMatrix Double           -- matrix      -> (SpVector Double, Double) -- initial guess of (eigenvector, eigenvalue)@@ -272,22 +288,131 @@   -- * LU factorization-+-- ** Doolittle algorithm {- Doolittle algorithm for factoring A' = P A, where P is a permutation matrix such that A' has a nonzero as its (0, 0) entry -} +-- | LU factors+lu :: SpMatrix Double -> (SpMatrix Double, SpMatrix Double)+lu aa = (lfin, ufin) where+  (ixf,lf,uf) = execState (modifyUntil q (luUpd aa)) (luInit aa)+  lfin = lf+  ufin = uUpd aa (ixf, lf, uf)+  q (i, _, _) = i == (nrows aa - 1) --- lu aa | isSquareSM aa = undefined---       | otherwise = error "LU factorization not currently defined for rectangular matrices" where---           n = nrows aa---           l0 = eye n---           aa0 = +-- | First iteration of LU+luInit ::+  (Num t, Fractional a) => SpMatrix a -> (t, SpMatrix a, SpMatrix a)+luInit aa = (1, l0, u0) where+  n = nrows aa+  l0 = insertCol (eye n) ((1/u00) .* extractSubCol aa 0 (1,n - 1)) 0  -- initial L+  u0 = insertRow (zeroSM n n) (extractRow aa 0) 0               -- initial U+  u00 = u0 @@ (0,0)  -- make sure this is non-zero by applying permutation --- luStep aa i l u +-- | LU update step+luUpd :: SpMatrix Double+     -> (Int, SpMatrix Double, SpMatrix Double)+     -> (Int, SpMatrix Double, SpMatrix Double)+luUpd aa (i, l, u) = (i', l', u') where+  n = nrows aa  +  u' = uUpdSparse aa (i, l, u)  -- update U+  l' = lUpdSparse aa (i, l, u') -- update L+  i' = i + 1     -- increment i  +uUpd' ::+  Num a =>+  ([(Int, a)] -> [(Int, a)]) ->+  SpMatrix a ->+  (Rows, SpMatrix a, SpMatrix a) ->+  SpMatrix a+uUpd' ff amat (ix, lmat, umat) = insertRow umat uv ix where+  n = nrows amat+  colsix = [ix .. n - 1]+  us = ff $ zip colsix $ map (solveForUij amat lmat umat ix) colsix+  uv = fromListSV n us++uUpd :: Num a => SpMatrix a -> (Rows, SpMatrix a, SpMatrix a) -> SpMatrix a+uUpd = uUpd' id++-- update U while sparsifying+uUpdSparse ::+  SpMatrix Double -> (Rows, SpMatrix Double, SpMatrix Double) -> SpMatrix Double+uUpdSparse = uUpd' (filter (isNz . snd))+++++-- solve for element Uij+solveForUij ::+  Num a => SpMatrix a -> SpMatrix a -> SpMatrix a -> IxRow -> IxCol -> a+solveForUij amat lmat umat i j = a - p where+  a = amat @@! (i, j)+  p = contractSub lmat umat i j (i - 1)+++-- solve for element Lij+solveForLij ::+  SpMatrix Double -> SpMatrix Double -> SpMatrix Double -> IxRow -> IxCol -> Double+solveForLij amat lmat umat i j | isNz uii = (a - p)/uii+                               | otherwise = error $ unwords ["solveForLij : U",show (i,i)," is close to 0. Permute rows in order to have a nonzero diagonal of U"]+  where+   a = amat @@! (i, j)+   uii = umat @@! (i-1, i-1) -- NB this must be /= 0+   p = contractSub lmat umat i j (i - 1)+++++lUpd' :: ([(Rows, Double)] -> [(Int, Double)])+     -> SpMatrix Double+     -> (Rows, SpMatrix Double, SpMatrix Double)+     -> SpMatrix Double+lUpd' ff amat (ix, lmat, umat) = insertCol lmat lv ix where+  n = nrows amat+  rowsix = [ix + 1 .. n - 1]+  ls = ff $ zip rowsix $ map (\i -> solveForLij amat lmat umat i ix) rowsix+  lv = fromListSV n ls++lUpd :: SpMatrix Double -> (Rows, SpMatrix Double, SpMatrix Double) -> SpMatrix Double+lUpd = lUpd' id++lUpdSparse ::+  SpMatrix Double -> (Rows, SpMatrix Double, SpMatrix Double) -> SpMatrix Double+lUpdSparse = lUpd' (filter (isNz . snd))++++++++++++++++++++++-- test data++aa0 = fromListDenseSM 3 [2, -4, -4, -1, 6, -2, -2, 3, 8] :: SpMatrix Double+++++++ -- Produces the permutation matrix necessary to have a nonzero in position (iref, jref). This is used in the LU factorization-permutAA :: Num b => SpMatrix a -> IxRow -> IxCol -> Maybe (SpMatrix b)-permutAA (SM (nro,_) mm) iref jref+permutAA :: Num b => IxRow -> IxCol -> SpMatrix a -> Maybe (SpMatrix b)+permutAA iref jref (SM (nro,_) mm)    | isJust (lookupIM2 iref jref mm) = Nothing -- eye nro   | otherwise = Just $ permutationSM nro [head u] where       u = IM.keys (ifilterIM2 ff mm)@@ -302,19 +427,106 @@   +++-- * Incomplete LU++-- | used for Incomplete LU : remove entries in `m` corresponding to zero entries in `m2`+ripHoles :: SpMatrix t -> SpMatrix a -> SpMatrix a+ripHoles (SM d m) m2 = SM d $ ifilterIM2 f (dat m2) where+  f i j _ = isJust (lookupSM m2 i j)++++++++++++++ -- * Iterative linear solvers   --- -- | residual of candidate solution x0 of a linear system--- residual :: Num a => SpMatrix a -> SpVector a -> SpVector a -> SpVector a--- residual aa b x0 = b ^-^ (aa #> x0) --- converged :: SpMatrix Double -> SpVector Double -> SpVector Double -> Bool--- converged aa b x0 = normSq (residual aa b x0) <= eps+-- ** CGNE +cgneStep :: SpMatrix Double -> CGNE -> CGNE+cgneStep aa (CGNE x r p) = CGNE x1 r1 p1 where+  alphai = r `dot` r / (p `dot` p)+  x1 = x ^+^ (alphai .* p)+  r1 = r ^-^ (alphai .* (aa #> p))+  beta = r1 `dot` r1 / (r `dot` r)+  p1 = transposeSM aa #> r ^+^ (beta .* p) +data CGNE =+  CGNE {_xCgne , _rCgne, _pCgne :: SpVector Double} deriving Eq +cgne :: SpMatrix Double -> SpVector Double -> SpVector Double -> CGNE+cgne aa b x0 = execState (untilConverged _xCgne (cgneStep aa)) cgneInit where+  r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution+  p0 = transposeSM aa #> r0+  cgneInit = CGNE x0 r0 p0+++-- ** TFQMR++-- | one step of TFQMR+tfqmrStep :: SpMatrix Double -> SpVector Double -> TFQMR -> TFQMR+tfqmrStep aa r0hat (TFQMR x w u v d m tau theta eta rho alpha) =+  TFQMR x1 w1 u1 v1 d1 (m+1) tau1 theta1 eta1 rho1 alpha1+  where+  -- alpham = alpha+  w1 = w ^-^ (alpha .* (aa #> u))+  d1 = u ^+^ ((theta**2/alpha*eta) .* d)+  theta1 = norm2 w1 / tau+  c = recip $ sqrt (1 + theta1**2)+  tau1 = tau * theta1 * c+  eta1 = c**2 * alpha+  x1 = x^+^ (eta1 .* d1)+  (alpha1, u1, rho1, v1)+    | even m = let+                   alpha' = rho / (v `dot` r0hat)+                   u' = u ^-^ (alpha' .* v)+               in+                   (alpha', u', rho, v)+    | otherwise = let+                   rho' = w1 `dot` r0hat+                   beta = rho'/rho+                   u' = w1 ^+^ (beta .* u)+                   v' = (aa #> u') ^+^ (beta .* (aa #> u ^+^ (beta .* v)) )+                  in (alpha, u', rho', v')++tfqmr :: SpMatrix Double -> SpVector Double -> SpVector Double -> TFQMR  +tfqmr aa b x0 = execState (untilConverged _xTfq (tfqmrStep aa r0)) tfqmrInit where+  n = dim b+  r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution+  w0 = r0+  u0 = r0+  v0 = aa #> u0+  d0 = zeroSV n+  r0hat = r0+  rho0 = r0hat `dot` r0+  alpha0 = rho0 / (v0 `dot` r0hat)+  m = 0+  tau0 = norm2 r0+  theta0 = 0+  eta0 = 0+  tfqmrInit = TFQMR x0 w0 u0 v0 d0 m tau0 theta0 eta0 rho0 alpha0++data TFQMR =+  TFQMR { _xTfq, _wTfq, _uTfq, _vTfq, _dTfq :: SpVector Double,+          _mTfq :: Int,+          _tauTfq, _thetaTfq, _etaTfq, _rhoTfq, _alphaTfq :: Double}+  deriving Eq+++ -- ** BCG  -- | one step of BCG@@ -449,7 +661,7 @@  -- * Linear solver interface -data LinSolveMethod = BCG_ | CGS_ | BICGSTAB_ deriving (Eq, Show) +data LinSolveMethod = CGNE_ | TFQMR_ | BCG_ | CGS_ | BICGSTAB_ deriving (Eq, Show)   -- -- | Linear solve with _random_ starting vector -- linSolveM ::@@ -472,9 +684,11 @@       solve aa' b' | isDiagonalSM aa' = reciprocal aa' #> b' -- diagonal solve is easy                    | otherwise = solveWith aa' b'        solveWith aa' b' = case method of-                                BCG_ -> _xBcg (bcg aa' b' x0)-                                CGS_ ->  _xBicgstab (bicgstab aa' b' x0 x0)-                                BICGSTAB_ -> _x (cgs aa' b' x0 x0)+        CGNE_ -> _xCgne (cgne aa' b' x0)+        TFQMR_ -> _xTfq (tfqmr aa' b' x0)+        BCG_ -> _xBcg (bcg aa' b' x0)+        CGS_ ->  _xBicgstab (bicgstab aa' b' x0 x0)+        BICGSTAB_ -> _x (cgs aa' b' x0 x0)       x0 = mkSpVectorD n $ replicate n 0.1        (m, n) = dim aa       nb     = dim b@@ -584,7 +798,7 @@  -- | iterate until convergence is verified or we run out of a fixed iteration budget untilConverged :: MonadState a m => (a -> SpVector Double) -> (a -> a) -> m a-untilConverged fproj = modifyInspectN 100 (normDiffConverged fproj)+untilConverged fproj = modifyInspectN 200 (normDiffConverged fproj)  -- | convergence check (FIXME) normDiffConverged :: (Foldable t, Functor t) =>@@ -669,103 +883,3 @@   ii <- replicateM nsp (MWC.uniformR (0, n-1) g :: IO Int)   return $ fromListSV n $ zip ii aav ------- * Pretty printing------showNonZero :: (Show a, Num a, Eq a) => a -> String-showNonZero x  = if x == 0 then " " else show x---toDenseRow :: Num a => SpMatrix a -> IM.Key -> [a]-toDenseRow (SM (_,ncol) im) irow =-  fmap (\icol -> im `lookupWD_IM` (irow,icol)) [0..ncol-1]--toDenseRowClip :: (Show a, Num a) => SpMatrix a -> IM.Key -> Int -> String-toDenseRowClip sm irow ncomax-  | ncols sm > ncomax = unwords (map show h) ++  " ... " ++ show t-  | otherwise = show dr-     where dr = toDenseRow sm irow-           h = take (ncomax - 2) dr-           t = last dr--newline :: IO ()-newline = putStrLn ""--printDenseSM :: (Show t, Num t) => SpMatrix t -> IO ()-printDenseSM sm = do-  newline-  putStrLn $ sizeStr sm-  newline-  printDenseSM' sm 5 5-  newline-  where    -    printDenseSM' :: (Show t, Num t) => SpMatrix t -> Int -> Int -> IO ()-    printDenseSM' sm'@(SM (nr,_) _) nromax ncomax = mapM_ putStrLn rr_' where-      rr_ = map (\i -> toDenseRowClip sm' i ncomax) [0..nr - 1]-      rr_' | nrows sm > nromax = take (nromax - 2) rr_ ++ [" ... "] ++[last rr_]-           | otherwise = rr_---toDenseListClip :: (Show a, Num a) => SpVector a -> Int -> String-toDenseListClip sv ncomax-  | dim sv > ncomax = unwords (map show h) ++  " ... " ++ show t-  | otherwise = show dr-     where dr = toDenseListSV sv-           h = take (ncomax - 2) dr-           t = last dr--printDenseSV :: (Show t, Num t) => SpVector t -> IO ()-printDenseSV sv = do-  newline-  printDenseSV' sv 5-  newline where-    printDenseSV' v nco = putStrLn rr_' where-      rr_ = toDenseListClip v nco :: String-      rr_' | dim sv > nco = unwords [take (nco - 2) rr_ , " ... " , [last rr_]]-           | otherwise = rr_---- ** Pretty printer typeclass-class PrintDense a where-  prd :: a -> IO ()--instance (Show a, Num a) => PrintDense (SpVector a) where-  prd = printDenseSV--instance (Show a, Num a) => PrintDense (SpMatrix a) where-  prd = printDenseSM----  ------------- * Type synonyms for SpMatrix--------------  

@@ -35,23 +35,33 @@   --- | ========= IntMap-of-IntMap (IM2) stuff+-- * set-like brackets  --- insert an element+-- unionWithKeyIM2 f im1 im2 = undefined where++++++-- * Insertion++-- | Insert an element insertIM2 ::   IM.Key -> IM.Key -> a -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a) insertIM2 i j x imm = IM.insert i ro imm where   ro = maybe (IM.singleton j x) (IM.insert j x) (IM.lookup i imm) {-# inline insertIM2 #-}   --- lookup a key+-- * Lookup++-- |Lookup a key lookupIM2 ::   IM.Key -> IM.Key -> IM.IntMap (IM.IntMap a) -> Maybe a lookupIM2 i j imm = IM.lookup i imm >>= IM.lookup j {-# inline lookupIM2 #-}   --- populate an IM2 from a list of (row index, column index, value)  +-- |Ppopulate an IM2 from a list of (row index, column index, value)   fromListIM2 ::   Foldable t =>      t (IM.Key, IM.Key, a) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)@@ -59,15 +69,16 @@   ins t (i,j,x) = insertIM2 i j x t  --- | folding+-- * folding --- indexed fold over an IM2+-- |Indexed left fold over an IM2, with general accumulator ifoldlIM2' :: (IM.Key -> IM.Key -> a -> b -> b) -> b -> IM.IntMap (IM.IntMap a) -> b ifoldlIM2' f empty mm = IM.foldlWithKey' accRow empty mm where   accRow acc i r = IM.foldlWithKey' (accElem i) acc r   accElem i acc j x = f i j x acc {-# inline ifoldlIM2' #-} +-- |Indexed left fold over an IM2 ifoldlIM2 ::   (IM.Key -> IM.Key -> t -> IM.IntMap a -> IM.IntMap a) ->   IM.IntMap (IM.IntMap t) ->  @@ -77,6 +88,7 @@   accElem i acc j x   = f i j x acc {-# inline ifoldlIM2 #-}   +-- |Left fold over an IM2, with general accumulator foldlIM2 :: (a -> b -> b) -> b -> IM.IntMap (IM.IntMap a) -> b foldlIM2 f empty mm = IM.foldl accRow empty mm where   accRow acc r = IM.foldl accElem acc r@@ -84,7 +96,7 @@ {-# inline foldlIM2 #-}  --- transposeIM2 : inner indices become outer ones and vice versa. No loss of information because both inner and outer IntMaps are nubbed.+-- | Inner indices become outer ones and vice versa. No loss of information because both inner and outer IntMaps are nubbed. transposeIM2 :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a) transposeIM2 = ifoldlIM2 (flip insertIM2) {-# inline transposeIM2 #-}@@ -98,9 +110,9 @@   --- | filtering+-- * filtering --- map over outer IM and filter all inner IM's+-- |Map over outer IM and filter all inner IM's ifilterIM2 ::   (IM.Key -> IM.Key -> a -> Bool) ->   IM.IntMap (IM.IntMap a) ->@@ -109,9 +121,7 @@   IM.mapWithKey (\irow row -> IM.filterWithKey (f irow) row)  {-# inline ifilterIM2 #-} --- specialized filtering function---- keep only sub-diagonal elements+-- |Specialized filtering : keep only sub-diagonal elements filterSubdiag :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a) filterSubdiag = ifilterIM2 (\i j _ -> i>j) @@ -119,7 +129,7 @@ countSubdiagonalNZ im =   IM.size $ IM.filter (not . IM.null) (filterSubdiag im) --- list of (row, col) indices of (nonzero) subdiagonal elements+-- |List of (row, col) indices of (nonzero) subdiagonal elements subdiagIndices :: IM.IntMap (IM.IntMap a) -> [(IM.Key, IM.Key)] subdiagIndices im = concatMap rpairs $ IM.toList (IM.map IM.keys im') where   im' = filterSubdiag im@@ -137,15 +147,14 @@     --- | mapping---- map over IM2+-- * mapping +-- |Map over IM2 mapIM2 :: (a -> b) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap b) mapIM2 = IM.map . IM.map   -- imapIM2 (\_ _ x -> f x)  --- indexed map over IM2+-- |Indexed map over IM2 imapIM2 ::   (IM.Key -> IM.Key -> a -> b) ->   IM.IntMap (IM.IntMap a) ->@@ -155,8 +164,7 @@   --- mapping keys-+-- |Mapping keys mapKeysIM2 ::   (IM.Key -> IM.Key) -> (IM.Key -> IM.Key) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a) mapKeysIM2 fi fj im = IM.map adjCols adjRows where@@ -173,7 +181,4 @@   ----- sparsification : 

@@ -28,7 +28,7 @@  spec :: Spec spec = do-  describe "Math.Linear.Sparse : library" $ do+  describe "Numeric.LinearAlgebra.Sparse : library" $ do     -- prop "subtraction is cancellative" $ \(x :: SpVector Double) ->     --   x ^-^ x `shouldBe` zero     it "dot : inner product" $@@ -43,58 +43,44 @@       (m1 `matMat` m2) `shouldBe` m1m2     it "eye : identity matrix" $       infoSM (eye 10) `shouldBe` SMInfo 10 0.1+    it "insertCol : insert a column in a SpMatrix" $+      insertCol (eye 3) (fromListDenseSV 3 [2,2,2]) 0 `shouldBe` (fromListSM (3,3) [(0,0,2),(1,0,2),(1,1,1),(2,0,2),(2,2,1)])+    it "insertRow : insert a row in a SpMatrix" $+      insertRow (eye 3) (fromListDenseSV 3 [2,2,2]) 1 `shouldBe` (fromListSM (3,3) [(0,0,1), (1,0,2), (1,1,2), (1,2,2), (2,2,1)])+    it "extractCol -> insertCol : identity" $+      insertCol (eye 3) (extractCol (eye 3) 1) 1 `shouldBe` eye 3+    it "extractRow -> insertRow : identity" $+      insertRow (eye 3) (extractRow (eye 3) 1) 1 `shouldBe` eye 3           it "countSubdiagonalNZ : # of nonzero elements below the diagonal" $       countSubdiagonalNZSM m3 `shouldBe` 1+    it "permutPairsSM : permutation matrices are orthogonal" $ do+      let pm0 = permutPairsSM 3 [(0,2), (1,2)]+      pm0 ##^ pm0 `shouldBe` eye 3+      pm0 #^# pm0 `shouldBe` eye 3              it "modifyInspectN : early termination by iteration count" $       execState (modifyInspectN 2 ((< eps) . diffSqL) (/2)) 1 `shouldBe` 1/8     it "modifyInspectN : termination by value convergence" $       execState (modifyInspectN (2^16) ((< eps) . head) (/2)) 1 < eps `shouldBe` True -  describe "Math.Linear.Sparse : Linear solvers" $ do+  describe "Numeric.LinearAlgebra.Sparse : Linear solvers" $ do+    -- it "TFQMR (2 x 2 dense)" $+    --   normSq (_xTfq (tfqmr aa0 b0 x0) ^-^ x0true) <= eps `shouldBe` True     it "BCG (2 x 2 dense)" $       normSq (_xBcg (bcg aa0 b0 x0) ^-^ x0true) <= eps `shouldBe` True     it "BiCGSTAB (2 x 2 dense)" $        normSq (aa0 <\> b0 ^-^ x0true) <= eps `shouldBe` True     it "CGS (2 x 2 dense)" $        normSq (_x (cgs aa0 b0 x0 x0) ^-^ x0true) <= eps `shouldBe` True-  describe "Math.Linear.Sparse : QR decomposition" $ do    +  describe "Numeric.LinearAlgebra.Sparse : QR decomposition" $ do         it "QR (4 x 4 sparse)" $       checkQr tm4 `shouldBe` True     it "QR (3 x 3 dense)" $        checkQr tm2 `shouldBe` True-    -  -- let n = 10-  --     nsp = 3-  -- describe ("random sparse linear system of size " ++ show n ++ " and sparsity " ++ show (fromIntegral nsp/fromIntegral n)) $ it "<\\>" $ do-  --   aa <- randSpMat n nsp-  --   xtrue <- randSpVec n nsp-  --   b <- randSpVec n nsp    -  --   let b = aa #> xtrue-  --   printDenseSM aa-  --   normSq (aa <\> b ^-^ xtrue) <= eps `shouldBe` True-  -- --     normSq (_xBicgstab (bicgstab aa b x0 x0) ^-^ x) <= eps `shouldBe` True+  describe "Numeric.LinearAlgebra.Sparse : LU decomposition" $ do+    it "LU (3 x 3 dense)" $+      checkLu tm2 `shouldBe` True   --- -- run N iterations ---- runNBiC :: Int -> SpMatrix Double -> SpVector Double -> BICGSTAB-runNBiC n aa b = map _xBicgstab $ runAppendN' (bicgstabStep aa x0) n bicgsInit where-   x0 = mkSpVectorD nd $ replicate nd 0.9-   nd = dim r0-   r0 = b ^-^ (aa #> x0)    -   p0 = r0-   bicgsInit = BICGSTAB x0 r0 p0---- runNCGS :: Int -> SpMatrix Double -> SpVector Double -> CGS-runNCGS n aa b = map _x $ runAppendN' (cgsStep aa x0) n cgsInit where-  x0 = mkSpVectorD nd $ replicate nd 0.1-  nd = dim r0-  r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution-  p0 = r0-  u0 = r0-  cgsInit = CGS x0 r0 p0 u0  -- {-  example 0 : 2x2 linear system@@ -193,17 +179,7 @@   --- `ndim` iterations -solveRandomN ndim nsp niter = do-  aa0 <- randSpMat ndim (nsp ^ 2)-  let aa = aa0 ^+^ eye ndim-  xtrue <- randSpVec ndim nsp-  let b = aa #> xtrue-      xhatB = head $ runNBiC niter aa b-      xhatC = head $ runNCGS niter aa b-  printDenseSM aa    -  return (normSq (xhatB ^-^ xtrue), normSq (xhatC ^-^ xtrue))  -- @@ -270,6 +246,14 @@   +{- LU -}++checkLu a = lup == a where+  (l, u) = lu a+  lup = l ## u+++ {- eigenvalues -}  @@ -332,3 +316,40 @@ --  tm4 = sparsifySM $ fromListDenseSM 4 [1,0,0,0,2,5,0,10,3,6,8,11,4,7,9,12]++++++++-- -- run N iterations ++-- -- runNBiC :: Int -> SpMatrix Double -> SpVector Double -> BICGSTAB+-- runNBiC n aa b = map _xBicgstab $ runAppendN' (bicgstabStep aa x0) n bicgsInit where+--    x0 = mkSpVectorD nd $ replicate nd 0.9+--    nd = dim r0+--    r0 = b ^-^ (aa #> x0)    +--    p0 = r0+--    bicgsInit = BICGSTAB x0 r0 p0++-- -- runNCGS :: Int -> SpMatrix Double -> SpVector Double -> CGS+-- runNCGS n aa b = map _x $ runAppendN' (cgsStep aa x0) n cgsInit where+--   x0 = mkSpVectorD nd $ replicate nd 0.1+--   nd = dim r0+--   r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution+--   p0 = r0+--   u0 = r0+--   cgsInit = CGS x0 r0 p0 u0++++-- solveRandomN ndim nsp niter = do+--   aa0 <- randSpMat ndim (nsp ^ 2)+--   let aa = aa0 ^+^ eye ndim+--   xtrue <- randSpVec ndim nsp+--   let b = aa #> xtrue+--       xhatB = head $ runNBiC niter aa b+--       xhatC = head $ runNCGS niter aa b+--   -- printDenseSM aa    +--   return (normSq (xhatB ^-^ xtrue), normSq (xhatC ^-^ xtrue))