sparse-linear-algebra 0.2.0.8 → 0.2.0.9

raw patch · 9 files changed

+316/−233 lines, 9 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• README.md +1/−1
• sparse-linear-algebra.cabal +2/−2
• src/Data/Sparse/Common.hs +1/−2
• src/Data/Sparse/IntMap2/IntMap2.hs +187/−0
• src/Data/Sparse/SpMatrix.hs +20/−1
• src/Data/Sparse/SpVector.hs +1/−1
• src/Numeric/LinearAlgebra/Sparse.hs +48/−3
• src/Numeric/LinearAlgebra/Sparse/IntMap.hs +0/−184
• test/LibSpec.hs +56/−39

@@ -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 )

@@ -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

@@ -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 

@@ -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++++

@@ -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

@@ -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 

@@ -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++ 

@@ -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----

@@ -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-}