summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--examples/benchmarks.hs66
-rw-r--r--lib/Data/Packed/Internal/Matrix.hs13
-rw-r--r--lib/Data/Packed/Internal/Vector.hs36
-rw-r--r--lib/Data/Packed/ST.hs2
-rw-r--r--lib/Data/Packed/Vector.hs8
5 files changed, 75 insertions, 50 deletions
diff --git a/examples/benchmarks.hs b/examples/benchmarks.hs
index 6490b25..0e95851 100644
--- a/examples/benchmarks.hs
+++ b/examples/benchmarks.hs
@@ -18,7 +18,7 @@ time act = do
18 18
19-------------------------------------------------------------------------------- 19--------------------------------------------------------------------------------
20 20
21main = sequence_ [bench1,bench2,bench3,bench4, 21main = sequence_ [bench1,bench2,bench4,
22 bench5 1000000 3, 22 bench5 1000000 3,
23 bench5 100000 50] 23 bench5 100000 50]
24 24
@@ -26,16 +26,25 @@ w :: Vector Double
26w = constant 1 5000000 26w = constant 1 5000000
27w2 = 1 * w 27w2 = 1 * w
28 28
29v = flatten $ ident 500 :: Vector Double
30
31
29bench1 = do 32bench1 = do
33 time $ print$ vectorMax (w+w2) -- evaluate it
30 putStrLn "Sum of a vector with 5M doubles:" 34 putStrLn "Sum of a vector with 5M doubles:"
31 print$ vectorMax (w+w2) -- evaluate it 35 print $ vectorMax v -- evaluate it
32 time $ printf " BLAS: %.2f: " $ sumVB w 36-- time $ printf " BLAS: %.2f: " $ sumVB w
33 time $ printf " Haskell: %.2f: " $ sumVH w 37 time $ printf " Haskell: %.2f: " $ sumVH w
34 time $ printf " BLAS: %.2f: " $ sumVB w 38 time $ printf " BLAS: %.2f: " $ w <.> w2
35 time $ printf " Haskell: %.2f: " $ sumVH w 39 time $ printf " Haskell: %.2f: " $ sumVH w
36 time $ printf " innerH: %.2f: " $ innerH w w2 40 time $ printf " innerH: %.2f: " $ innerH w w2
37 41 time $ printf "foldVector: %.2f: " $ sumVector w
38sumVB v = constant 1 (dim v) <.> v 42 let getPos k s = if k `mod` 500 < 200 && w@>k > 0 then k:s else s
43 putStrLn "foldLoop for element selection:"
44 time $ print $ (`divMod` 500) $ maximum $ foldLoop getPos [] (dim w)
45 putStrLn "constant 5M:"
46 time $ print $ constant (1::Double) 5000001 @> 7
47 time $ print $ constant i 5000001 @> 7
39 48
40sumVH v = go (d - 1) 0 49sumVH v = go (d - 1) 0
41 where 50 where
@@ -51,8 +60,10 @@ innerH u v = go (d - 1) 0
51 go 0 s = s + (u @> 0) * (v @> 0) 60 go 0 s = s + (u @> 0) * (v @> 0)
52 go !j !s = go (j - 1) (s + (u @> j) * (v @> j)) 61 go !j !s = go (j - 1) (s + (u @> j) * (v @> j))
53 62
54-- These functions are much faster if the library 63
55-- is configured with -funsafe 64-- sumVector = foldVectorG (\k v s -> v k + s) 0.0
65sumVector = foldVector (+) 0.0
66
56 67
57-------------------------------------------------------------------------------- 68--------------------------------------------------------------------------------
58 69
@@ -89,39 +100,6 @@ manymult n r = foldl1' (<>) (map r angles)
89 100
90-------------------------------------------------------------------------------- 101--------------------------------------------------------------------------------
91 102
92bench3 = do
93 putStrLn "-------------------------------------------------------"
94 putStrLn "foldVector"
95 let v = flatten $ ident 500 :: Vector Double
96 print $ vectorMax v -- evaluate it
97
98 putStrLn "sum, dim=5M:"
99 -- time $ print $ foldLoop (\k s -> w@>k + s) 0.0 (dim w)
100 time $ print $ sumVector w
101
102 putStrLn "sum, dim=0.25M:"
103 --time $ print $ foldLoop (\k s -> v@>k + s) 0.0 (dim v)
104 time $ print $ sumVector v
105
106 let getPos k s = if k `mod` 500 < 200 && v@>k > 0 then k:s else s
107 putStrLn "foldLoop for element selection, dim=0.25M:"
108 time $ print $ (`divMod` 500) $ maximum $ foldLoop getPos [] (dim v)
109
110foldLoop f s d = go (d - 1) s
111 where
112 go 0 s = f (0::Int) s
113 go !j !s = go (j - 1) (f j s)
114
115foldVector f s v = foldLoop g s (dim v)
116 where g !k !s = f k (v@>) s
117 {-# INLINE g #-} -- Thanks Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479)
118
119sumVector = foldVector (\k v s -> v k + s) 0.0
120
121-- foldVector is slower if used in two places unless we use the above INLINE
122-- this does not happen with foldLoop
123--------------------------------------------------------------------------------
124
125bench4 = do 103bench4 = do
126 putStrLn "-------------------------------------------------------" 104 putStrLn "-------------------------------------------------------"
127 putStrLn "1000x1000 inverse" 105 putStrLn "1000x1000 inverse"
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index 9def473..09f081a 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -217,14 +217,17 @@ class (Storable a, Floating a) => Element a where
217 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix 217 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix
218 -> Matrix a -> Matrix a 218 -> Matrix a -> Matrix a
219 transdata :: Int -> Vector a -> Int -> Vector a 219 transdata :: Int -> Vector a -> Int -> Vector a
220 constantD :: a -> Int -> Vector a
220 221
221instance Element Double where 222instance Element Double where
222 subMatrixD = subMatrixR 223 subMatrixD = subMatrixR
223 transdata = transdata' 224 transdata = transdata'
225 constantD = constant'
224 226
225instance Element (Complex Double) where 227instance Element (Complex Double) where
226 subMatrixD = subMatrixC 228 subMatrixD = subMatrixC
227 transdata = transdata' 229 transdata = transdata'
230 constantD = constant'
228 231
229------------------------------------------------------------------- 232-------------------------------------------------------------------
230 233
@@ -256,6 +259,16 @@ transdata' c1 v c2 =
256 259
257---------------------------------------------------------------------- 260----------------------------------------------------------------------
258 261
262constant' v n = unsafePerformIO $ do
263 w <- createVector n
264 withForeignPtr (fptr w) $ \p -> do
265 let go (-1) = return ()
266 go !k = pokeElemOff p k v >> go (k-1)
267 go (n-1)
268 return w
269
270----------------------------------------------------------------------
271
259-- | extraction of a submatrix from a real matrix 272-- | extraction of a submatrix from a real matrix
260subMatrixR :: (Int,Int) -> (Int,Int) -> Matrix Double -> Matrix Double 273subMatrixR :: (Int,Int) -> (Int,Int) -> Matrix Double -> Matrix Double
261subMatrixR (r0,c0) (rt,ct) x' = unsafePerformIO $ do 274subMatrixR (r0,c0) (rt,ct) x' = unsafePerformIO $ do
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index f590919..dd9b9b6 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -1,4 +1,4 @@
1{-# LANGUAGE MagicHash, CPP, UnboxedTuples #-} 1{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns #-}
2----------------------------------------------------------------------------- 2-----------------------------------------------------------------------------
3-- | 3-- |
4-- Module : Data.Packed.Internal.Vector 4-- Module : Data.Packed.Internal.Vector
@@ -182,7 +182,7 @@ asComplex v = V { dim = dim v `div` 2, fptr = castForeignPtr (fptr v) }
182 182
183-- | map on Vectors 183-- | map on Vectors
184liftVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b 184liftVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b
185liftVector f = fromList . map f . toList 185liftVector = mapVector
186 186
187-- | zipWith for Vectors 187-- | zipWith for Vectors
188liftVector2 :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c 188liftVector2 :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c
@@ -196,3 +196,35 @@ cloneVector (v@V {dim=n}) = do
196 let f _ s _ d = copyArray d s n >> return 0 196 let f _ s _ d = copyArray d s n >> return 0
197 app2 f vec v vec r "cloneVector" 197 app2 f vec v vec r "cloneVector"
198 return r 198 return r
199
200------------------------------------------------------------------
201
202mapVector f v = unsafePerformIO $ do
203 w <- createVector (dim v)
204 withForeignPtr (fptr v) $ \p ->
205 withForeignPtr (fptr w) $ \q -> do
206 let go (-1) = return ()
207 go !k = do x <- peekElemOff p k
208 pokeElemOff q k (f x)
209 go (k-1)
210 go (dim v -1)
211 return w
212{-# INLINE mapVector #-}
213
214foldVector f x v = unsafePerformIO $
215 withForeignPtr (fptr (v::Vector Double)) $ \p -> do
216 let go (-1) s = return s
217 go !k !s = do y <- peekElemOff p k
218 go (k-1::Int) (f y s)
219 go (dim v -1) x
220{-# INLINE foldVector #-}
221
222foldLoop f s0 d = go (d - 1) s0
223 where
224 go 0 s = f (0::Int) s
225 go !j !s = go (j - 1) (f j s)
226
227foldVectorG f s0 v = foldLoop g s0 (dim v)
228 where g !k !s = f k (at' v) s
229 {-# INLINE g #-} -- Thanks to Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479)
230{-# INLINE foldVectorG #-}
diff --git a/lib/Data/Packed/ST.hs b/lib/Data/Packed/ST.hs
index 9c9c3b3..a47da05 100644
--- a/lib/Data/Packed/ST.hs
+++ b/lib/Data/Packed/ST.hs
@@ -93,7 +93,7 @@ writeVector = safeIndexV unsafeWriteVector
93newUndefinedVector :: Element t => Int -> ST s (STVector s t) 93newUndefinedVector :: Element t => Int -> ST s (STVector s t)
94newUndefinedVector = unsafeIOToST . fmap STVector . createVector 94newUndefinedVector = unsafeIOToST . fmap STVector . createVector
95 95
96{-# NOINLINE newVector #-} 96{-# INLINE newVector #-}
97newVector :: Element t => t -> Int -> ST s (STVector s t) 97newVector :: Element t => t -> Int -> ST s (STVector s t)
98newVector x n = do 98newVector x n = do
99 v <- newUndefinedVector n 99 v <- newUndefinedVector n
diff --git a/lib/Data/Packed/Vector.hs b/lib/Data/Packed/Vector.hs
index 0bbbc34..b85f0bd 100644
--- a/lib/Data/Packed/Vector.hs
+++ b/lib/Data/Packed/Vector.hs
@@ -19,12 +19,13 @@ module Data.Packed.Vector (
19 subVector, join, 19 subVector, join,
20 constant, linspace, 20 constant, linspace,
21 vectorMax, vectorMin, vectorMaxIndex, vectorMinIndex, 21 vectorMax, vectorMin, vectorMaxIndex, vectorMinIndex,
22 liftVector, liftVector2 22 liftVector, liftVector2,
23 foldLoop, foldVector, foldVectorG, mapVector
23) where 24) where
24 25
25import Data.Packed.Internal 26import Data.Packed.Internal
26import Numeric.GSL.Vector 27import Numeric.GSL.Vector
27import Data.Packed.ST 28-- import Data.Packed.ST
28 29
29{- | Creates a real vector containing a range of values: 30{- | Creates a real vector containing a range of values:
30 31
@@ -55,4 +56,5 @@ vectorMinIndex = round . toScalarR MinIdx
557 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@ 567 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@
56-} 57-}
57constant :: Element a => a -> Int -> Vector a 58constant :: Element a => a -> Int -> Vector a
58constant x n = runSTVector (newVector x n) 59-- constant x n = runSTVector (newVector x n)
60constant = constantD -- about 2x faster