From e1b4cc06a5f98e576524b37ad0d9132f0678d722 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Fri, 14 Nov 2008 11:01:14 +0000 Subject: constantD --- examples/benchmarks.hs | 66 +++++++++++++------------------------- lib/Data/Packed/Internal/Matrix.hs | 13 ++++++++ lib/Data/Packed/Internal/Vector.hs | 36 +++++++++++++++++++-- lib/Data/Packed/ST.hs | 2 +- lib/Data/Packed/Vector.hs | 8 +++-- 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 -------------------------------------------------------------------------------- -main = sequence_ [bench1,bench2,bench3,bench4, +main = sequence_ [bench1,bench2,bench4, bench5 1000000 3, bench5 100000 50] @@ -26,16 +26,25 @@ w :: Vector Double w = constant 1 5000000 w2 = 1 * w +v = flatten $ ident 500 :: Vector Double + + bench1 = do + time $ print$ vectorMax (w+w2) -- evaluate it putStrLn "Sum of a vector with 5M doubles:" - print$ vectorMax (w+w2) -- evaluate it - time $ printf " BLAS: %.2f: " $ sumVB w - time $ printf " Haskell: %.2f: " $ sumVH w - time $ printf " BLAS: %.2f: " $ sumVB w - time $ printf " Haskell: %.2f: " $ sumVH w - time $ printf " innerH: %.2f: " $ innerH w w2 - -sumVB v = constant 1 (dim v) <.> v + print $ vectorMax v -- evaluate it +-- time $ printf " BLAS: %.2f: " $ sumVB w + time $ printf " Haskell: %.2f: " $ sumVH w + time $ printf " BLAS: %.2f: " $ w <.> w2 + time $ printf " Haskell: %.2f: " $ sumVH w + time $ printf " innerH: %.2f: " $ innerH w w2 + time $ printf "foldVector: %.2f: " $ sumVector w + let getPos k s = if k `mod` 500 < 200 && w@>k > 0 then k:s else s + putStrLn "foldLoop for element selection:" + time $ print $ (`divMod` 500) $ maximum $ foldLoop getPos [] (dim w) + putStrLn "constant 5M:" + time $ print $ constant (1::Double) 5000001 @> 7 + time $ print $ constant i 5000001 @> 7 sumVH v = go (d - 1) 0 where @@ -51,8 +60,10 @@ innerH u v = go (d - 1) 0 go 0 s = s + (u @> 0) * (v @> 0) go !j !s = go (j - 1) (s + (u @> j) * (v @> j)) --- These functions are much faster if the library --- is configured with -funsafe + +-- sumVector = foldVectorG (\k v s -> v k + s) 0.0 +sumVector = foldVector (+) 0.0 + -------------------------------------------------------------------------------- @@ -89,39 +100,6 @@ manymult n r = foldl1' (<>) (map r angles) -------------------------------------------------------------------------------- -bench3 = do - putStrLn "-------------------------------------------------------" - putStrLn "foldVector" - let v = flatten $ ident 500 :: Vector Double - print $ vectorMax v -- evaluate it - - putStrLn "sum, dim=5M:" - -- time $ print $ foldLoop (\k s -> w@>k + s) 0.0 (dim w) - time $ print $ sumVector w - - putStrLn "sum, dim=0.25M:" - --time $ print $ foldLoop (\k s -> v@>k + s) 0.0 (dim v) - time $ print $ sumVector v - - let getPos k s = if k `mod` 500 < 200 && v@>k > 0 then k:s else s - putStrLn "foldLoop for element selection, dim=0.25M:" - time $ print $ (`divMod` 500) $ maximum $ foldLoop getPos [] (dim v) - -foldLoop f s d = go (d - 1) s - where - go 0 s = f (0::Int) s - go !j !s = go (j - 1) (f j s) - -foldVector f s v = foldLoop g s (dim v) - where g !k !s = f k (v@>) s - {-# INLINE g #-} -- Thanks Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479) - -sumVector = foldVector (\k v s -> v k + s) 0.0 - --- foldVector is slower if used in two places unless we use the above INLINE --- this does not happen with foldLoop --------------------------------------------------------------------------------- - bench4 = do putStrLn "-------------------------------------------------------" 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 -> (Int,Int) -- ^ (rt,ct) dimensions of submatrix -> Matrix a -> Matrix a transdata :: Int -> Vector a -> Int -> Vector a + constantD :: a -> Int -> Vector a instance Element Double where subMatrixD = subMatrixR transdata = transdata' + constantD = constant' instance Element (Complex Double) where subMatrixD = subMatrixC transdata = transdata' + constantD = constant' ------------------------------------------------------------------- @@ -256,6 +259,16 @@ transdata' c1 v c2 = ---------------------------------------------------------------------- +constant' v n = unsafePerformIO $ do + w <- createVector n + withForeignPtr (fptr w) $ \p -> do + let go (-1) = return () + go !k = pokeElemOff p k v >> go (k-1) + go (n-1) + return w + +---------------------------------------------------------------------- + -- | extraction of a submatrix from a real matrix subMatrixR :: (Int,Int) -> (Int,Int) -> Matrix Double -> Matrix Double subMatrixR (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 @@ -{-# LANGUAGE MagicHash, CPP, UnboxedTuples #-} +{-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns #-} ----------------------------------------------------------------------------- -- | -- Module : Data.Packed.Internal.Vector @@ -182,7 +182,7 @@ asComplex v = V { dim = dim v `div` 2, fptr = castForeignPtr (fptr v) } -- | map on Vectors liftVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b -liftVector f = fromList . map f . toList +liftVector = mapVector -- | zipWith for Vectors liftVector2 :: (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 let f _ s _ d = copyArray d s n >> return 0 app2 f vec v vec r "cloneVector" return r + +------------------------------------------------------------------ + +mapVector f v = unsafePerformIO $ do + w <- createVector (dim v) + withForeignPtr (fptr v) $ \p -> + withForeignPtr (fptr w) $ \q -> do + let go (-1) = return () + go !k = do x <- peekElemOff p k + pokeElemOff q k (f x) + go (k-1) + go (dim v -1) + return w +{-# INLINE mapVector #-} + +foldVector f x v = unsafePerformIO $ + withForeignPtr (fptr (v::Vector Double)) $ \p -> do + let go (-1) s = return s + go !k !s = do y <- peekElemOff p k + go (k-1::Int) (f y s) + go (dim v -1) x +{-# INLINE foldVector #-} + +foldLoop f s0 d = go (d - 1) s0 + where + go 0 s = f (0::Int) s + go !j !s = go (j - 1) (f j s) + +foldVectorG f s0 v = foldLoop g s0 (dim v) + where g !k !s = f k (at' v) s + {-# INLINE g #-} -- Thanks to Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479) +{-# 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 newUndefinedVector :: Element t => Int -> ST s (STVector s t) newUndefinedVector = unsafeIOToST . fmap STVector . createVector -{-# NOINLINE newVector #-} +{-# INLINE newVector #-} newVector :: Element t => t -> Int -> ST s (STVector s t) newVector x n = do 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 ( subVector, join, constant, linspace, vectorMax, vectorMin, vectorMaxIndex, vectorMinIndex, - liftVector, liftVector2 + liftVector, liftVector2, + foldLoop, foldVector, foldVectorG, mapVector ) where import Data.Packed.Internal import Numeric.GSL.Vector -import Data.Packed.ST +-- import Data.Packed.ST {- | Creates a real vector containing a range of values: @@ -55,4 +56,5 @@ vectorMinIndex = round . toScalarR MinIdx 7 |> [2.0,2.0,2.0,2.0,2.0,2.0,2.0]@ -} constant :: Element a => a -> Int -> Vector a -constant x n = runSTVector (newVector x n) +-- constant x n = runSTVector (newVector x n) +constant = constantD -- about 2x faster -- cgit v1.2.3