summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--examples/speed.hs46
-rw-r--r--examples/tests.hs18
-rw-r--r--lib/Data/Packed/Internal/Vector.hs6
-rw-r--r--lib/Data/Packed/Vector.hs6
4 files changed, 59 insertions, 17 deletions
diff --git a/examples/speed.hs b/examples/speed.hs
index 865e3e1..22d7220 100644
--- a/examples/speed.hs
+++ b/examples/speed.hs
@@ -1,21 +1,42 @@
1{- speed tests 1{- speed tests
2 2
3$ ghc --make -O speed 3GNU-Octave (see speed.m in this folder):
4In my machine:
5$ ./speed 5 100000 1
6(3><3)
7 [ -1.7877285611885504e-2, 0.0, -0.9998401885597121
8 , 0.0, 1.0, 0.0
9 , 0.9998401885597168, 0.0, -1.7877285611891697e-2 ]
100.29 CPU seconds
11 4
12GNU-Octave:
13./speed.m 5./speed.m
14 -0.017877255967426 0.000000000000000 -0.999840189089781 6 -0.017877255967426 0.000000000000000 -0.999840189089781
15 0.000000000000000 1.000000000000000 0.000000000000000 7 0.000000000000000 1.000000000000000 0.000000000000000
16 0.999840189089763 0.000000000000000 -0.017877255967417 8 0.999840189089763 0.000000000000000 -0.017877255967417
179.69 seconds 99.69 seconds
18 10
11Mathematica:
12
13rot[a_]:={{ Cos[a], 0, Sin[a]},
14 { 0, 1, 0},
15 { -Sin[a],0,Cos[a]}}//N
16
17test := Timing[
18 n = 100000;
19 angles = Range[0.,n-1]/(n-1);
20 Fold[Dot,IdentityMatrix[3],rot/@angles] // MatrixForm
21]
22
232.08013 Second
24 {{\(-0.017877255967432837`\), "0.`", \(-0.9998401890898042`\)},
25 {"0.`", "1.`", "0.`"},
26 {"0.9998401890898042`", "0.`", \(-0.017877255967432837`\)}}
27
28$ ghc --make -O speed
29
30$ ./speed 5 100000 1
31(3><3)
32 [ -1.7877255967425523e-2, 0.0, -0.9998401890897632
33 , 0.0, 1.0, 0.0
34 , 0.999840189089781, 0.0, -1.7877255967416586e-2 ]
350.33 CPU seconds
36
37cos 50000 = -0.0178772559665563
38sin 50000 = -0.999840189089790
39
19-} 40-}
20 41
21import Numeric.LinearAlgebra 42import Numeric.LinearAlgebra
@@ -31,7 +52,7 @@ timing act = do
31 t0 <- getCPUTime 52 t0 <- getCPUTime
32 act 53 act
33 t1 <- getCPUTime 54 t1 <- getCPUTime
34 printf "%.2f CPU seconds\n" $ (fromIntegral ((t1 - t0) `div` (10^10)) / 100 :: Double) 55 printf "%.2f CPU seconds\n" $ (fromIntegral ((t1 - t0) `div` (10^10)) / 100 :: Double) :: IO ()
35 56
36op a b = trans $ (trans a) <> (trans b) 57op a b = trans $ (trans a) <> (trans b)
37 58
@@ -52,9 +73,8 @@ rot a = (3><3) [ c,0,s
52 where c = cos a 73 where c = cos a
53 s = sin a 74 s = sin a
54 75
55 76fun n r = foldl1' (<>) (map r angles)
56fun n r = foldl1' (<>) (map r [0,delta..1]) where delta = 1 /(fromIntegral n-1) 77 where angles = toList $ linspace n (0,1)
57
58 78
59main = do 79main = do
60 args <- getArgs 80 args <- getArgs
diff --git a/examples/tests.hs b/examples/tests.hs
index 7b85a9b..90437bb 100644
--- a/examples/tests.hs
+++ b/examples/tests.hs
@@ -13,6 +13,7 @@ import Test.QuickCheck hiding (test)
13import Test.HUnit hiding ((~:),test) 13import Test.HUnit hiding ((~:),test)
14import System.Random(randomRs,mkStdGen) 14import System.Random(randomRs,mkStdGen)
15import System.Info 15import System.Info
16import Data.List(foldl1')
16 17
17type RM = Matrix Double 18type RM = Matrix Double
18type CM = Matrix (Complex Double) 19type CM = Matrix (Complex Double)
@@ -310,6 +311,20 @@ mulF a b = multiply' ColumnMajor a b
310 311
311--------------------------------------------------------------------- 312---------------------------------------------------------------------
312 313
314rot :: Double -> Matrix Double
315rot a = (3><3) [ c,0,s
316 , 0,1,0
317 ,-s,0,c ]
318 where c = cos a
319 s = sin a
320
321fun n = foldl1' (<>) (map rot angles)
322 where angles = toList $ linspace n (0,1)
323
324rotTest = fun (10^5) :~12~: rot 5E4
325
326---------------------------------------------------------------------
327
313tests = do 328tests = do
314 putStrLn "--------- internal -----" 329 putStrLn "--------- internal -----"
315 quickCheck ((\m -> m == trans m).sym :: Sym Double -> Bool) 330 quickCheck ((\m -> m == trans m).sym :: Sym Double -> Bool)
@@ -322,6 +337,9 @@ tests = do
322 quickCheck $ \m -> m == asFortran (m :: CM) 337 quickCheck $ \m -> m == asFortran (m :: CM)
323 quickCheck $ \m -> m == (asC . asFortran) (m :: RM) 338 quickCheck $ \m -> m == (asC . asFortran) (m :: RM)
324 quickCheck $ \m -> m == (asC . asFortran) (m :: CM) 339 quickCheck $ \m -> m == (asC . asFortran) (m :: CM)
340 runTestTT $ TestList
341 [ test "1E5 rots" rotTest
342 ]
325 putStrLn "--------- multiply ----" 343 putStrLn "--------- multiply ----"
326 quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == mulF m1 (m2 :: RM) 344 quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == mulF m1 (m2 :: RM)
327 quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == mulF m1 (m2 :: CM) 345 quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == mulF m1 (m2 :: CM)
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index 76bd4d1..082e09d 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -80,6 +80,8 @@ fromList l = unsafePerformIO $ do
80 f // vec v // check "fromList" [] 80 f // vec v // check "fromList" []
81 return v 81 return v
82 82
83safeRead v f = unsafePerformIO $ withForeignPtr (fptr v) $ const $ f (ptr v)
84
83{- | extracts the Vector elements to a list 85{- | extracts the Vector elements to a list
84 86
85@> toList (linspace 5 (1,10)) 87@> toList (linspace 5 (1,10))
@@ -87,7 +89,7 @@ fromList l = unsafePerformIO $ do
87 89
88-} 90-}
89toList :: Storable a => Vector a -> [a] 91toList :: Storable a => Vector a -> [a]
90toList v = unsafePerformIO $ peekArray (dim v) (ptr v) 92toList v = safeRead v $ peekArray (dim v)
91 93
92-- | an alternative to 'fromList' with explicit dimension, used also in the instances for Show (Vector a). 94-- | an alternative to 'fromList' with explicit dimension, used also in the instances for Show (Vector a).
93(|>) :: (Storable a) => Int -> [a] -> Vector a 95(|>) :: (Storable a) => Int -> [a] -> Vector a
@@ -96,7 +98,7 @@ n |> l = if length l == n then fromList l else error "|> with wrong size"
96 98
97-- | access to Vector elements without range checking 99-- | access to Vector elements without range checking
98at' :: Storable a => Vector a -> Int -> a 100at' :: Storable a => Vector a -> Int -> a
99at' v n = unsafePerformIO $ peekElemOff (ptr v) n 101at' v n = safeRead v $ flip peekElemOff n
100 102
101-- | access to Vector elements with range checking. 103-- | access to Vector elements with range checking.
102at :: Storable a => Vector a -> Int -> a 104at :: Storable a => Vector a -> Int -> a
diff --git a/lib/Data/Packed/Vector.hs b/lib/Data/Packed/Vector.hs
index 4c331ef..60fa7c1 100644
--- a/lib/Data/Packed/Vector.hs
+++ b/lib/Data/Packed/Vector.hs
@@ -32,8 +32,10 @@ import Numeric.GSL.Vector
325 |> [-3.0,-0.5,2.0,4.5,7.0]@ 325 |> [-3.0,-0.5,2.0,4.5,7.0]@
33-} 33-}
34linspace :: Int -> (Double, Double) -> Vector Double 34linspace :: Int -> (Double, Double) -> Vector Double
35linspace n (a,b) = fromList [a, a+delta .. b] 35linspace n (a,b) = add a $ scale s $ fromList [0 .. fromIntegral n-1]
36 where delta = (b-a)/(fromIntegral n -1) 36 where scale = vectorMapValR Scale
37 add = vectorMapValR AddConstant
38 s = (b-a)/fromIntegral (n-1)
37 39
38vectorMax :: Vector Double -> Double 40vectorMax :: Vector Double -> Double
39vectorMax = toScalarR Max 41vectorMax = toScalarR Max