accelerate-fft (empty) → 0.13.0.0
raw patch · 8 files changed
+801/−0 lines, 8 filesdep +acceleratedep +accelerate-cudadep +basesetup-changed
Dependencies added: accelerate, accelerate-cuda, base, cuda, cufft
Files
- Data/Array/Accelerate/Math/Complex.hs +79/−0
- Data/Array/Accelerate/Math/DFT.hs +110/−0
- Data/Array/Accelerate/Math/DFT/Centre.hs +103/−0
- Data/Array/Accelerate/Math/DFT/Roots.hs +50/−0
- Data/Array/Accelerate/Math/FFT.hs +374/−0
- LICENSE +23/−0
- Setup.hs +2/−0
- accelerate-fft.cabal +60/−0
+ Data/Array/Accelerate/Math/Complex.hs view
@@ -0,0 +1,79 @@+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE IncoherentInstances #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeSynonymInstances #-}+{-# OPTIONS -fno-warn-orphans #-}++module Data.Array.Accelerate.Math.Complex+ where++import Prelude+import Data.Function+import Data.Array.Accelerate as A++type Complex a = (a, a)++instance (Elt a, IsFloating a) => Num (Exp (Complex a)) where+ c1 + c2 = lift ( on (+) real c1 c2, on (+) imag c1 c2 )+ c1 - c2 = lift ( on (-) real c1 c2, on (-) imag c1 c2 )+ c1 * c2 = let (x, y) = unlift c1+ (x', y') = unlift c2 :: Complex (Exp a)+ in lift (x*x'-y*y', x*y'+y*x')++ negate c = lift ( negate (real c), negate (imag c) )+ abs z = lift ( magnitude z, constant 0 )+ signum z = let r = magnitude z+ (x, y) = unlift z+ in r ==* 0 ? (constant (0,0), lift (x/r, y/r))++ fromInteger n+ = lift (constant (fromInteger n), constant 0)+++instance (Elt a, IsFloating a) => Fractional (Exp (Complex a)) where+ c1 / c2+ = let (a,b) = unlift c1+ (c,d) = unlift c2 :: Complex (Exp a)+ den = c^(2 :: Int) + d^(2 :: Int)+ re = (a * c + b * d) / den+ im = (b * c - a * d) / den+ in+ lift (re, im)++ fromRational x+ = lift (constant (fromRational x), constant 0)+++-- | Non-negative magnitude of a complex number+--+magnitude :: (Elt a, IsFloating a) => Exp (Complex a) -> Exp a+magnitude c =+ let (r, i) = unlift c+ in sqrt (r*r + i*i)++-- | The phase of a complex number, in the range (-pi, pi]. If the magnitude is+-- zero, then so is the phase.+--+phase :: (Elt a, IsFloating a) => Exp (Complex a) -> Exp a+phase c =+ let (x, y) = unlift c+ in atan2 y x+++-- | Return the real part of a complex number+--+real :: Elt a => Exp (Complex a) -> Exp a+real = A.fst++-- | Return the imaginary part of a complex number+--+imag :: Elt a => Exp (Complex a) -> Exp a+imag = A.snd++-- | Return the complex conjugate of a complex number, defined as+--+-- > conj(Z) = X - iY+--+conj :: (Elt a, IsNum a) => Exp (Complex a) -> Exp (Complex a)+conj z = lift (real z, - imag z)+
+ Data/Array/Accelerate/Math/DFT.hs view
@@ -0,0 +1,110 @@+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}+-- |+-- Module : Data.Array.Accelerate.Math.DFT+-- Copyright : [2012] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Manuel M T Chakravarty <chak@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--+-- Compute the Discrete Fourier Transform (DFT) along the lower order dimension+-- of an array.+--+-- This uses a naïve algorithm which takes O(n^2) time. However, you can+-- transform an array with an arbitrary extent, unlike with FFT which requires+-- each dimension to be a power of two.+--+-- The `dft` and `idft` functions compute the roots of unity as needed. If you+-- need to transform several arrays with the same extent than it is faster to+-- compute the roots once using `rootsOfUnity` or `inverseRootsOfUnity`+-- respectively, then call `dftG` directly.+--+-- You can also compute single values of the transform using `dftGS`+--+module Data.Array.Accelerate.Math.DFT (++ dft, idft, dftG, dftGS,++) where++import Prelude as P hiding ((!!))+import Data.Array.Accelerate as A+import Data.Array.Accelerate.Math.DFT.Roots+import Data.Array.Accelerate.Math.Complex+++-- | Compute the DFT along the low order dimension of an array+--+dft :: (Shape sh, Slice sh, Elt e, IsFloating e)+ => Acc (Array (sh:.Int) (Complex e))+ -> Acc (Array (sh:.Int) (Complex e))+dft v = dftG (rootsOfUnity (shape v)) v+++-- | Compute the inverse DFT along the low order dimension of an array+--+idft :: (Shape sh, Slice sh, Elt e, IsFloating e)+ => Acc (Array (sh:.Int) (Complex e))+ -> Acc (Array (sh:.Int) (Complex e))+idft v+ = let sh = shape v+ n = indexHead sh+ roots = inverseRootsOfUnity sh+ scale = lift (A.fromIntegral n, constant 0)+ in+ A.map (/scale) $ dftG roots v+++-- | Generic function for computation of forward and inverse DFT. This function+-- is also useful if you transform many arrays of the same extent, and don't+-- want to recompute the roots for each one.+--+-- The extent of the input and roots must match.+--+dftG :: forall sh e. (Shape sh, Slice sh, Elt e, IsFloating e)+ => Acc (Array (sh:.Int) (Complex e)) -- ^ roots of unity+ -> Acc (Array (sh:.Int) (Complex e)) -- ^ input array+ -> Acc (Array (sh:.Int) (Complex e))+dftG roots arr+ = A.fold (+) (constant (0,0))+ $ A.zipWith (*) arr' roots'+ where+ base = shape arr+ l = indexHead base+ extend = lift (base :. shapeSize base)++ -- Extend the entirety of the input arrays into a higher dimension, reading+ -- roots from the appropriate places and then reduce along this axis.+ --+ -- In the calculation for 'roots'', 'i' is the index into the extended+ -- dimension, with corresponding base index 'ix' which we are attempting to+ -- calculate the single DFT value of. The rest proceeds as per 'dftGS'.+ --+ arr' = A.generate extend (\ix' -> let i = indexHead ix' in arr !! i)+ roots' = A.generate extend (\ix' -> let ix :. i = unlift ix'+ sh :. n = unlift (fromIndex base i) :: Exp sh :. Exp Int+ k = indexHead ix+ in+ roots ! lift (sh :. (k*n) `mod` l))+++-- | Compute a single value of the DFT.+--+dftGS :: forall sh e. (Shape sh, Slice sh, Elt e, IsFloating e)+ => Exp (sh :. Int) -- ^ index of the value we want+ -> Acc (Array (sh:.Int) (Complex e)) -- ^ roots of unity+ -> Acc (Array (sh:.Int) (Complex e)) -- ^ input array+ -> Acc (Scalar (Complex e))+dftGS ix roots arr+ = let k = indexHead ix+ l = indexHead (shape arr)++ -- all the roots we need to multiply with+ roots' = A.generate (shape arr)+ (\ix' -> let sh :. n = unlift ix' :: Exp sh :. Exp Int+ in roots ! lift (sh :. (k*n) `mod` l))+ in+ A.foldAll (+) (constant (0,0)) $ A.zipWith (*) arr roots'+
+ Data/Array/Accelerate/Math/DFT/Centre.hs view
@@ -0,0 +1,103 @@+{-# LANGUAGE TypeOperators #-}+-- |+-- Module : Data.Array.Accelerate.Math.DFT.Centre+-- Copyright : [2012..2013] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell, Robert Clifton-Everest+-- License : BSD3+--+-- Maintainer : Manuel M T Chakravarty <chak@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--+-- These transforms allow the centering of the frequency domain of a DFT such+-- that the the zero frequency is in the middle. The centering transform, when+-- performed on the input of a DFT, will cause zero frequency to be centred in+-- the middle. The shifting transform however takes the output of a DFT to+-- give the same result. Therefore the relationship between the two is:+--+-- > fft(center(X)) = shift(fft(X))+--+module Data.Array.Accelerate.Math.DFT.Centre (++ centre1D, centre2D, centre3D,+ shift1D, shift2D, shift3D,++) where++import Prelude as P+import Data.Array.Accelerate as A+import Data.Array.Accelerate.Math.Complex+++-- | Apply the centring transform to a vector+--+centre1D :: (Elt e, IsFloating e)+ => Acc (Array DIM1 (Complex e))+ -> Acc (Array DIM1 (Complex e))+centre1D arr+ = A.generate (shape arr)+ (\ix -> let Z :. x = unlift ix :: Z :. Exp Int+ in lift (-1 ** A.fromIntegral x, A.constant 0) * arr!ix)++-- | Apply the centring transform to a matrix+--+centre2D :: (Elt e, IsFloating e)+ => Acc (Array DIM2 (Complex e))+ -> Acc (Array DIM2 (Complex e))+centre2D arr+ = A.generate (shape arr)+ (\ix -> let Z :. y :. x = unlift ix :: Z :. Exp Int :. Exp Int+ in lift (-1 ** A.fromIntegral (y + x), A.constant 0) * arr!ix)++-- | Apply the centring transform to a 3D array+--+centre3D :: (Elt e, IsFloating e)+ => Acc (Array DIM3 (Complex e))+ -> Acc (Array DIM3 (Complex e))+centre3D arr+ = A.generate (shape arr)+ (\ix -> let Z :. z :. y :. x = unlift ix :: Z :. Exp Int :. Exp Int :. Exp Int+ in lift (-1 ** A.fromIntegral (z + y + x), A.constant 0) * arr!ix)+++-- | Apply the shifting transform to a vector+--+shift1D :: Elt e => Acc (Vector e) -> Acc (Vector e)+shift1D arr+ = A.backpermute (A.shape arr) p arr+ where+ p ix+ = let Z:.x = unlift ix :: Z :. Exp Int+ in index1 (x <* mw ? (x + mw, x - mw))+ Z:.w = unlift (A.shape arr)+ mw = w `div` 2+++-- | Apply the shifting transform to a 2D array+--+shift2D :: Elt e => Acc (Array DIM2 e) -> Acc (Array DIM2 e)+shift2D arr+ = A.backpermute (A.shape arr) p arr+ where+ p ix+ = let Z:.y:.x = unlift ix :: Z :. Exp Int :. Exp Int+ in index2 (y <* mh ? (y + mh, y - mh))+ (x <* mw ? (x + mw, x - mw))+ Z:.h:.w = unlift (A.shape arr)+ (mh,mw) = (h `div` 2, w `div` 2)+++-- | Apply the shifting transform to a 3D array+--+shift3D :: Elt e => Acc (Array DIM3 e) -> Acc (Array DIM3 e)+shift3D arr+ = A.backpermute (A.shape arr) p arr+ where+ p ix+ = let Z:.z:.y:.x = unlift ix :: Z :. Exp Int :. Exp Int :. Exp Int+ in index3 (z <* md ? (z + md, z - md))+ (y <* mh ? (y + mh, y - mh))+ (x <* mw ? (x + mw, x - mw))+ Z:.h:.w:.d = unlift (A.shape arr)+ (mh,mw,md) = (h `div` 2, w `div` 2, d `div` 2)+ index3 i j k = lift (Z:.i:.j:.k)+
+ Data/Array/Accelerate/Math/DFT/Roots.hs view
@@ -0,0 +1,50 @@+{-# LANGUAGE TypeOperators #-}+-- |+-- Module : Data.Array.Accelerate.Math.DFT.Roots+-- Copyright : [2012] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Manuel M T Chakravarty <chak@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--+module Data.Array.Accelerate.Math.DFT.Roots (++ rootsOfUnity, inverseRootsOfUnity,++) where++import Prelude as P+import Data.Array.Accelerate as A+import Data.Array.Accelerate.Math.Complex+++-- | Calculate the roots of unity for the forward transform+--+rootsOfUnity+ :: (Elt e, IsFloating e, Shape sh, Slice sh)+ => Exp (sh :. Int)+ -> Acc (Array (sh:.Int) (Complex e))+rootsOfUnity sh =+ let n = A.fromIntegral (A.indexHead sh)+ in+ A.generate sh (\ix -> let i = A.fromIntegral (A.indexHead ix)+ k = 2 * pi * i / n+ in+ A.lift ( cos k, -sin k ))+++-- | Calculate the roots of unity for an inverse transform+--+inverseRootsOfUnity+ :: (Elt e, IsFloating e, Shape sh, Slice sh)+ => Exp (sh :. Int)+ -> Acc (Array (sh:.Int) (Complex e))+inverseRootsOfUnity sh =+ let n = A.fromIntegral (A.indexHead sh)+ in+ A.generate sh (\ix -> let i = A.fromIntegral (A.indexHead ix)+ k = 2 * pi * i / n+ in+ A.lift ( cos k, sin k ))+
+ Data/Array/Accelerate/Math/FFT.hs view
@@ -0,0 +1,374 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE EmptyDataDecls #-}+{-# LANGUAGE ForeignFunctionInterface #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ViewPatterns #-}+-- |+-- Module : Data.Array.Accelerate.Math.FFT+-- Copyright : [2012..2013] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell, Robert Clifton-Everest+-- License : BSD3+--+-- Maintainer : Manuel M T Chakravarty <chak@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--+-- Computation of a Discrete Fourier Transform using the Cooley-Tuckey+-- algorithm. The time complexity is O(n log n) in the size of the input.+--+-- This uses a naïve divide-and-conquer algorithm whose absolute performance is+-- appalling.+--+module Data.Array.Accelerate.Math.FFT (++ Mode(..),+ fft1D, fft1D',+ fft2D, fft2D',+ fft3D, fft3D',+ fft++) where++import Prelude as P+import Data.Array.Accelerate as A+import Data.Array.Accelerate.Array.Sugar ( showShape )+import Data.Array.Accelerate.Math.Complex++#ifdef ACCELERATE_CUDA_BACKEND+import Data.Array.Accelerate.CUDA.Foreign+import Data.Array.Accelerate.Array.Sugar as S ( shapeToList, shape, EltRepr )+import Data.Array.Accelerate.Type++import Data.Functor+import Foreign.CUDA.FFT+import qualified Foreign.CUDA.Driver as CUDA hiding (free)+#endif++import Data.Bits++data Mode = Forward | Reverse | Inverse+ deriving (Eq, Show)++isPow2 :: Int -> Bool+isPow2 x = x .&. (x-1) == 0++signOfMode :: Num a => Mode -> a+signOfMode m+ = case m of+ Forward -> -1+ Reverse -> 1+ Inverse -> 1+++-- Vector Transform+-- ----------------+--+-- Discrete Fourier Transform of a vector. Array dimensions must be powers of+-- two else error.+--+fft1D :: (Elt e, IsFloating e)+ => Mode+ -> Vector (Complex e)+ -> Acc (Vector (Complex e))+fft1D mode vec+ = let Z :. len = arrayShape vec+ in+ fft1D' mode len (use vec)++fft1D' :: forall e. (Elt e, IsFloating e)+ => Mode+ -> Int+ -> Acc (Vector (Complex e))+ -> Acc (Vector (Complex e))+fft1D' mode len vec+ = let sign = signOfMode mode :: e+ scale = P.fromIntegral len+#ifdef ACCELERATE_CUDA_BACKEND+ sh = (Z:.len)+ vec' = cudaFFT mode sh fft' vec+#else+ vec' = fft' vec+#endif+ fft' a = fft sign Z len a+ in+ if P.not (isPow2 len)+ then error $ unlines+ [ "Data.Array.Accelerate.FFT: fft1D"+ , " Array dimensions must be powers of two, but are: " ++ showShape (Z:.len) ]++ else case mode of+ Inverse -> A.map (/scale) vec'+ _ -> vec'+++-- Matrix Transform+-- ----------------+--+-- Discrete Fourier Transform of a matrix. Array dimensions must be powers of+-- two else error.+--+fft2D :: (Elt e, IsFloating e)+ => Mode+ -> Array DIM2 (Complex e)+ -> Acc (Array DIM2 (Complex e))+fft2D mode arr+ = let Z :. height :. width = arrayShape arr+ in+ fft2D' mode width height (use arr)+++fft2D' :: forall e. (Elt e, IsFloating e)+ => Mode+ -> Int -- ^ width+ -> Int -- ^ height+ -> Acc (Array DIM2 (Complex e))+ -> Acc (Array DIM2 (Complex e))+fft2D' mode width height arr+ = let sign = signOfMode mode :: e+ scale = P.fromIntegral (width * height)+#ifdef ACCELERATE_CUDA_BACKEND+ sh = (Z:.width:.height)+ arr' = cudaFFT mode sh fft' arr+#else+ arr' = fft' arr+#endif+ fft' a = A.transpose . fft sign (Z:.width) height+ >-> A.transpose . fft sign (Z:.height) width+ $ a+ in+ if P.not (isPow2 width && isPow2 height)+ then error $ unlines+ [ "Data.Array.Accelerate.FFT: fft2D"+ , " Array dimensions must be powers of two, but are: " ++ showShape (Z:.height:.width) ]++ else case mode of+ Inverse -> A.map (/scale) arr'+ _ -> arr'+++-- Cube Transform+-- --------------+--+-- Discrete Fourier Transform of a 3D array. Array dimensions must be power of+-- two else error.+--+fft3D :: (Elt e, IsFloating e)+ => Mode+ -> Array DIM3 (Complex e)+ -> Acc (Array DIM3 (Complex e))+fft3D mode arr+ = let Z :. depth :. height :. width = arrayShape arr+ in+ fft3D' mode width height depth (use arr)+++fft3D' :: forall e. (Elt e, IsFloating e)+ => Mode+ -> Int -- ^ width+ -> Int -- ^ height+ -> Int -- ^ depth+ -> Acc (Array DIM3 (Complex e))+ -> Acc (Array DIM3 (Complex e))+fft3D' mode width height depth arr+ = let sign = signOfMode mode :: e+ scale = P.fromIntegral (width * height)+#ifdef ACCELERATE_CUDA_BACKEND+ sh = (Z:.width:.height:.depth)+ arr' = cudaFFT mode sh fft' arr+#else+ arr' = fft' arr+#endif+ fft' a = rotate3D . fft sign (Z:.width :.depth) height+ >-> rotate3D . fft sign (Z:.height:.width) depth+ >-> rotate3D . fft sign (Z:.depth :.height) width+ $ a+ in+ if P.not (isPow2 width && isPow2 height && isPow2 depth)+ then error $ unlines+ [ "Data.Array.Accelerate.FFT: fft3D"+ , " Array dimensions must be powers of two, but are: " ++ showShape (Z:.depth:.height:.width) ]++ else case mode of+ Inverse -> A.map (/scale) arr'+ _ -> arr'++++rotate3D :: Elt e => Acc (Array DIM3 e) -> Acc (Array DIM3 e)+rotate3D arr+ = backpermute (swap (A.shape arr)) swap arr+ where+ swap :: Exp DIM3 -> Exp DIM3+ swap ix =+ let Z :. m :. k :. l = unlift ix :: Z :. Exp Int :. Exp Int :. Exp Int+ in lift $ Z :. k :. l :. m+++-- Rank-generalised Cooley-Tuckey DFT+--+-- We require the innermost dimension be passed as a Haskell value because we+-- can't do divide-and-conquer recursion directly in the meta-language.+--+fft :: forall sh e. (Slice sh, Shape sh, IsFloating e, Elt e)+ => e+ -> sh+ -> Int+ -> Acc (Array (sh:.Int) (Complex e))+ -> Acc (Array (sh:.Int) (Complex e))+fft sign sh sz arr = go sz 0 1+ where+ go :: Int -> Int -> Int -> Acc (Array (sh:.Int) (Complex e))+ go len offset stride+ | len == 2+ = A.generate (constant (sh :. len)) swivel++ | otherwise+ = combine+ (go (len `div` 2) offset (stride * 2))+ (go (len `div` 2) (offset + stride) (stride * 2))++ where+ len' = the (unit (constant len))+ offset' = the (unit (constant offset))+ stride' = the (unit (constant stride))++ swivel ix =+ let sh' :. sz' = unlift ix :: Exp sh :. Exp Int+ in+ sz' ==* 0 ? ( (arr ! lift (sh' :. offset')) + (arr ! lift (sh' :. offset' + stride'))+ {- ==* 1-} , (arr ! lift (sh' :. offset')) - (arr ! lift (sh' :. offset' + stride')) )++ combine evens odds =+ let odds' = A.generate (A.shape odds) (\ix -> twiddle len' (indexHead ix) * odds!ix)+ in+ append (A.zipWith (+) evens odds') (A.zipWith (-) evens odds')++ twiddle n' i' =+ let n = A.fromIntegral n'+ i = A.fromIntegral i'+ k = 2*pi*i/n+ in+ lift ( cos k, A.constant sign * sin k )++#ifdef ACCELERATE_CUDA_BACKEND+-- FFT using the CUFFT library to enable high performance for the CUDA backend of+-- Accelerate. The implementation works on all arrays of rank less than or equal+-- to 3. The result is un-normalised.+--+cudaFFT :: forall e sh. (Shape sh, Elt e, IsFloating e)+ => Mode+ -> sh+ -> (Acc (Array sh (Complex e)) -> Acc (Array sh (Complex e)))+ -> Acc (Array sh (Complex e))+ -> Acc (Array sh (Complex e))+cudaFFT mode sh p arr = deinterleave sh (foreignAcc ff pure (interleave arr))+ where+ ff = cudaAcc foreignFFT+ -- Unfortunately the pure version of the function needs to be wrapped in+ -- interleave and deinterleave to match how the foreign version works.+ --+ -- RCE: Do the interleaving and deinterleaving in foreignFFT+ --+ -- TLM: The interleaving might get fused into other parts of the+ -- computation and thus be okay. We should really support multi types+ -- such as float2 instead.+ --+ pure = interleave . p . deinterleave sh+ sign = signOfMode mode :: Int++ foreignFFT :: Array DIM1 e -> CIO (Array DIM1 e)+ foreignFFT arr' = do+ -- Create the plan+ -- TODO: Cache this.+ --+ hndl <- liftIO $+ case shapeToList sh of+ [width] -> plan1D width types 1+ [height, width] -> plan2D height width types+ [depth, height, width] -> plan3D depth height width types+ _ -> error "Accelerate-fft cannot use CUFFT for arrays of dimensions higher than 3"++ output <- allocateArray (S.shape arr')+ iptr <- floatingDevicePtr arr'+ optr <- floatingDevicePtr output++ --Execute+ liftIO $ execute hndl iptr optr++ liftIO $ destroy hndl++ return output++ types+ = case (floatingType :: FloatingType e) of+ TypeFloat{} -> C2C+ TypeDouble{} -> Z2Z+ TypeCFloat{} -> C2C+ TypeCDouble{} -> Z2Z++ execute :: Handle -> CUDA.DevicePtr e -> CUDA.DevicePtr e -> IO ()+ execute hndl iptr optr+ = case (floatingType :: FloatingType e) of+ TypeFloat{} -> execC2C hndl iptr optr sign+ TypeDouble{} -> execZ2Z hndl iptr optr sign+ TypeCFloat{} -> execC2C hndl (CUDA.castDevPtr iptr) (CUDA.castDevPtr optr) sign+ TypeCDouble{} -> execZ2Z hndl (CUDA.castDevPtr iptr) (CUDA.castDevPtr optr) sign++ floatingDevicePtr :: Vector e -> CIO (CUDA.DevicePtr e)+ floatingDevicePtr v+ = case (floatingType :: FloatingType e) of+ TypeFloat{} -> singleDevicePtr v+ TypeDouble{} -> singleDevicePtr v+ TypeCFloat{} -> CUDA.castDevPtr <$> singleDevicePtr v+ TypeCDouble{} -> CUDA.castDevPtr <$> singleDevicePtr v++ singleDevicePtr :: DevicePtrs (EltRepr e) ~ ((),CUDA.DevicePtr b) => Vector e -> CIO (CUDA.DevicePtr b)+ singleDevicePtr v = P.snd <$> devicePtrsOfArray v+#endif++-- Append two arrays. Doesn't do proper bounds checking or intersection...+--+append+ :: forall sh e. (Slice sh, Shape sh, Elt e)+ => Acc (Array (sh:.Int) e)+ -> Acc (Array (sh:.Int) e)+ -> Acc (Array (sh:.Int) e)+append xs ys+ = let sh :. n = unlift (A.shape xs) :: Exp sh :. Exp Int+ _ :. m = unlift (A.shape ys) :: Exp sh :. Exp Int+ in+ generate (lift (sh :. n+m))+ (\ix -> let sz :. i = unlift ix :: Exp sh :. Exp Int+ in i <* n ? (xs ! lift (sz:.i), ys ! lift (sz:.i-n) ))+++#ifdef ACCELERATE_CUDA_BACKEND+{-# RULES+ "interleave/deinterleave" forall sh x. deinterleave sh (interleave x) = x;+ "deinterleave/interleave" forall sh x. interleave (deinterleave sh x) = x+ #-}++-- Interleave the real and imaginary components in a complex array and produce a+-- flattened vector. This allows us to mimic the float2 structure used by CUFFT+-- to store complex numbers.+--+interleave :: (Shape sh, Elt e) => Acc (Array sh (Complex e)) -> Acc (Vector e)+interleave arr = generate sh swizzle+ where+ sh = index1 (2 * A.size arr)+ swizzle ix =+ let i = indexHead ix+ v = arr A.!! (i `div` 2)+ in+ i `mod` 2 ==* 0 ? (real v, imag v)++-- Deinterleave a vector into a complex array. Assumes the array is even in length.+--+deinterleave :: (Shape sh, Elt e) => sh -> Acc (Vector e) -> Acc (Array sh (Complex e))+deinterleave (constant -> sh) arr =+ generate sh (\ix -> let i = toIndex sh ix * 2+ in lift (arr A.!! i, arr A.!! (i+1)))+#endif+
+ LICENSE view
@@ -0,0 +1,23 @@+Copyright (c) [2007..2012] The Accelerate Team. All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:+ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.+ * Redistributions in binary form must reproduce the above copyright+ notice, this list of conditions and the following disclaimer in the+ documentation and/or other materials provided with the distribution.+ * Neither the names of the contributors nor of their affiliations may+ be used to endorse or promote products derived from this software+ without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ''AS IS'' AND ANY+EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE+DISCLAIMED. IN NO EVENT SHALL COPYRIGHT HOLDERS BE LIABLE FOR ANY+DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ accelerate-fft.cabal view
@@ -0,0 +1,60 @@+Name: accelerate-fft+Version: 0.13.0.0+Cabal-version: >= 1.6+Tested-with: GHC >= 7.4+Build-type: Simple++Synopsis: FFT using the Accelerate library+Description:+ Rank-polymorphic discrete Fourier transform (DFT), computed with a fast+ Fourier transform (FFT) algorithm using the Accelerate library+ .+ Refer to the main /Accelerate/ package for more information:+ <http://hackage.haskell.org/package/accelerate>+ .++License: BSD3+License-file: LICENSE+Author: Manuel M T Chakravarty,+ Gabriele Keller,+ Trevor L. McDonell+Maintainer: Manuel M T Chakravarty <chak@cse.unsw.edu.au>+Homepage: https://github.com/AccelerateHS/accelerate-fft+Bug-reports: https://github.com/AccelerateHS/accelerate/issues++Category: Compilers/Interpreters, Concurrency, Data, Parallelism+Stability: Experimental++Flag cuda+ Description: Enable support for using CUFFT via the CUDA backend's+ FFI+ Default: True++Library+ Build-depends: accelerate == 0.13.*,+ base == 4.*++ Exposed-modules: Data.Array.Accelerate.Math.Complex+ Data.Array.Accelerate.Math.FFT+ Data.Array.Accelerate.Math.DFT+ Data.Array.Accelerate.Math.DFT.Centre+ Data.Array.Accelerate.Math.DFT.Roots++ ghc-options: -O2 -Wall -funbox-strict-fields++ if flag(cuda)+ CPP-options: -DACCELERATE_CUDA_BACKEND+ Build-depends: accelerate-cuda == 0.13.*,+ cuda >= 0.5 && < 0.6,+ cufft >= 0.1 && < 0.2++ -- Don't add the extensions list here. Instead, place individual LANGUAGE+ -- pragmas in the files that require a specific extension. This means the+ -- project loads in GHCi, and avoids extension clashes.+ --+ -- Extensions:++Source-repository head+ Type: git+ Location: git://github.com/AccelerateHS/accelerate-fft.git+