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