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 @@
 {- speed tests
 
-$ ghc --make -O speed
-In my machine:
-$ ./speed 5 100000 1
-(3><3)
- [ -1.7877285611885504e-2, 0.0,    -0.9998401885597121
- ,                    0.0, 1.0,                    0.0
- ,     0.9998401885597168, 0.0, -1.7877285611891697e-2 ]
-0.29 CPU seconds
+GNU-Octave (see speed.m in this folder):
 
-GNU-Octave:
 ./speed.m
   -0.017877255967426   0.000000000000000  -0.999840189089781
    0.000000000000000   1.000000000000000   0.000000000000000
    0.999840189089763   0.000000000000000  -0.017877255967417
 9.69 seconds
 
+Mathematica:
+
+rot[a_]:={{ Cos[a], 0, Sin[a]},
+          { 0,      1,      0},
+          { -Sin[a],0,Cos[a]}}//N
+
+test := Timing[
+    n = 100000;
+    angles = Range[0.,n-1]/(n-1);
+    Fold[Dot,IdentityMatrix[3],rot/@angles] // MatrixForm
+]
+
+2.08013 Second
+     {{\(-0.017877255967432837`\), "0.`", \(-0.9998401890898042`\)},
+      {"0.`", "1.`", "0.`"},
+      {"0.9998401890898042`", "0.`", \(-0.017877255967432837`\)}}
+
+$ ghc --make -O speed
+
+$ ./speed 5 100000 1
+(3><3)
+ [ -1.7877255967425523e-2, 0.0,    -0.9998401890897632
+ ,                    0.0, 1.0,                    0.0
+ ,      0.999840189089781, 0.0, -1.7877255967416586e-2 ]
+0.33 CPU seconds
+
+cos 50000 = -0.0178772559665563
+sin 50000 = -0.999840189089790
+
 -}
 
 import Numeric.LinearAlgebra
@@ -31,7 +52,7 @@ timing act = do
     t0 <- getCPUTime
     act
     t1 <- getCPUTime
-    printf "%.2f CPU seconds\n" $ (fromIntegral ((t1 - t0) `div` (10^10)) / 100 :: Double)
+    printf "%.2f CPU seconds\n" $ (fromIntegral ((t1 - t0) `div` (10^10)) / 100 :: Double) :: IO ()
 
 op a b = trans $ (trans a) <> (trans b)
 
@@ -52,9 +73,8 @@ rot a = (3><3) [ c,0,s
     where c = cos a
           s = sin a
 
-
-fun n r = foldl1' (<>) (map r [0,delta..1]) where delta = 1 /(fromIntegral n-1)
-
+fun n r = foldl1' (<>) (map r angles)
+    where angles = toList $ linspace n (0,1)
 
 main = do
     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)
 import Test.HUnit hiding ((~:),test)
 import System.Random(randomRs,mkStdGen)
 import System.Info
+import Data.List(foldl1')
 
 type RM = Matrix Double
 type CM = Matrix (Complex Double)
@@ -310,6 +311,20 @@ mulF a b = multiply' ColumnMajor a b
 
 ---------------------------------------------------------------------
 
+rot :: Double -> Matrix Double
+rot a = (3><3) [ c,0,s
+               , 0,1,0
+               ,-s,0,c ]
+    where c = cos a
+          s = sin a
+
+fun n = foldl1' (<>) (map rot angles)
+    where angles = toList $ linspace n (0,1)
+
+rotTest = fun (10^5) :~12~: rot 5E4
+
+---------------------------------------------------------------------
+
 tests = do
     putStrLn "--------- internal -----"
     quickCheck ((\m -> m == trans m).sym :: Sym Double -> Bool)
@@ -322,6 +337,9 @@ tests = do
     quickCheck $ \m -> m == asFortran (m :: CM)
     quickCheck $ \m -> m == (asC . asFortran) (m :: RM)
     quickCheck $ \m -> m == (asC . asFortran) (m :: CM)
+    runTestTT $ TestList
+     [ test "1E5 rots" rotTest
+     ]
     putStrLn "--------- multiply ----"
     quickCheck $ \(PairM m1 m2) -> mulC m1 m2 == mulF m1 (m2 :: RM)
     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
     f // vec v // check "fromList" []
     return v
 
+safeRead v f = unsafePerformIO $ withForeignPtr (fptr v) $ const $ f (ptr v)
+
 {- | extracts the Vector elements to a list
 
 @> toList (linspace 5 (1,10))
@@ -87,7 +89,7 @@ fromList l = unsafePerformIO $ do
 
 -}
 toList :: Storable a => Vector a -> [a]
-toList v = unsafePerformIO $ peekArray (dim v) (ptr v)
+toList v = safeRead v $ peekArray (dim v)
 
 -- | an alternative to 'fromList' with explicit dimension, used also in the instances for Show (Vector a).
 (|>) :: (Storable a) => Int -> [a] -> Vector a
@@ -96,7 +98,7 @@ n |> l = if length l == n then fromList l else error "|> with wrong size"
 
 -- | access to Vector elements without range checking
 at' :: Storable a => Vector a -> Int -> a
-at' v n = unsafePerformIO $ peekElemOff (ptr v) n
+at' v n = safeRead v $ flip peekElemOff n
 
 -- | access to Vector elements with range checking.
 at :: 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
 5 |> [-3.0,-0.5,2.0,4.5,7.0]@
 -}
 linspace :: Int -> (Double, Double) -> Vector Double
-linspace n (a,b) = fromList [a, a+delta .. b]
-    where delta = (b-a)/(fromIntegral n -1)
+linspace n (a,b) = add a $ scale s  $ fromList [0 .. fromIntegral n-1]
+    where scale = vectorMapValR Scale
+          add   = vectorMapValR AddConstant
+          s = (b-a)/fromIntegral (n-1)
 
 vectorMax :: Vector Double -> Double
 vectorMax = toScalarR Max