foldVector benchmark

3 files changed, 30 insertions, 40 deletions

diff --git a/examples/benchmarks.hs b/examples/benchmarks.hs
index 398d1ce..073a95d 100644
--- a/examples/benchmarks.hs
+++ b/examples/benchmarks.hs
@@ -70,20 +70,34 @@ rot a = (3><3) [ c,0,s
 
 manymult n r = foldl1' (<>) (map r angles)
     where angles = toList $ linspace n (0,1)
+          -- angles = map (k*) [0..n']
+          -- n' = fromIntegral n - 1
+          -- k  = recip n'
 
 --------------------------------------------------------------------------------
 
 bench3 = do
     putStrLn "-------------------------------------------------------"
-    putStrLn "foldVector:"
-    let v = flatten $ ident 555 :: Vector Double
+    putStrLn "foldVector"
+    let v = flatten $ ident 500 :: Vector Double
     print $ vectorMax v  -- evaluate it
-    let getPos k s = if k `rem` 555 < 200 && v@>k > 0 then k:s else s
-    time $ print $ sum $ foldVector getPos [] v
-    time $ print $ foldVector (\k s -> w@>k + s) 0.0 w
-
-foldVector f s v = go (d - 1) s
+    let getPos k s = if k `mod` 500 < 200 && v@>k > 0 then k:s else s
+    putStrLn "indices extraction, dim=0.25M:"
+    time $ print $ (`divMod` 500) $ maximum $ foldLoop getPos [] (dim v)
+    putStrLn "sum, dim=30M:"
+    --time $ print $ foldLoop (\k s -> w@>k + s) 0.0 (dim w)
+    time $ print $ foldVector (\k v s -> v k + s) 0.0 w
+    putStrLn "sum, dim=0.25M:"
+    --time $ print $ foldVector (\k v s -> v k + s) 0.0 v
+    time $ print $ foldLoop (\k s -> v@>k + s) 0.0 (dim v)
+
+-- foldVector is slower if it is used in two places. (!?)
+-- this does not happen with foldLoop
+
+foldLoop f s d = go (d - 1) s
      where
-       d = dim v
        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
diff --git a/hmatrix.cabal b/hmatrix.cabal
index c2ee18a..75e2984 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -25,6 +25,10 @@ flag mkl
     description:    Link with Intel's MKL optimized libraries.
     default:        False
 
+flag gsl
+    description:    Link with GSL unoptimized blas.
+    default:        False
+
 flag unsafe
     description: Compile the library with bound checking disabled.
     default: False
@@ -108,7 +112,10 @@ library
       else
         extra-libraries:   gsl mkl_lapack mkl_intel mkl_sequential mkl_core
     else
-      extra-libraries:  gsl blas lapack
+      if flag(gsl)
+        extra-libraries:  gsl gslcblas lapack
+      else
+        extra-libraries:  gsl blas lapack
 
     cc-options:         -O4
     ghc-prof-options:   -auto-all
diff --git a/lib/Numeric/GSL/Matrix.hs b/lib/Numeric/GSL/Matrix.hs
index cb3fa8f..f62bb82 100644
--- a/lib/Numeric/GSL/Matrix.hs
+++ b/lib/Numeric/GSL/Matrix.hs
@@ -92,23 +92,6 @@ foreign import ccall "gsl-aux.h eigensystemC" c_eigH :: TCMVCM
 
 {- | Singular value decomposition of a real matrix, using /gsl_linalg_SV_decomp_mod/:
 
-@\> let (u,s,v) = svdg $ 'fromLists' [[1,2,3],[-4,1,7]]
-\
-\> u
-0.310 -0.951
-0.951  0.310
-\
-\> s
-8.497 2.792
-\
-\> v
--0.411 -0.785
- 0.185 -0.570
- 0.893 -0.243
-\
-\> u \<\> 'diag' s \<\> 'trans' v
- 1. 2. 3.
--4. 1. 7.@
 
 -}
 svdg :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
@@ -128,20 +111,6 @@ foreign import ccall "gsl-aux.h svd" c_svd :: TMMVM
 
 {- | QR decomposition of a real matrix using /gsl_linalg_QR_decomp/ and /gsl_linalg_QR_unpack/.
 
-@\> let (q,r) = qr $ 'fromLists' [[1,3,5,7],[2,0,-2,4]]
-\
-\> q
--0.447 -0.894
--0.894  0.447
-\
-\> r
--2.236 -1.342 -0.447 -6.708
-    0. -2.683 -5.367 -4.472
-\
-\> q \<\> r
-1.000 3.000  5.000 7.000
-2.000    0. -2.000 4.000@
-
 -}
 qr :: Matrix Double -> (Matrix Double, Matrix Double)
 qr = qr' . cmat