accelerate-fft 1.1.0.0 → 1.2.0.0
raw patch · 32 files changed
+2301/−1596 lines, 32 filesdep +accelerate-fftdep +containersdep +hashabledep −storable-complexdep ~acceleratedep ~accelerate-llvmdep ~accelerate-llvm-native
Dependencies added: accelerate-fft, containers, hashable, hedgehog, lens-accelerate, mtl, tasty, tasty-hedgehog, unordered-containers
Dependencies removed: storable-complex
Dependency ranges changed: accelerate, accelerate-llvm, accelerate-llvm-native, accelerate-llvm-ptx, base, cufft
Files
- CHANGELOG.md +21/−2
- Data/Array/Accelerate/Math/DFT.hs +0/−113
- Data/Array/Accelerate/Math/DFT/Centre.hs +0/−105
- Data/Array/Accelerate/Math/DFT/Roots.hs +0/−53
- Data/Array/Accelerate/Math/FFT.hs +0/−298
- Data/Array/Accelerate/Math/FFT/LLVM/Native.hs +0/−255
- Data/Array/Accelerate/Math/FFT/LLVM/PTX.hs +0/−279
- Data/Array/Accelerate/Math/FFT/Mode.hs +0/−28
- Data/Array/Accelerate/Math/FFT/Twine.hs +0/−88
- accelerate-fft.cabal +86/−33
- cubits/twine_f32.cu +0/−63
- cubits/twine_f32.ptx +0/−108
- cubits/twine_f64.cu +0/−63
- cubits/twine_f64.ptx +0/−108
- src/Data/Array/Accelerate/Math/DFT.hs +113/−0
- src/Data/Array/Accelerate/Math/DFT/Centre.hs +108/−0
- src/Data/Array/Accelerate/Math/DFT/Roots.hs +53/−0
- src/Data/Array/Accelerate/Math/FFT.hs +258/−0
- src/Data/Array/Accelerate/Math/FFT/Adhoc.hs +603/−0
- src/Data/Array/Accelerate/Math/FFT/LLVM/Native.hs +110/−0
- src/Data/Array/Accelerate/Math/FFT/LLVM/Native/Base.hs +119/−0
- src/Data/Array/Accelerate/Math/FFT/LLVM/Native/Ix.hs +57/−0
- src/Data/Array/Accelerate/Math/FFT/LLVM/PTX.hs +170/−0
- src/Data/Array/Accelerate/Math/FFT/LLVM/PTX/Base.hs +82/−0
- src/Data/Array/Accelerate/Math/FFT/LLVM/PTX/Plans.hs +91/−0
- src/Data/Array/Accelerate/Math/FFT/Mode.hs +28/−0
- src/Data/Array/Accelerate/Math/FFT/Type.hs +37/−0
- test/Test/Base.hs +60/−0
- test/Test/FFT.hs +237/−0
- test/Test/ShowType.hs +30/−0
- test/TestNative.hs +19/−0
- test/TestPTX.hs +19/−0
CHANGELOG.md view
@@ -6,14 +6,33 @@ project adheres to the [Haskell Package Versioning Policy (PVP)](https://pvp.haskell.org) +## [1.2.0.0] - 2018-04-03+### Changed+ * update for AoS representation of complex numbers+ * improve pure FFT implementation++### Contributors++Special thanks to those who contributed patches as part of this release:++ * Trevor L. McDonell (@tmcdonell)+ * Rinat Striungis (@Haskell-mouse)++ ## [1.1.0.0] - 2017-09-21- * [#5]: fix to ignore `sh` parameter in inverse mode+### Changed+ * build against FFTW and cuFFT foreign implementations by default++### Fixed+ * fix to ignore `sh` parameter in inverse mode ([#5])++### Removed * Drop support for (deprecated) `accelerate-cuda` backend- * Build against FFTW and CUFFT foreign implementations by default ## 1.0.0.0 - 2017-03-31 +[1.2.0.0]: https://github.com/AccelerateHS/accelerate-fft/compare/1.1.0.0...1.2.0.0 [1.1.0.0]: https://github.com/AccelerateHS/accelerate-fft/compare/1.0.0.0...1.1.0.0 [#5]: https://github.com/AccelerateHS/accelerate-fft/pull/5
− Data/Array/Accelerate/Math/DFT.hs
@@ -1,113 +0,0 @@-{-# LANGUAGE ConstraintKinds #-}-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TypeOperators #-}--- |--- Module : Data.Array.Accelerate.Math.DFT--- Copyright : [2012..2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell--- [2013..2017] Robert Clifton-Everest--- License : BSD3------ Maintainer : Trevor L. McDonell <tmcdonell@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.Data.Complex----- | Compute the DFT along the low order dimension of an array----dft :: (Shape sh, Slice sh, A.RealFloat e, A.FromIntegral Int 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, A.RealFloat e, A.FromIntegral Int 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 :+ 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, A.RealFloat 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 (+) 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, A.RealFloat 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 (+) 0 $ A.zipWith (*) arr roots'-
− Data/Array/Accelerate/Math/DFT/Centre.hs
@@ -1,105 +0,0 @@-{-# LANGUAGE ConstraintKinds #-}-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE TypeOperators #-}--- |--- Module : Data.Array.Accelerate.Math.DFT.Centre--- Copyright : [2012..2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell--- [2013..2017] Robert Clifton-Everest--- License : BSD3------ Maintainer : Trevor L. McDonell <tmcdonell@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.Data.Complex----- | Apply the centring transform to a vector----centre1D :: (A.RealFloat e, A.FromIntegral Int 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) :+ 0) * arr!ix)---- | Apply the centring transform to a matrix----centre2D :: (A.RealFloat e, A.FromIntegral Int 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)) :+ 0) * arr!ix)---- | Apply the centring transform to a 3D array----centre3D :: (A.RealFloat e, A.FromIntegral Int 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)) :+ 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 A.< 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 A.< mh ? (y + mh, y - mh))- (x A.< 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 A.< md ? (z + md, z - md))- (y A.< mh ? (y + mh, y - mh))- (x A.< 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)-
− Data/Array/Accelerate/Math/DFT/Roots.hs
@@ -1,53 +0,0 @@-{-# LANGUAGE ConstraintKinds #-}-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE TypeOperators #-}--- |--- Module : Data.Array.Accelerate.Math.DFT.Roots--- Copyright : [2012..2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell--- [2013..2017] Robert Clifton-Everest--- License : BSD3------ Maintainer : Trevor L. McDonell <tmcdonell@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.Data.Complex----- | Calculate the roots of unity for the forward transform----rootsOfUnity- :: (Shape sh, Slice sh, A.Floating e, A.FromIntegral Int e)- => 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- :: (Shape sh, Slice sh, A.Floating e, A.FromIntegral Int e)- => 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
@@ -1,298 +0,0 @@-{-# LANGUAGE CPP #-}-{-# LANGUAGE ConstraintKinds #-}-{-# LANGUAGE EmptyDataDecls #-}-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE ForeignFunctionInterface #-}-{-# LANGUAGE GADTs #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE TypeOperators #-}-{-# LANGUAGE ViewPatterns #-}--- |--- Module : Data.Array.Accelerate.Math.FFT--- Copyright : [2012..2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell--- [2013..2017] Robert Clifton-Everest--- License : BSD3------ Maintainer : Trevor L. McDonell <tmcdonell@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.------ The base (default) implementation uses a naïve divide-and-conquer algorithm--- whose absolute performance is appalling. It also requires that you know on--- the Haskell side the size of the data being transformed, and that this is--- a power-of-two in each dimension.------ For performance, compile against the foreign library bindings (using any--- number of '-fllvm-ptx', and '-fllvm-cpu' for the accelerate-llvm-ptx, and--- accelerate-llvm-native backends, respectively), which have none of the above--- restrictions.-----module Data.Array.Accelerate.Math.FFT (-- Mode(..),- FFTElt,- fft1D, fft1D',- fft2D, fft2D',- fft3D, fft3D',- fft--) where--import Data.Array.Accelerate as A-import Data.Array.Accelerate.Array.Sugar ( showShape, shapeToList )-import Data.Array.Accelerate.Data.Complex-import Data.Array.Accelerate.Math.FFT.Mode--#ifdef ACCELERATE_LLVM_NATIVE_BACKEND-import qualified Data.Array.Accelerate.Math.FFT.LLVM.Native as Native-#endif-#ifdef ACCELERATE_LLVM_PTX_BACKEND-import qualified Data.Array.Accelerate.Math.FFT.LLVM.PTX as PTX-#endif--import Data.Bits-import Text.Printf-import Prelude as P----- The type of supported FFT elements; namely 'Float' and 'Double'.----type FFTElt e = (P.Num e, A.RealFloat e, A.FromIntegral Int e, A.IsFloating e)----- Vector Transform--- -------------------- | Discrete Fourier Transform of a vector.------ The default implementation requires the array dimension to be a power of two--- (else error).----fft1D :: FFTElt e- => Mode- -> Array DIM1 (Complex e)- -> Acc (Array DIM1 (Complex e))-fft1D mode vec- = fft1D' mode (arrayShape vec) (use vec)----- | Discrete Fourier Transform of a vector.------ The default implementation requires the array dimension to be a power of two.--- The FFI-backed implementations ignore the Haskell-side size parameter (second--- argument).----fft1D' :: forall e. FFTElt e- => Mode- -> DIM1- -> Acc (Array DIM1 (Complex e))- -> Acc (Array DIM1 (Complex e))-fft1D' mode (Z :. len) arr- = let sign = signOfMode mode :: e- scale = A.fromIntegral (A.length arr)- go =-#ifdef ACCELERATE_LLVM_NATIVE_BACKEND- foreignAcc (Native.fft1D mode) $-#endif-#ifdef ACCELERATE_LLVM_PTX_BACKEND- foreignAcc (PTX.fft1D mode) $-#endif- fft sign Z len- in- case mode of- Inverse -> A.map (/scale) (go arr)- _ -> go arr----- Matrix Transform--- -------------------- | Discrete Fourier Transform of a matrix.------ The default implementation requires the array dimensions to be powers of two--- (else error).----fft2D :: FFTElt e- => Mode- -> Array DIM2 (Complex e)- -> Acc (Array DIM2 (Complex e))-fft2D mode arr- = fft2D' mode (arrayShape arr) (use arr)----- | Discrete Fourier Transform of a matrix.------ The default implementation requires the array dimensions to be powers of two.--- The FFI-backed implementations ignore the Haskell-side size parameter (second--- argument).----fft2D' :: forall e. FFTElt e- => Mode- -> DIM2- -> Acc (Array DIM2 (Complex e))- -> Acc (Array DIM2 (Complex e))-fft2D' mode (Z :. height :. width) arr- = let sign = signOfMode mode :: e- scale = A.fromIntegral (A.size arr)- go =-#ifdef ACCELERATE_LLVM_NATIVE_BACKEND- foreignAcc (Native.fft2D mode) $-#endif-#ifdef ACCELERATE_LLVM_PTX_BACKEND- foreignAcc (PTX.fft2D mode) $-#endif- fft'-- fft' a = A.transpose . fft sign (Z:.height) width- >-> A.transpose . fft sign (Z:.width) height- $ a- in- case mode of- Inverse -> A.map (/scale) (go arr)- _ -> go arr----- Cube Transform--- ------------------ | Discrete Fourier Transform of a 3D array.------ The default implementation requires the array dimensions to be powers of two--- (else error).----fft3D :: FFTElt e- => Mode- -> Array DIM3 (Complex e)- -> Acc (Array DIM3 (Complex e))-fft3D mode arr- = fft3D' mode (arrayShape arr) (use arr)----- | Discrete Fourier Transform of a 3D array.------ The default implementation requires the array dimensions to be powers of two.--- The FFI-backed implementations ignore the Haskell-side size parameter (second--- argument).----fft3D' :: forall e. FFTElt e- => Mode- -> DIM3- -> Acc (Array DIM3 (Complex e))- -> Acc (Array DIM3 (Complex e))-fft3D' mode (Z :. depth :. height :. width) arr- = let sign = signOfMode mode :: e- scale = A.fromIntegral (A.size arr)- go =-#ifdef ACCELERATE_LLVM_NATIVE_BACKEND- foreignAcc (Native.fft3D mode) $-#endif-#ifdef ACCELERATE_LLVM_PTX_BACKEND- foreignAcc (PTX.fft3D mode) $-#endif- fft'-- fft' a = rotate3D . fft sign (Z:.depth :.height) width- >-> rotate3D . fft sign (Z:.height:.width) depth- >-> rotate3D . fft sign (Z:.width :.depth) height- $ a- in- case mode of- Inverse -> A.map (/scale) (go arr)- _ -> go arr---rotate3D :: Elt e => Acc (Array DIM3 e) -> Acc (Array DIM3 e)-rotate3D arr = backpermute sh rot arr- where- sh :: Exp DIM3- sh =- let Z :. z :. y :. x = unlift (shape arr) :: Z :. Exp Int :. Exp Int :. Exp Int- in index3 y x z- --- rot :: Exp DIM3 -> Exp DIM3- rot ix =- let Z :. z :. y :. x = unlift ix :: Z :. Exp Int :. Exp Int :. Exp Int- in index3 x z y----- 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, A.RealFloat e, A.FromIntegral Int e)- => e- -> sh- -> Int- -> Acc (Array (sh:.Int) (Complex e))- -> Acc (Array (sh:.Int) (Complex e))-fft sign sh sz arr- | P.any (P.not . isPow2) (shapeToList (sh:.sz))- = error $ printf "fft: array dimensions must be powers-of-two, but are: %s" (showShape (sh:.sz))- --- | otherwise- = go sz 0 1- where- go :: Int -> Int -> Int -> Acc (Array (sh:.Int) (Complex e))- go len offset stride- | len P.== 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' A.== 0 ? ( (arr ! lift (sh' :. offset')) + (arr ! lift (sh' :. offset' + stride'))- {- A.== 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 )----- Append two arrays. This is a specialised version of (A.++) which does not do--- 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 A.< n ? (xs ! lift (sz:.i), ys ! lift (sz:.i-n) ))---isPow2 :: Int -> Bool-isPow2 0 = True-isPow2 1 = False-isPow2 x = x .&. (x-1) P.== 0---
− Data/Array/Accelerate/Math/FFT/LLVM/Native.hs
@@ -1,255 +0,0 @@-{-# LANGUAGE BangPatterns #-}-{-# LANGUAGE GADTs #-}-{-# LANGUAGE PatternGuards #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TemplateHaskell #-}-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE TypeOperators #-}--- |--- Module : Data.Array.Accelerate.Math.FFT.LLVM.Native--- Copyright : [2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell--- License : BSD3------ Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>--- Stability : experimental--- Portability : non-portable (GHC extensions)-----module Data.Array.Accelerate.Math.FFT.LLVM.Native (-- fft1D,- fft2D,- fft3D,--) where--import Data.Array.Accelerate.Math.FFT.Mode--import Data.Array.Accelerate as A-import Data.Array.Accelerate.Type as A-import Data.Array.Accelerate.Array.Sugar as S-import Data.Array.Accelerate.Error as A-import Data.Array.Accelerate.Array.Data as A-import Data.Array.Accelerate.Array.Unique as A-import Data.Array.Accelerate.Data.Complex as A--import Data.Array.Accelerate.LLVM.Native.Foreign--import Data.Ix ( Ix )-import Data.Array.CArray ( CArray )-import qualified Data.Array.CArray as C--import Math.FFT.Base ( FFTWReal, Sign(..), Flag, measure, destroyInput )-import qualified Math.FFT as FFT--import Foreign.Ptr-import Foreign.Storable-import Foreign.Storable.Complex ()--import Data.Bits-import Text.Printf-import Prelude as P---fft1D :: forall e. (Elt e, IsFloating e)- => Mode- -> ForeignAcc (Vector (Complex e) -> Vector (Complex e))-fft1D mode- = ForeignAcc (nameOf mode (undefined::DIM1))- $ case floatingType :: FloatingType e of- TypeFloat{} -> liftIO . liftAtoC go- TypeDouble{} -> liftIO . liftAtoC go- TypeCFloat{} -> liftIO . liftAtoC go- TypeCDouble{} -> liftIO . liftAtoC go- where- go :: FFTWReal r => CArray Int (Complex r) -> CArray Int (Complex r)- go = FFT.dftGU (signOf mode) flags [0]---fft2D :: forall e. (Elt e, IsFloating e)- => Mode- -> ForeignAcc (Array DIM2 (Complex e) -> Array DIM2 (Complex e))-fft2D mode- = ForeignAcc (nameOf mode (undefined::DIM2))- $ case floatingType :: FloatingType e of- TypeFloat{} -> liftIO . liftAtoC go- TypeDouble{} -> liftIO . liftAtoC go- TypeCFloat{} -> liftIO . liftAtoC go- TypeCDouble{} -> liftIO . liftAtoC go- where- go :: FFTWReal r => CArray (Int,Int) (Complex r) -> CArray (Int,Int) (Complex r)- go = FFT.dftGU (signOf mode) flags [0,1]---fft3D :: forall e. (Elt e, IsFloating e)- => Mode- -> ForeignAcc (Array DIM3 (Complex e) -> Array DIM3 (Complex e))-fft3D mode- = ForeignAcc (nameOf mode (undefined::DIM3))- $ case floatingType :: FloatingType e of- TypeFloat{} -> liftIO . liftAtoC go- TypeDouble{} -> liftIO . liftAtoC go- TypeCFloat{} -> liftIO . liftAtoC go- TypeCDouble{} -> liftIO . liftAtoC go- where- go :: FFTWReal r => CArray (Int,Int,Int) (Complex r) -> CArray (Int,Int,Int) (Complex r)- go = FFT.dftGU (signOf mode) flags [0,1,2]---signOf :: Mode -> Sign-signOf Forward = DFTForward-signOf _ = DFTBackward--flags :: Flag-flags = measure .|. destroyInput--nameOf :: forall sh. Shape sh => Mode -> sh -> String-nameOf Forward _ = printf "FFTW.dft%dD" (rank (undefined::sh))-nameOf _ _ = printf "FFTW.idft%dD" (rank (undefined::sh))----- | Lift an operation on CArray into an operation on Accelerate arrays----liftAtoC- :: (IxShapeRepr (EltRepr ix) ~ EltRepr sh, Shape sh, Ix ix, Elt ix, Elt e, IsFloating e, Storable e', ArrayPtrs e ~ Ptr e')- => (CArray ix (Complex e') -> CArray ix (Complex e'))- -> Array sh (Complex e)- -> IO (Array sh (Complex e))-liftAtoC f a = c2a . f =<< a2c a----- | Convert a multidimensional Accelerate array of complex numbers into--- a packed CArray of complex numbers suitable for use by FFTW.----a2c :: forall ix sh e e'. (IxShapeRepr (EltRepr ix) ~ EltRepr sh, Ix ix, Elt ix, Shape sh, IsFloating e, Storable e', ArrayPtrs e ~ Ptr e')- => Array sh (Complex e)- -> IO (CArray ix (Complex e'))-a2c arr- | FloatingDict <- floatingDict (floatingType :: FloatingType e)- = let- (lo,hi) = shapeToRange (arrayShape arr)- bnds = (fromIxShapeRepr lo, fromIxShapeRepr hi)- n = S.size (arrayShape arr)- in- C.createCArray bnds $ \p_cs ->- withComplexArrayPtrs arr $ \p_re p_im ->- let- -- TLM: Should we execute this in parallel using the worker threads of- -- the current target? (Native)- go !i | i P.>= n = return ()- go !i = do- re <- peekElemOff p_re i- im <- peekElemOff p_im i- pokeElemOff p_cs i (re :+ im)- go (i+1)- in- go 0----- | Convert a packed CArray of complex numbers into an unzipped (SoA)--- multidimensional Accelerate array of complex numbers.----c2a :: forall ix sh e e'. (IxShapeRepr (EltRepr ix) ~ EltRepr sh, Ix ix, Elt ix, Shape sh, Elt e, IsFloating e, Storable e', ArrayPtrs e ~ Ptr e')- => CArray ix (Complex e')- -> IO (Array sh (Complex e))-c2a carr- | FloatingDict <- floatingDict (floatingType :: FloatingType e)- = let- (lo,hi) = C.bounds carr- n = C.rangeSize (lo,hi)- sh = rangeToShape (toIxShapeRepr lo, toIxShapeRepr hi)- in do- arr <- allocateArray sh- C.withCArray carr $ \p_cs -> do- withComplexArrayPtrs arr $ \p_re p_im -> do- let- -- TLM: Should we execute this in parallel using the worker threads- -- of the current target? (Native)- go !i | i P.>= n = return ()- go !i = do- re :+ im <- peekElemOff p_cs i- pokeElemOff p_re i re- pokeElemOff p_im i im- go (i+1)- --- go 0- return arr----- Converting between Accelerate multidimensional shapes/indices and those used--- by the CArray package (Data.Ix)-----type family IxShapeRepr e where- IxShapeRepr () = ()- IxShapeRepr Int = ((),Int)- IxShapeRepr (t,h) = (IxShapeRepr t, h)--fromIxShapeRepr- :: forall ix sh. (IxShapeRepr (EltRepr ix) ~ EltRepr sh, Shape sh, Elt ix)- => sh- -> ix-fromIxShapeRepr = liftToElt (go (eltType (undefined::ix)))- where- go :: forall ix'. TupleType ix' -> IxShapeRepr ix' -> ix'- go UnitTuple () = ()- go (PairTuple tt _) (t, h) = (go tt t, h)- go (SingleTuple (NumScalarType (IntegralNumType TypeInt{}))) ((),h) = h- go _ _- = $internalError "fromIxShapeRepr" "expected Int dimensions"--toIxShapeRepr- :: forall ix sh. (IxShapeRepr (EltRepr ix) ~ EltRepr sh, Shape sh, Elt ix)- => ix- -> sh-toIxShapeRepr = liftToElt (go (eltType (undefined::ix)))- where- go :: forall ix'. TupleType ix' -> ix' -> IxShapeRepr ix'- go UnitTuple () = ()- go (SingleTuple (NumScalarType (IntegralNumType TypeInt{}))) h = ((), h)- go (PairTuple tt _) (t, h) = (go tt t, h)- go _ _- = error "toIxShapeRepr: not a valid Data.Ix index"----- Dig out the underlying pointers of the Accelerate SoA data structure-----withComplexArrayPtrs- :: forall sh e a. IsFloating e- => Array sh (Complex e)- -> (ArrayPtrs e -> ArrayPtrs e -> IO a)- -> IO a-withComplexArrayPtrs (Array _ adata) k- | AD_Pair (AD_Pair AD_Unit ad1) ad2 <- adata- = case floatingType :: FloatingType e of- TypeFloat{} -> withArrayData arrayElt ad1 $ \p1 -> withArrayData arrayElt ad2 $ \p2 -> k p1 p2- TypeDouble{} -> withArrayData arrayElt ad1 $ \p1 -> withArrayData arrayElt ad2 $ \p2 -> k p1 p2- TypeCFloat{} -> withArrayData arrayElt ad1 $ \p1 -> withArrayData arrayElt ad2 $ \p2 -> k p1 p2- TypeCDouble{} -> withArrayData arrayElt ad1 $ \p1 -> withArrayData arrayElt ad2 $ \p2 -> k p1 p2---- withScalarArrayPtrs--- :: forall sh e a. IsFloating e--- => Array sh e--- -> (ArrayPtrs e -> IO a)--- -> IO a--- withScalarArrayPtrs (Array _ adata) =--- case floatingType :: FloatingType e of--- TypeFloat{} -> withArrayData arrayElt adata--- TypeDouble{} -> withArrayData arrayElt adata--- TypeCFloat{} -> withArrayData arrayElt adata--- TypeCDouble{} -> withArrayData arrayElt adata--withArrayData- :: (ArrayPtrs e ~ Ptr a)- => ArrayEltR e- -> ArrayData e- -> (Ptr a -> IO b)- -> IO b-withArrayData ArrayEltRfloat (AD_Float ua) = withUniqueArrayPtr ua-withArrayData ArrayEltRdouble (AD_Double ua) = withUniqueArrayPtr ua-withArrayData ArrayEltRcfloat (AD_CFloat ua) = withUniqueArrayPtr ua-withArrayData ArrayEltRcdouble (AD_CDouble ua) = withUniqueArrayPtr ua-withArrayData _ _ =- $internalError "withArrayData" "expected array of [C]Float or [C]Double"-
− Data/Array/Accelerate/Math/FFT/LLVM/PTX.hs
@@ -1,279 +0,0 @@-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE GADTs #-}-{-# LANGUAGE PatternGuards #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TupleSections #-}-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE TypeOperators #-}-{-# LANGUAGE ViewPatterns #-}--- |--- Module : Data.Array.Accelerate.Math.FFT.LLVM.PTX--- Copyright : [2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell--- License : BSD3------ Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>--- Stability : experimental--- Portability : non-portable (GHC extensions)-----module Data.Array.Accelerate.Math.FFT.LLVM.PTX (-- fft1D,- fft2D,- fft3D,--) where--import Data.Array.Accelerate.Math.FFT.Mode-import Data.Array.Accelerate.Math.FFT.Twine--import Data.Array.Accelerate.Array.Data-import Data.Array.Accelerate.Lifetime-import Data.Array.Accelerate.Array.Sugar-import Data.Array.Accelerate.Data.Complex-import Data.Array.Accelerate.Type--import Data.Array.Accelerate.LLVM.PTX.Foreign--import Foreign.CUDA.Ptr ( DevicePtr )-import Foreign.Ptr-import Foreign.Storable-import Foreign.CUDA.Analysis-import qualified Foreign.CUDA.FFT as FFT-import qualified Foreign.CUDA.Driver as CUDA hiding ( device )-import qualified Foreign.CUDA.Driver.Context as CUDA ( device )--import Control.Concurrent.MVar-import Control.Exception-import Control.Monad-import Data.Maybe-import Data.Typeable-import System.IO.Unsafe---fft1D :: IsFloating e- => Mode- -> ForeignAcc (Vector (Complex e) -> (Vector (Complex e)))-fft1D mode = ForeignAcc "fft1D" $ liftAtoC (cuFFT mode)--fft2D :: IsFloating e- => Mode- -> ForeignAcc (Array DIM2 (Complex e) -> (Array DIM2 (Complex e)))-fft2D mode = ForeignAcc "fft2D" $ liftAtoC (cuFFT mode)--fft3D :: IsFloating e- => Mode- -> ForeignAcc (Array DIM3 (Complex e) -> (Array DIM3 (Complex e)))-fft3D mode = ForeignAcc "fft3D" $ liftAtoC (cuFFT mode)---liftAtoC- :: forall sh e. (Shape sh, IsFloating e)- => (Stream -> Array (sh:.Int) e -> LLVM PTX (Array (sh:.Int) e))- -> Stream- -> Array (sh:.Int) (Complex e)- -> LLVM PTX (Array (sh:.Int) (Complex e))-liftAtoC f s =- case floatingType :: FloatingType e of- TypeFloat{} -> c2a s <=< f s <=< a2c s- TypeDouble{} -> c2a s <=< f s <=< a2c s- TypeCFloat{} -> c2a s <=< f s <=< a2c s- TypeCDouble{} -> c2a s <=< f s <=< a2c s----- | Call the cuFFT library to execute the FFT (inplace)----cuFFT :: forall sh e. (Shape sh, IsFloating e)- => Mode- -> Stream- -> Array (sh:.Int) e- -> LLVM PTX (Array (sh:.Int) e)-cuFFT mode stream arr =- withScalarArrayPtr arr stream $ \d_arr -> liftIO $- withLifetime stream $ \st -> do- let sh :. sz = shape arr- p <- plan (sh :. sz `quot` 2) (undefined::e) -- recall this is an array of packed (Vec2 e)- FFT.setStream p st- case floatingType :: FloatingType e of- TypeFloat{} -> FFT.execC2C p d_arr d_arr (signOfMode mode) >> return arr- TypeDouble{} -> FFT.execZ2Z p d_arr d_arr (signOfMode mode) >> return arr- TypeCFloat{} -> FFT.execC2C p d_arr d_arr (signOfMode mode) >> return arr- TypeCDouble{} -> FFT.execZ2Z p d_arr d_arr (signOfMode mode) >> return arr----- | Convert an unzipped Accelerate array of complex numbers into a (new) packed--- array suitable for use with CUFFT.----a2c :: forall sh e. (Shape sh, Elt e, IsFloating e, Storable (DevicePtrs e))- => Stream- -> Array (sh:.Int) (Complex e)- -> LLVM PTX (Array (sh:.Int) e) -- this is really a packed array of (Vec2 e) type-a2c stream arr | FloatingDict <- floatingDict (floatingType :: FloatingType e) = do- let- sh :. sz = shape arr- n = size sh * sz- --- cs <- allocateRemote (sh :. 2*sz)- withComplexArrayPtrs arr stream $ \d_re d_im -> do- withScalarArrayPtr cs stream $ \d_cs -> liftIO $ do- withLifetime stream $ \st -> do- mdl <- twine (sizeOf (undefined::e))- pack <- CUDA.getFun mdl "interleave"- dev <- CUDA.device- prp <- CUDA.props dev- regs <- CUDA.requires pack CUDA.NumRegs- let- blockSize = 256- sharedMem = 0- maxBlocks = maxResidentBlocks prp blockSize regs sharedMem- numBlocks = maxBlocks `min` ((n + blockSize - 1) `div` blockSize)- --- CUDA.launchKernel pack (numBlocks,1,1) (blockSize,1,1) sharedMem (Just st)- [ CUDA.VArg d_cs, CUDA.VArg d_re, CUDA.VArg d_im, CUDA.IArg (fromIntegral n) ]- return cs---- | Convert a packed array of complex numbers into a (new) unzipped Accelerate--- array.----c2a :: forall sh e. (Shape sh, Elt e, IsFloating e, Storable (DevicePtrs e))- => Stream- -> Array (sh:.Int) e- -> LLVM PTX (Array (sh:.Int) (Complex e))-c2a stream cs | FloatingDict <- floatingDict (floatingType :: FloatingType e) = do- let- sh :. sz2 = shape cs- sz = sz2 `quot` 2- n = size sh * sz- --- arr <- allocateRemote (sh :. sz)- withComplexArrayPtrs arr stream $ \d_re d_im -> do- withScalarArrayPtr cs stream $ \d_cs -> liftIO $ do- withLifetime stream $ \st -> do- mdl <- twine (sizeOf (undefined::e))- unpack <- CUDA.getFun mdl "deinterleave"- dev <- CUDA.device- prp <- CUDA.props dev- regs <- CUDA.requires unpack CUDA.NumRegs- let- blockSize = 256- sharedMem = 0- maxBlocks = maxResidentBlocks prp blockSize regs sharedMem- numBlocks = maxBlocks `min` ((n + blockSize - 1) `div` blockSize)- --- CUDA.launchKernel unpack (numBlocks,1,1) (blockSize,1,1) sharedMem (Just st)- [ CUDA.VArg d_re, CUDA.VArg d_im, CUDA.VArg d_cs, CUDA.IArg (fromIntegral n) ]- return arr----- | Generate an execute plan for a given type and size of FFT. These plans are--- cached so that subsequent invocations are quicker.----plan :: forall sh e. (Shape sh, IsFloating e) => sh -> e -> IO FFT.Handle-plan (shapeToList -> sh) _ =- modifyMVar fft_plans $ \ps ->- case lookup (ty, sh) ps of- Just p -> return (ps, p)- Nothing -> do- p <- case sh of- [w] -> FFT.plan1D w ty 1- [w,h] -> FFT.plan2D h w ty- [w,h,d] -> FFT.plan3D d h w ty- _ -> error "cuFFT only supports 1D, 2D, and 3D transforms"- return (((ty,sh),p) : ps, p)- where- ty = case floatingType :: FloatingType e of- TypeFloat{} -> FFT.C2C- TypeDouble{} -> FFT.Z2Z- TypeCFloat{} -> FFT.C2C- TypeCDouble{} -> FFT.Z2Z----- | Load the module to convert between SoA and AoS representation for the given--- type. This is cached for subsequent reuse.----twine :: Int -> IO CUDA.Module-twine bitsize = do- ctx <- fromMaybe (error "could not determine current CUDA context") `fmap` CUDA.get- modifyMVar ptx_twine_modules $ \ms -> do- case lookup (bitsize,ctx) ms of- Just m -> return (ms, m)- Nothing -> do- m <- CUDA.loadData $ case bitsize of- 4 -> ptx_twine_f32- 8 -> ptx_twine_f64- _ -> error "cuFFT only supports Float and Double"- return (((bitsize,ctx), m) : ms, m)----- | Dig out the two device pointers for an unzipped array of complex numbers.----withComplexArrayPtrs- :: forall sh e a. IsFloating e- => Array sh (Complex e)- -> Stream- -> (DevicePtrs e -> DevicePtrs e -> LLVM PTX a)- -> LLVM PTX a-withComplexArrayPtrs (Array _ adata) st k- | AD_Pair (AD_Pair AD_Unit ad1) ad2 <- adata- = case floatingType :: FloatingType e of- TypeFloat{} -> withArrayData arrayElt ad1 st $ \p1 -> withArrayData arrayElt ad2 st $ \p2 -> k p1 p2- TypeDouble{} -> withArrayData arrayElt ad1 st $ \p1 -> withArrayData arrayElt ad2 st $ \p2 -> k p1 p2- TypeCDouble{} -> withArrayData arrayElt ad1 st $ \p1 -> withArrayData arrayElt ad2 st $ \p2 -> k p1 p2- TypeCFloat{} -> withArrayData arrayElt ad1 st $ \p1 -> withArrayData arrayElt ad2 st $ \p2 -> k p1 p2---- | Dig out the device pointer for a scalar array----withScalarArrayPtr- :: forall sh e a. IsFloating e- => Array sh e- -> Stream- -> (DevicePtrs e -> LLVM PTX a)- -> LLVM PTX a-withScalarArrayPtr (Array _ ad) st k- = case floatingType :: FloatingType e of- TypeFloat{} -> withArrayData arrayElt ad st $ \p -> k p- TypeDouble{} -> withArrayData arrayElt ad st $ \p -> k p- TypeCDouble{} -> withArrayData arrayElt ad st $ \p -> k p- TypeCFloat{} -> withArrayData arrayElt ad st $ \p -> k p--withArrayData- :: (Typeable e, Typeable a, ArrayElt e, Storable a, ArrayPtrs e ~ Ptr a)- => ArrayEltR e- -> ArrayData e- -> Stream- -> (DevicePtr a -> LLVM PTX b)- -> LLVM PTX b-withArrayData _ ad s k =- withDevicePtr ad $ \p -> do- r <- k p- e <- checkpoint s- return (Just e,r)--type family DevicePtrs e :: *--type instance DevicePtrs Float = DevicePtr Float-type instance DevicePtrs Double = DevicePtr Double-type instance DevicePtrs CFloat = DevicePtr Float-type instance DevicePtrs CDouble = DevicePtr Double----- Cache the FFT planning step for faster repeat evaluations.-{-# NOINLINE fft_plans #-}-fft_plans :: MVar [((FFT.Type, [Int]), FFT.Handle)]-fft_plans = unsafePerformIO $ do- mv <- newMVar []- _ <- mkWeakMVar mv- $ withMVar mv- $ mapM_ (\(_,p) -> FFT.destroy p)- return mv---- Cache the functions which convert between SoA and AoS format.-{-# NOINLINE ptx_twine_modules #-}-ptx_twine_modules :: MVar [((Int, CUDA.Context), CUDA.Module)]-ptx_twine_modules = unsafePerformIO $ do- mv <- newMVar []- _ <- mkWeakMVar mv- $ withMVar mv- $ mapM_ (\((_,ctx),mdl) -> bracket_ (CUDA.push ctx) CUDA.pop (CUDA.unload mdl))- return mv-
− Data/Array/Accelerate/Math/FFT/Mode.hs
@@ -1,28 +0,0 @@--- |--- Module : Data.Array.Accelerate.Math.FFT.Mode--- Copyright : [2012..2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell--- [2013..2017] Robert Clifton-Everest--- License : BSD3------ Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>--- Stability : experimental--- Portability : non-portable (GHC extensions)-----module Data.Array.Accelerate.Math.FFT.Mode- where---data Mode- = Forward -- ^ Forward DFT- | Reverse -- ^ Inverse DFT, un-normalised- | Inverse -- ^ Inverse DFT, normalised- deriving (Eq, Show)--signOfMode :: Num a => Mode -> a-signOfMode m- = case m of- Forward -> -1- Reverse -> 1- Inverse -> 1-
− Data/Array/Accelerate/Math/FFT/Twine.hs
@@ -1,88 +0,0 @@-{-# LANGUAGE CPP #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TemplateHaskell #-}--- |--- Module : Data.Array.Accelerate.Math.FFT.Twine--- Copyright : [2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell--- License : BSD3------ Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>--- Stability : experimental--- Portability : non-portable (GHC extensions)--------module Data.Array.Accelerate.Math.FFT.Twine- where--import Data.Array.Accelerate as A-import Data.Array.Accelerate.Data.Complex--#ifdef ACCELERATE_LLVM_PTX_BACKEND-import Data.FileEmbed-import Data.ByteString ( ByteString )-#endif----- Interleave the real and imaginary components in a complex array and produce a--- flattened vector. This allows us to mimic the array-of-struct representation--- commonly used by FFT libraries to store complex numbers (CUFFT, FFTW).------ We would really prefer to implement this with a zipWith of the two arrays,--- but we can't represent the packed structure in Accelerate.----{-# NOINLINE interleave #-}-interleave :: Elt e => Acc (Vector (Complex e)) -> Acc (Vector e)-interleave arr = generate sh swizzle- where- reals = A.map real arr- imags = A.map imag arr- --- sh = index1 (2 * A.size arr)- swizzle ix =- let i = indexHead ix- (j,k) = i `quotRem` 2- in- k A.== 0 ? ( reals A.!! j, imags A.!! j )----- Deinterleave a vector into a complex array. Requires the array to have an--- even number of elements.----{-# NOINLINE deinterleave #-}-deinterleave :: forall e. Elt e => Acc (Vector e) -> Acc (Vector (Complex e))-deinterleave arr = generate sh swizzle- where- sh = index1 (A.size arr `quot` 2)- swizzle ix =- let i = indexHead ix `quot` 2- in lift ( arr A.!! i :+ arr A.!! (i+1) ) :: Exp (Complex e)---{-# RULES- "interleave/deinterleave" forall x. deinterleave (interleave x) = x;- "deinterleave/interleave" forall x. interleave (deinterleave x) = x- #-}---#ifdef ACCELERATE_LLVM_PTX_BACKEND---- Embedded PTX code for interleave and deinterleave for 32- and 64-bit floating--- point numbers respectively. These can be loaded and executed by the CUDA--- driver at runtime as required.------ The PTX code was compiled for SM-2.0 and 64-bit address space (the default--- settings of nvcc-7.5), but the code is simple enough that the CUDA device--- driver should be able to compile it for the actual target architecture--- without issue. This has been confirmed with respect to SM, but I don't have--- a 32-bit machine available to test that aspect with.-----ptx_twine_f32 :: ByteString-ptx_twine_f32 = $(makeRelativeToProject "cubits/twine_f32.ptx" >>= embedFile)--ptx_twine_f64 :: ByteString-ptx_twine_f64 = $(makeRelativeToProject "cubits/twine_f64.ptx" >>= embedFile)--#endif-
accelerate-fft.cabal view
@@ -1,6 +1,6 @@ Name: accelerate-fft-Version: 1.1.0.0-Cabal-version: >= 1.6+Version: 1.2.0.0+Cabal-version: >= 1.8 Tested-with: GHC >= 7.10 Build-type: Simple @@ -31,25 +31,22 @@ extra-source-files: README.md CHANGELOG.md- cubits/twine_f32.ptx- cubits/twine_f64.ptx- cubits/twine_f32.cu- cubits/twine_f64.cu -Flag llvm-ptx+flag llvm-ptx Description: Use CUFFT-based implementation in the LLVM.PTX backend Default: True -Flag llvm-cpu+flag llvm-cpu Description: Use FFTW-based implementation in the LLVM.Native backend Default: True -Library+library build-depends:- base >= 4.7 && < 4.11- , accelerate >= 1.0 && < 1.2+ base >= 4.7 && < 4.12+ , accelerate >= 1.2 && < 1.3 , bytestring >= 0.9+ , lens-accelerate >= 0.2 exposed-modules: Data.Array.Accelerate.Math.FFT@@ -58,45 +55,43 @@ Data.Array.Accelerate.Math.DFT.Roots other-modules:+ Data.Array.Accelerate.Math.FFT.Adhoc Data.Array.Accelerate.Math.FFT.Mode- Data.Array.Accelerate.Math.FFT.Twine+ Data.Array.Accelerate.Math.FFT.Type + hs-source-dirs: src ghc-options: -O2 -Wall -funbox-strict-fields - -- if flag(cuda)- -- cpp-options: -DACCELERATE_CUDA_BACKEND- -- build-depends:- -- accelerate-cuda >= 0.16- -- , cuda >= 0.5- -- , cufft >= 0.1.2- -- , file-embed >= 0.0.10- --- -- other-modules:- -- Data.Array.Accelerate.Math.FFT.CUDA- if flag(llvm-cpu) cpp-options: -DACCELERATE_LLVM_NATIVE_BACKEND build-depends:- accelerate-llvm >= 1.0 && < 1.2- , accelerate-llvm-native >= 1.0 && < 1.2+ accelerate-llvm >= 1.0 && < 1.3+ , accelerate-llvm-native >= 1.0 && < 1.3 , carray >= 0.1.5 , fft >= 0.1.8- , storable-complex >= 0.2 other-modules: Data.Array.Accelerate.Math.FFT.LLVM.Native+ Data.Array.Accelerate.Math.FFT.LLVM.Native.Base+ Data.Array.Accelerate.Math.FFT.LLVM.Native.Ix if flag(llvm-ptx) cpp-options: -DACCELERATE_LLVM_PTX_BACKEND build-depends:- accelerate-llvm >= 1.0 && < 1.2- , accelerate-llvm-ptx >= 1.0 && < 1.2+ accelerate-llvm >= 1.0 && < 1.3+ , accelerate-llvm-ptx >= 1.0 && < 1.3+ , containers >= 0.5+ , hashable >= 1.0+ , unordered-containers >= 0.2 , cuda >= 0.5- , cufft >= 0.1.2+ , cufft >= 0.9 , file-embed >= 0.0.10+ , mtl >= 2.2 other-modules: Data.Array.Accelerate.Math.FFT.LLVM.PTX+ Data.Array.Accelerate.Math.FFT.LLVM.PTX.Base+ Data.Array.Accelerate.Math.FFT.LLVM.PTX.Plans -- Don't add the extensions list here. Instead, place individual LANGUAGE -- pragmas in the files that require a specific extension. This means the@@ -104,14 +99,72 @@ -- -- Extensions: -Source-repository head++test-suite test-llvm-native+ type: exitcode-stdio-1.0+ hs-source-dirs: test+ main-is: TestNative.hs+ ghc-options: -main-is TestNative++ if !flag(llvm-cpu)+ buildable: False++ build-depends:+ base >= 4.7 && < 4.12+ , accelerate+ , accelerate-fft+ , accelerate-llvm-native+ , hedgehog >= 0.5+ , tasty >= 0.11+ , tasty-hedgehog >= 0.1++ ghc-options:+ -Wall+ -threaded+ -rtsopts++ other-modules:+ Test.FFT+ Test.Base+ Test.ShowType+++test-suite test-llvm-ptx+ type: exitcode-stdio-1.0+ hs-source-dirs: test+ main-is: TestPTX.hs+ ghc-options: -main-is TestPTX++ if !flag(llvm-ptx)+ buildable: False++ build-depends:+ base >= 4.7 && < 4.12+ , accelerate+ , accelerate-fft+ , accelerate-llvm-ptx+ , hedgehog >= 0.5+ , tasty >= 0.11+ , tasty-hedgehog >= 0.1++ ghc-options:+ -Wall+ -threaded+ -rtsopts++ other-modules:+ Test.FFT+ Test.Base+ Test.ShowType+++source-repository head Type: git Location: git://github.com/AccelerateHS/accelerate-fft.git -Source-repository this+source-repository this Type: git- Tag: 1.1.0.0+ Tag: 1.2.0.0 Location: git://github.com/AccelerateHS/accelerate-fft.git -- vim: nospell-
− cubits/twine_f32.cu
@@ -1,63 +0,0 @@-/*- * Module : Twine- * Copyright : [2016] Trevor L. McDonell- * License : BSD3- *- * Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>- * Stability : experimental- * Portability : non-portable (GHC extensions)- *- * Convert between Accelerate's Struct-of-Array representation of complex- * numbers and the Array-of-Struct representation necessary for CUBLAS.- *- */--#include <cuda.h>-#include <cuComplex.h>--#ifdef __cplusplus-extern "C" {-#endif--__global__ void interleave-(- cuFloatComplex * __restrict__ cplx,- const float * __restrict__ real,- const float * __restrict__ imag,- const int size-)-{- const int gridSize = blockDim.x * gridDim.x;- int ix;-- for (ix = blockDim.x * blockIdx.x + threadIdx.x; ix < size; ix += gridSize) {- const float re = real[ix];- const float im = imag[ix];-- cplx[ix] = make_cuFloatComplex(re, im);- }-}--__global__ void deinterleave-(- float * __restrict__ real,- float * __restrict__ imag,- const cuFloatComplex * __restrict__ cplx,- const int size-)-{- const int gridSize = blockDim.x * gridDim.x;- int ix;-- for (ix = blockDim.x * blockIdx.x + threadIdx.x; ix < size; ix += gridSize) {- const cuFloatComplex c = cplx[ix];-- real[ix] = cuCrealf(c);- imag[ix] = cuCimagf(c);- }-}--#ifdef __cplusplus-}-#endif-
− cubits/twine_f32.ptx
@@ -1,108 +0,0 @@-//-// Generated by NVIDIA NVVM Compiler-//-// Compiler Build ID: CL-20633761-// Cuda compilation tools, release 7.5, V7.5.26-// Based on LLVM 3.4svn-//--.version 4.3-.target sm_20-.address_size 64-- // .globl interleave--.visible .entry interleave(- .param .u64 interleave_param_0,- .param .u64 interleave_param_1,- .param .u64 interleave_param_2,- .param .u32 interleave_param_3-)-{- .reg .pred %p<3>;- .reg .f32 %f<3>;- .reg .b32 %r<11>;- .reg .b64 %rd<12>;--- ld.param.u64 %rd4, [interleave_param_0];- ld.param.u64 %rd5, [interleave_param_1];- ld.param.u64 %rd6, [interleave_param_2];- ld.param.u32 %r5, [interleave_param_3];- cvta.to.global.u64 %rd1, %rd4;- cvta.to.global.u64 %rd2, %rd6;- cvta.to.global.u64 %rd3, %rd5;- mov.u32 %r6, %nctaid.x;- mov.u32 %r7, %ntid.x;- mul.lo.s32 %r1, %r6, %r7;- mov.u32 %r8, %ctaid.x;- mov.u32 %r9, %tid.x;- mad.lo.s32 %r10, %r8, %r7, %r9;- setp.ge.s32 %p1, %r10, %r5;- @%p1 bra BB0_2;--BB0_1:- mul.wide.s32 %rd7, %r10, 4;- add.s64 %rd8, %rd3, %rd7;- add.s64 %rd9, %rd2, %rd7;- mul.wide.s32 %rd10, %r10, 8;- add.s64 %rd11, %rd1, %rd10;- ld.global.f32 %f1, [%rd9];- ld.global.f32 %f2, [%rd8];- st.global.v2.f32 [%rd11], {%f2, %f1};- add.s32 %r10, %r10, %r1;- setp.lt.s32 %p2, %r10, %r5;- @%p2 bra BB0_1;--BB0_2:- ret;-}-- // .globl deinterleave-.visible .entry deinterleave(- .param .u64 deinterleave_param_0,- .param .u64 deinterleave_param_1,- .param .u64 deinterleave_param_2,- .param .u32 deinterleave_param_3-)-{- .reg .pred %p<3>;- .reg .f32 %f<5>;- .reg .b32 %r<11>;- .reg .b64 %rd<12>;--- ld.param.u64 %rd4, [deinterleave_param_0];- ld.param.u64 %rd5, [deinterleave_param_1];- ld.param.u64 %rd6, [deinterleave_param_2];- ld.param.u32 %r5, [deinterleave_param_3];- cvta.to.global.u64 %rd1, %rd5;- cvta.to.global.u64 %rd2, %rd4;- cvta.to.global.u64 %rd3, %rd6;- mov.u32 %r6, %nctaid.x;- mov.u32 %r7, %ntid.x;- mul.lo.s32 %r1, %r6, %r7;- mov.u32 %r8, %ctaid.x;- mov.u32 %r9, %tid.x;- mad.lo.s32 %r10, %r8, %r7, %r9;- setp.ge.s32 %p1, %r10, %r5;- @%p1 bra BB1_2;--BB1_1:- mul.wide.s32 %rd7, %r10, 8;- add.s64 %rd8, %rd3, %rd7;- ld.global.v2.f32 {%f1, %f2}, [%rd8];- mul.wide.s32 %rd9, %r10, 4;- add.s64 %rd10, %rd2, %rd9;- st.global.f32 [%rd10], %f1;- add.s64 %rd11, %rd1, %rd9;- st.global.f32 [%rd11], %f2;- add.s32 %r10, %r10, %r1;- setp.lt.s32 %p2, %r10, %r5;- @%p2 bra BB1_1;--BB1_2:- ret;-}--
− cubits/twine_f64.cu
@@ -1,63 +0,0 @@-/*- * Module : Twine- * Copyright : [2016] Trevor L. McDonell- * License : BSD3- *- * Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>- * Stability : experimental- * Portability : non-portable (GHC extensions)- *- * Convert between Accelerate's Struct-of-Array representation of complex- * numbers and the Array-of-Struct representation necessary for CUBLAS.- *- */--#include <cuda.h>-#include <cuComplex.h>--#ifdef __cplusplus-extern "C" {-#endif--__global__ void interleave-(- cuDoubleComplex * __restrict__ cplx,- const double * __restrict__ real,- const double * __restrict__ imag,- const int size-)-{- const int gridSize = blockDim.x * gridDim.x;- int ix;-- for (ix = blockDim.x * blockIdx.x + threadIdx.x; ix < size; ix += gridSize) {- const double re = real[ix];- const double im = imag[ix];-- cplx[ix] = make_cuDoubleComplex(re, im);- }-}--__global__ void deinterleave-(- double * __restrict__ real,- double * __restrict__ imag,- const cuDoubleComplex * __restrict__ cplx,- const int size-)-{- const int gridSize = blockDim.x * gridDim.x;- int ix;-- for (ix = blockDim.x * blockIdx.x + threadIdx.x; ix < size; ix += gridSize) {- const cuDoubleComplex c = cplx[ix];-- real[ix] = cuCreal(c);- imag[ix] = cuCimag(c);- }-}--#ifdef __cplusplus-}-#endif-
− cubits/twine_f64.ptx
@@ -1,108 +0,0 @@-//-// Generated by NVIDIA NVVM Compiler-//-// Compiler Build ID: CL-20633761-// Cuda compilation tools, release 7.5, V7.5.26-// Based on LLVM 3.4svn-//--.version 4.3-.target sm_20-.address_size 64-- // .globl interleave--.visible .entry interleave(- .param .u64 interleave_param_0,- .param .u64 interleave_param_1,- .param .u64 interleave_param_2,- .param .u32 interleave_param_3-)-{- .reg .pred %p<3>;- .reg .b32 %r<11>;- .reg .f64 %fd<3>;- .reg .b64 %rd<12>;--- ld.param.u64 %rd4, [interleave_param_0];- ld.param.u64 %rd5, [interleave_param_1];- ld.param.u64 %rd6, [interleave_param_2];- ld.param.u32 %r5, [interleave_param_3];- cvta.to.global.u64 %rd1, %rd4;- cvta.to.global.u64 %rd2, %rd6;- cvta.to.global.u64 %rd3, %rd5;- mov.u32 %r6, %nctaid.x;- mov.u32 %r7, %ntid.x;- mul.lo.s32 %r1, %r6, %r7;- mov.u32 %r8, %ctaid.x;- mov.u32 %r9, %tid.x;- mad.lo.s32 %r10, %r8, %r7, %r9;- setp.ge.s32 %p1, %r10, %r5;- @%p1 bra BB0_2;--BB0_1:- mul.wide.s32 %rd7, %r10, 8;- add.s64 %rd8, %rd3, %rd7;- add.s64 %rd9, %rd2, %rd7;- mul.wide.s32 %rd10, %r10, 16;- add.s64 %rd11, %rd1, %rd10;- ld.global.f64 %fd1, [%rd9];- ld.global.f64 %fd2, [%rd8];- st.global.v2.f64 [%rd11], {%fd2, %fd1};- add.s32 %r10, %r10, %r1;- setp.lt.s32 %p2, %r10, %r5;- @%p2 bra BB0_1;--BB0_2:- ret;-}-- // .globl deinterleave-.visible .entry deinterleave(- .param .u64 deinterleave_param_0,- .param .u64 deinterleave_param_1,- .param .u64 deinterleave_param_2,- .param .u32 deinterleave_param_3-)-{- .reg .pred %p<3>;- .reg .b32 %r<11>;- .reg .f64 %fd<5>;- .reg .b64 %rd<12>;--- ld.param.u64 %rd4, [deinterleave_param_0];- ld.param.u64 %rd5, [deinterleave_param_1];- ld.param.u64 %rd6, [deinterleave_param_2];- ld.param.u32 %r5, [deinterleave_param_3];- cvta.to.global.u64 %rd1, %rd5;- cvta.to.global.u64 %rd2, %rd4;- cvta.to.global.u64 %rd3, %rd6;- mov.u32 %r6, %nctaid.x;- mov.u32 %r7, %ntid.x;- mul.lo.s32 %r1, %r6, %r7;- mov.u32 %r8, %ctaid.x;- mov.u32 %r9, %tid.x;- mad.lo.s32 %r10, %r8, %r7, %r9;- setp.ge.s32 %p1, %r10, %r5;- @%p1 bra BB1_2;--BB1_1:- mul.wide.s32 %rd7, %r10, 16;- add.s64 %rd8, %rd3, %rd7;- ld.global.v2.f64 {%fd1, %fd2}, [%rd8];- mul.wide.s32 %rd9, %r10, 8;- add.s64 %rd10, %rd2, %rd9;- st.global.f64 [%rd10], %fd1;- add.s64 %rd11, %rd1, %rd9;- st.global.f64 [%rd11], %fd2;- add.s32 %r10, %r10, %r1;- setp.lt.s32 %p2, %r10, %r5;- @%p2 bra BB1_1;--BB1_2:- ret;-}--
+ src/Data/Array/Accelerate/Math/DFT.hs view
@@ -0,0 +1,113 @@+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}+-- |+-- Module : Data.Array.Accelerate.Math.DFT+-- Copyright : [2012..2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- [2013..2017] Robert Clifton-Everest+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@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.Data.Complex+++-- | Compute the DFT along the low order dimension of an array+--+dft :: (Shape sh, Slice sh, Elt (Complex e), A.RealFloat e, A.FromIntegral Int 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 (Complex e), A.RealFloat e, A.FromIntegral Int 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 :+ 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 (Complex e), A.RealFloat 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 (+) 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 (Complex e), A.RealFloat 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 (+) 0 $ A.zipWith (*) arr roots'+
+ src/Data/Array/Accelerate/Math/DFT/Centre.hs view
@@ -0,0 +1,108 @@+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE TypeOperators #-}+-- |+-- Module : Data.Array.Accelerate.Math.DFT.Centre+-- Copyright : [2012..2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- [2013..2017] Robert Clifton-Everest+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@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 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.Data.Complex+++-- | Apply the centring transform to a vector+--+centre1D+ :: (Elt (Complex e), A.RealFloat e, A.FromIntegral Int 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) :+ 0) * arr!ix)++-- | Apply the centring transform to a matrix+--+centre2D+ :: (Elt (Complex e), A.RealFloat e, A.FromIntegral Int 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)) :+ 0) * arr!ix)++-- | Apply the centring transform to a 3D array+--+centre3D+ :: (Elt (Complex e), A.RealFloat e, A.FromIntegral Int 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)) :+ 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 A.< 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 A.< mh ? (y + mh, y - mh))+ (x A.< 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 A.< md ? (z + md, z - md))+ (y A.< mh ? (y + mh, y - mh))+ (x A.< 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)+
+ src/Data/Array/Accelerate/Math/DFT/Roots.hs view
@@ -0,0 +1,53 @@+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE TypeOperators #-}+-- |+-- Module : Data.Array.Accelerate.Math.DFT.Roots+-- Copyright : [2012..2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- [2013..2017] Robert Clifton-Everest+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@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.Data.Complex+++-- | Calculate the roots of unity for the forward transform+--+rootsOfUnity+ :: (Shape sh, Slice sh, Elt (Complex e), A.Floating e, A.FromIntegral Int e)+ => 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+ :: (Shape sh, Slice sh, Elt (Complex e), A.Floating e, A.FromIntegral Int e)+ => 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 ))+
+ src/Data/Array/Accelerate/Math/FFT.hs view
@@ -0,0 +1,258 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE EmptyDataDecls #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ViewPatterns #-}+-- |+-- Module : Data.Array.Accelerate.Math.FFT+-- Copyright : [2012..2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- [2013..2017] Robert Clifton-Everest+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--+-- For performance, compile against the foreign library bindings (using any+-- number of '-fllvm-ptx', and '-fllvm-cpu' for the accelerate-llvm-ptx, and+-- accelerate-llvm-native backends, respectively).+--++module Data.Array.Accelerate.Math.FFT (++ Mode(..),+ Numeric,+ fft,++ fft1D,+ fft2D,+ fft3D,++) where++import Data.Array.Accelerate as A+import Data.Array.Accelerate.Data.Complex+import Data.Array.Accelerate.Math.FFT.Type+import Data.Array.Accelerate.Math.FFT.Mode+import qualified Data.Array.Accelerate.Array.Sugar as A ( rank )+import qualified Data.Array.Accelerate.Math.FFT.Adhoc as Adhoc++#ifdef ACCELERATE_LLVM_NATIVE_BACKEND+import qualified Data.Array.Accelerate.Math.FFT.LLVM.Native as Native+#endif+#ifdef ACCELERATE_LLVM_PTX_BACKEND+import qualified Data.Array.Accelerate.Math.FFT.LLVM.PTX as PTX+#endif++import Prelude as P+++-- | Discrete Fourier Transform along the innermost dimension of an array.+--+-- Notes for FFI implementations:+--+-- * fftw supports arrays of dimension 1-5+-- * cuFFT supports arrays of dimension 1-3+--+-- The pure implementation will be used otherwise.+--+fft :: forall sh e. (Shape sh, Slice sh, Numeric e)+ => Mode+ -> Acc (Array (sh:.Int) (Complex e))+ -> Acc (Array (sh:.Int) (Complex e))+fft mode arr+ = let+ scale = A.fromIntegral (indexHead (shape arr))+ rank = A.rank (undefined :: sh:.Int)+ go =+#ifdef ACCELERATE_LLVM_NATIVE_BACKEND+ (if rank P.<= 5 then foreignAcc (Native.fft mode) else id) $+#endif+#ifdef ACCELERATE_LLVM_PTX_BACKEND+ (if rank P.<= 3 then foreignAcc (PTX.fft mode) else id) $+#endif+ Adhoc.fft mode+ in+ case mode of+ Inverse -> A.map (/scale) (go arr)+ _ -> go arr+++-- Vector Transform+-- ----------------++-- | Discrete Fourier Transform of a vector.+--+fft1D :: forall e. Numeric e+ => Mode+ -> Acc (Array DIM1 (Complex e))+ -> Acc (Array DIM1 (Complex e))+fft1D mode arr+ = let+ scale = A.fromIntegral (A.length arr)+ go =+#ifdef ACCELERATE_LLVM_NATIVE_BACKEND+ foreignAcc (Native.fft1D mode) $+#endif+#ifdef ACCELERATE_LLVM_PTX_BACKEND+ foreignAcc (PTX.fft1D mode) $+#endif+ Adhoc.fft mode+ in+ case mode of+ Inverse -> A.map (/scale) (go arr)+ _ -> go arr+++-- Matrix Transform+-- ----------------++-- | Discrete Fourier Transform of a matrix.+--+fft2D :: forall e. Numeric e+ => Mode+ -> Acc (Array DIM2 (Complex e))+ -> Acc (Array DIM2 (Complex e))+fft2D mode arr+ = let+ scale = A.fromIntegral (A.size arr)+ go =+#ifdef ACCELERATE_LLVM_NATIVE_BACKEND+ foreignAcc (Native.fft2D mode) $+#endif+#ifdef ACCELERATE_LLVM_PTX_BACKEND+ foreignAcc (PTX.fft2D mode) $+#endif+ fft'++ fft' a = A.transpose . Adhoc.fft mode+ >-> A.transpose . Adhoc.fft mode+ $ a+ in+ case mode of+ Inverse -> A.map (/scale) (go arr)+ _ -> go arr+++-- Cube Transform+-- --------------++-- | Discrete Fourier Transform of a 3D array.+--+fft3D :: forall e. Numeric e+ => Mode+ -> Acc (Array DIM3 (Complex e))+ -> Acc (Array DIM3 (Complex e))+fft3D mode arr+ = let scale = A.fromIntegral (A.size arr)+ go =+#ifdef ACCELERATE_LLVM_NATIVE_BACKEND+ foreignAcc (Native.fft3D mode) $+#endif+#ifdef ACCELERATE_LLVM_PTX_BACKEND+ foreignAcc (PTX.fft3D mode) $+#endif+ fft'++ fft' a = rotate3D . Adhoc.fft mode+ >-> rotate3D . Adhoc.fft mode+ >-> rotate3D . Adhoc.fft mode+ $ a+ in+ case mode of+ Inverse -> A.map (/scale) (go arr)+ _ -> go arr+++rotate3D :: Elt e => Acc (Array DIM3 e) -> Acc (Array DIM3 e)+rotate3D arr = backpermute sh rot arr+ where+ sh :: Exp DIM3+ sh =+ let Z :. z :. y :. x = unlift (shape arr) :: Z :. Exp Int :. Exp Int :. Exp Int+ in index3 y x z+ --+ rot :: Exp DIM3 -> Exp DIM3+ rot ix =+ let Z :. z :. y :. x = unlift ix :: Z :. Exp Int :. Exp Int :. Exp Int+ in index3 x z y++{--+-- 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, A.RealFloat e, A.FromIntegral Int e)+ => e+ -> sh+ -> Int+ -> Acc (Array (sh:.Int) (Complex e))+ -> Acc (Array (sh:.Int) (Complex e))+fft sign sh sz arr+ | P.any (P.not . isPow2) (shapeToList (sh:.sz))+ = error $ printf "fft: array dimensions must be powers-of-two, but are: %s" (showShape (sh:.sz))+ --+ | otherwise+ = go sz 0 1+ where+ go :: Int -> Int -> Int -> Acc (Array (sh:.Int) (Complex e))+ go len offset stride+ | len P.== 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' A.== 0 ? ( (arr ! lift (sh' :. offset')) + (arr ! lift (sh' :. offset' + stride'))+ {- A.== 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 )+++-- Append two arrays. This is a specialised version of (A.++) which does not do+-- 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 A.< n ? (xs ! lift (sz:.i), ys ! lift (sz:.i-n) ))++isPow2 :: Int -> Bool+isPow2 0 = True+isPow2 1 = False+isPow2 x = x .&. (x-1) P.== 0+--}+
+ src/Data/Array/Accelerate/Math/FFT/Adhoc.hs view
@@ -0,0 +1,603 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE RebindableSyntax #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ViewPatterns #-}+-- |+-- Module : Data.Array.Accelerate.Math.FFT.Adhoc+-- Copyright : [2017] Henning Thielemann+-- [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--+-- Implementation of ad-hoc FFT stolen from the accelerate-fourier by Henning+-- Thielemann (BSD3 licensed), and updated to work with current Accelerate. That+-- package contains other more sophisticated algorithms as well.+--++module Data.Array.Accelerate.Math.FFT.Adhoc ( fft )+ where++import Data.Array.Accelerate hiding ( transpose )+import Data.Array.Accelerate.Data.Bits+import Data.Array.Accelerate.Data.Complex+import Data.Array.Accelerate.Control.Lens.Shape++import Data.Array.Accelerate.Math.FFT.Mode+import Data.Array.Accelerate.Math.FFT.Type+++fft :: (Shape sh, Slice sh, Numeric e)+ => Mode+ -> Acc (Array (sh:.Int) (Complex e))+ -> Acc (Array (sh:.Int) (Complex e))+fft mode arr =+ let len = indexHead (shape arr)+ (pow2, smooth5) = is2or5smooth len+ in+ if len <= 1 then arr else+ if pow2 then ditSplitRadixLoop mode arr else+ if smooth5 then dit235 mode arr+ else transformChirp235 mode arr+++-- Implementations+-- ---------------++is2or5smooth :: Exp Int -> (Exp Bool, Exp Bool)+is2or5smooth len =+ let maxPowerOfTwo = len .&. negate len+ lenOdd = len `quot` maxPowerOfTwo+ in+ ( 1 == lenOdd+ , 1 == divideMaxPower 5 (divideMaxPower 3 lenOdd)+ )++divideMaxPower :: Exp Int -> Exp Int -> Exp Int+divideMaxPower fac =+ while (\n -> n `rem` fac == 0)+ (\n -> n `quot` fac)+++-- -- | Split-radix for power-of-two sizes+-- --+-- ditSplitRadix+-- :: (Shape sh, Slice sh, Numeric e)+-- => Mode+-- -> Acc (Array (sh:.Int) (Complex e))+-- -> Acc (Array (sh:.Int) (Complex e))+-- ditSplitRadix mode arr =+-- if indexHead (shape arr) <= 1+-- then arr+-- else ditSplitRadixLoop mode arr++ditSplitRadixLoop+ :: forall sh e. (Shape sh, Slice sh, Numeric e)+ => Mode+ -> Acc (Array (sh:.Int) (Complex e))+ -> Acc (Array (sh:.Int) (Complex e))+ditSplitRadixLoop mode arr =+ let+ twiddleSR (fromIntegral -> n4) k (fromIntegral -> j) =+ let w = pi * k * j / (2*n4)+ in lift (cos w :+ signOfMode mode * sin w)++ twiddle len4 k =+ generate (index1 len4) (twiddleSR len4 k . indexHead)++ step (unlift -> (us,zs)) =+ let+ k = indexHead (shape zs)+ tw1 = twiddle k 1+ tw3 = twiddle k 3+ --+ im = lift (0 :+ signOfMode mode)+ twidZeven = zipWithExtrude1 (*) tw1 (sieveV 2 0 zs)+ twidZodd = zipWithExtrude1 (*) tw3 (sieveV 2 1 zs)+ zsum = zipWith (+) twidZeven twidZodd+ zdiff = map (im *) (zipWith (-) twidZeven twidZodd)+ zcomplete = zsum ++ zdiff+ _ :. n :. _ = unlift (shape zcomplete) :: Exp sh :. Exp Int :. Exp Int+ in+ lift ( zipWith (+) us zcomplete ++ zipWith (-) us zcomplete+ , dropV n us+ )++ rebase s = lift (transform2 (-1) (afst s), asnd s)++ reorder (unlift -> (xs,ys)) =+ let evens = sieve 2 0 xs+ odds = sieve 2 1 xs+ in+ lift (evens ++^ ys, twist 2 odds)++ initial =+ let sh :. n = unlift (shape arr) :: Exp sh :. Exp Int+ in lift ( reshape (lift (sh :. 1 :. n)) arr+ , fill (lift (sh :. 0 :. n `quot` 2)) 0+ )+ in+ headV+ $ afst+ $ awhile (\s -> unit (indexHead (indexTail (shape (asnd s))) > 0)) step+ $ rebase+ $ awhile (\s -> unit (indexHead (shape (asnd s)) > 1)) reorder+ $ initial+++-- | Decimation in time for sizes that are composites of the factors 2,3 and 5.+-- These sizes are known as 5-smooth numbers or the Hamming sequence.+--+-- <http://oeis.org/A051037>+--+dit235+ :: forall sh e. (Shape sh, Slice sh, Numeric e)+ => Mode+ -> Acc (Array (sh:.Int) (Complex e))+ -> Acc (Array (sh:.Int) (Complex e))+dit235 mode arr =+ let+ merge :: forall sh' a. (Shape sh', Slice sh', Elt a)+ => Acc (Array (sh':.Int:.Int) a)+ -> Acc (Array (sh':.Int) a)+ merge xs =+ let sh :. m :. n = unlift (shape xs) :: Exp sh' :. Exp Int :. Exp Int+ in backpermute+ (lift (sh :. m*n))+ (\(unlift -> ix :. k :: Exp sh' :. Exp Int) ->+ let (q,r) = k `quotRem` m+ in lift (ix :. r :. q))+ xs++ step fac xs =+ let sh :. count :. len = unlift (shape xs) :: Exp sh :. Exp Int :. Exp Int+ twiddled = transpose+ $ zipWithExtrude2 (*) (twiddleFactors fac len)+ $ reshape (lift (sh :. count `quot` fac :. fac :. len)) xs+ in+ merge $ if fac == 5 then transform5 cache5 twiddled else+ if fac == 4 then transform4 cache4 twiddled else+ if fac == 3 then transform3 cache3 twiddled+ else transform2 cache2 twiddled++ initial :: Acc (Array (sh:.Int:.Int) (Complex e), Vector Int)+ initial =+ let sh :. n = unlift (shape arr) :: Exp sh :. Exp Int+ in lift ( reshape (lift (sh :. 1 :. n)) arr+ , fill (index1 0) 0+ )++ twiddleFactors :: Exp Int -> Exp Int -> Acc (Matrix (Complex e))+ twiddleFactors m n =+ generate (index2 m n)+ (\(unlift -> Z :. j :. i) -> twiddle (m*n) j i)++ cisrat :: Exp Int -> Exp Int -> Exp (Complex e)+ cisrat d n =+ let w = 2*pi * fromIntegral n / fromIntegral d+ in lift (cos w :+ signOfMode mode * sin w)++ twiddle :: Exp Int -> Exp Int -> Exp Int -> Exp (Complex e)+ twiddle n k j = cisrat n ((k*j) `rem` n)++ cache2 :: Exp (Complex e)+ cache2 = -1++ cache3 :: Exp (Complex e, Complex e)+ cache3 =+ let sqrt3d2 = sqrt 3 / 2+ mhalf = -1/2+ s = signOfMode mode+ u = s * sqrt3d2+ in+ lift (mhalf :+ u, mhalf :+ (-u))++ cache4 :: Exp (Complex e, Complex e, Complex e)+ cache4 =+ let s = signOfMode mode+ in lift (0 :+ s, (-1) :+ (-0), 0 :+ (-s))++ cache5 :: Exp (Complex e, Complex e, Complex e, Complex e)+ cache5 =+ let z = cisrat 5+ in lift (z 1, z 2, z 3, z 4)+ in+ headV+ $ afst+ $ awhile+ (\s -> unit (length (asnd s) > 0))+ (\s -> let (xs,fs) = unlift s+ f = fs !! 0+ in+ lift (step f xs, tail fs))+ $ awhile+ (\s -> unit (indexHead (shape (afst s)) > 1))+ (\s -> let (xs,fs) = unlift s+ len = indexHead (shape xs)+ divides k n = n `rem` k == 0+ f = if divides 3 len then 3 else+ if divides 4 len then 4 else+ if divides 5 len then 5+ else 2+ in+ lift (twist f xs, unit f `cons` fs))+ $ initial+++-- | Transformation of arbitrary length base on Bluestein on a 5-smooth size.+--+transformChirp235+ :: (Shape sh, Slice sh, Numeric e)+ => Mode+ -> Acc (Array (sh:.Int) (Complex e))+ -> Acc (Array (sh:.Int) (Complex e))+transformChirp235 mode arr =+ let n = indexHead (shape arr)+ f = ceiling5Smooth (2*n)+ in+ transformChirp mode f (dit235 Forward) (dit235 Inverse) arr+++transformChirp+ :: (Shape sh, Slice sh, Numeric e)+ => Mode+ -> Exp Int+ -> (forall sh'. (Shape sh', Slice sh') => Acc (Array (sh':.Int) (Complex e)) -> Acc (Array (sh':.Int) (Complex e)))+ -> (forall sh'. (Shape sh', Slice sh') => Acc (Array (sh':.Int) (Complex e)) -> Acc (Array (sh':.Int) (Complex e)))+ -> Acc (Array (sh:.Int) (Complex e))+ -> Acc (Array (sh:.Int) (Complex e))+transformChirp mode p analysis synthesis arr =+ let sz :. n = unlift (shape arr)+ --+ chirp =+ generate (index1 p) $ \ix ->+ let k = unindex1 ix+ sk = fromIntegral (if p > 2*k then k else k-p)+ w = pi * sk * sk / fromIntegral n+ in+ lift $ cos w :+ signOfMode mode * sin w+ --+ spectrum = analysis+ $ map conjugate chirp+ `consV`+ reshape (lift (Z :. shapeSize sz :. p))+ (pad p 0 (zipWithExtrude1 (*) chirp arr))+ scaleDown xs =+ let scale x (unlift -> r :+ i) = lift (x*r :+ x*i)+ len = indexHead (shape xs)+ in map (scale (recip (fromIntegral len))) xs+ in+ if n <= 1+ then arr+ else take n+ $ scaleDown+ $ zipWithExtrude1 (*) chirp+ $ synthesis+ $ zipWithExtrude1 (*) (headV spectrum)+ $ reshape (lift (sz:.p)) (tailV spectrum)+++ceiling5Smooth :: Exp Int -> Exp Int+ceiling5Smooth n =+ let (i2,i3,i5) = unlift (snd (ceiling5Smooth' (fromIntegral n :: Exp Double)))+ in pow i2 2 * pow i3 3 * pow i5 5++ceiling5Smooth'+ :: (RealFloat a, Ord a, FromIntegral Int a)+ => Exp a+ -> Exp (a, (Int,Int,Int))+ceiling5Smooth' n =+ let d3 = ceiling (logBase 3 n)+ d5 = ceiling (logBase 5 n)+ --+ argmin x y = if fst x < fst y then x else y+ in+ the $ fold1All argmin+ $ generate (index2 d5 d3) -- this is probably quite small!+ (\(unlift -> Z :. i5 :. i3) ->+ let+ p53 = 5 ** fromIntegral i5 * 3 ** fromIntegral i3+ i2 = 0 `max` ceiling (logBase 2 (n/p53))+ in+ lift ( p53 * 2 ** fromIntegral i2+ , (i2,i3,i5)+ ))++-- Utilities+-- ---------++pow :: Exp Int -> Exp Int -> Exp Int+pow x k+ = snd+ $ while (\ip -> fst ip < k)+ (\ip -> lift (fst ip + 1, snd ip * x))+ (lift (0,1))++pad :: (Shape sh, Slice sh, Elt e)+ => Exp Int+ -> Exp e+ -> Acc (Array (sh:.Int) e)+ -> Acc (Array (sh:.Int) e)+pad n x xs =+ let sz = indexTail (shape xs)+ sh = lift (sz :. n)+ in+ xs ++ fill sh x++cons :: forall sh e. (Shape sh, Slice sh, Elt e)+ => Acc (Array sh e)+ -> Acc (Array (sh:.Int) e)+ -> Acc (Array (sh:.Int) e)+cons x xs =+ let x' = reshape (lift (shape x :. 1)) x+ in x' ++ xs++consV :: forall sh e. (Shape sh, Slice sh, Elt e)+ => Acc (Array (sh:.Int) e)+ -> Acc (Array (sh:.Int:.Int) e)+ -> Acc (Array (sh:.Int:.Int) e)+consV x xs =+ let sh :. n = unlift (shape x) :: Exp sh :. Exp Int+ in reshape (lift (sh :. 1 :. n)) x ++^ xs++headV :: (Shape sh, Slice sh, Elt e)+ => Acc (Array (sh:.Int:.Int) e)+ -> Acc (Array (sh:.Int) e)+headV xs = slice xs (lift (Any :. (0 :: Exp Int) :. All))++tailV :: forall sh e. (Shape sh, Slice sh, Elt e)+ => Acc (Array (sh:.Int:.Int) e)+ -> Acc (Array (sh:.Int:.Int) e)+tailV = tailOn _2++dropV :: forall sh e. (Shape sh, Slice sh, Elt e)+ => Exp Int+ -> Acc (Array (sh:.Int:.Int) e)+ -> Acc (Array (sh:.Int:.Int) e)+dropV = dropOn _2++sieve+ :: forall sh e. (Shape sh, Slice sh, Elt e)+ => Exp Int+ -> Exp Int+ -> Acc (Array (sh:.Int) e)+ -> Acc (Array (sh:.Int) e)+sieve fac start xs =+ let sh :. n = unlift (shape xs) :: Exp sh :. Exp Int+ in+ backpermute+ (lift (sh :. n `quot` fac))+ (\(unlift -> ix :. j :: Exp sh :. Exp Int) -> lift (ix :. fac*j + start))+ xs++sieveV+ :: forall sh e. (Shape sh, Slice sh, Elt e)+ => Exp Int+ -> Exp Int+ -> Acc (Array (sh:.Int:.Int) e)+ -> Acc (Array (sh:.Int:.Int) e)+sieveV fac start xs =+ let sh :. m :. n = unlift (shape xs) :: Exp sh :. Exp Int :. Exp Int+ in+ backpermute+ (lift (sh :. m `quot` fac :. n))+ (\(unlift -> ix :. j :. i :: Exp sh :. Exp Int :. Exp Int) -> lift (ix :. fac*j+start :. i))+ xs++twist :: forall sh e. (Shape sh, Slice sh, Elt e)+ => Exp Int+ -> Acc (Array (sh:.Int:.Int) e)+ -> Acc (Array (sh:.Int:.Int) e)+twist fac xs =+ let sh :. m :. n = unlift (shape xs) :: Exp sh :. Exp Int :. Exp Int+ in+ backpermute+ (lift (sh :. fac*m :. n `quot` fac))+ (\(unlift -> ix :. j :. i :: Exp sh :. Exp Int :. Exp Int) -> lift (ix :. j `quot` fac :. fac*i + j `rem` fac))+ xs+++infixr 5 ++^+(++^) :: forall sh e. (Slice sh, Shape sh, Elt e)+ => Acc (Array (sh:.Int:.Int) e)+ -> Acc (Array (sh:.Int:.Int) e)+ -> Acc (Array (sh:.Int:.Int) e)+(++^) = concatOn _2++zipWithExtrude1+ :: (Shape sh, Slice sh, Elt a, Elt b, Elt c)+ => (Exp a -> Exp b -> Exp c)+ -> Acc (Array DIM1 a)+ -> Acc (Array (sh:.Int) b)+ -> Acc (Array (sh:.Int) c)+zipWithExtrude1 f xs ys =+ zipWith f (replicate (lift (indexTail (shape ys) :. All)) xs) ys++zipWithExtrude2+ :: (Shape sh, Slice sh, Elt a, Elt b, Elt c)+ => (Exp a -> Exp b -> Exp c)+ -> Acc (Array DIM2 a)+ -> Acc (Array (sh:.Int:.Int) b)+ -> Acc (Array (sh:.Int:.Int) c)+zipWithExtrude2 f xs ys =+ zipWith f (replicate (lift (indexTail (indexTail (shape ys)) :. All :. All)) xs) ys++transpose+ :: forall sh e. (Shape sh, Slice sh, Elt e)+ => Acc (Array (sh:.Int:.Int) e)+ -> Acc (Array (sh:.Int:.Int) e)+transpose = transposeOn _1 _2++transform2+ :: (Shape sh, Slice sh, Num e)+ => Exp e+ -> Acc (Array (sh:.Int) e)+ -> Acc (Array (sh:.Int) e)+transform2 v xs =+ generate+ (lift (indexTail (shape xs) :. 2))+ (\(unlift -> ix :. k :: Exp sh :. Exp Int) ->+ let x0 = xs ! lift (ix :. 0)+ x1 = xs ! lift (ix :. 1)+ in+ if k == 0 then x0+x1+ else x0+v*x1)++transform3+ :: forall sh e. (Shape sh, Slice sh, Num e)+ => Exp (e,e)+ -> Acc (Array (sh:.Int) e)+ -> Acc (Array (sh:.Int) e)+transform3 (unlift -> (z1,z2)) xs =+ generate+ (lift (indexTail (shape xs) :. 3))+ (\(unlift -> ix :. k :: Exp sh :. Exp Int) ->+ let+ x0 = xs ! lift (ix :. 0)+ x1 = xs ! lift (ix :. 1)+ x2 = xs ! lift (ix :. 2)+ --+ ((s,_), (zx1,zx2)) = sumAndConvolve2 (x1,x2) (z1,z2)+ in+ if k == 0 then x0 + s else+ if k == 1 then x0 + zx1+ {- k == 2 -} else x0 + zx2)++transform4+ :: forall sh e. (Shape sh, Slice sh, Num e)+ => Exp (e,e,e)+ -> Acc (Array (sh:.Int) e)+ -> Acc (Array (sh:.Int) e)+transform4 (unlift -> (z1,z2,z3)) xs =+ generate+ (lift (indexTail (shape xs) :. 4))+ (\(unlift -> ix :. k :: Exp sh :. Exp Int) ->+ let+ x0 = xs ! lift (ix :. 0)+ x1 = xs ! lift (ix :. 1)+ x2 = xs ! lift (ix :. 2)+ x3 = xs ! lift (ix :. 3)+ --+ x02a = x0+x2+ x02b = x0+z2*x2+ x13a = x1+x3+ x13b = x1+z2*x3+ in+ if k == 0 then x02a + x13a else+ if k == 1 then x02b + z1 * x13b else+ if k == 2 then x02a + z2 * x13a+ {- k == 3 -} else x02b + z3 * x13b)++-- Use Rader's trick for mapping the transform to a convolution and apply+-- Karatsuba's trick at two levels (i.e. total three times) to that convolution.+--+-- 0 0 0 0 0+-- 0 1 2 3 4+-- 0 2 4 1 3+-- 0 3 1 4 2+-- 0 4 3 2 1+--+-- Permutation.T: 0 1 2 4 3+--+-- 0 0 0 0 0+-- 0 1 2 4 3+-- 0 2 4 3 1+-- 0 4 3 1 2+-- 0 3 1 2 4+--+transform5+ :: forall sh e. (Shape sh, Slice sh, Num e)+ => Exp (e,e,e,e)+ -> Acc (Array (sh:.Int) e)+ -> Acc (Array (sh:.Int) e)+transform5 (unlift -> (z1,z2,z3,z4)) xs =+ generate+ (lift (indexTail (shape xs) :. 5))+ (\(unlift -> ix :. k :: Exp sh :. Exp Int) ->+ let+ x0 = xs ! lift (ix :. 0)+ x1 = xs ! lift (ix :. 1)+ x2 = xs ! lift (ix :. 2)+ x3 = xs ! lift (ix :. 3)+ x4 = xs ! lift (ix :. 4)+ --+ ((s,_), (d1,d2,d4,d3)) = sumAndConvolve4 (x1,x3,x4,x2) (z1,z2,z4,z3)+ in+ if k == 0 then x0 + s else+ if k == 1 then x0 + d1 else+ if k == 2 then x0 + d2 else+ if k == 3 then x0 + d3+ {- k == 4 -} else x0 + d4)+++-- Some small size convolutions using the Karatsuba trick.+--+-- This does not use Toom-3 multiplication, because this requires division by+-- 2 and 6, and thus 'Fractional' constraints.+--+sumAndConvolve2+ :: Num e+ => (Exp e, Exp e)+ -> (Exp e, Exp e)+ -> ((Exp e, Exp e), (Exp e, Exp e))+sumAndConvolve2 (a0,a1) (b0,b1) =+ let sa01 = a0+a1+ sb01 = b0+b1+ ab0ab1 = a0*b0+a1*b1+ in+ ((sa01, sb01), (ab0ab1, sa01*sb01-ab0ab1))++-- sumAndConvolve3+-- :: Num e+-- => (Exp e, Exp e, Exp e)+-- -> (Exp e, Exp e, Exp e)+-- -> ((Exp e, Exp e), (Exp e, Exp e, Exp e))+-- sumAndConvolve3 (a0,a1,a2) (b0,b1,b2) =+-- let ab0 = a0*b0+-- dab12 = a1*b1 - a2*b2+-- sa01 = a0+a1; sb01 = b0+b1; tab01 = sa01*sb01 - ab0+-- sa02 = a0+a2; sb02 = b0+b2; tab02 = sa02*sb02 - ab0+-- sa012 = sa01+a2+-- sb012 = sb01+b2+-- --+-- d0 = sa012*sb012 - tab01 - tab02+-- d1 = tab01 - dab12+-- d2 = tab02 + dab12+-- in+-- ((sa012, sb012), (d0, d1, d2))++sumAndConvolve4+ :: Num e+ => (Exp e, Exp e, Exp e, Exp e)+ -> (Exp e, Exp e, Exp e, Exp e)+ -> ((Exp e, Exp e), (Exp e, Exp e, Exp e, Exp e))+sumAndConvolve4 (a0,a1,a2,a3) (b0,b1,b2,b3) =+ let ab0 = a0*b0+ ab1 = a1*b1+ sa01 = a0+a1; sb01 = b0+b1+ ab01 = sa01*sb01 - (ab0+ab1)+ ab2 = a2*b2+ ab3 = a3*b3+ sa23 = a2+a3; sb23 = b2+b3+ ab23 = sa23*sb23 - (ab2+ab3)+ c0 = ab0 + ab2 - (ab1 + ab3)+ c1 = ab01 + ab23+ ab02 = (a0+a2)*(b0+b2)+ ab13 = (a1+a3)*(b1+b3)+ sa0123 = sa01+sa23+ sb0123 = sb01+sb23+ ab0123 = sa0123*sb0123 - (ab02+ab13)+ --+ d0 = ab13 + c0+ d1 = c1+ d2 = ab02 - c0+ d3 = ab0123 - c1+ in+ ((sa0123, sb0123), (d0, d1, d2, d3))+
+ src/Data/Array/Accelerate/Math/FFT/LLVM/Native.hs view
@@ -0,0 +1,110 @@+{-# LANGUAGE GADTs #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+-- |+-- Module : Data.Array.Accelerate.Math.FFT.LLVM.Native+-- Copyright : [2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module Data.Array.Accelerate.Math.FFT.LLVM.Native (++ fft,+ fft1D,+ fft2D,+ fft3D,++) where++import Data.Array.Accelerate.Math.FFT.Mode+import Data.Array.Accelerate.Math.FFT.Type+import Data.Array.Accelerate.Math.FFT.LLVM.Native.Ix+import Data.Array.Accelerate.Math.FFT.LLVM.Native.Base++import Data.Array.Accelerate+import Data.Array.Accelerate.Analysis.Match+import Data.Array.Accelerate.Array.Sugar+import Data.Array.Accelerate.Data.Complex+import Data.Array.Accelerate.Error++import Data.Array.Accelerate.LLVM.Native.Foreign++import Data.Array.CArray ( CArray )+import Math.FFT.Base ( FFTWReal )+import Prelude as P+import qualified Math.FFT as FFT+++fft :: forall sh e. (Shape sh, Numeric e)+ => Mode+ -> ForeignAcc (Array sh (Complex e) -> Array sh (Complex e))+fft mode+ = ForeignAcc (nameOf mode (undefined::sh))+ $ case numericR::NumericR e of+ NumericRfloat32 -> go+ NumericRfloat64 -> go+ where+ go :: FFTWReal e => Array sh (Complex e) -> LLVM Native (Array sh (Complex e))+ go | Just Refl <- matchShapeType (undefined::sh) (undefined::DIM1) = liftCtoA (FFT.dftGU (signOf mode) flags [0] `ix` (undefined :: (Int)))+ | Just Refl <- matchShapeType (undefined::sh) (undefined::DIM2) = liftCtoA (FFT.dftGU (signOf mode) flags [1] `ix` (undefined :: (Int,Int)))+ | Just Refl <- matchShapeType (undefined::sh) (undefined::DIM3) = liftCtoA (FFT.dftGU (signOf mode) flags [2] `ix` (undefined :: (Int,Int,Int)))+ | Just Refl <- matchShapeType (undefined::sh) (undefined::DIM4) = liftCtoA (FFT.dftGU (signOf mode) flags [3] `ix` (undefined :: (Int,Int,Int,Int)))+ | Just Refl <- matchShapeType (undefined::sh) (undefined::DIM5) = liftCtoA (FFT.dftGU (signOf mode) flags [4] `ix` (undefined :: (Int,Int,Int,Int,Int)))+ | otherwise = $internalError "fft" "only for 1D..5D inner-dimension transforms"+ --+ ix :: (a i r -> a i r) -> i -> (a i r -> a i r)+ ix f _ = f+++fft1D :: forall e. Numeric e+ => Mode+ -> ForeignAcc (Array DIM1 (Complex e) -> Array DIM1 (Complex e))+fft1D mode+ = ForeignAcc (nameOf mode (undefined::DIM1))+ $ case numericR::NumericR e of+ NumericRfloat32 -> liftCtoA go+ NumericRfloat64 -> liftCtoA go+ where+ go :: FFTWReal r => CArray Int (Complex r) -> CArray Int (Complex r)+ go = FFT.dftGU (signOf mode) flags [0]++fft2D :: forall e. Numeric e+ => Mode+ -> ForeignAcc (Array DIM2 (Complex e) -> Array DIM2 (Complex e))+fft2D mode+ = ForeignAcc (nameOf mode (undefined::DIM2))+ $ case numericR::NumericR e of+ NumericRfloat32 -> liftCtoA go+ NumericRfloat64 -> liftCtoA go+ where+ go :: FFTWReal r => CArray (Int,Int) (Complex r) -> CArray (Int,Int) (Complex r)+ go = FFT.dftGU (signOf mode) flags [0,1]++fft3D :: forall e. Numeric e+ => Mode+ -> ForeignAcc (Array DIM3 (Complex e) -> Array DIM3 (Complex e))+fft3D mode+ = ForeignAcc (nameOf mode (undefined::DIM3))+ $ case numericR::NumericR e of+ NumericRfloat32 -> liftCtoA go+ NumericRfloat64 -> liftCtoA go+ where+ go :: FFTWReal r => CArray (Int,Int,Int) (Complex r) -> CArray (Int,Int,Int) (Complex r)+ go = FFT.dftGU (signOf mode) flags [0,1,2]++{-# INLINE liftCtoA #-}+liftCtoA+ :: forall ix sh e. (IxShapeRepr (EltRepr ix) ~ EltRepr sh, Shape sh, Elt ix, Numeric e)+ => (CArray ix (Complex e) -> CArray ix (Complex e))+ -> Array sh (Complex e)+ -> LLVM Native (Array sh (Complex e))+liftCtoA f a =+ liftIO $ withCArray a (fromCArray . f)+
+ src/Data/Array/Accelerate/Math/FFT/LLVM/Native/Base.hs view
@@ -0,0 +1,119 @@+{-# LANGUAGE GADTs #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}+-- |+-- Module : Data.Array.Accelerate.Math.FFT.LLVM.Native.Base+-- Copyright : [2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module Data.Array.Accelerate.Math.FFT.LLVM.Native.Base+ where++import Data.Array.Accelerate.Analysis.Match+import Data.Array.Accelerate.Array.Data+import Data.Array.Accelerate.Array.Sugar+import Data.Array.Accelerate.Array.Unique+import Data.Array.Accelerate.Data.Complex+import Data.Array.Accelerate.Lifetime++import Data.Array.Accelerate.Math.FFT.Mode+import Data.Array.Accelerate.Math.FFT.Type++import Data.Array.Accelerate.Math.FFT.LLVM.Native.Ix++import Data.Array.CArray.Base ( CArray(..) )+import Math.FFT.Base ( Sign(..), Flag, measure, preserveInput )++import Data.Bits+import Data.Typeable+import Foreign.ForeignPtr+import Text.Printf+import Prelude as P+++signOf :: Mode -> Sign+signOf Forward = DFTForward+signOf _ = DFTBackward++flags :: Flag+flags = measure .|. preserveInput++nameOf :: forall sh. Shape sh => Mode -> sh -> String+nameOf Forward _ = printf "FFTW.dft%dD" (rank (undefined::sh))+nameOf _ _ = printf "FFTW.idft%dD" (rank (undefined::sh))+++-- /O(1)/ Convert a CArray to an Accelerate array+--+{-# INLINE fromCArray #-}+fromCArray+ :: forall ix sh e. (IxShapeRepr (EltRepr ix) ~ EltRepr sh, Shape sh, Elt ix, Numeric e)+ => CArray ix (Complex e)+ -> IO (Array sh (Complex e))+fromCArray (CArray lo hi _ fp) = do+ --+ sh <- return $ rangeToShape (toIxShapeRepr lo, toIxShapeRepr hi) :: IO sh+ ua <- newUniqueArray (castForeignPtr fp :: ForeignPtr e)+ --+ case numericR::NumericR e of+ NumericRfloat32 -> return $ Array (fromElt sh) (AD_V2 (AD_Float ua))+ NumericRfloat64 -> return $ Array (fromElt sh) (AD_V2 (AD_Double ua))++-- /O(1)/ Use an Accelerate array as a CArray+--+{-# INLINE withCArray #-}+withCArray+ :: forall ix sh e a. (IxShapeRepr (EltRepr ix) ~ EltRepr sh, Shape sh, Elt ix, Numeric e)+ => Array sh (Complex e)+ -> (CArray ix (Complex e) -> IO a)+ -> IO a+withCArray arr f =+ let+ sh = shape arr+ (lo, hi) = shapeToRange sh+ wrap fp = CArray (fromIxShapeRepr lo) (fromIxShapeRepr hi) (size sh) (castForeignPtr fp)+ in+ withArray arr (f . wrap)+++-- Use underlying array pointers+--+{-# INLINE withArray #-}+withArray+ :: forall sh e a. Numeric e+ => Array sh (Complex e)+ -> (ForeignPtr e -> IO a)+ -> IO a+withArray (Array _ adata) = withArrayData (numericR::NumericR e) adata++{-# INLINE withArrayData #-}+withArrayData+ :: NumericR e+ -> ArrayData (EltRepr (Complex e))+ -> (ForeignPtr e -> IO a)+ -> IO a+withArrayData NumericRfloat32 (AD_V2 (AD_Float ua)) = withLifetime (uniqueArrayData ua)+withArrayData NumericRfloat64 (AD_V2 (AD_Double ua)) = withLifetime (uniqueArrayData ua)+++-- Match shape surface types+--+{-# INLINE matchShapeType #-}+matchShapeType+ :: forall sh sh'. (Shape sh, Shape sh')+ => sh+ -> sh'+ -> Maybe (sh :~: sh')+matchShapeType _ _+ | Just Refl <- matchTupleType (eltType (undefined::sh)) (eltType (undefined::sh'))+ = gcast Refl++matchShapeType _ _+ = Nothing+
+ src/Data/Array/Accelerate/Math/FFT/LLVM/Native/Ix.hs view
@@ -0,0 +1,57 @@+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TypeFamilies #-}+-- |+-- Module : Data.Array.Accelerate.Math.FFT.LLVM.Native.Ix+-- Copyright : [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module Data.Array.Accelerate.Math.FFT.LLVM.Native.Ix+ where++import Data.Array.Accelerate.Array.Sugar+import Data.Array.Accelerate.Error+import Data.Array.Accelerate.Type+++-- Converting between Accelerate multidimensional shapes/indices and those used+-- by the CArray package (Data.Ix)+--+type family IxShapeRepr e where+ IxShapeRepr () = ()+ IxShapeRepr Int = ((),Int)+ IxShapeRepr (t,h) = (IxShapeRepr t, h)++{-# INLINE fromIxShapeRepr #-}+fromIxShapeRepr+ :: forall ix sh. (IxShapeRepr (EltRepr ix) ~ EltRepr sh, Shape sh, Elt ix)+ => sh+ -> ix+fromIxShapeRepr = liftToElt (go (eltType (undefined::ix)))+ where+ go :: forall ix'. TupleType ix' -> IxShapeRepr ix' -> ix'+ go TypeRunit () = ()+ go (TypeRpair tt _) (t, h) = (go tt t, h)+ go (TypeRscalar (SingleScalarType (NumSingleType (IntegralNumType TypeInt{})))) ((),h) = h+ go _ _+ = $internalError "fromIxShapeRepr" "expected Int dimensions"++{-# INLINE toIxShapeRepr #-}+toIxShapeRepr+ :: forall ix sh. (IxShapeRepr (EltRepr ix) ~ EltRepr sh, Shape sh, Elt ix)+ => ix+ -> sh+toIxShapeRepr = liftToElt (go (eltType (undefined::ix)))+ where+ go :: forall ix'. TupleType ix' -> ix' -> IxShapeRepr ix'+ go TypeRunit () = ()+ go (TypeRscalar (SingleScalarType (NumSingleType (IntegralNumType TypeInt{})))) h = ((), h)+ go (TypeRpair tt _) (t, h) = (go tt t, h)+ go _ _+ = $internalError "toIxShapeRepr" "not a valid Data.Ix index"+
+ src/Data/Array/Accelerate/Math/FFT/LLVM/PTX.hs view
@@ -0,0 +1,170 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TupleSections #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ViewPatterns #-}+-- |+-- Module : Data.Array.Accelerate.Math.FFT.LLVM.PTX+-- Copyright : [2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module Data.Array.Accelerate.Math.FFT.LLVM.PTX (++ fft,+ fft1D,+ fft2D,+ fft3D,++) where++import Data.Array.Accelerate.Math.FFT.Mode+import Data.Array.Accelerate.Math.FFT.Type+import Data.Array.Accelerate.Math.FFT.LLVM.PTX.Base+import Data.Array.Accelerate.Math.FFT.LLVM.PTX.Plans++import Data.Array.Accelerate.Array.Sugar+import Data.Array.Accelerate.Data.Complex+import Data.Array.Accelerate.Error+import Data.Array.Accelerate.Lifetime++import Data.Array.Accelerate.LLVM.PTX.Foreign++import Foreign.CUDA.Ptr ( DevicePtr, castDevPtr )+import qualified Foreign.CUDA.FFT as FFT++import Data.Hashable+import Data.Proxy+import Data.Typeable+import System.IO.Unsafe+++fft :: forall sh e. (Shape sh, Numeric e)+ => Mode+ -> ForeignAcc (Array (sh:.Int) (Complex e) -> Array (sh:.Int) (Complex e))+fft mode+ | Just Refl <- matchShapeType (undefined::sh) (undefined::DIM0) = fft1D mode+ | Just Refl <- matchShapeType (undefined::sh) (undefined::DIM1) = ForeignAcc "cuda.fft2.many" $ fft' fft2DMany_plans mode+ | Just Refl <- matchShapeType (undefined::sh) (undefined::DIM2) = ForeignAcc "cuda.fft3.many" $ fft' fft3DMany_plans mode+ | otherwise = $internalError "fft" "only for 1D..3D inner-dimension transforms"++fft1D :: Numeric e+ => Mode+ -> ForeignAcc (Vector (Complex e) -> Vector (Complex e))+fft1D mode = ForeignAcc "cuda.fft1d" $ fft' fft1D_plans mode++fft2D :: Numeric e+ => Mode+ -> ForeignAcc (Array DIM2 (Complex e) -> Array DIM2 (Complex e))+fft2D mode = ForeignAcc "cuda.fft2d" $ fft' fft2D_plans mode++fft3D :: Numeric e+ => Mode+ -> ForeignAcc (Array DIM3 (Complex e) -> Array DIM3 (Complex e))+fft3D mode = ForeignAcc "cuda.fft3d" $ fft' fft3D_plans mode+++-- Internals+-- ---------++{-# INLINEABLE fft' #-}+fft' :: forall sh e. (Shape sh, Numeric e)+ => Plans (sh, FFT.Type)+ -> Mode+ -> Stream+ -> Array sh (Complex e)+ -> LLVM PTX (Array sh (Complex e))+fft' plans mode stream =+ let+ go :: Numeric e => Array sh (Complex e) -> LLVM PTX (Array sh (Complex e))+ go ain = do+ let+ sh = shape ain+ t = fftType (Proxy::Proxy e)+ --+ aout <- allocateRemote sh+ withArray ain stream $ \d_in -> do+ withArray aout stream $ \d_out -> do+ withPlan plans (sh,t) $ \h -> do+ liftIO $ cuFFT (Proxy::Proxy e) h mode stream (castDevPtr d_in) (castDevPtr d_out)+ return aout+ in+ case numericR::NumericR e of+ NumericRfloat32 -> go+ NumericRfloat64 -> go+++-- Execute the FFT+--+{-# INLINE cuFFT #-}+cuFFT :: forall e. Numeric e+ => Proxy e+ -> FFT.Handle+ -> Mode+ -> Stream+ -> DevicePtr (Complex e)+ -> DevicePtr (Complex e)+ -> IO ()+cuFFT _ p mode stream d_in d_out =+ withLifetime stream $ \s -> do+ FFT.setStream p s+ case numericR::NumericR e of+ NumericRfloat32 -> FFT.execC2C p (fftMode mode) d_in d_out+ NumericRfloat64 -> FFT.execZ2Z p (fftMode mode) d_in d_out++fftType :: forall e. Numeric e => Proxy e -> FFT.Type+fftType _ =+ case numericR::NumericR e of+ NumericRfloat32 -> FFT.C2C+ NumericRfloat64 -> FFT.Z2Z++fftMode :: Mode -> FFT.Mode+fftMode Forward = FFT.Forward+fftMode _ = FFT.Inverse+++-- Plan caches+-- -----------++{-# NOINLINE fft1D_plans #-}+fft1D_plans :: Plans (DIM1, FFT.Type)+fft1D_plans+ = unsafePerformIO+ $ createPlan (\(Z:.n, t) -> FFT.plan1D n t 1)+ (\(Z:.n, t) -> fromEnum t `hashWithSalt` n)++{-# NOINLINE fft2D_plans #-}+fft2D_plans :: Plans (DIM2, FFT.Type)+fft2D_plans+ = unsafePerformIO+ $ createPlan (\(Z:.h:.w, t) -> FFT.plan2D h w t)+ (\(Z:.h:.w, t) -> fromEnum t `hashWithSalt` h `hashWithSalt` w)++{-# NOINLINE fft3D_plans #-}+fft3D_plans :: Plans (DIM3, FFT.Type)+fft3D_plans+ = unsafePerformIO+ $ createPlan (\(Z:.d:.h:.w, t) -> FFT.plan3D d h w t)+ (\(Z:.d:.h:.w, t) -> fromEnum t `hashWithSalt` d `hashWithSalt` h `hashWithSalt` w)++{-# NOINLINE fft2DMany_plans #-}+fft2DMany_plans :: Plans (DIM2, FFT.Type)+fft2DMany_plans+ = unsafePerformIO+ $ createPlan (\(Z:.h:.w, t) -> FFT.planMany [h,w] Nothing Nothing t 1)+ (\(Z:.h:.w, t) -> fromEnum t `hashWithSalt` h `hashWithSalt` w)++{-# NOINLINE fft3DMany_plans #-}+fft3DMany_plans :: Plans (DIM3, FFT.Type)+fft3DMany_plans+ = unsafePerformIO+ $ createPlan (\(Z:.d:.h:.w, t) -> FFT.planMany [d,h,w] Nothing Nothing t 1)+ (\(Z:.d:.h:.w, t) -> fromEnum t `hashWithSalt` d `hashWithSalt` h `hashWithSalt` w)+
+ src/Data/Array/Accelerate/Math/FFT/LLVM/PTX/Base.hs view
@@ -0,0 +1,82 @@+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+-- |+-- Module : Data.Array.Accelerate.Math.FFT.LLVM.PTX.Base+-- Copyright : [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module Data.Array.Accelerate.Math.FFT.LLVM.PTX.Base+ where++import Data.Array.Accelerate.Math.FFT.Type++import Data.Array.Accelerate.Analysis.Match+import Data.Array.Accelerate.Array.Data+import Data.Array.Accelerate.Array.Sugar+import Data.Array.Accelerate.Data.Complex+import Data.Array.Accelerate.Lifetime++import Data.Array.Accelerate.LLVM.PTX.Foreign++import Foreign.CUDA.Ptr ( DevicePtr )++import Data.Typeable+++{-# INLINE withArray #-}+withArray+ :: forall sh e b. Numeric e+ => Array sh (Complex e)+ -> Stream+ -> (DevicePtr e -> LLVM PTX b)+ -> LLVM PTX b+withArray (Array _ adata) = withArrayData (numericR::NumericR e) adata++{-# INLINE withArrayData #-}+withArrayData+ :: NumericR e+ -> ArrayData (EltRepr (Complex e))+ -> Stream+ -> (DevicePtr e -> LLVM PTX b)+ -> LLVM PTX b+withArrayData NumericRfloat32 (AD_V2 ad) s k =+ withDevicePtr ad $ \p -> do+ r <- k p+ e <- checkpoint s+ return (Just e,r)+withArrayData NumericRfloat64 (AD_V2 ad) s k =+ withDevicePtr ad $ \p -> do+ r <- k p+ e <- checkpoint s+ return (Just e, r)++{-# INLINE withLifetime' #-}+withLifetime' :: Lifetime a -> (a -> LLVM PTX b) -> LLVM PTX b+withLifetime' l k = do+ r <- k (unsafeGetValue l)+ liftIO $ touchLifetime l+ return r+++-- Match shape surface types+--+{-# INLINE matchShapeType #-}+matchShapeType+ :: forall sh sh'. (Shape sh, Shape sh')+ => sh+ -> sh'+ -> Maybe (sh :~: sh')+matchShapeType _ _+ | Just Refl <- matchTupleType (eltType (undefined::sh)) (eltType (undefined::sh'))+ = gcast Refl++matchShapeType _ _+ = Nothing+
+ src/Data/Array/Accelerate/Math/FFT/LLVM/PTX/Plans.hs view
@@ -0,0 +1,91 @@+{-# LANGUAGE MagicHash #-}+{-# LANGUAGE RecordWildCards #-}+-- |+-- Module : Data.Array.Accelerate.Math.FFT.LLVM.PTX.Plans+-- Copyright : [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module Data.Array.Accelerate.Math.FFT.LLVM.PTX.Plans (++ Plans,+ createPlan,+ withPlan,++) where++import Data.Array.Accelerate.Lifetime+import Data.Array.Accelerate.LLVM.PTX+import Data.Array.Accelerate.LLVM.PTX.Foreign++import Data.Array.Accelerate.Math.FFT.LLVM.PTX.Base++import Control.Concurrent.MVar+import Control.Monad.State+import Data.HashMap.Strict+import qualified Data.HashMap.Strict as Map++import qualified Foreign.CUDA.Driver.Context as CUDA+import qualified Foreign.CUDA.FFT as FFT++import GHC.Ptr+import GHC.Base+import Prelude hiding ( lookup )+++data Plans a = Plans+ { plans :: {-# UNPACK #-} !(MVar ( HashMap (Int, Int) (Lifetime FFT.Handle)))+ , create :: a -> IO FFT.Handle+ , hash :: a -> Int+ }+++-- Create a new plan cache+--+{-# INLINE createPlan #-}+createPlan :: (a -> IO FFT.Handle) -> (a -> Int) -> IO (Plans a)+createPlan via mix =+ Plans <$> newMVar Map.empty <*> pure via <*> pure mix+++-- Execute an operation with a cuFFT handle appropriate for the current+-- execution context.+--+-- Initial creation of the context is an atomic operation, but subsequently+-- multiple threads may use the context concurrently.+--+-- TLM: check that plans can be used concurrently+--+-- <http://docs.nvidia.com/cuda/cufft/index.html#thread-safety>+--+{-# INLINE withPlan #-}+withPlan :: Plans a -> a -> (FFT.Handle -> LLVM PTX b) -> LLVM PTX b+withPlan Plans{..} a k = do+ lc <- gets (deviceContext . ptxContext)+ h <- liftIO $+ withLifetime lc $ \ctx ->+ modifyMVar plans $ \pm ->+ let key = (toKey ctx, hash a) in+ case Map.lookup key pm of+ -- handle does not exist yet; create it and add to the global+ -- state for reuse+ Nothing -> do+ h <- create a+ l <- newLifetime h+ addFinalizer lc $ modifyMVar plans (\pm' -> return (Map.delete key pm', ()))+ addFinalizer l $ FFT.destroy h+ return ( Map.insert key l pm, l )++ -- return existing handle+ Just h -> return (pm, h)+ --+ withLifetime' h k++{-# INLINE toKey #-}+toKey :: CUDA.Context -> Int+toKey (CUDA.Context (Ptr addr#)) = I# (addr2Int# addr#)+
+ src/Data/Array/Accelerate/Math/FFT/Mode.hs view
@@ -0,0 +1,28 @@+-- |+-- Module : Data.Array.Accelerate.Math.FFT.Mode+-- Copyright : [2012..2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- [2013..2017] Robert Clifton-Everest+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module Data.Array.Accelerate.Math.FFT.Mode+ where+++data Mode+ = Forward -- ^ Forward DFT+ | Reverse -- ^ Inverse DFT, un-normalised+ | Inverse -- ^ Inverse DFT, normalised+ deriving (Eq, Show)++signOfMode :: Num a => Mode -> a+signOfMode m+ = case m of+ Forward -> -1+ Reverse -> 1+ Inverse -> 1+
+ src/Data/Array/Accelerate/Math/FFT/Type.hs view
@@ -0,0 +1,37 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE RebindableSyntax #-}+-- |+-- Module : Data.Array.Accelerate.Math.FFT.Type+-- Copyright : [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module Data.Array.Accelerate.Math.FFT.Type+ where++import Data.Array.Accelerate as A+import Data.Array.Accelerate.Data.Complex as A+++-- For explicit dictionary reification, to discover the concrete type the+-- operation should be performed at.+--+data NumericR a where+ NumericRfloat32 :: NumericR Float+ NumericRfloat64 :: NumericR Double++class (RealFloat a, FromIntegral Int a, Elt (Complex a)) => Numeric a where+ numericR :: NumericR a++instance Numeric Float where+ numericR = NumericRfloat32++instance Numeric Double where+ numericR = NumericRfloat64+
+ test/Test/Base.hs view
@@ -0,0 +1,60 @@+{-# LANGUAGE RankNTypes #-}+-- |+-- Module : Test.Base+-- Copyright : [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module Test.Base+ where++import Data.Array.Accelerate ( Z(..), (:.)(..), DIM1, DIM2, DIM3, Shape, Elt, Acc, Array )+import Data.Array.Accelerate.Array.Sugar ( fromList, size )+import Data.Array.Accelerate.Trafo ( Afunction, AfunctionR )+import Data.Array.Accelerate.Data.Complex+import Data.Array.Accelerate.Math.FFT++import Hedgehog+import qualified Hedgehog.Gen as Gen+import qualified Hedgehog.Range as Range++import Prelude as P+++type RunN = forall f. Afunction f => f -> AfunctionR f++type Transform sh e = Mode -> Acc (Array sh e) -> Acc (Array sh e)+++f32 :: Gen Float+f32 = Gen.realFloat (Range.linearFracFrom 0 (-1) 1)++f64 :: Gen Double+f64 = Gen.realFloat (Range.linearFracFrom 0 (-1) 1)++complex :: Gen a -> Gen (Complex a)+complex f = (:+) <$> f <*> f++dim1 :: Gen DIM1+dim1 = (Z :.) <$> Gen.int (Range.linear 1 1024)++dim2 :: Gen DIM2+dim2 = do+ x <- Gen.int (Range.linear 1 128)+ y <- Gen.int (Range.linear 1 48)+ return (Z :. y :. x)++dim3 :: Gen DIM3+dim3 = do+ x <- Gen.int (Range.linear 1 64)+ y <- Gen.int (Range.linear 1 32)+ z <- Gen.int (Range.linear 1 16)+ return (Z :. z :. y :. x)++array :: (Shape sh, Elt e) => sh -> Gen e -> Gen (Array sh e)+array sh gen = fromList sh <$> Gen.list (Range.singleton (size sh)) gen+
+ test/Test/FFT.hs view
@@ -0,0 +1,237 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ViewPatterns #-}+-- |+-- Module : Test.FFT+-- Copyright : [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module Test.FFT ( testFFT )+ where++import Test.Base+import Test.ShowType++import Data.Array.Accelerate as A hiding ( RealFloat, Eq, reverse )+import Data.Array.Accelerate.Data.Complex+import Data.Array.Accelerate.Math.FFT+import Data.Array.Accelerate.Test.Similar++import Hedgehog+import qualified Hedgehog.Gen as Gen++import Test.Tasty+import Test.Tasty.Hedgehog++import Data.Proxy+import Prelude as P hiding ( reverse )+++testFFT :: RunN -> TestTree+testFFT runN =+ testGroup "FFT"+ [ testFFT' f32 runN+ , testFFT' f64 runN+ ]++testFFT'+ :: forall e. (Numeric e, Similar e, RealFloat e, Show (ArgType e))+ => Gen e+ -> RunN+ -> TestTree+testFFT' e runN =+ testGroup (showType (Proxy::Proxy e))+ [ testGroup "DIM1"+ [ testProperty "homogeneity" $ test_homogeneity runN fft1D dim1 e+ , testProperty "additivity" $ test_additivity runN fft1D dim1 e+ , testProperty "inverse" $ test_inverse runN fft1D dim1 e+ , testProperty "reverse" $ test_reverse runN fft1D dim1 e+ , testProperty "conjugate" $ test_conjugate runN fft1D dim1 e+ , testProperty "isometry" $ test_isometry runN fft1D dim1 e+ , testProperty "unitarity" $ test_unitarity runN fft1D dim1 e+ ]+ , testGroup "DIM2"+ [ testProperty "homogeneity" $ test_homogeneity runN fft2D dim2 e+ , testProperty "additivity" $ test_additivity runN fft2D dim2 e+ , testProperty "inverse" $ test_inverse runN fft2D dim2 e+ , testProperty "reverse" $ test_reverse runN fft dim2 e+ , testProperty "conjugate" $ test_conjugate runN fft dim2 e+ , testProperty "isometry" $ test_isometry runN fft dim2 e+ , testProperty "unitarity" $ test_unitarity runN fft dim2 e+ ]+ , testGroup "DIM3"+ [ testProperty "homogeneity" $ test_homogeneity runN fft3D dim3 e+ , testProperty "additivity" $ test_additivity runN fft3D dim3 e+ , testProperty "inverse" $ test_inverse runN fft3D dim3 e+ , testProperty "reverse" $ test_reverse runN fft dim3 e+ , testProperty "conjugate" $ test_conjugate runN fft dim3 e+ , testProperty "isometry" $ test_isometry runN fft dim3 e+ , testProperty "unitarity" $ test_unitarity runN fft dim3 e+ ]+ ]+++mode :: Gen Mode+mode = Gen.element [Forward, Reverse, Inverse]++reverse+ :: forall sh e. (Shape sh, Slice sh, Elt e)+ => Acc (Array (sh:.Int) e)+ -> Acc (Array (sh:.Int) e)+reverse arr =+ let sh = A.shape arr+ n = A.indexHead sh+ in+ A.backpermute sh (\(A.unlift -> ix:.k :: Exp sh :. Exp Int) -> A.lift (ix :. (-k) `mod` n)) arr++norm2+ :: (Numeric e, Shape sh)+ => Acc (Array (sh:.Int) (Complex e))+ -> Acc (Array sh e)+norm2 = A.map sqrt . A.sum . A.map (\c -> real c * real c + imag c * imag c)++dotc :: (Numeric e, Shape sh)+ => Acc (Array (sh:.Int) (Complex e))+ -> Acc (Array (sh:.Int) (Complex e))+ -> Acc (Array sh (Complex e))+dotc xs ys = A.sum $ A.zipWith (*) xs (A.map conjugate ys)++scalar :: Elt e => e -> Scalar e+scalar x = fromFunction Z (const x)+++test_homogeneity+ :: (Numeric e, Similar e, Shape sh, Eq sh)+ => RunN+ -> Transform sh (Complex e)+ -> Gen sh+ -> Gen e+ -> Property+test_homogeneity runN transform dim e =+ property $ do+ sign <- forAll mode+ sh <- forAll dim+ arr <- forAll (array sh (complex e))+ x <- forAll (complex e)+ --+ let !go1 = runN (\u -> transform sign . A.map (the u *))+ !go2 = runN (\u -> A.map (the u *) . transform sign)+ --+ go1 (scalar x) arr ~~~ go2 (scalar x) arr++test_additivity+ :: (Numeric e, Similar e, Shape sh, Eq sh)+ => RunN+ -> Transform sh (Complex e)+ -> Gen sh+ -> Gen e+ -> Property+test_additivity runN transform dim e =+ property $ do+ sign <- forAll mode+ sh <- forAll dim+ xs <- forAll (array sh (complex e))+ ys <- forAll (array sh (complex e))+ --+ let !go1 = runN (\u v -> transform sign (A.zipWith (+) u v))+ !go2 = runN (\u v -> A.zipWith (+) (transform sign u) (transform sign v))+ --+ go1 xs ys ~~~ go2 xs ys++test_inverse+ :: (Numeric e, Similar e, Shape sh, Eq sh)+ => RunN+ -> Transform sh (Complex e)+ -> Gen sh+ -> Gen e+ -> Property+test_inverse runN transform dim e =+ property $ do+ sh <- forAll dim+ xs <- forAll (array sh (complex e))+ --+ let !go = runN (transform Inverse . transform Forward)+ xs ~~~ go xs++test_reverse+ :: (Numeric e, Similar e, Shape sh, Slice sh, Eq sh)+ => RunN+ -> Transform (sh:.Int) (Complex e)+ -> Gen (sh:.Int)+ -> Gen e+ -> Property+test_reverse runN transform dim e =+ property $ do+ sign <- forAll mode+ sh <- forAll dim+ xs <- forAll (array sh (complex e))+ --+ let !go1 = runN (reverse . transform sign)+ !go2 = runN (transform sign . reverse)+ --+ go1 xs ~~~ go2 xs++test_conjugate+ :: (Numeric e, Similar e, Shape sh, Slice sh, Eq sh)+ => RunN+ -> Transform (sh:.Int) (Complex e)+ -> Gen (sh:.Int)+ -> Gen e+ -> Property+test_conjugate runN transform dim e =+ property $ do+ sign <- forAll mode+ sh <- forAll dim+ xs <- forAll (array sh (complex e))+ --+ let !go1 = runN (A.map conjugate . transform sign)+ !go2 = runN (transform sign . A.map conjugate . reverse)+ --+ go1 xs ~~~ go2 xs++test_isometry+ :: forall sh e. (Numeric e, Similar e, Shape sh, Slice sh, Eq sh, P.Floating e)+ => RunN+ -> Transform (sh:.Int) (Complex e)+ -> Gen (sh:.Int)+ -> Gen e+ -> Property+test_isometry runN transform dim e =+ property $ do+ sign <- forAll (Gen.element [Forward, Reverse])+ sh@(_:.n) <- forAll dim+ xs <- forAll (array sh (complex e))+ --+ let !go1 = runN (norm2 . transform sign)+ !go2 = runN (\u -> A.map (the u *) . norm2)+ --+ go1 xs ~~~ go2 (scalar (sqrt (P.fromIntegral n))) xs++test_unitarity+ :: forall sh e. (Numeric e, Similar e, RealFloat e, Shape sh, Slice sh, Eq sh)+ => RunN+ -> Transform (sh:.Int) (Complex e)+ -> Gen (sh:.Int)+ -> Gen e+ -> Property+test_unitarity runN transform dim e =+ property $ do+ sign <- forAll (Gen.element [Forward, Reverse])+ sh@(_:.n) <- forAll dim+ xs <- forAll (array sh (complex e))+ ys <- forAll (array sh (complex e))+ --+ let !go1 = runN (\u v -> dotc (transform sign u) (transform sign v))+ !go2 = runN (\m u v -> A.map (the m *) (dotc u v))+ --+ go1 xs ys ~~~ go2 (scalar (P.fromIntegral n)) xs ys+
+ test/Test/ShowType.hs view
@@ -0,0 +1,30 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE PolyKinds #-}+{-# LANGUAGE ScopedTypeVariables #-}+-- |+-- Module : Test.ShowType+-- Copyright : [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module Test.ShowType+ where++import Data.Complex++data ArgType (a :: *) = AT++showType :: forall proxy a. Show (ArgType a) => proxy a -> String+showType _ = show (AT :: ArgType a)++instance Show (ArgType a) => Show (ArgType (Complex a)) where+ show _ = "Complex " ++ show (AT :: ArgType a)++instance Show (ArgType Float) where show _ = "Float"+instance Show (ArgType Double) where show _ = "Double"+
+ test/TestNative.hs view
@@ -0,0 +1,19 @@+-- |+-- Module : TestNative+-- Copyright : [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module TestNative where++import Test.FFT+import Test.Tasty+import Data.Array.Accelerate.LLVM.Native as CPU++main :: IO ()+main = defaultMain (testFFT CPU.runN)+
+ test/TestPTX.hs view
@@ -0,0 +1,19 @@+-- |+-- Module : TestPTX+-- Copyright : [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module TestPTX where++import Test.FFT+import Test.Tasty+import Data.Array.Accelerate.LLVM.PTX as PTX++main :: IO ()+main = defaultMain (testFFT PTX.runN)+