diff --git a/Data/Array/Accelerate/Math/Complex.hs b/Data/Array/Accelerate/Math/Complex.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Math/Complex.hs @@ -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) + diff --git a/Data/Array/Accelerate/Math/DFT.hs b/Data/Array/Accelerate/Math/DFT.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Math/DFT.hs @@ -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 +-- 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' + diff --git a/Data/Array/Accelerate/Math/DFT/Centre.hs b/Data/Array/Accelerate/Math/DFT/Centre.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Math/DFT/Centre.hs @@ -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 +-- 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) + diff --git a/Data/Array/Accelerate/Math/DFT/Roots.hs b/Data/Array/Accelerate/Math/DFT/Roots.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Math/DFT/Roots.hs @@ -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 +-- 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 )) + diff --git a/Data/Array/Accelerate/Math/FFT.hs b/Data/Array/Accelerate/Math/FFT.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Math/FFT.hs @@ -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 +-- 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 + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/accelerate-fft.cabal b/accelerate-fft.cabal new file mode 100644 --- /dev/null +++ b/accelerate-fft.cabal @@ -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: + + . + +License: BSD3 +License-file: LICENSE +Author: Manuel M T Chakravarty, + Gabriele Keller, + Trevor L. McDonell +Maintainer: Manuel M T Chakravarty +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 +