summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/Data/Packed/Convert.hs101
-rw-r--r--lib/Data/Packed/Internal/Matrix.hs24
-rw-r--r--lib/Data/Packed/Internal/Vector.hs75
-rw-r--r--lib/Data/Packed/ST.hs4
-rw-r--r--lib/Numeric/GSL/Internal.hs8
-rw-r--r--lib/Numeric/GSL/Minimization.hs4
-rw-r--r--lib/Numeric/GSL/ODE.hs4
-rw-r--r--lib/Numeric/GSL/Root.hs4
-rw-r--r--lib/Numeric/LinearAlgebra/Instances.hs14
-rw-r--r--lib/Numeric/LinearAlgebra/Tests.hs14
-rw-r--r--lib/Numeric/LinearAlgebra/Tests/Properties.hs4
11 files changed, 92 insertions, 164 deletions
diff --git a/lib/Data/Packed/Convert.hs b/lib/Data/Packed/Convert.hs
deleted file mode 100644
index 8f0eef6..0000000
--- a/lib/Data/Packed/Convert.hs
+++ /dev/null
@@ -1,101 +0,0 @@
1{-# OPTIONS -XTypeOperators -XRank2Types -XFlexibleContexts #-}
2
3-----------------------------------------------------------------------------
4-- |
5-- Module : Data.Packed.Convert
6-- Copyright : (c) Alberto Ruiz 2008
7-- License : GPL-style
8--
9-- Maintainer : Alberto Ruiz <aruiz@um.es>
10-- Stability : provisional
11-- Portability : portable
12--
13-- Conversion of Vectors and Matrices to and from the standard Haskell arrays.
14-- (provisional)
15--
16-----------------------------------------------------------------------------
17
18module Data.Packed.Convert (
19 arrayFromVector, vectorFromArray,
20 mArrayFromVector, vectorFromMArray,
21 vectorToStorableArray, storableArrayToVector,
22 arrayFromMatrix, matrixFromArray,
23 mArrayFromMatrix, matrixFromMArray,
24-- matrixToStorableArray, storableArrayToMatrix
25) where
26
27import Data.Packed.Internal
28import Data.Array.Storable
29import Foreign
30import Control.Monad.ST
31import Data.Array.ST
32import Data.Array.Unboxed
33
34-- | Creates a StorableArray indexed from 0 to dim -1.
35-- (Memory is efficiently copied, so you can then freely modify the obtained array)
36vectorToStorableArray :: Storable t => Vector t -> IO (StorableArray Int t)
37vectorToStorableArray v = do
38 r <- cloneVector v
39 unsafeForeignPtrToStorableArray (fptr r) (0,dim r -1)
40
41-- | Creates a Vector from a StorableArray.
42-- (Memory is efficiently copied, so posterior changes in the array will not affect the result)
43storableArrayToVector :: Storable t => StorableArray Int t -> IO (Vector t)
44storableArrayToVector s = do
45 (a,b) <- getBounds s
46 let n = (b-a+1)
47 r <- createVector n
48 withStorableArray s $ \p -> do
49 let f _ d = copyArray d p n >> return 0
50 app1 f vec r "storableArrayToVector"
51 return r
52
53
54unsafeVectorToStorableArray :: Storable t => Vector t -> IO (StorableArray Int t)
55unsafeVectorToStorableArray v = unsafeForeignPtrToStorableArray (fptr v) (0,dim v -1)
56
57--unsafeStorableArrayToVector :: Storable t => StorableArray Int t -> IO (Vector t)
58--unsafeStorableArrayToVector s = undefined -- the foreign ptr of Storable Array is not available?
59
60-----------------------------------------------------------------
61-- provisional, we need Unboxed arrays for Complex Double
62
63
64unsafeFreeze' :: (MArray a e m, Ix i) => a i e -> m (Array i e)
65unsafeFreeze' = unsafeFreeze
66
67-- | creates an immutable Array from an hmatrix Vector (to do: unboxed)
68arrayFromVector :: (Storable t) => Vector t -> Array Int t
69arrayFromVector x = runSTArray (mArrayFromVector x)
70
71-- | creates a mutable array from an hmatrix Vector (to do: unboxed)
72mArrayFromVector :: (MArray b t (ST s), Storable t) => Vector t -> ST s (b Int t)
73mArrayFromVector v = unsafeThaw =<< unsafeIOToST ( unsafeFreeze' =<< (vectorToStorableArray $ v))
74
75
76-- (creates an hmatrix Vector from an immutable array (to do: unboxed))
77vectorFromArray :: (Storable t) => Array Int t -> Vector t
78vectorFromArray a = unsafePerformIO $ storableArrayToVector =<< unsafeThaw a
79
80-- | creates a mutable Array from an hmatrix Vector for manipulation with runSTUArray (to do: unboxed)
81vectorFromMArray :: (Storable t, MArray a t (ST s)) => a Int t -> ST s (Vector t)
82vectorFromMArray x = fmap vectorFromArray (unsafeFreeze' x)
83
84--------------------------------------------------------------------
85-- provisional
86
87matrixFromArray :: UArray (Int, Int) Double -> Matrix Double
88matrixFromArray m = reshape c . fromList . elems $ m
89 where ((r1,c1),(r2,c2)) = bounds m
90 _r = r2-r1+1
91 c = c2-c1+1
92
93arrayFromMatrix :: Matrix Double -> UArray (Int, Int) Double
94arrayFromMatrix m = listArray ((0,0),(rows m -1, cols m -1)) (toList $ flatten m)
95
96
97mArrayFromMatrix :: (MArray b Double m) => Matrix Double -> m (b (Int, Int) Double)
98mArrayFromMatrix = unsafeThaw . arrayFromMatrix
99
100matrixFromMArray :: (MArray a Double (ST s)) => a (Int,Int) Double -> ST s (Matrix Double)
101matrixFromMArray x = fmap matrixFromArray (unsafeFreeze x)
diff --git a/lib/Data/Packed/Internal/Matrix.hs b/lib/Data/Packed/Internal/Matrix.hs
index a6b6215..b8e7ab0 100644
--- a/lib/Data/Packed/Internal/Matrix.hs
+++ b/lib/Data/Packed/Internal/Matrix.hs
@@ -18,7 +18,7 @@
18module Data.Packed.Internal.Matrix( 18module Data.Packed.Internal.Matrix(
19 Matrix(..), rows, cols, 19 Matrix(..), rows, cols,
20 MatrixOrder(..), orderOf, 20 MatrixOrder(..), orderOf,
21 createMatrix, withMatrix, mat, 21 createMatrix, mat,
22 cmat, fmat, 22 cmat, fmat,
23 toLists, flatten, reshape, 23 toLists, flatten, reshape,
24 Element(..), 24 Element(..),
@@ -115,11 +115,11 @@ fmat m@MF{} = m
115fmat MC {irows = r, icols = c, cdat = d } = MF {irows = r, icols = c, fdat = transdata c d r} 115fmat MC {irows = r, icols = c, cdat = d } = MF {irows = r, icols = c, fdat = transdata c d r}
116 116
117-- C-Haskell matrix adapter 117-- C-Haskell matrix adapter
118mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r 118-- mat :: Adapt (CInt -> CInt -> Ptr t -> r) (Matrix t) r
119mat = withMatrix
120 119
121withMatrix a f = 120mat :: (Storable t) => Matrix t -> (((CInt -> CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b
122 withForeignPtr (fptr (xdat a)) $ \p -> do 121mat a f =
122 unsafeWith (xdat a) $ \p -> do
123 let m g = do 123 let m g = do
124 g (fi (rows a)) (fi (cols a)) p 124 g (fi (rows a)) (fi (cols a)) p
125 f m 125 f m
@@ -273,8 +273,8 @@ transdata' c1 v c2 =
273 then v 273 then v
274 else unsafePerformIO $ do 274 else unsafePerformIO $ do
275 w <- createVector (r2*c2) 275 w <- createVector (r2*c2)
276 withForeignPtr (fptr v) $ \p -> 276 unsafeWith v $ \p ->
277 withForeignPtr (fptr w) $ \q -> do 277 unsafeWith w $ \q -> do
278 let go (-1) _ = return () 278 let go (-1) _ = return ()
279 go !i (-1) = go (i-1) (c1-1) 279 go !i (-1) = go (i-1) (c1-1)
280 go !i !j = do x <- peekElemOff p (i*c1+j) 280 go !i !j = do x <- peekElemOff p (i*c1+j)
@@ -300,8 +300,8 @@ transdataAux fun c1 d c2 =
300 then d 300 then d
301 else unsafePerformIO $ do 301 else unsafePerformIO $ do
302 v <- createVector (dim d) 302 v <- createVector (dim d)
303 withForeignPtr (fptr d) $ \pd -> 303 unsafeWith d $ \pd ->
304 withForeignPtr (fptr v) $ \pv -> 304 unsafeWith v $ \pv ->
305 fun (fi r1) (fi c1) pd (fi r2) (fi c2) pv // check "transdataAux" 305 fun (fi r1) (fi c1) pd (fi r2) (fi c2) pv // check "transdataAux"
306 return v 306 return v
307 where r1 = dim d `div` c1 307 where r1 = dim d `div` c1
@@ -314,7 +314,7 @@ foreign import ccall "transC" ctransC :: TCMCM
314 314
315constant' v n = unsafePerformIO $ do 315constant' v n = unsafePerformIO $ do
316 w <- createVector n 316 w <- createVector n
317 withForeignPtr (fptr w) $ \p -> do 317 unsafeWith w $ \p -> do
318 let go (-1) = return () 318 let go (-1) = return ()
319 go !k = pokeElemOff p k v >> go (k-1) 319 go !k = pokeElemOff p k v >> go (k-1)
320 go (n-1) 320 go (n-1)
@@ -352,8 +352,8 @@ subMatrix (r0,c0) (rt,ct) m
352 352
353subMatrix'' (r0,c0) (rt,ct) c v = unsafePerformIO $ do 353subMatrix'' (r0,c0) (rt,ct) c v = unsafePerformIO $ do
354 w <- createVector (rt*ct) 354 w <- createVector (rt*ct)
355 withForeignPtr (fptr v) $ \p -> 355 unsafeWith v $ \p ->
356 withForeignPtr (fptr w) $ \q -> do 356 unsafeWith w $ \q -> do
357 let go (-1) _ = return () 357 let go (-1) _ = return ()
358 go !i (-1) = go (i-1) (ct-1) 358 go !i (-1) = go (i-1) (ct-1)
359 go !i !j = do x <- peekElemOff p ((i+r0)*c+j+c0) 359 go !i !j = do x <- peekElemOff p ((i+r0)*c+j+c0)
diff --git a/lib/Data/Packed/Internal/Vector.hs b/lib/Data/Packed/Internal/Vector.hs
index 51c0d92..6bbd207 100644
--- a/lib/Data/Packed/Internal/Vector.hs
+++ b/lib/Data/Packed/Internal/Vector.hs
@@ -12,20 +12,20 @@
12-- Vector implementation 12-- Vector implementation
13-- 13--
14----------------------------------------------------------------------------- 14-----------------------------------------------------------------------------
15-- #hide
16 15
17module Data.Packed.Internal.Vector ( 16module Data.Packed.Internal.Vector (
18 Vector(..), dim, 17 Vector, dim,
19 fromList, toList, (|>), 18 fromList, toList, (|>),
20 join, (@>), safe, at, at', subVector, 19 join, (@>), safe, at, at', subVector,
21 mapVector, zipVector, 20 mapVector, zipVector,
22 foldVector, foldVectorG, foldLoop, 21 foldVector, foldVectorG, foldLoop,
23 createVector, withVector, vec, 22 createVector, vec,
24 asComplex, asReal, 23 asComplex, asReal,
25 fwriteVector, freadVector, fprintfVector, fscanfVector, 24 fwriteVector, freadVector, fprintfVector, fscanfVector,
26 cloneVector, 25 cloneVector,
27 unsafeToForeignPtr, 26 unsafeToForeignPtr,
28 unsafeFromForeignPtr 27 unsafeFromForeignPtr,
28 unsafeWith
29) where 29) where
30 30
31import Data.Packed.Internal.Common 31import Data.Packed.Internal.Common
@@ -47,39 +47,45 @@ import GHC.Base
47import GHC.IOBase 47import GHC.IOBase
48#endif 48#endif
49 49
50-- | A one-dimensional array of objects stored in a contiguous memory block. 50-- | One-dimensional array of objects stored in a contiguous memory block.
51data Vector t = 51data Vector t =
52 V { idim :: {-# UNPACK #-} !Int -- ^ number of elements 52 V { ioff :: {-# UNPACK #-} !Int -- ^ offset of first element
53 , fptr :: {-# UNPACK #-} !(ForeignPtr t) -- ^ foreign pointer to the memory block 53 , idim :: {-# UNPACK #-} !Int -- ^ number of elements
54 , fptr :: {-# UNPACK #-} !(ForeignPtr t) -- ^ foreign pointer to the memory block
54 } 55 }
55 56
56unsafeToForeignPtr :: Vector a -> (ForeignPtr a, Int, Int) 57unsafeToForeignPtr :: Vector a -> (ForeignPtr a, Int, Int)
57unsafeToForeignPtr v = (fptr v, 0, idim v) 58unsafeToForeignPtr v = (fptr v, ioff v, idim v)
58 59
59-- | Same convention as in Roman Leshchinskiy's vector package. 60-- | Same convention as in Roman Leshchinskiy's vector package.
60unsafeFromForeignPtr :: ForeignPtr a -> Int -> Int -> Vector a 61unsafeFromForeignPtr :: ForeignPtr a -> Int -> Int -> Vector a
61unsafeFromForeignPtr fp i n | i == 0 = V {idim = n, fptr = fp} 62unsafeFromForeignPtr fp i n | n > 0 = V {ioff = i, idim = n, fptr = fp}
62 | otherwise = error "unsafeFromForeignPtr with nonzero offset" 63 | otherwise = error "unsafeFromForeignPtr with dim < 1"
64
65unsafeWith (V i _ fp) m = withForeignPtr fp $ \p -> m (p `advancePtr` i)
66{-# INLINE unsafeWith #-}
67
63 68
64-- | Number of elements 69-- | Number of elements
65dim :: Vector t -> Int 70dim :: Vector t -> Int
66dim = idim 71dim = idim
67 72
68-- C-Haskell vector adapter 73-- C-Haskell vector adapter
69vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r 74-- vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r
70vec = withVector 75vec :: (Storable t) => Vector t -> (((CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b
71 76vec x f = unsafeWith x $ \p -> do
72withVector (V n fp) f = withForeignPtr fp $ \p -> do
73 let v g = do 77 let v g = do
74 g (fi n) p 78 g (fi $ dim x) p
75 f v 79 f v
80{-# INLINE vec #-}
81
76 82
77-- allocates memory for a new vector 83-- allocates memory for a new vector
78createVector :: Storable a => Int -> IO (Vector a) 84createVector :: Storable a => Int -> IO (Vector a)
79createVector n = do 85createVector n = do
80 when (n <= 0) $ error ("trying to createVector of dim "++show n) 86 when (n <= 0) $ error ("trying to createVector of dim "++show n)
81 fp <- doMalloc undefined 87 fp <- doMalloc undefined
82 return $ V n fp 88 return $ V 0 n fp
83 where 89 where
84 -- 90 --
85 -- Use the much cheaper Haskell heap allocated storage 91 -- Use the much cheaper Haskell heap allocated storage
@@ -102,11 +108,10 @@ createVector n = do
102fromList :: Storable a => [a] -> Vector a 108fromList :: Storable a => [a] -> Vector a
103fromList l = unsafePerformIO $ do 109fromList l = unsafePerformIO $ do
104 v <- createVector (length l) 110 v <- createVector (length l)
105 let f _ p = pokeArray p l >> return 0 111 unsafeWith v $ \ p -> pokeArray p l
106 app1 f vec v "fromList"
107 return v 112 return v
108 113
109safeRead v = inlinePerformIO . withForeignPtr (fptr v) 114safeRead v = inlinePerformIO . unsafeWith v
110{-# INLINE safeRead #-} 115{-# INLINE safeRead #-}
111 116
112inlinePerformIO :: IO a -> a 117inlinePerformIO :: IO a -> a
@@ -171,7 +176,12 @@ subVector :: Storable t => Int -- ^ index of the starting element
171 -> Int -- ^ number of elements to extract 176 -> Int -- ^ number of elements to extract
172 -> Vector t -- ^ source 177 -> Vector t -- ^ source
173 -> Vector t -- ^ result 178 -> Vector t -- ^ result
174subVector k l (v@V {idim=n}) 179
180subVector k l v@V{idim = n, ioff = i}
181 | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range"
182 | otherwise = v {idim = l, ioff = i+k}
183
184subVectorCopy k l (v@V {idim=n})
175 | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range" 185 | k<0 || k >= n || k+l > n || l < 0 = error "subVector out of range"
176 | otherwise = unsafePerformIO $ do 186 | otherwise = unsafePerformIO $ do
177 r <- createVector l 187 r <- createVector l
@@ -201,23 +211,24 @@ join [] = error "joining zero vectors"
201join [v] = v 211join [v] = v
202join as = unsafePerformIO $ do 212join as = unsafePerformIO $ do
203 let tot = sum (map dim as) 213 let tot = sum (map dim as)
204 r@V {fptr = p} <- createVector tot 214 r <- createVector tot
205 withForeignPtr p $ \ptr -> 215 unsafeWith r $ \ptr ->
206 joiner as tot ptr 216 joiner as tot ptr
207 return r 217 return r
208 where joiner [] _ _ = return () 218 where joiner [] _ _ = return ()
209 joiner (V {idim = n, fptr = b} : cs) _ p = do 219 joiner (v:cs) _ p = do
210 withForeignPtr b $ \pb -> copyArray p pb n 220 let n = dim v
221 unsafeWith v $ \pb -> copyArray p pb n
211 joiner cs 0 (advancePtr p n) 222 joiner cs 0 (advancePtr p n)
212 223
213 224
214-- | transforms a complex vector into a real vector with alternating real and imaginary parts 225-- | transforms a complex vector into a real vector with alternating real and imaginary parts
215asReal :: Vector (Complex Double) -> Vector Double 226asReal :: Vector (Complex Double) -> Vector Double
216asReal v = V { idim = 2*dim v, fptr = castForeignPtr (fptr v) } 227asReal v = V { ioff = 2*ioff v, idim = 2*dim v, fptr = castForeignPtr (fptr v) }
217 228
218-- | transforms a real vector into a complex vector with alternating real and imaginary parts 229-- | transforms a real vector into a complex vector with alternating real and imaginary parts
219asComplex :: Vector Double -> Vector (Complex Double) 230asComplex :: Vector Double -> Vector (Complex Double)
220asComplex v = V { idim = dim v `div` 2, fptr = castForeignPtr (fptr v) } 231asComplex v = V { ioff = ioff v `div` 2, idim = dim v `div` 2, fptr = castForeignPtr (fptr v) }
221 232
222---------------------------------------------------------------- 233----------------------------------------------------------------
223 234
@@ -234,8 +245,8 @@ cloneVector (v@V {idim=n}) = do
234mapVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b 245mapVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b
235mapVector f v = unsafePerformIO $ do 246mapVector f v = unsafePerformIO $ do
236 w <- createVector (dim v) 247 w <- createVector (dim v)
237 withForeignPtr (fptr v) $ \p -> 248 unsafeWith v $ \p ->
238 withForeignPtr (fptr w) $ \q -> do 249 unsafeWith w $ \q -> do
239 let go (-1) = return () 250 let go (-1) = return ()
240 go !k = do x <- peekElemOff p k 251 go !k = do x <- peekElemOff p k
241 pokeElemOff q k (f x) 252 pokeElemOff q k (f x)
@@ -249,9 +260,9 @@ zipVector :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a ->
249zipVector f u v = unsafePerformIO $ do 260zipVector f u v = unsafePerformIO $ do
250 let n = min (dim u) (dim v) 261 let n = min (dim u) (dim v)
251 w <- createVector n 262 w <- createVector n
252 withForeignPtr (fptr u) $ \pu -> 263 unsafeWith u $ \pu ->
253 withForeignPtr (fptr v) $ \pv -> 264 unsafeWith v $ \pv ->
254 withForeignPtr (fptr w) $ \pw -> do 265 unsafeWith w $ \pw -> do
255 let go (-1) = return () 266 let go (-1) = return ()
256 go !k = do x <- peekElemOff pu k 267 go !k = do x <- peekElemOff pu k
257 y <- peekElemOff pv k 268 y <- peekElemOff pv k
@@ -262,7 +273,7 @@ zipVector f u v = unsafePerformIO $ do
262{-# INLINE zipVector #-} 273{-# INLINE zipVector #-}
263 274
264foldVector f x v = unsafePerformIO $ 275foldVector f x v = unsafePerformIO $
265 withForeignPtr (fptr (v::Vector Double)) $ \p -> do 276 unsafeWith (v::Vector Double) $ \p -> do
266 let go (-1) s = return s 277 let go (-1) s = return s
267 go !k !s = do y <- peekElemOff p k 278 go !k !s = do y <- peekElemOff p k
268 go (k-1::Int) (f y s) 279 go (k-1::Int) (f y s)
diff --git a/lib/Data/Packed/ST.hs b/lib/Data/Packed/ST.hs
index a47da05..48e35b4 100644
--- a/lib/Data/Packed/ST.hs
+++ b/lib/Data/Packed/ST.hs
@@ -37,11 +37,11 @@ import Foreign
37 37
38{-# INLINE ioReadV #-} 38{-# INLINE ioReadV #-}
39ioReadV :: Storable t => Vector t -> Int -> IO t 39ioReadV :: Storable t => Vector t -> Int -> IO t
40ioReadV v k = withForeignPtr (fptr v) $ \s -> peekElemOff s k 40ioReadV v k = unsafeWith v $ \s -> peekElemOff s k
41 41
42{-# INLINE ioWriteV #-} 42{-# INLINE ioWriteV #-}
43ioWriteV :: Storable t => Vector t -> Int -> t -> IO () 43ioWriteV :: Storable t => Vector t -> Int -> t -> IO ()
44ioWriteV v k x = withForeignPtr (fptr v) $ \s -> pokeElemOff s k x 44ioWriteV v k x = unsafeWith v $ \s -> pokeElemOff s k x
45 45
46newtype STVector s t = STVector (Vector t) 46newtype STVector s t = STVector (Vector t)
47 47
diff --git a/lib/Numeric/GSL/Internal.hs b/lib/Numeric/GSL/Internal.hs
index 37bcc1b..2ccc3ef 100644
--- a/lib/Numeric/GSL/Internal.hs
+++ b/lib/Numeric/GSL/Internal.hs
@@ -35,12 +35,12 @@ foreign import ccall "wrapper"
35 35
36aux_vTov :: (Vector Double -> Vector Double) -> TVV 36aux_vTov :: (Vector Double -> Vector Double) -> TVV
37aux_vTov f n p nr r = g where 37aux_vTov f n p nr r = g where
38 V {fptr = pr} = f x 38 v = f x
39 x = createV (fromIntegral n) copy "aux_vTov" 39 x = createV (fromIntegral n) copy "aux_vTov"
40 copy n' q = do 40 copy n' q = do
41 copyArray q p (fromIntegral n') 41 copyArray q p (fromIntegral n')
42 return 0 42 return 0
43 g = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral nr) 43 g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral nr)
44 return 0 44 return 0
45 45
46foreign import ccall "wrapper" 46foreign import ccall "wrapper"
@@ -51,12 +51,12 @@ foreign import ccall "wrapper"
51 51
52aux_vTom :: (Vector Double -> Matrix Double) -> TVM 52aux_vTom :: (Vector Double -> Matrix Double) -> TVM
53aux_vTom f n p rr cr r = g where 53aux_vTom f n p rr cr r = g where
54 V {fptr = pr} = flatten $ f x 54 v = flatten $ f x
55 x = createV (fromIntegral n) copy "aux_vTov" 55 x = createV (fromIntegral n) copy "aux_vTov"
56 copy n' q = do 56 copy n' q = do
57 copyArray q p (fromIntegral n') 57 copyArray q p (fromIntegral n')
58 return 0 58 return 0
59 g = do withForeignPtr pr $ \p' -> copyArray r p' (fromIntegral $ rr*cr) 59 g = do unsafeWith v $ \p' -> copyArray r p' (fromIntegral $ rr*cr)
60 return 0 60 return 0
61 61
62createV n fun msg = unsafePerformIO $ do 62createV n fun msg = unsafePerformIO $ do
diff --git a/lib/Numeric/GSL/Minimization.hs b/lib/Numeric/GSL/Minimization.hs
index 3240ddf..3bbb2b3 100644
--- a/lib/Numeric/GSL/Minimization.hs
+++ b/lib/Numeric/GSL/Minimization.hs
@@ -109,7 +109,7 @@ ww2 w1 o1 w2 o2 f = w1 o1 $ \a1 -> w2 o2 $ \a2 -> f a1 a2
109minimizeV method eps maxit szv f xiv = unsafePerformIO $ do 109minimizeV method eps maxit szv f xiv = unsafePerformIO $ do
110 let n = dim xiv 110 let n = dim xiv
111 fp <- mkVecfun (iv f) 111 fp <- mkVecfun (iv f)
112 rawpath <- ww2 withVector xiv withVector szv $ \xiv' szv' -> 112 rawpath <- ww2 vec xiv vec szv $ \xiv' szv' ->
113 createMIO maxit (n+3) 113 createMIO maxit (n+3)
114 (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv') 114 (c_minimize (fi (fromEnum method)) fp eps (fi maxit) // xiv' // szv')
115 "minimize" 115 "minimize"
@@ -166,7 +166,7 @@ minimizeVD method eps maxit istep tol f df xiv = unsafePerformIO $ do
166 df' = (checkdim1 n . df) 166 df' = (checkdim1 n . df)
167 fp <- mkVecfun (iv f') 167 fp <- mkVecfun (iv f')
168 dfp <- mkVecVecfun (aux_vTov df') 168 dfp <- mkVecVecfun (aux_vTov df')
169 rawpath <- withVector xiv $ \xiv' -> 169 rawpath <- vec xiv $ \xiv' ->
170 createMIO maxit (n+2) 170 createMIO maxit (n+2)
171 (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv') 171 (c_minimizeD (fi (fromEnum method)) fp dfp istep tol eps (fi maxit) // xiv')
172 "minimizeD" 172 "minimizeD"
diff --git a/lib/Numeric/GSL/ODE.hs b/lib/Numeric/GSL/ODE.hs
index 9d890fc..f6f11f9 100644
--- a/lib/Numeric/GSL/ODE.hs
+++ b/lib/Numeric/GSL/ODE.hs
@@ -83,8 +83,8 @@ odeSolveV method h epsAbs epsRel f mbjac xiv ts = unsafePerformIO $ do
83 jp <- case mbjac of 83 jp <- case mbjac of
84 Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t)) 84 Just jac -> mkDoubleVecMatfun (\t -> aux_vTom (checkdim2 n . jac t))
85 Nothing -> return nullFunPtr 85 Nothing -> return nullFunPtr
86 sol <- withVector xiv $ \xiv' -> 86 sol <- vec xiv $ \xiv' ->
87 withVector (checkTimes ts) $ \ts' -> 87 vec (checkTimes ts) $ \ts' ->
88 createMIO (dim ts) n 88 createMIO (dim ts) n
89 (ode_c (fi (fromEnum method)) h epsAbs epsRel fp jp // xiv' // ts' ) 89 (ode_c (fi (fromEnum method)) h epsAbs epsRel fp jp // xiv' // ts' )
90 "ode" 90 "ode"
diff --git a/lib/Numeric/GSL/Root.hs b/lib/Numeric/GSL/Root.hs
index 840e8ee..ddb090d 100644
--- a/lib/Numeric/GSL/Root.hs
+++ b/lib/Numeric/GSL/Root.hs
@@ -79,7 +79,7 @@ rootGen m f xi epsabs maxit = unsafePerformIO $ do
79 let xiv = fromList xi 79 let xiv = fromList xi
80 n = dim xiv 80 n = dim xiv
81 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList)) 81 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
82 rawpath <- withVector xiv $ \xiv' -> 82 rawpath <- vec xiv $ \xiv' ->
83 createMIO maxit (2*n+1) 83 createMIO maxit (2*n+1)
84 (c_root m fp epsabs (fi maxit) // xiv') 84 (c_root m fp epsabs (fi maxit) // xiv')
85 "root" 85 "root"
@@ -117,7 +117,7 @@ rootJGen m f jac xi epsabs maxit = unsafePerformIO $ do
117 n = dim xiv 117 n = dim xiv
118 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList)) 118 fp <- mkVecVecfun (aux_vTov (checkdim1 n . fromList . f . toList))
119 jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList)) 119 jp <- mkVecMatfun (aux_vTom (checkdim2 n . fromLists . jac . toList))
120 rawpath <- withVector xiv $ \xiv' -> 120 rawpath <- vec xiv $ \xiv' ->
121 createMIO maxit (2*n+1) 121 createMIO maxit (2*n+1)
122 (c_rootj m fp jp epsabs (fi maxit) // xiv') 122 (c_rootj m fp jp epsabs (fi maxit) // xiv')
123 "root" 123 "root"
diff --git a/lib/Numeric/LinearAlgebra/Instances.hs b/lib/Numeric/LinearAlgebra/Instances.hs
index 67496f2..b17da5a 100644
--- a/lib/Numeric/LinearAlgebra/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Instances.hs
@@ -25,7 +25,7 @@ import Data.Packed.Matrix
25import Complex 25import Complex
26import Data.List(transpose,intersperse) 26import Data.List(transpose,intersperse)
27import Foreign(Storable) 27import Foreign(Storable)
28import Data.Monoid 28-- import Data.Monoid
29import Data.Packed.Internal.Vector 29import Data.Packed.Internal.Vector
30-- import Control.Parallel.Strategies 30-- import Control.Parallel.Strategies
31 31
@@ -182,12 +182,12 @@ instance (Linear Vector a, Floating (Vector a), Fractional (Matrix a)) => Floati
182 182
183--------------------------------------------------------------- 183---------------------------------------------------------------
184 184
185instance (Storable a, Num (Vector a)) => Monoid (Vector a) where 185-- instance (Storable a, Num (Vector a)) => Monoid (Vector a) where
186 mempty = 0 { idim = 0 } 186-- mempty = 0 { idim = 0 }
187 mappend a b = mconcat [a,b] 187-- mappend a b = mconcat [a,b]
188 mconcat = j . filter ((>0).dim) 188-- mconcat = j . filter ((>0).dim)
189 where j [] = mempty 189-- where j [] = mempty
190 j l = join l 190-- j l = join l
191 191
192--------------------------------------------------------------- 192---------------------------------------------------------------
193 193
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index 016b9a1..3f7c847 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -180,6 +180,9 @@ runTests n = do
180 test (multProp1 . cConsist) 180 test (multProp1 . cConsist)
181 test (multProp2 . rConsist) 181 test (multProp2 . rConsist)
182 test (multProp2 . cConsist) 182 test (multProp2 . cConsist)
183 putStrLn "------ sub-trans"
184 test (subProp . rM)
185 test (subProp . cM)
183 putStrLn "------ lu" 186 putStrLn "------ lu"
184 test (luProp . rM) 187 test (luProp . rM)
185 test (luProp . cM) 188 test (luProp . cM)
@@ -305,6 +308,7 @@ makeUnitary v | realPart n > 1 = v / scalar n
305-- | Performance measurements. 308-- | Performance measurements.
306runBenchmarks :: IO () 309runBenchmarks :: IO ()
307runBenchmarks = do 310runBenchmarks = do
311 subBench
308 multBench 312 multBench
309 svdBench 313 svdBench
310 eigBench 314 eigBench
@@ -335,6 +339,16 @@ multb n = foldl1' (<>) (replicate (10^6) (ident n :: Matrix Double))
335 339
336-------------------------------- 340--------------------------------
337 341
342subBench = do
343 putStrLn ""
344 let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (dim v -1) v))
345 time "0.1M subVector " (g (constant 1 (1+10^5) :: Vector Double) @> 0)
346 let f = foldl1' (.) (replicate (10^5) (fromRows.toRows))
347 time "subVector-join 3" (f (ident 3 :: Matrix Double) @@>(0,0))
348 time "subVector-join 10" (f (ident 10 :: Matrix Double) @@>(0,0))
349
350--------------------------------
351
338multBench = do 352multBench = do
339 let a = ident 1000 :: Matrix Double 353 let a = ident 1000 :: Matrix Double
340 let b = ident 2000 :: Matrix Double 354 let b = ident 2000 :: Matrix Double
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
index 93175f2..6d176fa 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -38,6 +38,7 @@ module Numeric.LinearAlgebra.Tests.Properties (
38 cholProp, 38 cholProp,
39 expmDiagProp, 39 expmDiagProp,
40 multProp1, multProp2, 40 multProp1, multProp2,
41 subProp,
41 linearSolveProp, linearSolveProp2 42 linearSolveProp, linearSolveProp2
42) where 43) where
43 44
@@ -243,3 +244,6 @@ linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)
243 where q = min (rows a) (cols a) 244 where q = min (rows a) (cols a)
244 b = a <> x 245 b = a <> x
245 wc = rank a == q 246 wc = rank a == q
247
248subProp m = m == (trans . fromColumns . toRows) m
249
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-09 13:35:04 +0000