diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -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 diff --git a/Data/Array/Accelerate/Math/DFT.hs b/Data/Array/Accelerate/Math/DFT.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Math/DFT.hs +++ /dev/null @@ -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 --- 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' - diff --git a/Data/Array/Accelerate/Math/DFT/Centre.hs b/Data/Array/Accelerate/Math/DFT/Centre.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Math/DFT/Centre.hs +++ /dev/null @@ -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 --- 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) - diff --git a/Data/Array/Accelerate/Math/DFT/Roots.hs b/Data/Array/Accelerate/Math/DFT/Roots.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Math/DFT/Roots.hs +++ /dev/null @@ -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 --- 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 )) - diff --git a/Data/Array/Accelerate/Math/FFT.hs b/Data/Array/Accelerate/Math/FFT.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Math/FFT.hs +++ /dev/null @@ -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 --- 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 - - - diff --git a/Data/Array/Accelerate/Math/FFT/LLVM/Native.hs b/Data/Array/Accelerate/Math/FFT/LLVM/Native.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Math/FFT/LLVM/Native.hs +++ /dev/null @@ -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 --- 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" - diff --git a/Data/Array/Accelerate/Math/FFT/LLVM/PTX.hs b/Data/Array/Accelerate/Math/FFT/LLVM/PTX.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Math/FFT/LLVM/PTX.hs +++ /dev/null @@ -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 --- 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 - diff --git a/Data/Array/Accelerate/Math/FFT/Mode.hs b/Data/Array/Accelerate/Math/FFT/Mode.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Math/FFT/Mode.hs +++ /dev/null @@ -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 --- 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 - diff --git a/Data/Array/Accelerate/Math/FFT/Twine.hs b/Data/Array/Accelerate/Math/FFT/Twine.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Math/FFT/Twine.hs +++ /dev/null @@ -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 --- 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 - diff --git a/accelerate-fft.cabal b/accelerate-fft.cabal --- a/accelerate-fft.cabal +++ b/accelerate-fft.cabal @@ -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 - diff --git a/cubits/twine_f32.cu b/cubits/twine_f32.cu deleted file mode 100644 --- a/cubits/twine_f32.cu +++ /dev/null @@ -1,63 +0,0 @@ -/* - * Module : Twine - * Copyright : [2016] Trevor L. McDonell - * License : BSD3 - * - * Maintainer : Trevor L. McDonell - * 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 -#include - -#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 - diff --git a/cubits/twine_f32.ptx b/cubits/twine_f32.ptx deleted file mode 100644 --- a/cubits/twine_f32.ptx +++ /dev/null @@ -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; -} - - diff --git a/cubits/twine_f64.cu b/cubits/twine_f64.cu deleted file mode 100644 --- a/cubits/twine_f64.cu +++ /dev/null @@ -1,63 +0,0 @@ -/* - * Module : Twine - * Copyright : [2016] Trevor L. McDonell - * License : BSD3 - * - * Maintainer : Trevor L. McDonell - * 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 -#include - -#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 - diff --git a/cubits/twine_f64.ptx b/cubits/twine_f64.ptx deleted file mode 100644 --- a/cubits/twine_f64.ptx +++ /dev/null @@ -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; -} - - diff --git a/src/Data/Array/Accelerate/Math/DFT.hs b/src/Data/Array/Accelerate/Math/DFT.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Math/DFT.hs @@ -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 +-- 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' + diff --git a/src/Data/Array/Accelerate/Math/DFT/Centre.hs b/src/Data/Array/Accelerate/Math/DFT/Centre.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Math/DFT/Centre.hs @@ -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 +-- 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) + diff --git a/src/Data/Array/Accelerate/Math/DFT/Roots.hs b/src/Data/Array/Accelerate/Math/DFT/Roots.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Math/DFT/Roots.hs @@ -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 +-- 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 )) + diff --git a/src/Data/Array/Accelerate/Math/FFT.hs b/src/Data/Array/Accelerate/Math/FFT.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Math/FFT.hs @@ -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 +-- 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 +--} + diff --git a/src/Data/Array/Accelerate/Math/FFT/Adhoc.hs b/src/Data/Array/Accelerate/Math/FFT/Adhoc.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Math/FFT/Adhoc.hs @@ -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 +-- 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. +-- +-- +-- +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)) + diff --git a/src/Data/Array/Accelerate/Math/FFT/LLVM/Native.hs b/src/Data/Array/Accelerate/Math/FFT/LLVM/Native.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Math/FFT/LLVM/Native.hs @@ -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 +-- 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) + diff --git a/src/Data/Array/Accelerate/Math/FFT/LLVM/Native/Base.hs b/src/Data/Array/Accelerate/Math/FFT/LLVM/Native/Base.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Math/FFT/LLVM/Native/Base.hs @@ -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 +-- 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 + diff --git a/src/Data/Array/Accelerate/Math/FFT/LLVM/Native/Ix.hs b/src/Data/Array/Accelerate/Math/FFT/LLVM/Native/Ix.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Math/FFT/LLVM/Native/Ix.hs @@ -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 +-- 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" + diff --git a/src/Data/Array/Accelerate/Math/FFT/LLVM/PTX.hs b/src/Data/Array/Accelerate/Math/FFT/LLVM/PTX.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Math/FFT/LLVM/PTX.hs @@ -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 +-- 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) + diff --git a/src/Data/Array/Accelerate/Math/FFT/LLVM/PTX/Base.hs b/src/Data/Array/Accelerate/Math/FFT/LLVM/PTX/Base.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Math/FFT/LLVM/PTX/Base.hs @@ -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 +-- 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 + diff --git a/src/Data/Array/Accelerate/Math/FFT/LLVM/PTX/Plans.hs b/src/Data/Array/Accelerate/Math/FFT/LLVM/PTX/Plans.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Math/FFT/LLVM/PTX/Plans.hs @@ -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 +-- 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 +-- +-- +-- +{-# 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#) + diff --git a/src/Data/Array/Accelerate/Math/FFT/Mode.hs b/src/Data/Array/Accelerate/Math/FFT/Mode.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Math/FFT/Mode.hs @@ -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 +-- 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 + diff --git a/src/Data/Array/Accelerate/Math/FFT/Type.hs b/src/Data/Array/Accelerate/Math/FFT/Type.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Math/FFT/Type.hs @@ -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 +-- 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 + diff --git a/test/Test/Base.hs b/test/Test/Base.hs new file mode 100644 --- /dev/null +++ b/test/Test/Base.hs @@ -0,0 +1,60 @@ +{-# LANGUAGE RankNTypes #-} +-- | +-- Module : Test.Base +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- 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 + diff --git a/test/Test/FFT.hs b/test/Test/FFT.hs new file mode 100644 --- /dev/null +++ b/test/Test/FFT.hs @@ -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 +-- 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 + diff --git a/test/Test/ShowType.hs b/test/Test/ShowType.hs new file mode 100644 --- /dev/null +++ b/test/Test/ShowType.hs @@ -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 +-- 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" + diff --git a/test/TestNative.hs b/test/TestNative.hs new file mode 100644 --- /dev/null +++ b/test/TestNative.hs @@ -0,0 +1,19 @@ +-- | +-- Module : TestNative +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- 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) + diff --git a/test/TestPTX.hs b/test/TestPTX.hs new file mode 100644 --- /dev/null +++ b/test/TestPTX.hs @@ -0,0 +1,19 @@ +-- | +-- Module : TestPTX +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- 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) +