sparse-linear-algebra 0.2.1.1 → 0.2.2.0

raw patch · 10 files changed

+491/−339 lines, 10 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• README.md +38/−9
• sparse-linear-algebra.cabal +14/−3
• src/Data/Sparse/Common.hs +1/−1
• src/Data/Sparse/IntMap2/IntMap2.hs +0/−194
• src/Data/Sparse/Internal/CSR.hs +34/−0
• src/Data/Sparse/Internal/IntMap2.hs +194/−0
• src/Data/Sparse/SpMatrix.hs +11/−3
• src/Data/Sparse/SpVector.hs +17/−8
• src/Numeric/LinearAlgebra/Sparse.hs +115/−80
• test/LibSpec.hs +67/−41

@@ -8,16 +8,18 @@  Contents : -* Iterative linear solvers+* Iterative linear solvers (`linSolve`) -    * BiConjugate Gradient (`bcg`)+    * Generalized Minimal Residual (GMRES) (non-Hermitian systems)  -    * Conjugate Gradient Squared (`cgs`)+    * BiConjugate Gradient (BCG) -    * BiConjugate Gradient Stabilized (`bicgstab`) (non-Hermitian systems)+    * Conjugate Gradient Squared (CGS) -    * Transpose-Free Quasi-Minimal Residual (`tfqmr`)+    * BiConjugate Gradient Stabilized (BiCGSTAB) (non-Hermitian systems) +    * Transpose-Free Quasi-Minimal Residual (TFQMR)+ * Direct linear solvers      * LU-based (`luSolve`)@@ -64,7 +66,11 @@     fromListSV :: Int -> [(Int, a)] -> SpVector a  can be used to create sparse vectors.+Alternatively, the user can copy the contents of a list to a (dense) SpVector using +    fromListDenseSV :: Int -> [a] -> SpVector a++ ### Displaying sparse data  Both sparse vectors and matrices can be pretty-printed using `prd`:@@ -80,7 +86,7 @@  ### Matrix operations -Matrix factorizations are available as `lu` and `qr` respectively, and are straightforward to verify by using the matrix product `##`  :+There are a few common matrix factorizations available; in the following example we compute the LU factorization of a matrix and verify it with the matrix-matrix product `##`  :      > (l, u) = lu amat     > prd $ l ## u@@ -100,12 +106,33 @@     [4.0,3.0,2.0]     [0.0,0.0,5.0] +A matrix is transposed using `transposeSM`. +Sometimes we need to compute matrix-matrix transpose products, which is why the library offers the infix operators `#^#` (M^T N) and `##^` (M N^T):++    > amat' = amat #^# amat+    > prd amat'+    ( 3 rows, 3 columns ) , 9 NZ ( sparsity 1.0 )++    [20.0,12.0,8.0]+    [12.0,9.0,6.0]+    [8.0,6.0,29.0]++    > l = chol amat'+    > prd $ l ##^ l+    ( 3 rows, 3 columns ) , 9 NZ ( sparsity 1.0 )++    [20.000000000000004,12.0,8.0]+    [12.0,9.0,10.8]+    [8.0,10.8,29.0]++In the above example we have also shown the Cholesky decomposition (M = L L^T where L is a lower-triangular matrix), which is only possible for symmetric positive-definite matrices.+ ### Linear systems -Large sparse linear systems are best solved with iterative methods. `sparse-linear-algebra` provides a selection of these via the `linSolve` function, or alternatively `<\>` (which uses BiCGSTAB as default) :+Large sparse linear systems are best solved with iterative methods. `sparse-linear-algebra` provides a selection of these via the `linSolve` function, or alternatively `<\>` (which uses GMRES as default solver method) : -    > b = fromListSV 3 [(0,3),(1,2),(2,5)]+    > b = fromListDenseSV 3 [3,2,5]     > x = amat <\> b     > prd x     ( 3 elements ) ,  3 NZ ( sparsity 1.0 )@@ -140,7 +167,7 @@  * define appropriate data structures, and how they relate to those properties (sparse vectors and matrices, defined internally via `Data.IntMap`, are made instances of the VectorSpace and Additive classes respectively). This allows to decouple the algorithms from the actual implementation of the backend, -* implement the algorithms, following 1:1 the textbook [1] +* implement the algorithms, following 1:1 the textbook [1, 2]    ## License@@ -157,3 +184,5 @@ ## References  [1] : Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., 2000++[2] : L. N. Trefethen, D. Bau, Numerical Linear Algebra, SIAM, 1997

@@ -1,7 +1,17 @@ name:                sparse-linear-algebra-version:             0.2.1.1+version:             0.2.2.0 synopsis:            Numerical computation in native Haskell -description:         Currently the library provides iterative linear solvers, matrix decompositions, eigenvalue computations and related utilities. Please see README.md for details+description:+  /Overview/+  .+  The @sparse-linear-algebra@ library provides iterative linear solvers, matrix decompositions, eigenvalue computations and related utilities. The user interface is provided by the top-level module "Numeric.LinearAlgebra.Sparse":+  .+  @+  import           Numeric.LinearAlgebra.Sparse+  @+  .+  Please refer to the README file for usage examples.+  . homepage:            https://github.com/ocramz/sparse-linear-algebra license:             GPL-3 license-file:        LICENSE@@ -20,10 +30,11 @@   hs-source-dirs:      src   exposed-modules:     Numeric.LinearAlgebra.Sparse                        Numeric.LinearAlgebra.Class-                       Data.Sparse.IntMap2.IntMap2                        Data.Sparse.SpVector                        Data.Sparse.SpMatrix                        Data.Sparse.Common+  other-modules:       Data.Sparse.Internal.IntMap2+                       Data.Sparse.Internal.CSR                        Data.Sparse.Utils                        Data.Sparse.Types                        Numeric.Eps

@@ -24,7 +24,7 @@  import Data.Sparse.Utils as X import Data.Sparse.Types as X-import Data.Sparse.IntMap2.IntMap2 -- as X+import Data.Sparse.Internal.IntMap2 -- as X import Data.Sparse.SpMatrix as X import Data.Sparse.SpVector as X 

@@ -1,194 +0,0 @@-module Data.Sparse.IntMap2.IntMap2 where--import Numeric.LinearAlgebra.Class--import qualified Data.IntMap.Strict as IM-import Data.Sparse.Types--import Data.Maybe----instance Set IM.IntMap where-  liftU2 = IM.unionWith-  {-# INLINE liftU2 #-}-  liftI2 = IM.intersectionWith-  {-# INLINE liftI2 #-}--instance Additive IM.IntMap where-  zero = IM.empty-  {-# INLINE zero #-}-  (^+^) = liftU2 (+)-  {-# INLINE (^+^) #-}---instance VectorSpace IM.IntMap where-  n .* im = IM.map (* n) im-  -instance Hilbert IM.IntMap where-   a `dot` b = sum $ liftI2 (*) a b-              --instance Normed IM.IntMap where-  norm p v | p==1 = norm1 v-           | p==2 = norm2 v-           | otherwise = normP p v------ set-like brackets---------- * Insertion---- | Insert an element-insertIM2 ::-  IM.Key -> IM.Key -> a -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)-insertIM2 i j x imm = IM.insert i ro imm where-  ro = maybe (IM.singleton j x) (IM.insert j x) (IM.lookup i imm)-{-# inline insertIM2 #-}  ---- * Lookup---- |Lookup a key-lookupIM2 ::-  IM.Key -> IM.Key -> IM.IntMap (IM.IntMap a) -> Maybe a-lookupIM2 i j imm = IM.lookup i imm >>= IM.lookup j-{-# inline lookupIM2 #-}  ---- | Lookup with default 0-lookupWD_IM :: Num a => IM.IntMap (IM.IntMap a) -> (IxRow, IxCol) -> a-lookupWD_IM im (i,j) = fromMaybe 0 (IM.lookup i im >>= IM.lookup j)------- |Populate an IM2 from a list of (row index, column index, value)  -fromListIM2 ::-  Foldable t =>-     t (IM.Key, IM.Key, a) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)-fromListIM2 iix sm = foldl ins sm iix where-  ins t (i,j,x) = insertIM2 i j x t----- * Folding---- |Indexed left fold over an IM2, with general accumulator-ifoldlIM2' :: (IM.Key -> IM.Key -> a -> b -> b) -> b -> IM.IntMap (IM.IntMap a) -> b-ifoldlIM2' f empty mm = IM.foldlWithKey' accRow empty mm where-  accRow acc i r = IM.foldlWithKey' (accElem i) acc r-  accElem i acc j x = f i j x acc-{-# inline ifoldlIM2' #-}---- |Indexed left fold over an IM2-ifoldlIM2 ::-  (IM.Key -> IM.Key -> t -> IM.IntMap a -> IM.IntMap a) ->-  IM.IntMap (IM.IntMap t) ->  -  IM.IntMap a-ifoldlIM2 f m         = IM.foldlWithKey' accRow IM.empty m where-  accRow    acc i row = IM.foldlWithKey' (accElem i) acc row-  accElem i acc j x   = f i j x acc-{-# inline ifoldlIM2 #-}  ---- |Left fold over an IM2, with general accumulator-foldlIM2 :: (a -> b -> b) -> b -> IM.IntMap (IM.IntMap a) -> b-foldlIM2 f empty mm = IM.foldl accRow empty mm where-  accRow acc r = IM.foldl accElem acc r-  accElem acc x = f x acc-{-# inline foldlIM2 #-}----- | Inner indices become outer ones and vice versa. No loss of information because both inner and outer IntMaps are nubbed.-transposeIM2 :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)-transposeIM2 = ifoldlIM2 (flip insertIM2)-{-# inline transposeIM2 #-}---- specialized folds---- -- extract diagonal elements--- extractDiagonalIM2 :: IM.IntMap (IM.IntMap a) -> [a]--- extractDiagonalIM2 = ifoldlIM2' (\i j x xs -> if i==j then x : xs else xs) []------- * Filtering---- |Map over outer IM and filter all inner IM's-ifilterIM2 ::-  (IM.Key -> IM.Key -> a -> Bool) ->-  IM.IntMap (IM.IntMap a) ->-  IM.IntMap (IM.IntMap a)-ifilterIM2 f  =-  IM.mapWithKey (\irow row -> IM.filterWithKey (f irow) row) -{-# inline ifilterIM2 #-}---- |Specialized filtering : keep only sub-diagonal elements-filterSubdiag :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)-filterSubdiag = ifilterIM2 (\i j _ -> i>j)--countSubdiagonalNZ :: IM.IntMap (IM.IntMap a) -> Int-countSubdiagonalNZ im =-  IM.size $ IM.filter (not . IM.null) (filterSubdiag im)---- |List of (row, col) indices of (nonzero) subdiagonal elements-subdiagIndices :: IM.IntMap (IM.IntMap a) -> [(IM.Key, IM.Key)]-subdiagIndices im = concatMap rpairs $ IM.toList (IM.map IM.keys im') where-  im' = filterSubdiag im--rpairs :: (a, [b]) -> [(a, b)]-rpairs (i, jj@(_:_)) = zip (replicate (length jj) i) jj-rpairs (_, []) = []---- -- list of (row, col) indices of elements that satisfy a criterion--- indicesThatIM2 ::---   (IM.Key -> IM.IntMap a -> Bool) -> IM.IntMap (IM.IntMap a) -> [(IM.Key, IM.Key)]--- indicesThatIM2 f im = concatMap rpairs $ IM.toList (IM.map IM.keys im') where---   im' = IM.filterWithKey f im--  ----- * Mapping---- |Map over IM2-mapIM2 :: (a -> b) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap b)-mapIM2 = IM.map . IM.map   -- imapIM2 (\_ _ x -> f x)-------- |Indexed map over IM2-imapIM2 ::-  (IM.Key -> IM.Key -> a -> b) ->-  IM.IntMap (IM.IntMap a) ->-  IM.IntMap (IM.IntMap b)-imapIM2 f im = IM.mapWithKey ff im where-  ff j x = IM.mapWithKey (`f` j) x------ |Mapping keys-mapKeysIM2 ::-  (IM.Key -> IM.Key) -> (IM.Key -> IM.Key) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)-mapKeysIM2 fi fj im = IM.map adjCols adjRows where-  adjRows = IM.mapKeys fi im-  adjCols = IM.mapKeys fj ------- map over a single `column`--mapColumnIM2 :: (b -> b) -> IM.IntMap (IM.IntMap b) -> Int -> IM.IntMap (IM.IntMap b)-mapColumnIM2 f im jj = imapIM2 (\i j x -> if j == jj then f x else x) im----

@@ -0,0 +1,34 @@+module Data.Sparse.Internal.CSR where++import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as VU+import qualified Data.Vector.Unboxed.Mutable as VM++{-| Compressed Row Storage specification :++   http://netlib.org/utk/people/JackDongarra/etemplates/node373.html++   The compressed row storage (CRS) format puts the subsequent nonzeros of the matrix+   rows in contiguous memory locations. Assuming we have a nonsymmetric sparse matrix+   $A$, we create three vectors: one for floating point numbers (val) and the other+   two for integers (col_ind, row_ptr).++   The val vector stores the values of the nonzero elements of the matrix $A$ as+   they are traversed in a row-wise fashion.+ +   The col_ind vector stores the column indexes of the elements in the val vector,+   that is, if val(k)=a_{i,j}, then  col_ind(k)=j$.++   The row_ptr vector stores the locations in the val vector that start a row;+   that is, if  val(k)=a_{i,j}, then row_ptr(i) <= k < row_ptr(i+1)++-}++data CsrMatrix a =+  CsrMatrix { csrVal :: VU.Vector a,+              csrColInd :: VU.Vector Int,+              csrRowPtr :: VU.Vector Int,+              csrNnz :: {-# UNPACK #-} !Int,+              csrNrows :: {-# UNPACK #-} !Int,+              csrNcols :: {-# UNPACK #-} !Int } deriving Eq+

@@ -0,0 +1,194 @@+module Data.Sparse.Internal.IntMap2 where++import Numeric.LinearAlgebra.Class++import qualified Data.IntMap.Strict as IM+import Data.Sparse.Types++import Data.Maybe++++instance Set IM.IntMap where+  liftU2 = IM.unionWith+  {-# INLINE liftU2 #-}+  liftI2 = IM.intersectionWith+  {-# INLINE liftI2 #-}++instance Additive IM.IntMap where+  zero = IM.empty+  {-# INLINE zero #-}+  (^+^) = liftU2 (+)+  {-# INLINE (^+^) #-}+++instance VectorSpace IM.IntMap where+  n .* im = IM.map (* n) im+  +instance Hilbert IM.IntMap where+   a `dot` b = sum $ liftI2 (*) a b+              ++instance Normed IM.IntMap where+  norm p v | p==1 = norm1 v+           | p==2 = norm2 v+           | otherwise = normP p v++++-- set-like brackets++++++++-- * Insertion++-- | Insert an element+insertIM2 ::+  IM.Key -> IM.Key -> a -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+insertIM2 i j x imm = IM.insert i ro imm where+  ro = maybe (IM.singleton j x) (IM.insert j x) (IM.lookup i imm)+{-# inline insertIM2 #-}  ++-- * Lookup++-- |Lookup a key+lookupIM2 ::+  IM.Key -> IM.Key -> IM.IntMap (IM.IntMap a) -> Maybe a+lookupIM2 i j imm = IM.lookup i imm >>= IM.lookup j+{-# inline lookupIM2 #-}  ++-- | Lookup with default 0+lookupWD_IM :: Num a => IM.IntMap (IM.IntMap a) -> (IxRow, IxCol) -> a+lookupWD_IM im (i,j) = fromMaybe 0 (IM.lookup i im >>= IM.lookup j)+++++-- |Populate an IM2 from a list of (row index, column index, value)  +fromListIM2 ::+  Foldable t =>+     t (IM.Key, IM.Key, a) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+fromListIM2 iix sm = foldl ins sm iix where+  ins t (i,j,x) = insertIM2 i j x t+++-- * Folding++-- |Indexed left fold over an IM2, with general accumulator+ifoldlIM2' :: (IM.Key -> IM.Key -> a -> b -> b) -> b -> IM.IntMap (IM.IntMap a) -> b+ifoldlIM2' f empty mm = IM.foldlWithKey' accRow empty mm where+  accRow acc i r = IM.foldlWithKey' (accElem i) acc r+  accElem i acc j x = f i j x acc+{-# inline ifoldlIM2' #-}++-- |Indexed left fold over an IM2+ifoldlIM2 ::+  (IM.Key -> IM.Key -> t -> IM.IntMap a -> IM.IntMap a) ->+  IM.IntMap (IM.IntMap t) ->  +  IM.IntMap a+ifoldlIM2 f m         = IM.foldlWithKey' accRow IM.empty m where+  accRow    acc i row = IM.foldlWithKey' (accElem i) acc row+  accElem i acc j x   = f i j x acc+{-# inline ifoldlIM2 #-}  ++-- |Left fold over an IM2, with general accumulator+foldlIM2 :: (a -> b -> b) -> b -> IM.IntMap (IM.IntMap a) -> b+foldlIM2 f empty mm = IM.foldl accRow empty mm where+  accRow acc r = IM.foldl accElem acc r+  accElem acc x = f x acc+{-# inline foldlIM2 #-}+++-- | Inner indices become outer ones and vice versa. No loss of information because both inner and outer IntMaps are nubbed.+transposeIM2 :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+transposeIM2 = ifoldlIM2 (flip insertIM2)+{-# inline transposeIM2 #-}++-- specialized folds++-- -- extract diagonal elements+-- extractDiagonalIM2 :: IM.IntMap (IM.IntMap a) -> [a]+-- extractDiagonalIM2 = ifoldlIM2' (\i j x xs -> if i==j then x : xs else xs) []+++++-- * Filtering++-- |Map over outer IM and filter all inner IM's+ifilterIM2 ::+  (IM.Key -> IM.Key -> a -> Bool) ->+  IM.IntMap (IM.IntMap a) ->+  IM.IntMap (IM.IntMap a)+ifilterIM2 f  =+  IM.mapWithKey (\irow row -> IM.filterWithKey (f irow) row) +{-# inline ifilterIM2 #-}++-- |Specialized filtering : keep only sub-diagonal elements+filterSubdiag :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+filterSubdiag = ifilterIM2 (\i j _ -> i>j)++countSubdiagonalNZ :: IM.IntMap (IM.IntMap a) -> Int+countSubdiagonalNZ im =+  IM.size $ IM.filter (not . IM.null) (filterSubdiag im)++-- |List of (row, col) indices of (nonzero) subdiagonal elements+subdiagIndices :: IM.IntMap (IM.IntMap a) -> [(IM.Key, IM.Key)]+subdiagIndices im = concatMap rpairs $ IM.toList (IM.map IM.keys im') where+  im' = filterSubdiag im++rpairs :: (a, [b]) -> [(a, b)]+rpairs (i, jj@(_:_)) = zip (replicate (length jj) i) jj+rpairs (_, []) = []++-- -- list of (row, col) indices of elements that satisfy a criterion+-- indicesThatIM2 ::+--   (IM.Key -> IM.IntMap a -> Bool) -> IM.IntMap (IM.IntMap a) -> [(IM.Key, IM.Key)]+-- indicesThatIM2 f im = concatMap rpairs $ IM.toList (IM.map IM.keys im') where+--   im' = IM.filterWithKey f im++  +++-- * Mapping++-- |Map over IM2+mapIM2 :: (a -> b) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap b)+mapIM2 = IM.map . IM.map   -- imapIM2 (\_ _ x -> f x)++++++-- |Indexed map over IM2+imapIM2 ::+  (IM.Key -> IM.Key -> a -> b) ->+  IM.IntMap (IM.IntMap a) ->+  IM.IntMap (IM.IntMap b)+imapIM2 f im = IM.mapWithKey ff im where+  ff j x = IM.mapWithKey (`f` j) x++++-- |Mapping keys+mapKeysIM2 ::+  (IM.Key -> IM.Key) -> (IM.Key -> IM.Key) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+mapKeysIM2 fi fj im = IM.map adjCols adjRows where+  adjRows = IM.mapKeys fi im+  adjCols = IM.mapKeys fj +++++-- map over a single `column`++mapColumnIM2 :: (b -> b) -> IM.IntMap (IM.IntMap b) -> Int -> IM.IntMap (IM.IntMap b)+mapColumnIM2 f im jj = imapIM2 (\i j x -> if j == jj then f x else x) im++++

@@ -16,7 +16,7 @@ import Numeric.Eps import Numeric.LinearAlgebra.Class -import Data.Sparse.IntMap2.IntMap2+import Data.Sparse.Internal.IntMap2  import qualified Data.IntMap as IM @@ -58,12 +58,13 @@ instance Set SpMatrix where   liftU2 f2 (SM n1 x1) (SM n2 x2) = SM (maxTup n1 n2) ((liftU2.liftU2) f2 x1 x2)   liftI2 f2 (SM n1 x1) (SM n2 x2) = SM (minTup n1 n2) ((liftI2.liftI2) f2 x1 x2)-  ++-- | 'SpMatrix'es form a ring, in that they can be added and possess a zero element   instance Additive SpMatrix where   zero = SM (0,0) IM.empty   (^+^) = liftU2 (+) -+-- | 'SpMatrix'es are maps between finite-dimensional spaces instance FiniteDim SpMatrix where   type FDSize SpMatrix = (Rows, Cols)   dim = smDim@@ -75,6 +76,7 @@ instance Sparse SpMatrix a where   spy = spySM +-- | 'SpMatrix'es are sparse containers too, i.e. any specific component may be missing (so it is assumed to be 0) instance Num a => SpContainer SpMatrix a where   type ScIx SpMatrix = (Rows, Cols)   scInsert (i,j) = insertSpMatrix i j@@ -288,6 +290,12 @@ extractSubmatrix :: SpMatrix a -> (IxRow, IxRow) -> (IxCol, IxCol) -> SpMatrix a extractSubmatrix = extractSubmatrixSM id id ++takeRows :: IxRow -> SpMatrix a -> SpMatrix a+takeRows n mm = extractSubmatrix mm (0, n-1) (0, ncols mm - 1)++takeCols :: IxCol -> SpMatrix a -> SpMatrix a+takeCols n mm = extractSubmatrix mm (0, nrows mm - 1) (0, n - 1)   

@@ -12,7 +12,7 @@  import Data.Sparse.Utils import Data.Sparse.Types-import Data.Sparse.IntMap2.IntMap2+import Data.Sparse.Internal.IntMap2  import Numeric.Eps import Numeric.LinearAlgebra.Class@@ -58,11 +58,11 @@   (^+^) = liftU2 (+)  -                      +-- | 'SpVector's form a vector space because they can be multiplied by a scalar instance VectorSpace SpVector where   n .* v = scale n v -+-- | 'SpVector's are finite-dimensional vectors instance FiniteDim SpVector where   type FDSize SpVector = Int   dim = svDim  @@ -75,7 +75,7 @@   spy = spySV  -+-- | 'SpVector's are sparse containers too, i.e. any specific component may be missing (so it is assumed to be 0) instance Num a => SpContainer SpVector a where   type ScIx SpVector = Int   scInsert = insertSpVector@@ -84,14 +84,14 @@   -+-- | 'SpVector's form a Hilbert space, in that we can define an inner product over them instance Hilbert SpVector where   a `dot` b | dim a == dim b = dot (dat a) (dat b)             | otherwise =                      error $ "dot : sizes must coincide, instead we got " ++                            show (dim a, dim b) -+-- | Since 'SpVector's form a Hilbert space, we can define a norm for them  instance Normed SpVector where   norm p (SV _ v) = norm p v @@ -104,17 +104,22 @@  -- ** Creation --- | empty sparse vector (length n, no entries)+-- | Empty sparse vector (length n, no entries) zeroSV :: Int -> SpVector a zeroSV n = SV n IM.empty  --- | singleton sparse vector (length 1)+-- | Singleton sparse vector (length 1) singletonSV :: a -> SpVector a singletonSV x = SV 1 (IM.singleton 0 x)  +-- | Canonical basis vector in R^n+ei :: Num a => Int -> IM.Key -> SpVector a+ei n i = SV n (IM.insert (i - 1) 1 IM.empty) ++ -- | create a sparse vector from an association list while discarding all zero entries mkSpVector :: Epsilon a => Int -> IM.IntMap a -> SpVector a mkSpVector d im = SV d $ IM.filterWithKey (\k v -> isNz v && inBounds0 d k) im@@ -203,6 +208,10 @@ fromListSV :: Int -> [(Int, a)] -> SpVector a fromListSV d iix = SV d (IM.fromList (filter (inBounds0 d . fst) iix )) +-- fromListSV' d iix = SV d (F.foldl' insf IM.empty iix') where+--   insf im (i, x) = IM.insert i x im+--   iix' = filter (inBounds0 d . fst) iix   -- filtering forces list as only instance+  -- ** toList toListSV :: SpVector a -> [(IM.Key, a)] toListSV sv = IM.toList (dat sv)

@@ -3,7 +3,8 @@ module Numeric.LinearAlgebra.Sparse        (          -- * Matrix factorizations-         qr, lu,+         qr,+         lu,          chol,          -- * Condition number          conditionNumberSM,@@ -12,18 +13,15 @@          -- * Givens' rotation          givens,          -- * Arnoldi iteration-         arnoldi,+         arnoldi,           -- * Eigensolvers          eigsQR, eigRayleigh,          -- * Linear solvers-         linSolve, LinSolveMethod, (<\>),+         -- ** Iterative methods+         linSolve, LinSolveMethod(..), (<\>),+         pinv,          -- ** Direct methods          luSolve,-         -- ** Iterative methods-         cgne, tfqmr, bicgstab, cgs, bcg,-         _xCgne, _xTfq, _xBicgstab, _x, _xBcg,-         cgsStep, bicgstabStep,-         CGNE, TFQMR, BICGSTAB, CGS, BCG,          -- * Preconditioners          ilu0, mSsor,          -- * Matrix partitioning@@ -49,7 +47,7 @@ import Control.Monad.Primitive import Control.Monad (mapM_, forM_, replicateM) import Control.Monad.State.Strict-import Control.Monad.Writer+-- import Control.Monad.Writer -- import Control.Monad.Trans -- import Control.Monad.Trans.State (runStateT) -- import Control.Monad.Trans.Writer (runWriterT)@@ -85,7 +83,7 @@  -- |uses the R matrix from the QR factorization conditionNumberSM :: (Epsilon a, RealFloat a) => SpMatrix a -> a-conditionNumberSM m | isInfinite kappa = error "Infinite condition number : rank-deficient system"+conditionNumberSM m | nearZero lmin = error "Infinite condition number : rank-deficient system"                     | otherwise = kappa where   kappa = lmax / lmin   (_, r) = qr m@@ -162,8 +160,8 @@ {-# inline givens #-} givens :: (Floating a, Epsilon a, Ord a) => SpMatrix a -> IxRow -> IxCol -> SpMatrix a givens mm i j -  | isValidIxSM mm (i,j) && isSquareSM mm =-       sparsifySM $ fromListSM' [(i,i,c),(j,j,c),(j,i,-s),(i,j,s)] (eye (nrows mm))+  | isValidIxSM mm (i,j) && nrows mm >= ncols mm =+       fromListSM' [(i,i,c),(j,j,c),(j,i,-s),(i,j,s)] (eye (nrows mm))   | otherwise = error "givens : indices out of bounds"         where     (c, s, _) = givensCoef a b@@ -188,40 +186,33 @@   --- * QR decomposition  --- | Given a matrix A, returns a pair of matrices (Q, R) such that Q R = A, Q is orthogonal and R is upper triangular. Applies Givens rotation iteratively to zero out sub-diagonal elements-qr :: (Epsilon a, Floating a, Real a) => SpMatrix a -> (SpMatrix a, SpMatrix a)-qr mm = (transposeSM qmatt, rmat)  where-  qmatt = F.foldl' (#~#) ee $ gmats mm -- Q^T = (G_n * G_n-1 ... * G_1)-  rmat = qmatt #~# mm                  -- R = Q^T A-  ee = eye (nrows mm)-      --- | Givens matrices in order [G1, G2, .. , G_N ]-gmats :: (Epsilon a, Real a, Floating a) => SpMatrix a -> [SpMatrix a]-gmats mm = gm mm (subdiagIndicesSM mm) where- gm m ((i,j):is) = let g = givens m i j-                   in g : gm (g #~# m) is- gm _ [] = []+-- * QR decomposition  +-- | Given a matrix A, returns a pair of matrices (Q, R) such that Q R = A, where Q is orthogonal and R is upper triangular. Applies Givens rotation iteratively to zero out sub-diagonal elements.+qr :: (Epsilon a, Ord a, Floating a) => SpMatrix a -> (SpMatrix a, SpMatrix a)+qr mm = (transposeSM qt, r) where+  (qt, r, _) = execState (modifyUntil qf stepf) gminit+  qf (_, _, iis) = null iis+  stepf (qmatt, m, iis) = (qmatt', m', tail iis) where+    (i, j) = head iis+    g = givens m i j+    qmatt' = g #~# qmatt  -- update Q'+    m' = g #~# m          -- update R+  gminit = (eye (nrows mm), mm, subdiagIndicesSM mm) +    +    --- -- | QR algorithm, state transformer version--- gmatST0 (m, (i,j):is) = (m', is) where    -- WRONG, possible access to []---   g = givens m i j                        ---   m' = g #~# m--- gmatST0 (m, []) = (eye (nrows m), []) --- gmatST m = gmatST0 (m, subdiagIndicesSM m)     - -- * Eigenvalue algorithms  -- ** QR algorithm@@ -353,7 +344,7 @@ -- ** Doolittle algorithm {- Doolittle algorithm for factoring A' = P A, where P is a permutation matrix such that A' has a nonzero as its (0, 0) entry -} --- | Given a matrix A, returns a pair of matrices (L, U) such that L U = A+-- | Given a matrix A, returns a pair of matrices (L, U) where L is lower triangular and U is upper triangular such that L U = A lu :: (Epsilon a, Fractional a, Real a) => SpMatrix a -> (SpMatrix a, SpMatrix a) lu aa = (lf, ufin) where   (ixf, lf, uf) = execState (modifyUntil q luUpd) luInit@@ -429,61 +420,48 @@   --- * Incomplete LU --- | used for Incomplete LU : remove entries in `m` corresponding to zero entries in `m2`  -ilu0 :: (Epsilon a, Real a, Fractional a) => SpMatrix a -> (SpMatrix a, SpMatrix a)-ilu0 aa = (lh, uh) where-  (l, u) = lu aa-  lh = sparsifyLU l aa-  uh = sparsifyLU u aa-  sparsifyLU m m2 = ifilterSM f m where-    f i j _ = isJust (lookupSM m2 i j)  ---- -- * Arnoldi iteration --- | Given a matrix A and a positive integer `n`, this procedure finds the basis of an order `n` Krylov subspace (as the columns of matrix Q), along with an upper Hessenberg matrix H, such that A = Q^T H Q.+-- | Given a matrix A, a vector b and a positive integer `n`, this procedure finds the basis of an order `n` Krylov subspace (as the columns of matrix Q), along with an upper Hessenberg matrix H, such that A = Q^T H Q. -- At the i`th iteration, it finds (i + 1) coefficients (the i`th column of the Hessenberg matrix H) and the (i + 1)`th Krylov vector.+ arnoldi ::-  (Floating a, Eq a) => SpMatrix a -> Int -> (SpMatrix a, SpMatrix a)-arnoldi aa kn = (fromCols qvfin, hhfin) where-  (qvfin, hhfin, _) = execState (modifyUntil tf arnoldiStep) arnInit -  tf (_, _, ii) = ii == kn -- termination criterion+  (Epsilon a, Floating a, Eq a) => SpMatrix a -> SpVector a -> Int -> (SpMatrix a, SpMatrix a)+arnoldi aa b kn = (fromCols qvfin, fromListSM (nmax + 1, nmax) hhfin)+  where+  (qvfin, hhfin, nmax, _) = execState (modifyUntil tf arnoldiStep) arnInit +  tf (_, _, ii, fbreak) = ii == kn || fbreak -- termination criterion   (m, n) = dim aa-  arnInit = (qv1, hh1, 1) where      -      q0 = normalize 2 $ onesSV n -- starting basis vector-      aq0 = aa #> q0-      h11 = q0 `dot` aq0+  arnInit = (qv1, hh1, 1, False) where+      q0 = normalize 2 b   -- starting basis vector+      aq0 = aa #> q0       -- A q0+      h11 = q0 `dot` aq0                 q1nn = (aq0 ^-^ (h11 .* q0))-      hh1 = fromListSM (m + 1, n) [(0, 0, h11), (1, 0, h21)] where        +      hh1 = V.fromList [(0, 0, h11), (1, 0, h21)] where                 h21 = norm 2 q1nn-      q1 = normalize 2 q1nn+      q1 = normalize 2 q1nn       -- q1 `dot` q0 ~ 0       qv1 = V.fromList [q0, q1]-  arnoldiStep (qv, hh, i) = (qv', hh', i + 1)-   where+  arnoldiStep (qv, hh, i, _) = (qv', hh', i + 1, fb') where     qi = V.last qv     aqi = aa #> qi     hhcoli = fmap (`dot` aqi) qv -- H_{1, i}, H_{2, i}, .. , H_{m + 1, i}     zv = zeroSV m     qipnn =       aqi ^-^ (V.foldl' (^+^) zv (V.zipWith (.*) hhcoli qv)) -- unnormalized q_{i+1}-    qipnorm = singletonSV $ norm 2 qipnn      -- normalization factor H_{i+1, i}+    qipnorm = norm 2 qipnn      -- normalization factor H_{i+1, i}     qip = normalize 2 qipnn              -- q_{i + 1}-    hh' = insertCol hh (concatSV (fromVector hhcoli) qipnorm) i -- update H+    hh' = (V.++) hh (indexed2 $ V.snoc hhcoli qipnorm) where -- update H+      indexed2 v = V.zip3 ii jj v+      ii = V.fromList [0 .. n]    -- nth col of upper Hessenberg has `n+1` nz+      jj = V.replicate (n + 1) i  -- `n+1` replicas of `i`     qv' = V.snoc qv qip        -- append q_{i+1} to Krylov basis Q_i---  ---+    fb' | nearZero qipnorm = True  -- breakdown condition+        | otherwise = False   @@ -501,6 +479,25 @@   f = extractSuperDiag aa  +-- -- ** Jacobi preconditioner++-- -- | Returns the reciprocal of the diagonal +-- jacobiPreconditioner :: SpMatrix Double -> SpMatrix Double+-- jacobiPreconditioner = reciprocal . extractDiag+++-- ** Incomplete LU++-- | Used for Incomplete LU : remove entries in `m` corresponding to zero entries in `m2`+ilu0 :: (Epsilon a, Real a, Fractional a) => SpMatrix a -> (SpMatrix a, SpMatrix a)+ilu0 aa = (lh, uh) where+  (l, u) = lu aa+  lh = sparsifyLU l aa+  uh = sparsifyLU u aa+  sparsifyLU m m2 = ifilterSM f m where+    f i j _ = isJust (lookupSM m2 i j)++ -- ** SSOR  -- | `mSsor aa omega` : if `omega = 1` it returns the symmetric Gauss-Seidel preconditioner@@ -524,18 +521,19 @@ luSolve ::   (Fractional a, Eq a, Epsilon a) => SpMatrix a -> SpMatrix a -> SpVector a -> SpVector a luSolve ll uu b-  | isLowerTriSM ll && isUpperTriSM uu = lubwSolve uu (lufwSolve ll b)+  | isLowerTriSM ll && isUpperTriSM uu = triUpperSolve uu (triLowerSolve ll b)   | otherwise = error "luSolve : factors must be triangular matrices"  -lufwSolve ll b = sparsifySV v where+triLowerSolve :: (Epsilon a, Fractional a) => SpMatrix a -> SpVector a -> SpVector a+triLowerSolve ll b = sparsifySV v where   (v, _) = execState (modifyUntil q lStep) lInit where   q (_, i) = i == dim b-  lStep (ww, i) = (wwi, i + 1) where+  lStep (ww, i) = (ww', i + 1) where     lii = ll @@ (i, i)     bi = b @@ i     wi = (bi - r)/lii where       r = extractSubRow ll i (0, i-1) `dot` takeSV i ww-    wwi = insertSpVector i wi ww+    ww' = insertSpVector i wi ww   lInit = (ww0, 1) where     l00 = ll @@ (0, 0)     b0 = b @@ 0@@ -543,15 +541,16 @@     ww0 = insertSpVector 0 w0 $ zeroSV (dim b)    -- | NB in the computation of `xi` we must rebalance the subrow indices because `dropSV` does that too, in order to take the inner product with consistent index pairs-lubwSolve uu w = sparsifySV x where+triUpperSolve :: (Epsilon a, Fractional a) => SpMatrix a -> SpVector a -> SpVector a+triUpperSolve uu w = sparsifySV x where   (x, _) = execState (modifyUntil q uStep) uInit   q (_, i) = i == (- 1)-  uStep (xx, i) = (xxi, i - 1) where+  uStep (xx, i) = (xx', i - 1) where     uii = uu @@ (i, i)     wi = w @@ i     xi = (wi - r) / uii where         r = extractSubRow_RK uu i (i + 1, dim w - 1) `dot` dropSV (i + 1) xx-    xxi = insertSpVector i xi xx+    xx' = insertSpVector i xi xx   uInit = (xx0, i - 1) where     i = dim w - 1     u00 = uu @@ (i, i)@@ -569,12 +568,42 @@  -- ** GMRES +-- | Given a linear system `A x = b` where `A` is an (m x m) real-valued matrix, the GMRES method finds an approximate solution `xhat` such that the Euclidean norm of the residual `A xhat - b` is minimized. `xhat` is spanned by the order-`n` Krylov subspace of (A, b).+-- In this implementation:+-- 1) the Arnoldi iteration is carried out until numerical breakdown (therefore yielding _at_most_ `m` Krylov basis vectors)+-- 2) the resulting Hessenberg matrix is factorized in QR form+-- 3) the Krylov-subspace solution `yhat` is found by backsubstitution+-- 4) the approximate solution in the original space `xhat` is computed using the Krylov basis, `xhat = Q_n yhat`+--+-- Many optimizations are possible, for example interleaving the QR factorization with the Arnoldi process (and employing an updating QR factorization which only requires one Givens' rotation at every update). ++gmres :: (Epsilon a, Ord a, Floating a) => SpMatrix a -> SpVector a -> SpVector a+gmres aa b = qa' #> yhat where+  m = ncols aa+  (qa, ha) = arnoldi aa b m   -- at most m steps of Arnoldi (aa, b)+  -- b' = (transposeSe qa) #> b+  b' = norm2 b .* (ei mp1 1)  -- b rotated back to canonical basis by Q^T+     where mp1 = nrows ha     -- 1 + # Arnoldi iterations+  (qh, rh) = qr ha            -- QR fact.of H+  yhat = triUpperSolve rh' rhs' where+    rhs' = takeSV (dim b' - 1) (transposeSM qh #> b')+    rh' = takeRows (nrows rh - 1) rh+  qa' = takeCols (ncols qa - 1) qa+++++++ -- *** Left-preconditioning     ++ -- ** CGNE  cgneStep :: SpMatrix Double -> CGNE -> CGNE@@ -587,6 +616,10 @@  data CGNE =   CGNE {_xCgne , _rCgne, _pCgne :: SpVector Double} deriving Eq+instance Show CGNE where+    show (CGNE x r p) = "x = " ++ show x ++ "\n" +++                       "r = " ++ show r ++ "\n" +++                       "p = " ++ show p ++ "\n"  cgne :: SpMatrix Double -> SpVector Double -> SpVector Double -> CGNE cgne aa b x0 = execState (untilConverged _xCgne (cgneStep aa)) cgneInit where@@ -645,7 +678,8 @@           _mTfq :: Int,           _tauTfq, _thetaTfq, _etaTfq, _rhoTfq, _alphaTfq :: Double}   deriving Eq-+instance Show TFQMR where+    show (TFQMR x _ _ _ _ _ _ _ _ _ _) = "x = " ++ show x ++ "\n"   -- ** BCG@@ -782,7 +816,7 @@  -- * Linear solver interface -data LinSolveMethod = CGNE_ | TFQMR_ | BCG_ | CGS_ | BICGSTAB_ deriving (Eq, Show) +data LinSolveMethod = GMRES_ | CGNE_ | TFQMR_ | BCG_ | CGS_ | BICGSTAB_ deriving (Eq, Show)   -- -- | Linear solve with _random_ starting vector -- linSolveM ::@@ -805,18 +839,19 @@       solve aa' b' | isDiagonalSM aa' = reciprocal aa' #> b' -- diagonal solve is easy                    | otherwise = solveWith aa' b'        solveWith aa' b' = case method of+        GMRES_ -> gmres aa' b'         CGNE_ -> _xCgne (cgne aa' b' x0)         TFQMR_ -> _xTfq (tfqmr aa' b' x0)         BCG_ -> _xBcg (bcg aa' b' x0)-        CGS_ ->  _xBicgstab (bicgstab aa' b' x0 x0)-        BICGSTAB_ -> _x (cgs aa' b' x0 x0)+        BICGSTAB_ ->  _xBicgstab (bicgstab aa' b' x0 x0)+        CGS_ -> _x (cgs aa' b' x0 x0)       x0 = mkSpVectorD n $ replicate n 0.1        (m, n) = dim aa       nb     = dim b  -- | <\\> : linSolve using the BiCGSTAB method as default (<\>) :: SpMatrix Double -> SpVector Double -> SpVector Double      -(<\>) = linSolve BICGSTAB_ +(<\>) = linSolve GMRES_      

@@ -13,14 +13,13 @@ import Numeric.LinearAlgebra.Sparse -- import Numeric.LinearAlgebra.Class -import Control.Monad (liftM, liftM2)+import Control.Applicative (liftA2)+-- import Control.Monad (liftM, liftM2, replicateM) import Control.Monad.Primitive import Data.Foldable (foldrM)  import Data.Sparse.Common --import Control.Monad (replicateM) import Control.Monad.State.Strict (execState)  import qualified System.Random.MWC as MWC@@ -55,9 +54,9 @@     it "eye : identity matrix" $       infoSM (eye 10) `shouldBe` SMInfo 10 0.1     it "insertCol : insert a column in a SpMatrix" $-      insertCol (eye 3) (fromListDenseSV 3 [2,2,2]) 0 `shouldBe` (fromListSM (3,3) [(0,0,2),(1,0,2),(1,1,1),(2,0,2),(2,2,1)])+      insertCol (eye 3) (fromListDenseSV 3 [2,2,2]) 0 `shouldBe` fromListSM (3,3) [(0,0,2),(1,0,2),(1,1,1),(2,0,2),(2,2,1)]     it "insertRow : insert a row in a SpMatrix" $-      insertRow (eye 3) (fromListDenseSV 3 [2,2,2]) 1 `shouldBe` (fromListSM (3,3) [(0,0,1), (1,0,2), (1,1,2), (1,2,2), (2,2,1)])+      insertRow (eye 3) (fromListDenseSV 3 [2,2,2]) 1 `shouldBe` fromListSM (3,3) [(0,0,1), (1,0,2), (1,1,2), (1,2,2), (2,2,1)]     it "extractCol -> insertCol : identity" $       insertCol (eye 3) (extractCol (eye 3) 1) 1 `shouldBe` eye 3     it "extractRow -> insertRow : identity" $@@ -77,61 +76,76 @@   describe "Numeric.LinearAlgebra.Sparse : Iterative linear solvers" $ do     -- it "TFQMR (2 x 2 dense)" $     --   normSq (_xTfq (tfqmr aa0 b0 x0) ^-^ x0true) <= eps `shouldBe` True+    it "GMRES (2 x 2 dense)" $+      nearZero (normSq (linSolve GMRES_ aa0 b0 ^-^ x0true)) `shouldBe` True+    it "GMRES (3 x 3 sparse, s.p.d.)" $+      nearZero (normSq (linSolve GMRES_ aa2 b2 ^-^ x2)) `shouldBe` True+    it "GMRES (4 x 4 sparse)" $+      nearZero (normSq (linSolve GMRES_ aa1 b1 ^-^ x1)) `shouldBe` True           it "BCG (2 x 2 dense)" $-      nearZero (normSq (_xBcg (bcg aa0 b0 x0) ^-^ x0true)) `shouldBe` True-    it "BiCGSTAB (2 x 2 dense)" $ -      nearZero (normSq (aa0 <\> b0 ^-^ x0true)) `shouldBe` True+      nearZero (normSq (linSolve BCG_ aa0 b0 ^-^ x0true)) `shouldBe` True+    it "BCG (3 x 3 sparse, s.p.d.)" $+      nearZero (normSq (linSolve BCG_ aa2 b2 ^-^ x2)) `shouldBe` True      +    -- it "BiCGSTAB (2 x 2 dense)" $ +    --   nearZero (normSq (linSolve BICGSTAB_ aa0 b0 ^-^ x0true)) `shouldBe` True+    it "BiCGSTAB (3 x 3 sparse, s.p.d.)" $ +      nearZero (normSq (linSolve BICGSTAB_ aa2 b2 ^-^ x2)) `shouldBe` True           it "CGS (2 x 2 dense)" $ -      nearZero (normSq (_x (cgs aa0 b0 x0 x0) ^-^ x0true)) `shouldBe` True-  describe "Numeric.LinearAlgebra.Sparse : Direct linear solvers" $ do-    it "LU (unoptimized) (4 x 4 sparse)" $ +      nearZero (normSq (linSolve CGS_ aa0 b0 ^-^ x0true)) `shouldBe` True+    it "CGS (3 x 3 sparse, s.p.d.)" $ +      nearZero (normSq (linSolve CGS_ aa2 b2 ^-^ x2)) `shouldBe` True      +  describe "Numeric.LinearAlgebra.Sparse : Direct linear solvers" $ +    it "luSolve (4 x 4 sparse)" $        checkLuSolve aa1 b1 `shouldBe` True            describe "Numeric.LinearAlgebra.Sparse : QR decomposition" $ do    -    it "QR (4 x 4 sparse)" $+    it "qr (4 x 4 sparse)" $       checkQr tm4 `shouldBe` True-    it "QR (3 x 3 dense)" $ +    it "qr (3 x 3 dense)" $        checkQr tm2 `shouldBe` True   describe "Numeric.LinearAlgebra.Sparse : LU decomposition" $ do-    it "LU (4 x 4 dense)" $+    it "lu (4 x 4 dense)" $       checkLu tm6 `shouldBe` True-    it "LU (10 x 10 sparse)" $+    it "lu (10 x 10 sparse)" $       checkLu tm7 `shouldBe` True-  describe "Numeric.LinearAlgebra.Sparse : Cholesky decomposition (PSD matrices only)" $ +  describe "Numeric.LinearAlgebra.Sparse : Cholesky decomposition (PSD matrices)" $      it "chol (5 x 5 sparse)" $       checkChol tm7 `shouldBe` True-  describe "Numeric.LinearAlgebra.Sparse : Arnoldi iteration" $ do-    it "Arnoldi iteration (3 x 3 dense)" $-      checkArnoldi aa2 3 `shouldBe` True-    it "Arnoldi iteration (5 x 5 sparse)" $+  describe "Numeric.LinearAlgebra.Sparse : Arnoldi iteration, early breakdown detection" $ do      +    it "arnoldi (4 x 4 dense)" $+      checkArnoldi tm6 4 `shouldBe` True+    it "arnoldi (5 x 5 sparse)" $       checkArnoldi tm7 5 `shouldBe` True      - {- QR-}   checkQr :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool-checkQr a = c1 && c2 where+checkQr a = c1 && c2 && c3 where   (q, r) = qr a   c1 = nearZero $ normFrobenius ((q #~# r) ^-^ a)   c2 = isOrthogonalSM q+  c3 = isUpperTriSM r    {- LU -}  checkLu :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool-checkLu a = lup == a where+checkLu a = c1 && c2 where   (l, u) = lu a-  lup = l #~# u+  c1 = nearZero (normFrobenius ((l #~# u) ^-^ a))+  c2 = isUpperTriSM u && isLowerTriSM l    {- Cholesky -}  checkChol :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool-checkChol a = nearZero $ normFrobenius ((l ##^ l) ^-^ a) where+checkChol a = c1 && c2 where   l = chol a+  c1 = nearZero $ normFrobenius ((l ##^ l) ^-^ a)+  c2 = isLowerTriSM l   {- direct linear solver -}@@ -145,12 +159,13 @@    {- Arnoldi iteration -} checkArnoldi :: (Epsilon a, Floating a, Eq a) => SpMatrix a -> Int -> Bool-checkArnoldi aa kn = nearZero $ normFrobenius $ (aa #~# qvprev) ^-^ (qv #~# hh) where-  (qv, hh) = arnoldi aa kn-  (m, n) = dim qv-  qvprev = extractSubmatrix qv (0, m - 1) (0, n - 2)--+checkArnoldi aa kn = nearZero (normFrobenius $ lhs ^-^ rhs) where+  b = onesSV (nrows aa)+  (q, h) = arnoldi aa b kn+  (m, n) = dim q+  q' = extractSubmatrix q (0, m - 1) (0, n - 2) -- q' = all but one column of q+  rhs = q #~# h+  lhs = aa #~# q'   {-@@ -258,19 +273,20 @@ --   let --     diags = [-bw .. bw - 1] -randDiagMat :: PrimMonad m =>-     Rows -> Double -> Double -> Int -> m (SpMatrix Double)-randDiagMat n mu sig i = do-  x <- randArray n mu sig-  return $ mkSubDiagonal n i x+-- randDiagMat :: PrimMonad m =>+--      Rows -> Double -> Double -> Int -> m (SpMatrix Double)+-- randDiagMat n mu sig i = do+--   x <- randArray n mu sig+--   return $ mkSubDiagonal n i x  -go (m:ms) mat =-  m ^+^ go ms mat-go [] mat = mat+-- go (m:ms) mat =+--   m ^+^ go ms mat+-- go [] mat = mat -  -plusM x y = return $ x ^+^ y+-- plusM ::+--   (Additive f1, Applicative f, Num a) => f (f1 a) -> f (f1 a) -> f (f1 a)  +-- plusM = liftA2 (^+^)   @@ -397,6 +413,11 @@   b = mkSubDiagonal n 0 $ replicate n 2   c = mkSubDiagonal n (-1) $ replicate n (-1) +++++ -- -- run N iterations   -- -- runNBiC :: Int -> SpMatrix Double -> SpVector Double -> BICGSTAB@@ -434,3 +455,8 @@  tm8' :: SpMatrix Double tm8' = fromListSM (2,2) [(0,0,1), (1,0,1), (1,1,1)]++++tm9 :: SpMatrix Double+tm9 = fromListSM (4, 3) [(0,0,pi), (1,1, 3), (2,2,4), (3,2, 1), (3,1, 5)]