diff options
Diffstat (limited to 'packages/base/src/Data/Packed/Internal/Vector.hs')
-rw-r--r-- | packages/base/src/Data/Packed/Internal/Vector.hs | 471 |
1 files changed, 471 insertions, 0 deletions
diff --git a/packages/base/src/Data/Packed/Internal/Vector.hs b/packages/base/src/Data/Packed/Internal/Vector.hs new file mode 100644 index 0000000..d0bc143 --- /dev/null +++ b/packages/base/src/Data/Packed/Internal/Vector.hs | |||
@@ -0,0 +1,471 @@ | |||
1 | {-# LANGUAGE MagicHash, CPP, UnboxedTuples, BangPatterns, FlexibleContexts #-} | ||
2 | -- | | ||
3 | -- Module : Data.Packed.Internal.Vector | ||
4 | -- Copyright : (c) Alberto Ruiz 2007 | ||
5 | -- License : BSD3 | ||
6 | -- Maintainer : Alberto Ruiz | ||
7 | -- Stability : provisional | ||
8 | -- | ||
9 | -- Vector implementation | ||
10 | -- | ||
11 | -------------------------------------------------------------------------------- | ||
12 | |||
13 | module Data.Packed.Internal.Vector ( | ||
14 | Vector, dim, | ||
15 | fromList, toList, (|>), | ||
16 | vjoin, (@>), safe, at, at', subVector, takesV, | ||
17 | mapVector, mapVectorWithIndex, zipVectorWith, unzipVectorWith, | ||
18 | mapVectorM, mapVectorM_, mapVectorWithIndexM, mapVectorWithIndexM_, | ||
19 | foldVector, foldVectorG, foldLoop, foldVectorWithIndex, | ||
20 | createVector, vec, | ||
21 | asComplex, asReal, float2DoubleV, double2FloatV, | ||
22 | stepF, stepD, condF, condD, | ||
23 | conjugateQ, conjugateC, | ||
24 | cloneVector, | ||
25 | unsafeToForeignPtr, | ||
26 | unsafeFromForeignPtr, | ||
27 | unsafeWith | ||
28 | ) where | ||
29 | |||
30 | import Data.Packed.Internal.Common | ||
31 | import Data.Packed.Internal.Signatures | ||
32 | import Foreign.Marshal.Array(peekArray, copyArray, advancePtr) | ||
33 | import Foreign.ForeignPtr(ForeignPtr, castForeignPtr) | ||
34 | import Foreign.Ptr(Ptr) | ||
35 | import Foreign.Storable(Storable, peekElemOff, pokeElemOff, sizeOf) | ||
36 | import Foreign.C.Types | ||
37 | import Data.Complex | ||
38 | import Control.Monad(when) | ||
39 | import System.IO.Unsafe(unsafePerformIO) | ||
40 | |||
41 | #if __GLASGOW_HASKELL__ >= 605 | ||
42 | import GHC.ForeignPtr (mallocPlainForeignPtrBytes) | ||
43 | #else | ||
44 | import Foreign.ForeignPtr (mallocForeignPtrBytes) | ||
45 | #endif | ||
46 | |||
47 | import GHC.Base | ||
48 | #if __GLASGOW_HASKELL__ < 612 | ||
49 | import GHC.IOBase hiding (liftIO) | ||
50 | #endif | ||
51 | |||
52 | import qualified Data.Vector.Storable as Vector | ||
53 | import Data.Vector.Storable(Vector, | ||
54 | fromList, | ||
55 | unsafeToForeignPtr, | ||
56 | unsafeFromForeignPtr, | ||
57 | unsafeWith) | ||
58 | |||
59 | |||
60 | -- | Number of elements | ||
61 | dim :: (Storable t) => Vector t -> Int | ||
62 | dim = Vector.length | ||
63 | |||
64 | |||
65 | -- C-Haskell vector adapter | ||
66 | -- vec :: Adapt (CInt -> Ptr t -> r) (Vector t) r | ||
67 | vec :: (Storable t) => Vector t -> (((CInt -> Ptr t -> t1) -> t1) -> IO b) -> IO b | ||
68 | vec x f = unsafeWith x $ \p -> do | ||
69 | let v g = do | ||
70 | g (fi $ dim x) p | ||
71 | f v | ||
72 | {-# INLINE vec #-} | ||
73 | |||
74 | |||
75 | -- allocates memory for a new vector | ||
76 | createVector :: Storable a => Int -> IO (Vector a) | ||
77 | createVector n = do | ||
78 | when (n < 0) $ error ("trying to createVector of negative dim: "++show n) | ||
79 | fp <- doMalloc undefined | ||
80 | return $ unsafeFromForeignPtr fp 0 n | ||
81 | where | ||
82 | -- | ||
83 | -- Use the much cheaper Haskell heap allocated storage | ||
84 | -- for foreign pointer space we control | ||
85 | -- | ||
86 | doMalloc :: Storable b => b -> IO (ForeignPtr b) | ||
87 | doMalloc dummy = do | ||
88 | #if __GLASGOW_HASKELL__ >= 605 | ||
89 | mallocPlainForeignPtrBytes (n * sizeOf dummy) | ||
90 | #else | ||
91 | mallocForeignPtrBytes (n * sizeOf dummy) | ||
92 | #endif | ||
93 | |||
94 | {- | creates a Vector from a list: | ||
95 | |||
96 | @> fromList [2,3,5,7] | ||
97 | 4 |> [2.0,3.0,5.0,7.0]@ | ||
98 | |||
99 | -} | ||
100 | |||
101 | safeRead v = inlinePerformIO . unsafeWith v | ||
102 | {-# INLINE safeRead #-} | ||
103 | |||
104 | inlinePerformIO :: IO a -> a | ||
105 | inlinePerformIO (IO m) = case m realWorld# of (# _, r #) -> r | ||
106 | {-# INLINE inlinePerformIO #-} | ||
107 | |||
108 | {- | extracts the Vector elements to a list | ||
109 | |||
110 | >>> toList (linspace 5 (1,10)) | ||
111 | [1.0,3.25,5.5,7.75,10.0] | ||
112 | |||
113 | -} | ||
114 | toList :: Storable a => Vector a -> [a] | ||
115 | toList v = safeRead v $ peekArray (dim v) | ||
116 | |||
117 | {- | Create a vector from a list of elements and explicit dimension. The input | ||
118 | list is explicitly truncated if it is too long, so it may safely | ||
119 | be used, for instance, with infinite lists. | ||
120 | |||
121 | >>> 5 |> [1..] | ||
122 | fromList [1.0,2.0,3.0,4.0,5.0] | ||
123 | |||
124 | -} | ||
125 | (|>) :: (Storable a) => Int -> [a] -> Vector a | ||
126 | infixl 9 |> | ||
127 | n |> l = if length l' == n | ||
128 | then fromList l' | ||
129 | else error "list too short for |>" | ||
130 | where l' = take n l | ||
131 | |||
132 | |||
133 | -- | access to Vector elements without range checking | ||
134 | at' :: Storable a => Vector a -> Int -> a | ||
135 | at' v n = safeRead v $ flip peekElemOff n | ||
136 | {-# INLINE at' #-} | ||
137 | |||
138 | -- | ||
139 | -- turn off bounds checking with -funsafe at configure time. | ||
140 | -- ghc will optimise away the salways true case at compile time. | ||
141 | -- | ||
142 | #if defined(UNSAFE) | ||
143 | safe :: Bool | ||
144 | safe = False | ||
145 | #else | ||
146 | safe = True | ||
147 | #endif | ||
148 | |||
149 | -- | access to Vector elements with range checking. | ||
150 | at :: Storable a => Vector a -> Int -> a | ||
151 | at v n | ||
152 | | safe = if n >= 0 && n < dim v | ||
153 | then at' v n | ||
154 | else error "vector index out of range" | ||
155 | | otherwise = at' v n | ||
156 | {-# INLINE at #-} | ||
157 | |||
158 | {- | takes a number of consecutive elements from a Vector | ||
159 | |||
160 | >>> subVector 2 3 (fromList [1..10]) | ||
161 | fromList [3.0,4.0,5.0] | ||
162 | |||
163 | -} | ||
164 | subVector :: Storable t => Int -- ^ index of the starting element | ||
165 | -> Int -- ^ number of elements to extract | ||
166 | -> Vector t -- ^ source | ||
167 | -> Vector t -- ^ result | ||
168 | subVector = Vector.slice | ||
169 | |||
170 | |||
171 | {- | Reads a vector position: | ||
172 | |||
173 | >>> fromList [0..9] @> 7 | ||
174 | 7.0 | ||
175 | |||
176 | -} | ||
177 | (@>) :: Storable t => Vector t -> Int -> t | ||
178 | infixl 9 @> | ||
179 | (@>) = at | ||
180 | |||
181 | |||
182 | {- | concatenate a list of vectors | ||
183 | |||
184 | >>> vjoin [fromList [1..5::Double], konst 1 3] | ||
185 | fromList [1.0,2.0,3.0,4.0,5.0,1.0,1.0,1.0] | ||
186 | |||
187 | -} | ||
188 | vjoin :: Storable t => [Vector t] -> Vector t | ||
189 | vjoin [] = fromList [] | ||
190 | vjoin [v] = v | ||
191 | vjoin as = unsafePerformIO $ do | ||
192 | let tot = sum (map dim as) | ||
193 | r <- createVector tot | ||
194 | unsafeWith r $ \ptr -> | ||
195 | joiner as tot ptr | ||
196 | return r | ||
197 | where joiner [] _ _ = return () | ||
198 | joiner (v:cs) _ p = do | ||
199 | let n = dim v | ||
200 | unsafeWith v $ \pb -> copyArray p pb n | ||
201 | joiner cs 0 (advancePtr p n) | ||
202 | |||
203 | |||
204 | {- | Extract consecutive subvectors of the given sizes. | ||
205 | |||
206 | >>> takesV [3,4] (linspace 10 (1,10::Double)) | ||
207 | [fromList [1.0,2.0,3.0],fromList [4.0,5.0,6.0,7.0]] | ||
208 | |||
209 | -} | ||
210 | takesV :: Storable t => [Int] -> Vector t -> [Vector t] | ||
211 | takesV ms w | sum ms > dim w = error $ "takesV " ++ show ms ++ " on dim = " ++ (show $ dim w) | ||
212 | | otherwise = go ms w | ||
213 | where go [] _ = [] | ||
214 | go (n:ns) v = subVector 0 n v | ||
215 | : go ns (subVector n (dim v - n) v) | ||
216 | |||
217 | --------------------------------------------------------------- | ||
218 | |||
219 | -- | transforms a complex vector into a real vector with alternating real and imaginary parts | ||
220 | asReal :: (RealFloat a, Storable a) => Vector (Complex a) -> Vector a | ||
221 | asReal v = unsafeFromForeignPtr (castForeignPtr fp) (2*i) (2*n) | ||
222 | where (fp,i,n) = unsafeToForeignPtr v | ||
223 | |||
224 | -- | transforms a real vector into a complex vector with alternating real and imaginary parts | ||
225 | asComplex :: (RealFloat a, Storable a) => Vector a -> Vector (Complex a) | ||
226 | asComplex v = unsafeFromForeignPtr (castForeignPtr fp) (i `div` 2) (n `div` 2) | ||
227 | where (fp,i,n) = unsafeToForeignPtr v | ||
228 | |||
229 | --------------------------------------------------------------- | ||
230 | |||
231 | float2DoubleV :: Vector Float -> Vector Double | ||
232 | float2DoubleV v = unsafePerformIO $ do | ||
233 | r <- createVector (dim v) | ||
234 | app2 c_float2double vec v vec r "float2double" | ||
235 | return r | ||
236 | |||
237 | double2FloatV :: Vector Double -> Vector Float | ||
238 | double2FloatV v = unsafePerformIO $ do | ||
239 | r <- createVector (dim v) | ||
240 | app2 c_double2float vec v vec r "double2float2" | ||
241 | return r | ||
242 | |||
243 | |||
244 | foreign import ccall unsafe "float2double" c_float2double:: TFV | ||
245 | foreign import ccall unsafe "double2float" c_double2float:: TVF | ||
246 | |||
247 | --------------------------------------------------------------- | ||
248 | |||
249 | stepF :: Vector Float -> Vector Float | ||
250 | stepF v = unsafePerformIO $ do | ||
251 | r <- createVector (dim v) | ||
252 | app2 c_stepF vec v vec r "stepF" | ||
253 | return r | ||
254 | |||
255 | stepD :: Vector Double -> Vector Double | ||
256 | stepD v = unsafePerformIO $ do | ||
257 | r <- createVector (dim v) | ||
258 | app2 c_stepD vec v vec r "stepD" | ||
259 | return r | ||
260 | |||
261 | foreign import ccall unsafe "stepF" c_stepF :: TFF | ||
262 | foreign import ccall unsafe "stepD" c_stepD :: TVV | ||
263 | |||
264 | --------------------------------------------------------------- | ||
265 | |||
266 | condF :: Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float -> Vector Float | ||
267 | condF x y l e g = unsafePerformIO $ do | ||
268 | r <- createVector (dim x) | ||
269 | app6 c_condF vec x vec y vec l vec e vec g vec r "condF" | ||
270 | return r | ||
271 | |||
272 | condD :: Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double | ||
273 | condD x y l e g = unsafePerformIO $ do | ||
274 | r <- createVector (dim x) | ||
275 | app6 c_condD vec x vec y vec l vec e vec g vec r "condD" | ||
276 | return r | ||
277 | |||
278 | foreign import ccall unsafe "condF" c_condF :: CInt -> PF -> CInt -> PF -> CInt -> PF -> TFFF | ||
279 | foreign import ccall unsafe "condD" c_condD :: CInt -> PD -> CInt -> PD -> CInt -> PD -> TVVV | ||
280 | |||
281 | -------------------------------------------------------------------------------- | ||
282 | |||
283 | conjugateAux fun x = unsafePerformIO $ do | ||
284 | v <- createVector (dim x) | ||
285 | app2 fun vec x vec v "conjugateAux" | ||
286 | return v | ||
287 | |||
288 | conjugateQ :: Vector (Complex Float) -> Vector (Complex Float) | ||
289 | conjugateQ = conjugateAux c_conjugateQ | ||
290 | foreign import ccall unsafe "conjugateQ" c_conjugateQ :: TQVQV | ||
291 | |||
292 | conjugateC :: Vector (Complex Double) -> Vector (Complex Double) | ||
293 | conjugateC = conjugateAux c_conjugateC | ||
294 | foreign import ccall unsafe "conjugateC" c_conjugateC :: TCVCV | ||
295 | |||
296 | -------------------------------------------------------------------------------- | ||
297 | |||
298 | cloneVector :: Storable t => Vector t -> IO (Vector t) | ||
299 | cloneVector v = do | ||
300 | let n = dim v | ||
301 | r <- createVector n | ||
302 | let f _ s _ d = copyArray d s n >> return 0 | ||
303 | app2 f vec v vec r "cloneVector" | ||
304 | return r | ||
305 | |||
306 | ------------------------------------------------------------------ | ||
307 | |||
308 | -- | map on Vectors | ||
309 | mapVector :: (Storable a, Storable b) => (a-> b) -> Vector a -> Vector b | ||
310 | mapVector f v = unsafePerformIO $ do | ||
311 | w <- createVector (dim v) | ||
312 | unsafeWith v $ \p -> | ||
313 | unsafeWith w $ \q -> do | ||
314 | let go (-1) = return () | ||
315 | go !k = do x <- peekElemOff p k | ||
316 | pokeElemOff q k (f x) | ||
317 | go (k-1) | ||
318 | go (dim v -1) | ||
319 | return w | ||
320 | {-# INLINE mapVector #-} | ||
321 | |||
322 | -- | zipWith for Vectors | ||
323 | zipVectorWith :: (Storable a, Storable b, Storable c) => (a-> b -> c) -> Vector a -> Vector b -> Vector c | ||
324 | zipVectorWith f u v = unsafePerformIO $ do | ||
325 | let n = min (dim u) (dim v) | ||
326 | w <- createVector n | ||
327 | unsafeWith u $ \pu -> | ||
328 | unsafeWith v $ \pv -> | ||
329 | unsafeWith w $ \pw -> do | ||
330 | let go (-1) = return () | ||
331 | go !k = do x <- peekElemOff pu k | ||
332 | y <- peekElemOff pv k | ||
333 | pokeElemOff pw k (f x y) | ||
334 | go (k-1) | ||
335 | go (n -1) | ||
336 | return w | ||
337 | {-# INLINE zipVectorWith #-} | ||
338 | |||
339 | -- | unzipWith for Vectors | ||
340 | unzipVectorWith :: (Storable (a,b), Storable c, Storable d) | ||
341 | => ((a,b) -> (c,d)) -> Vector (a,b) -> (Vector c,Vector d) | ||
342 | unzipVectorWith f u = unsafePerformIO $ do | ||
343 | let n = dim u | ||
344 | v <- createVector n | ||
345 | w <- createVector n | ||
346 | unsafeWith u $ \pu -> | ||
347 | unsafeWith v $ \pv -> | ||
348 | unsafeWith w $ \pw -> do | ||
349 | let go (-1) = return () | ||
350 | go !k = do z <- peekElemOff pu k | ||
351 | let (x,y) = f z | ||
352 | pokeElemOff pv k x | ||
353 | pokeElemOff pw k y | ||
354 | go (k-1) | ||
355 | go (n-1) | ||
356 | return (v,w) | ||
357 | {-# INLINE unzipVectorWith #-} | ||
358 | |||
359 | foldVector :: Storable a => (a -> b -> b) -> b -> Vector a -> b | ||
360 | foldVector f x v = unsafePerformIO $ | ||
361 | unsafeWith v $ \p -> do | ||
362 | let go (-1) s = return s | ||
363 | go !k !s = do y <- peekElemOff p k | ||
364 | go (k-1::Int) (f y s) | ||
365 | go (dim v -1) x | ||
366 | {-# INLINE foldVector #-} | ||
367 | |||
368 | -- the zero-indexed index is passed to the folding function | ||
369 | foldVectorWithIndex :: Storable a => (Int -> a -> b -> b) -> b -> Vector a -> b | ||
370 | foldVectorWithIndex f x v = unsafePerformIO $ | ||
371 | unsafeWith v $ \p -> do | ||
372 | let go (-1) s = return s | ||
373 | go !k !s = do y <- peekElemOff p k | ||
374 | go (k-1::Int) (f k y s) | ||
375 | go (dim v -1) x | ||
376 | {-# INLINE foldVectorWithIndex #-} | ||
377 | |||
378 | foldLoop f s0 d = go (d - 1) s0 | ||
379 | where | ||
380 | go 0 s = f (0::Int) s | ||
381 | go !j !s = go (j - 1) (f j s) | ||
382 | |||
383 | foldVectorG f s0 v = foldLoop g s0 (dim v) | ||
384 | where g !k !s = f k (at' v) s | ||
385 | {-# INLINE g #-} -- Thanks to Ryan Ingram (http://permalink.gmane.org/gmane.comp.lang.haskell.cafe/46479) | ||
386 | {-# INLINE foldVectorG #-} | ||
387 | |||
388 | ------------------------------------------------------------------- | ||
389 | |||
390 | -- | monadic map over Vectors | ||
391 | -- the monad @m@ must be strict | ||
392 | mapVectorM :: (Storable a, Storable b, Monad m) => (a -> m b) -> Vector a -> m (Vector b) | ||
393 | mapVectorM f v = do | ||
394 | w <- return $! unsafePerformIO $! createVector (dim v) | ||
395 | mapVectorM' w 0 (dim v -1) | ||
396 | return w | ||
397 | where mapVectorM' w' !k !t | ||
398 | | k == t = do | ||
399 | x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k | ||
400 | y <- f x | ||
401 | return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y | ||
402 | | otherwise = do | ||
403 | x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k | ||
404 | y <- f x | ||
405 | _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y | ||
406 | mapVectorM' w' (k+1) t | ||
407 | {-# INLINE mapVectorM #-} | ||
408 | |||
409 | -- | monadic map over Vectors | ||
410 | mapVectorM_ :: (Storable a, Monad m) => (a -> m ()) -> Vector a -> m () | ||
411 | mapVectorM_ f v = do | ||
412 | mapVectorM' 0 (dim v -1) | ||
413 | where mapVectorM' !k !t | ||
414 | | k == t = do | ||
415 | x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k | ||
416 | f x | ||
417 | | otherwise = do | ||
418 | x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k | ||
419 | _ <- f x | ||
420 | mapVectorM' (k+1) t | ||
421 | {-# INLINE mapVectorM_ #-} | ||
422 | |||
423 | -- | monadic map over Vectors with the zero-indexed index passed to the mapping function | ||
424 | -- the monad @m@ must be strict | ||
425 | mapVectorWithIndexM :: (Storable a, Storable b, Monad m) => (Int -> a -> m b) -> Vector a -> m (Vector b) | ||
426 | mapVectorWithIndexM f v = do | ||
427 | w <- return $! unsafePerformIO $! createVector (dim v) | ||
428 | mapVectorM' w 0 (dim v -1) | ||
429 | return w | ||
430 | where mapVectorM' w' !k !t | ||
431 | | k == t = do | ||
432 | x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k | ||
433 | y <- f k x | ||
434 | return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y | ||
435 | | otherwise = do | ||
436 | x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k | ||
437 | y <- f k x | ||
438 | _ <- return $! inlinePerformIO $! unsafeWith w' $! \q -> pokeElemOff q k y | ||
439 | mapVectorM' w' (k+1) t | ||
440 | {-# INLINE mapVectorWithIndexM #-} | ||
441 | |||
442 | -- | monadic map over Vectors with the zero-indexed index passed to the mapping function | ||
443 | mapVectorWithIndexM_ :: (Storable a, Monad m) => (Int -> a -> m ()) -> Vector a -> m () | ||
444 | mapVectorWithIndexM_ f v = do | ||
445 | mapVectorM' 0 (dim v -1) | ||
446 | where mapVectorM' !k !t | ||
447 | | k == t = do | ||
448 | x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k | ||
449 | f k x | ||
450 | | otherwise = do | ||
451 | x <- return $! inlinePerformIO $! unsafeWith v $! \p -> peekElemOff p k | ||
452 | _ <- f k x | ||
453 | mapVectorM' (k+1) t | ||
454 | {-# INLINE mapVectorWithIndexM_ #-} | ||
455 | |||
456 | |||
457 | mapVectorWithIndex :: (Storable a, Storable b) => (Int -> a -> b) -> Vector a -> Vector b | ||
458 | --mapVectorWithIndex g = head . mapVectorWithIndexM (\a b -> [g a b]) | ||
459 | mapVectorWithIndex f v = unsafePerformIO $ do | ||
460 | w <- createVector (dim v) | ||
461 | unsafeWith v $ \p -> | ||
462 | unsafeWith w $ \q -> do | ||
463 | let go (-1) = return () | ||
464 | go !k = do x <- peekElemOff p k | ||
465 | pokeElemOff q k (f k x) | ||
466 | go (k-1) | ||
467 | go (dim v -1) | ||
468 | return w | ||
469 | {-# INLINE mapVectorWithIndex #-} | ||
470 | |||
471 | |||