sparse-linear-algebra 0.2.9.8 → 0.2.9.9
raw patch · 10 files changed
+30/−275 lines, 10 files
Files
- CHANGELOG.markdown +4/−6
- CONTRIBUTORS.md +5/−0
- README.md +8/−2
- sparse-linear-algebra.cabal +9/−6
- src/Data/Sparse/Common.hs +1/−1
- src/Data/Sparse/Internal/CSB.hs +0/−111
- src/Data/Sparse/Internal/CSC.hs +0/−125
- src/Data/Sparse/Internal/IntM.hs +1/−1
- src/Data/Sparse/Internal/TriMatrix.hs +1/−1
- src/Data/Sparse/Internal/Utils.hs +1/−22
CHANGELOG.markdown view
@@ -1,7 +1,9 @@ + 0.2.9.9+ Moved to IntMap.Strict (Gregory Schwartz)+ Stackage LTS bump to 10.4 (GHC 8.2)+ 0.2.9.7- ------- Improved pretty printer: * Fixed display precision (e.g. 2 decimal digits), fixed column width output for vectors and matrices@@ -9,14 +11,10 @@ * Small and large values (wrt fixed precision) switch to scientific notation 0.2.9.4- -------- Exceptions constructors are exported by Numeric.LinearAlgebra.Sparse 0.2.9.1- -------- * Uses classes from `vector-space` such as AdditiveGroup, VectorSpace and InnerSpace * QuickCheck tests for algebraic properties, such as matrix-vector products and soon more abstract ones e.g. positive semi-definite matrices
+ CONTRIBUTORS.md view
@@ -0,0 +1,5 @@+Original author : Marco Zocca (@ocramz)++Contributors :++Gregory Schwartz (@GregorySchwartz)
README.md view
@@ -3,6 +3,8 @@ Numerical computation in native Haskell [](https://hackage.haskell.org/package/sparse-linear-algebra) [](https://travis-ci.org/ocramz/sparse-linear-algebra)+[](http://stackage.org/lts/package/sparse-linear-algebra)+[](http://stackage.org/nightly/package/sparse-linear-algebra) This library provides common numerical analysis functionality, without requiring any external bindings. It aims to serve as an experimental platform for scientific computation in a purely functional setting. @@ -15,7 +17,7 @@ * The `vector`-based backend is being reworked. -* An `accelerate`-based backend is under development.+* An `accelerate`-based backend is under development [6, 7]. ## Contents@@ -233,4 +235,8 @@ [4] L.N. Trefethen, D. Bau, Numerical Linear Algebra, SIAM, 1997 -[5] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran 77, 2nd ed., 1992 +[5] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran 77, 2nd ed., 1992++[6] M. M. T. Chakravarty, et al., Accelerating Haskell array codes with multicore GPUs - DAMP'11++[7] [`accelerate`](http://hackage.haskell.org/package/accelerate)
sparse-linear-algebra.cabal view
@@ -1,10 +1,10 @@ name: sparse-linear-algebra-version: 0.2.9.8-synopsis: Numerical computation in native Haskell+version: 0.2.9.9+synopsis: Numerical computing in native Haskell 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":+ The @sparse-linear-algebra@ library provides iterative linear solvers, matrix decompositions, eigenvalue algorithms and related utilities. The user interface is provided by the top-level module "Numeric.LinearAlgebra.Sparse": . @ import Numeric.LinearAlgebra.Sparse@@ -22,6 +22,7 @@ build-type: Simple extra-source-files: README.md CHANGELOG.markdown+ CONTRIBUTORS.md data-dir: test/data data-files: e05r0000.mtx e05r0000_rhs1.mtx@@ -38,7 +39,9 @@ library default-language: Haskell2010- ghc-options: -Wall -O2+ ghc-options:+ -- -Wall+ -O2 -Wno-name-shadowing -Wno-unused-top-binds -Wno-missing-signatures -Wno-unused-imports if flag(dump) ghc-options: -ddump-simpl -ddump-stg -ddump-to-file@@ -60,8 +63,8 @@ -- Data.Sparse.Internal.SHVector Data.Sparse.Internal.SList Data.Sparse.Internal.TriMatrix- Data.Sparse.Internal.CSB- Data.Sparse.Internal.CSC + -- Data.Sparse.Internal.CSB+ -- Data.Sparse.Internal.CSC Data.Sparse.Utils Data.Sparse.PPrint Data.Sparse.Types
src/Data/Sparse/Common.hs view
@@ -40,7 +40,7 @@ import Numeric.Eps as X import Numeric.LinearAlgebra.Class as X -import qualified Data.IntMap as IM+import qualified Data.IntMap.Strict as IM import GHC.Exts import Data.Complex
− src/Data/Sparse/Internal/CSB.hs
@@ -1,111 +0,0 @@-{-# language DeriveFunctor #-}-module Data.Sparse.Internal.CSB where--import Control.Applicative--- import Control.Monad.Primitive-import Control.Monad.ST--import Data.Foldable (foldl')--- import Data.List (group, groupBy)--import qualified Data.Vector as V -import qualified Data.Vector.Unboxed as VU-import qualified Data.Vector.Mutable as VM-import qualified Data.Vector.Algorithms.Merge as VA (sortBy)--- import qualified Data.Vector.Generic as VG (convert)--import qualified Data.IntMap as IM--import Control.Monad-import Data.Maybe-import Data.Ord (comparing)---- import Data.Complex--- import Foreign.C.Types (CSChar, CInt, CShort, CLong, CLLong, CIntMax, CFloat, CDouble)--import Control.Concurrent--- import qualified Control.Monad.Par as Par--- import Control.DeepSeq---import Data.Sparse.Utils-import Data.Sparse.Types--- import Data.Sparse.Internal.CSRVector-import Data.Sparse.Internal.Utils--import Numeric.LinearAlgebra.Class---{--Specification of the "Compressed Sparse Block" matrix storage format, as published in [1]:--Let f(i, j) be the bijection from pairs of block indices to integers in the range 0,1,..., n^2/beta^2 − 1 that describes the ordering among blocks.-That is, f(i,j) < f(i',j') if and only if Aij appears before Ai'j' in val.--The row_ind and col_ind arrays store the row and column indices,-respectively, of the elements in the val array. These indices-are relative to the block containing the particular element, not the-entire matrix, and hence they range from 0 to β−1.--If val[k] stores the matrix element aiβ+r,jβ+c, which is located in the rth row-and cth column of the block Aij, then row_ind = r and col_ind = c.--The blk_ptr array stores the index of each block in the val array, which is analogous to the row_ptr array for CSR. If val[k] stores a matrix element falling in block Aij, then blk_ptr[ f(i,j) ] ≤ k < blk_ptr[ f(i,j)+1 ].--[1] A. Buluc, et al., Parallel Sparse Matrix-Vector and Matrix-Transpose-Vector Multiplication Using Compressed Sparse Blocks--}----- | CSB Block--- --- Invariants :--- 1) rowIx, colIx and val have same length--- 2) " have at most (bRows x bCols) NZ----data CsbMatrix a = CSB {- csbNrows, csbNcols, csbBeta :: {-# UNPACK #-} !Int,- csbVal :: V.Vector a,- csbBlkPtr, csbRowIx, csbColIx :: V.Vector Int- } deriving (Eq, Functor)---csbParams :: (Int, Int) -- ^ Matrix size- -> Int -- ^ Block parameter (i.e. block edge length)- -> (Int, Int) -- ^ # of blocks per side-csbParams (m,n) beta = (ceiling m', ceiling n')- where m' = fromIntegral m / fromIntegral beta- n' = fromIntegral n / fromIntegral beta---- | Which block (in row, col format) do matrix indices (i,j) fall in ?-bBin, bCoords :: Integral t => t -> t -> t -> (t, t)-bBin beta i j = (i `div` beta, j `div` beta)---- | Block coordinates-bCoords beta i j = (i `mod` beta, j `mod` beta)---- | Block index (row-major order)-blockIx :: (Int, Int) -> Int -> Int -> Int -> Int-blockIx dims beta i j = bx + by*nbx where- (bx, by) = bBin beta i j- (nbx, _) = csbParams dims beta----type Block i a = [(i,i,a)]--{-# inline consBlockElem #-}-consBlockElem ::- IM.Key -> (i, i, a) -> IM.IntMap (Block i a) -> IM.IntMap (Block i a)-consBlockElem ib x bb = IM.insert ib (x : blocki) bb where- blocki = fromMaybe [] (IM.lookup ib bb)---consBlocks :: Foldable t =>- (Int, Int) -> Int -> t (Int, Int, a) -> IM.IntMap (Block Int a)-consBlocks dims beta ijx = foldl' ins IM.empty ijx where- ins accb (i,j,x) = consBlockElem ib (i,j,x) accb where- ib = blockIx dims beta i j - ---
− src/Data/Sparse/Internal/CSC.hs
@@ -1,125 +0,0 @@-{-# language TypeFamilies, FlexibleInstances, MultiParamTypeClasses #-}-module Data.Sparse.Internal.CSC where--import Control.Monad (forM_, when)-import Control.Monad.Primitive (PrimMonad(..), PrimState(..))--import qualified Data.Graph as G-import qualified Data.Vector as V-import qualified Data.Vector.Mutable as VM---- import Data.Sparse.Types-import Data.Sparse.Internal.SVector-import Data.Sparse.Internal.Utils-import Numeric.LinearAlgebra.Class---data CscMatrix a =- CscM {- cscNrows :: {-# UNPACK #-} !Int,- cscNcols :: {-# UNPACK #-} !Int,- cscNz :: {-# UNPACK #-} !Int,- cscRowIx :: V.Vector Int,- cscColPtr :: V.Vector Int,- cscVal :: V.Vector a} deriving Eq--instance Functor CscMatrix where- fmap f (CscM m n nz cc rp x) = CscM m n nz cc rp (fmap f x)--instance Foldable CscMatrix where- foldr f z cm = foldr f z (cscVal cm)--instance Show a => Show (CscMatrix a) where- show m'@(CscM m n nz cix rp x) = szs where- szs = unwords ["CSC (",show m, "x", show n,"),",show nz, "NZ ( sparsity",show (spy m'),"), row indices:",show cix,", column pointers:", show rp,", data:",show x]--instance FiniteDim (CscMatrix a) where- type FDSize (CscMatrix a) = (Int, Int)- dim m = (cscNrows m, cscNcols m)--instance HasData (CscMatrix a) where- nnz = cscNz- -instance Sparse (CscMatrix a) where- spy m = fromIntegral (nnz m) / fromIntegral (cscNrows m * cscNcols m)-- ----- * Creation--- | O(N log N) : Copy a Vector containing (row index, column index, entry) into a CSC structure. Sorts the Vector by row columns ( O(log N) ), unzips row indices and data ( O(N) ) and generates the column pointer vector ( 2 O(N) passes )-toCSC :: Int -> Int -> V.Vector (Int, Int, a) -> CscMatrix a-toCSC m n ijxv = CscM m n nz rix crp x where- nz = V.length x- (rix, cix, x) = V.unzip3 (sortWith snd3 ijxv) -- sort by columns- crp = csPtrV (==) m cix---- | O(N) : Rebuilds the (row, column, entry) Vector from the CSC representation. -fromCSC :: CscMatrix a -> V.Vector (Int, Int, a)-fromCSC mc = V.zip3 ii jj xx where (ii,jj,xx) = fromCSC0 mc--fromCSC0 :: CscMatrix a -> (V.Vector Int, V.Vector Int, V.Vector a)-fromCSC0 mc = (rowIx, cols, cscVal mc) where- (_, n) = dim mc- rowIx = cscRowIx mc- l = length rowIx- cp = cscColPtr mc- cols = V.create $ do- colv <- VM.replicate l 0- forM_ [0 .. n-1] (\i -> go colv i 0)- return colv- go :: PrimMonad m => VM.MVector (PrimState m) Int -> Int -> Int -> m () - go vm irow j = when (j < nj) $ do- VM.write vm (j + jmin) irow- go vm irow (succ j) where- jmin = cp V.! irow- jmax = cp V.! (irow + 1)- nj = jmax - jmin----- ** Extract a column--- | O(1) : extract a column from the CSC matrix.-extractColCSC :: CscMatrix a -> Int -> SVector a-extractColCSC (CscM m _ _ rix cp x) jcol = SV m ixs vals where- jmin = cp V.! jcol- jmax = cp V.! (jcol + 1)- ixs = V.slice jmin (jmax - jmin) rix- vals = V.slice jmin (jmax - jmin) x---- ** Extract the diagonal element, if this exists (NB this only holds for a lower triangular matrix)-extractDiagSubdiagCSC :: CscMatrix a -> Int -> Maybe (a, SVector a)-extractDiagSubdiagCSC cscm j- | V.length ixs > 1 && j == V.head ixs = Just (diagEl, subdiagVec)- | otherwise = Nothing where- (SV m ixs vals) = extractColCSC cscm j- diagEl = V.head vals- subdiagVec = SV (m-1) (V.tail ixs) (V.tail vals) --------- | O(N log N) : Transpose CSC matrix-transposeCSC :: CscMatrix a -> CscMatrix a-transposeCSC mm = toCSC n m $ V.zip3 jj ii xx where- (m,n) = dim mm- (ii, jj, xx) = fromCSC0 mm----- * Helpers--cscToGraph :: CscMatrix a -> G.Graph-cscToGraph ll = G.buildG (0, m-1) $ V.toList (V.zip ill jll) -- graph of L- where- (m, _) = dim ll -- dimensions of L = bounds of G(L)- (ill, jll, _) = fromCSC0 ll------- -- example data---- -- row = np.array([0, 0, 1, 2, 2, 2])--- -- col = np.array([0, 2, 2, 0, 1, 2])--- -- data = np.array([1, 2, 3, 4, 5, 6])
src/Data/Sparse/Internal/IntM.hs view
@@ -7,7 +7,7 @@ import GHC.Exts import Data.Complex -import qualified Data.IntMap as IM+import qualified Data.IntMap.Strict as IM
src/Data/Sparse/Internal/TriMatrix.hs view
@@ -23,7 +23,7 @@ import Data.Sparse.Internal.SList import Numeric.LinearAlgebra.Class-import Data.Sparse.Internal.CSC+-- import Data.Sparse.Internal.CSC import Data.Sparse.Internal.SVector import qualified Data.Sparse.Internal.SVector.Mutable as SMV import Data.Sparse.SpMatrix (fromListSM, fromListDenseSM, insertSpMatrix, zeroSM, transposeSM, sparsifySM)
src/Data/Sparse/Internal/Utils.hs view
@@ -3,8 +3,8 @@ import Control.Monad (unless) import Control.Monad.State-import Control.Monad.ST + import Data.Ord (comparing) import qualified Data.Vector as V @@ -20,27 +20,6 @@ ----- | Given a number of rows(resp. columns) `n` and a _sorted_ Vector of Integers in increasing order (containing the row(col) indices of nonzero entries), return the cumulative vector of nonzero entries of length `n + 1` (the "row(col) pointer" of the CSR(CSC) format). NB: Fused count-and-accumulate--- E.g.:--- > csPtrV (==) 4 (V.fromList [1,1,2,3])--- [0,0,2,3,4]-csPtrV :: (a -> Int -> Bool) -> Int -> V.Vector a -> V.Vector Int-csPtrV eqf n xs = V.create createf where- createf :: ST s (VM.MVector s Int)- createf = do- let c = 0- vm <- VM.new (n + 1)- VM.write vm 0 0 -- write `0` at position 0- let loop v ll i count | i == n = return ()- | otherwise = do- let lp = V.length $ V.takeWhile (`eqf` i) ll- count' = count + lp- VM.write v (i + 1) count'- loop v (V.drop lp ll) (succ i) count'- loop vm xs 0 c- return vm -- csrPtrV' eqf n xs = V.create createf where