accelerate-fft 0.15.1.0 → 1.0.0.0
raw patch · 14 files changed
+1507/−278 lines, 14 filesdep +accelerate-llvmdep +accelerate-llvm-nativedep +accelerate-llvm-ptxdep ~acceleratedep ~accelerate-cudadep ~base
Dependencies added: accelerate-llvm, accelerate-llvm-native, accelerate-llvm-ptx, bytestring, carray, fft, file-embed, storable-complex
Dependency ranges changed: accelerate, accelerate-cuda, base
Files
- Data/Array/Accelerate/Math/DFT.hs +12/−10
- Data/Array/Accelerate/Math/DFT/Centre.hs +16/−16
- Data/Array/Accelerate/Math/DFT/Roots.hs +8/−6
- Data/Array/Accelerate/Math/FFT.hs +151/−226
- Data/Array/Accelerate/Math/FFT/CUDA.hs +249/−0
- Data/Array/Accelerate/Math/FFT/LLVM/Native.hs +255/−0
- Data/Array/Accelerate/Math/FFT/LLVM/PTX.hs +279/−0
- Data/Array/Accelerate/Math/FFT/Mode.hs +28/−0
- Data/Array/Accelerate/Math/FFT/Twine.hs +88/−0
- accelerate-fft.cabal +79/−20
- cubits/twine_f32.cu +63/−0
- cubits/twine_f32.ptx +108/−0
- cubits/twine_f64.cu +63/−0
- cubits/twine_f64.ptx +108/−0
Data/Array/Accelerate/Math/DFT.hs view
@@ -1,12 +1,14 @@+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeOperators #-} -- | -- Module : Data.Array.Accelerate.Math.DFT--- Copyright : [2012..2014] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell--- [2013..2014] Robert Clifton-Everest+-- Copyright : [2012..2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- [2013..2017] Robert Clifton-Everest -- License : BSD3 ----- Maintainer : Manuel M T Chakravarty <chak@cse.unsw.edu.au>+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au> -- Stability : experimental -- Portability : non-portable (GHC extensions) --@@ -38,7 +40,7 @@ -- | Compute the DFT along the low order dimension of an array ---dft :: (Shape sh, Slice sh, Elt e, IsFloating e)+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@@ -46,14 +48,14 @@ -- | Compute the inverse DFT along the low order dimension of an array ---idft :: (Shape sh, Slice sh, Elt e, IsFloating e)+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 :+ constant 0)+ scale = lift (A.fromIntegral n :+ 0) in A.map (/scale) $ dftG roots v @@ -64,12 +66,12 @@ -- -- The extent of the input and roots must match. ---dftG :: forall sh e. (Shape sh, Slice sh, Elt e, IsFloating e)+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 (+) (constant (0 :+ 0))+ = A.fold (+) 0 $ A.zipWith (*) arr' roots' where base = shape arr@@ -93,7 +95,7 @@ -- | Compute a single value of the DFT. ---dftGS :: forall sh e. (Shape sh, Slice sh, Elt e, IsFloating e)+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@@ -107,5 +109,5 @@ (\ix' -> let sh :. n = unlift ix' :: Exp sh :. Exp Int in roots ! lift (sh :. (k*n) `mod` l)) in- A.foldAll (+) (constant (0 :+ 0)) $ A.zipWith (*) arr roots'+ A.foldAll (+) 0 $ A.zipWith (*) arr roots'
Data/Array/Accelerate/Math/DFT/Centre.hs view
@@ -1,12 +1,13 @@+{-# LANGUAGE ConstraintKinds #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE TypeOperators #-} -- | -- Module : Data.Array.Accelerate.Math.DFT.Centre--- Copyright : [2012..2014] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell--- [2013..2014] Robert Clifton-Everest+-- Copyright : [2012..2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- [2013..2017] Robert Clifton-Everest -- License : BSD3 ----- Maintainer : Manuel M T Chakravarty <chak@cse.unsw.edu.au>+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au> -- Stability : experimental -- Portability : non-portable (GHC extensions) --@@ -32,33 +33,33 @@ -- | Apply the centring transform to a vector ---centre1D :: (Elt e, IsFloating e)+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) :+ A.constant 0) * arr!ix)+ in lift (((-1) ** A.fromIntegral x) :+ 0) * arr!ix) -- | Apply the centring transform to a matrix ---centre2D :: (Elt e, IsFloating e)+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)) :+ A.constant 0) * arr!ix)+ in lift (((-1) ** A.fromIntegral (y + x)) :+ 0) * arr!ix) -- | Apply the centring transform to a 3D array ---centre3D :: (Elt e, IsFloating e)+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)) :+ A.constant 0) * arr!ix)+ in lift (((-1) ** A.fromIntegral (z + y + x)) :+ 0) * arr!ix) -- | Apply the shifting transform to a vector@@ -69,7 +70,7 @@ where p ix = let Z:.x = unlift ix :: Z :. Exp Int- in index1 (x A.<* mw ? (x + mw, x - mw))+ in index1 (x A.< mw ? (x + mw, x - mw)) Z:.w = unlift (A.shape arr) mw = w `div` 2 @@ -82,8 +83,8 @@ 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))+ 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) @@ -96,10 +97,9 @@ 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))+ 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)- index3 i j k = lift (Z:.i:.j:.k)
Data/Array/Accelerate/Math/DFT/Roots.hs view
@@ -1,11 +1,13 @@-{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE TypeOperators #-} -- | -- Module : Data.Array.Accelerate.Math.DFT.Roots--- Copyright : [2012..2014] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell--- [2013..2014] Robert Clifton-Everest+-- Copyright : [2012..2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- [2013..2017] Robert Clifton-Everest -- License : BSD3 ----- Maintainer : Manuel M T Chakravarty <chak@cse.unsw.edu.au>+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au> -- Stability : experimental -- Portability : non-portable (GHC extensions) --@@ -23,7 +25,7 @@ -- | Calculate the roots of unity for the forward transform -- rootsOfUnity- :: (Elt e, IsFloating e, Shape sh, Slice sh)+ :: (Shape sh, Slice sh, A.Floating e, A.FromIntegral Int e) => Exp (sh :. Int) -> Acc (Array (sh:.Int) (Complex e)) rootsOfUnity sh =@@ -38,7 +40,7 @@ -- | Calculate the roots of unity for an inverse transform -- inverseRootsOfUnity- :: (Elt e, IsFloating e, Shape sh, Slice sh)+ :: (Shape sh, Slice sh, A.Floating e, A.FromIntegral Int e) => Exp (sh :. Int) -> Acc (Array (sh:.Int) (Complex e)) inverseRootsOfUnity sh =
Data/Array/Accelerate/Math/FFT.hs view
@@ -1,5 +1,7 @@ {-# LANGUAGE CPP #-}+{-# LANGUAGE ConstraintKinds #-} {-# LANGUAGE EmptyDataDecls #-}+{-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE ForeignFunctionInterface #-} {-# LANGUAGE GADTs #-} {-# LANGUAGE ScopedTypeVariables #-}@@ -8,23 +10,32 @@ {-# LANGUAGE ViewPatterns #-} -- | -- Module : Data.Array.Accelerate.Math.FFT--- Copyright : [2012..2013] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell--- [2013..2014] Robert Clifton-Everest+-- Copyright : [2012..2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- [2013..2017] Robert Clifton-Everest -- License : BSD3 ----- Maintainer : Manuel M T Chakravarty <chak@cse.unsw.edu.au>+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au> -- Stability : experimental -- Portability : non-portable (GHC extensions) -- -- Computation of a Discrete Fourier Transform using the Cooley-Tuckey -- algorithm. The time complexity is O(n log n) in the size of the input. ----- This uses a naïve divide-and-conquer algorithm whose absolute performance is--- appalling.+-- 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 '-fcuda', '-fllvm-gpu', and '-fllvm-cpu' for the accelerate-cuda,+-- 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',@@ -32,182 +43,193 @@ ) where -import Data.Array.Accelerate as A-import Data.Array.Accelerate.Array.Sugar ( showShape )+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 #ifdef ACCELERATE_CUDA_BACKEND-import Data.Array.Accelerate.CUDA.Foreign-import Data.Array.Accelerate.Array.Sugar as S ( shapeToList, shape, EltRepr )-import Data.Array.Accelerate.Type--import Data.Functor-import System.Mem.Weak-import System.IO.Unsafe-import Foreign.CUDA.FFT-import qualified Foreign.CUDA.Types as CUDA-import qualified Foreign.CUDA.Driver as CUDA+import qualified Data.Array.Accelerate.Math.FFT.CUDA as CUDA #endif import Data.Bits-import Prelude as P+import Text.Printf+import Prelude as P -data Mode = Forward | Reverse | Inverse- deriving (Eq, Show) -isPow2 :: Int -> Bool-isPow2 x = x .&. (x-1) == 0--signOfMode :: Num a => Mode -> a-signOfMode m- = case m of- Forward -> -1- Reverse -> 1- Inverse -> 1+-- 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. ----- Discrete Fourier Transform of a vector. Array dimensions must be powers of--- two else error.+-- The default implementation requires the array dimension to be a power of two+-- (else error). ---fft1D :: (Elt e, IsFloating e)+fft1D :: FFTElt e => Mode- -> Vector (Complex e)- -> Acc (Vector (Complex e))+ -> Array DIM1 (Complex e)+ -> Acc (Array DIM1 (Complex e)) fft1D mode vec- = let Z :. len = arrayShape vec- in- fft1D' mode len (use vec)+ = fft1D' mode (arrayShape vec) (use vec) -fft1D' :: forall e. (Elt e, IsFloating e)++-- | 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- -> Int- -> Acc (Vector (Complex e))- -> Acc (Vector (Complex e))-fft1D' mode len vec+ -> DIM1+ -> Acc (Array DIM1 (Complex e))+ -> Acc (Array DIM1 (Complex e))+fft1D' mode (Z :. len) arr = let sign = signOfMode mode :: e scale = P.fromIntegral len+ go =+#ifdef ACCELERATE_LLVM_NATIVE_BACKEND+ foreignAcc (Native.fft1D mode) $+#endif+#ifdef ACCELERATE_LLVM_PTX_BACKEND+ foreignAcc (PTX.fft1D mode) $+#endif #ifdef ACCELERATE_CUDA_BACKEND- sh = (Z:.len)- vec' = cudaFFT mode sh fft' vec-#else- vec' = fft' vec+ foreignAcc (CUDA.fft1D mode) $ #endif- fft' a = fft sign Z len a+ fft sign Z len in- if P.not (isPow2 len)- then error $ unlines- [ "Data.Array.Accelerate.FFT: fft1D"- , " Array dimensions must be powers of two, but are: " P.++ showShape (Z:.len) ]-- else case mode of- Inverse -> A.map (/scale) vec'- _ -> vec'+ case mode of+ Inverse -> A.map (/scale) (go arr)+ _ -> go arr -- Matrix Transform -- ----------------++-- | Discrete Fourier Transform of a matrix. ----- Discrete Fourier Transform of a matrix. Array dimensions must be powers of--- two else error.+-- The default implementation requires the array dimensions to be powers of two+-- (else error). ---fft2D :: (Elt e, IsFloating e)+fft2D :: FFTElt e => Mode -> Array DIM2 (Complex e) -> Acc (Array DIM2 (Complex e)) fft2D mode arr- = let Z :. height :. width = arrayShape arr- in- fft2D' mode width height (use arr)+ = fft2D' mode (arrayShape arr) (use arr) -fft2D' :: forall e. (Elt e, IsFloating e)+-- | 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- -> Int -- ^ width- -> Int -- ^ height+ -> DIM2 -> Acc (Array DIM2 (Complex e)) -> Acc (Array DIM2 (Complex e))-fft2D' mode width height arr+fft2D' mode (Z :. height :. width) arr = let sign = signOfMode mode :: e scale = P.fromIntegral (width * height)+ go =+#ifdef ACCELERATE_LLVM_NATIVE_BACKEND+ foreignAcc (Native.fft2D mode) $+#endif+#ifdef ACCELERATE_LLVM_PTX_BACKEND+ foreignAcc (PTX.fft2D mode) $+#endif #ifdef ACCELERATE_CUDA_BACKEND- sh = (Z:.height:.width)- arr' = cudaFFT mode sh fft' arr-#else- arr' = fft' arr+ foreignAcc (CUDA.fft2D mode) $ #endif- fft' a = A.transpose . fft sign (Z:.width) height- >-> A.transpose . fft sign (Z:.height) width+ fft'++ fft' a = A.transpose . fft sign (Z:.height) width+ >-> A.transpose . fft sign (Z:.width) height $ a in- if P.not (isPow2 width && isPow2 height)- then error $ unlines- [ "Data.Array.Accelerate.FFT: fft2D"- , " Array dimensions must be powers of two, but are: " P.++ showShape (Z:.height:.width) ]-- else case mode of- Inverse -> A.map (/scale) arr'- _ -> arr'+ case mode of+ Inverse -> A.map (/scale) (go arr)+ _ -> go arr -- Cube Transform -- --------------++-- | Discrete Fourier Transform of a 3D array. ----- Discrete Fourier Transform of a 3D array. Array dimensions must be power of--- two else error.+-- The default implementation requires the array dimensions to be powers of two+-- (else error). ---fft3D :: (Elt e, IsFloating e)+fft3D :: FFTElt e => Mode -> Array DIM3 (Complex e) -> Acc (Array DIM3 (Complex e)) fft3D mode arr- = let Z :. depth :. height :. width = arrayShape arr- in- fft3D' mode width height depth (use arr)+ = fft3D' mode (arrayShape arr) (use arr) -fft3D' :: forall e. (Elt e, IsFloating e)+-- | 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- -> Int -- ^ width- -> Int -- ^ height- -> Int -- ^ depth+ -> DIM3 -> Acc (Array DIM3 (Complex e)) -> Acc (Array DIM3 (Complex e))-fft3D' mode width height depth arr+fft3D' mode (Z :. depth :. height :. width) arr = let sign = signOfMode mode :: e- scale = P.fromIntegral (width * height)+ scale = P.fromIntegral (width * height * depth)+ go =+#ifdef ACCELERATE_LLVM_NATIVE_BACKEND+ foreignAcc (Native.fft3D mode) $+#endif+#ifdef ACCELERATE_LLVM_PTX_BACKEND+ foreignAcc (PTX.fft3D mode) $+#endif #ifdef ACCELERATE_CUDA_BACKEND- sh = (Z:.depth:.height:.width)- arr' = cudaFFT mode sh fft' arr-#else- arr' = fft' arr+ foreignAcc (CUDA.fft3D mode) $ #endif- fft' a = rotate3D . fft sign (Z:.width :.depth) height+ fft'++ fft' a = rotate3D . fft sign (Z:.depth :.height) width >-> rotate3D . fft sign (Z:.height:.width) depth- >-> rotate3D . fft sign (Z:.depth :.height) width+ >-> rotate3D . fft sign (Z:.width :.depth) height $ a in- if P.not (isPow2 width && isPow2 height && isPow2 depth)- then error $ unlines- [ "Data.Array.Accelerate.FFT: fft3D"- , " Array dimensions must be powers of two, but are: " P.++ showShape (Z:.depth:.height:.width) ]-- else case mode of- Inverse -> A.map (/scale) arr'- _ -> arr'-+ 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 (swap (A.shape arr)) swap arr+rotate3D arr = backpermute sh rot arr where- swap :: Exp DIM3 -> Exp DIM3- swap ix =- let Z :. m :. k :. l = unlift ix :: Z :. Exp Int :. Exp Int :. Exp Int- in lift $ Z :. k :. l :. m+ 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@@ -215,17 +237,22 @@ -- We require the innermost dimension be passed as a Haskell value because we -- can't do divide-and-conquer recursion directly in the meta-language. ---fft :: forall sh e. (Slice sh, Shape sh, IsFloating e, Elt e)+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 = go sz 0 1+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 == 2+ | len P.== 2 = A.generate (constant (sh :. len)) swivel | otherwise@@ -241,8 +268,8 @@ swivel ix = let sh' :. sz' = unlift ix :: Exp sh :. Exp Int in- sz' ==* 0 ? ( (arr ! lift (sh' :. offset')) + (arr ! lift (sh' :. offset' + stride'))- {- ==* 1-} , (arr ! lift (sh' :. offset')) - (arr ! lift (sh' :. offset' + stride')) )+ 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)@@ -256,87 +283,9 @@ in lift ( cos k :+ A.constant sign * sin k ) -#ifdef ACCELERATE_CUDA_BACKEND--- FFT using the CUFFT library to enable high performance for the CUDA backend of--- Accelerate. The implementation works on all arrays of rank less than or equal--- to 3. The result is un-normalised.----cudaFFT :: forall e sh. (Shape sh, Elt e, IsFloating e)- => Mode- -> sh- -> (Acc (Array sh (Complex e)) -> Acc (Array sh (Complex e)))- -> Acc (Array sh (Complex e))- -> Acc (Array sh (Complex e))-cudaFFT mode sh = cudaFFT'- where- -- Plan the FFT.- -- Doing this in unsafePerformIO so it is not reperformed every time the- -- AST is evaluated.- --- hndl = unsafePerformIO $ do- plan <- case shapeToList sh of- [width] -> plan1D width types 1- [width, height] -> plan2D height width types- [width, height, depth] -> plan3D depth height width types- _ -> error "Accelerate-fft cannot use CUFFT for arrays of dimensions higher than 3"- addFinalizer plan (destroy plan)- return plan - types = case (floatingType :: FloatingType e) of- TypeFloat{} -> C2C- TypeDouble{} -> Z2Z- TypeCFloat{} -> C2C- TypeCDouble{} -> Z2Z-- cudaFFT' p arr = deinterleave sh (foreignAcc ff pure (interleave arr))- where- ff = CUDAForeignAcc "foreignFFT" foreignFFT- -- Unfortunately the pure version of the function needs to be wrapped in- -- interleave and deinterleave to match how the foreign version works.- --- -- RCE: Do the interleaving and deinterleaving in foreignFFT- --- -- TLM: The interleaving might get fused into other parts of the- -- computation and thus be okay. We should really support multi types- -- such as float2 instead.- --- pure = interleave . p . deinterleave sh- sign = signOfMode mode :: Int-- foreignFFT :: CUDA.Stream -> Array DIM1 e -> CIO (Array DIM1 e)- foreignFFT stream arr' = do- output <- allocateArray (S.shape arr')- iptr <- floatingDevicePtr arr'- optr <- floatingDevicePtr output-- --Execute- liftIO $ do- setStream hndl stream- execute iptr optr-- return output-- execute :: CUDA.DevicePtr e -> CUDA.DevicePtr e -> IO ()- execute iptr optr- = case (floatingType :: FloatingType e) of- TypeFloat{} -> execC2C hndl iptr optr sign- TypeDouble{} -> execZ2Z hndl iptr optr sign- TypeCFloat{} -> execC2C hndl (CUDA.castDevPtr iptr) (CUDA.castDevPtr optr) sign- TypeCDouble{} -> execZ2Z hndl (CUDA.castDevPtr iptr) (CUDA.castDevPtr optr) sign-- floatingDevicePtr :: Vector e -> CIO (CUDA.DevicePtr e)- floatingDevicePtr v- = case (floatingType :: FloatingType e) of- TypeFloat{} -> singleDevicePtr v- TypeDouble{} -> singleDevicePtr v- TypeCFloat{} -> CUDA.castDevPtr <$> singleDevicePtr v- TypeCDouble{} -> CUDA.castDevPtr <$> singleDevicePtr v-- singleDevicePtr :: DevicePtrs (EltRepr e) ~ ((),CUDA.DevicePtr b) => Vector e -> CIO (CUDA.DevicePtr b)- singleDevicePtr v = P.snd <$> devicePtrsOfArray v-#endif---- Append two arrays. Doesn't do proper bounds checking or intersection...+-- Append 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)@@ -349,35 +298,11 @@ 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) ))+ in i A.< n ? (xs ! lift (sz:.i), ys ! lift (sz:.i-n) )) -#ifdef ACCELERATE_CUDA_BACKEND-{-# RULES- "interleave/deinterleave" forall sh x. deinterleave sh (interleave x) = x;- "deinterleave/interleave" forall sh x. interleave (deinterleave sh x) = x- #-}---- Interleave the real and imaginary components in a complex array and produce a--- flattened vector. This allows us to mimic the float2 structure used by CUFFT--- to store complex numbers.----{-# NOINLINE interleave #-}-interleave :: (Shape sh, Elt e) => Acc (Array sh (Complex e)) -> Acc (Vector e)-interleave arr = generate sh swizzle- where- sh = index1 (2 * A.size arr)- swizzle ix =- let i = indexHead ix- v = arr A.!! (i `div` 2)- in- i `mod` 2 ==* 0 ? (real v, imag v)+isPow2 :: Int -> Bool+isPow2 0 = True+isPow2 1 = False+isPow2 x = x .&. (x-1) P.== 0 --- Deinterleave a vector into a complex array. Assumes the array is even in length.----{-# NOINLINE deinterleave #-}-deinterleave :: (Shape sh, Elt e) => sh -> Acc (Vector e) -> Acc (Array sh (Complex e))-deinterleave (constant -> sh) arr =- generate sh (\ix -> let i = toIndex sh ix * 2- in lift (arr A.!! i :+ arr A.!! (i+1)))-#endif
+ Data/Array/Accelerate/Math/FFT/CUDA.hs view
@@ -0,0 +1,249 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ViewPatterns #-}+-- |+-- Module : Data.Array.Accelerate.Math.FFT.CUDA+-- Copyright : [2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--+--++module Data.Array.Accelerate.Math.FFT.CUDA (++ fft1D,+ fft2D,+ fft3D,++) where++import Data.Array.Accelerate.Math.FFT.Mode+import Data.Array.Accelerate.Math.FFT.Twine+import Data.Array.Accelerate.Data.Complex++import Data.Array.Accelerate.CUDA.Foreign+import Data.Array.Accelerate.Array.Sugar as S hiding ( allocateArray )+import Data.Array.Accelerate.Type++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 System.IO.Unsafe+++fft1D :: IsFloating e+ => Mode+ -> CUDAForeignAcc (Vector (Complex e) -> (Vector (Complex e)))+fft1D mode = CUDAForeignAcc "fft1D" $ liftAtoC (cuFFT mode)++fft2D :: IsFloating e+ => Mode+ -> CUDAForeignAcc (Array DIM2 (Complex e) -> (Array DIM2 (Complex e)))+fft2D mode = CUDAForeignAcc "fft2D" $ liftAtoC (cuFFT mode)++fft3D :: IsFloating e+ => Mode+ -> CUDAForeignAcc (Array DIM3 (Complex e) -> (Array DIM3 (Complex e)))+fft3D mode = CUDAForeignAcc "fft3D" $ liftAtoC (cuFFT mode)+++liftAtoC+ :: forall sh e. (Shape sh, IsFloating e)+ => (Stream -> Array (sh:.Int) e -> CIO (Array (sh:.Int) e))+ -> Stream+ -> Array (sh:.Int) (Complex e)+ -> CIO (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+ -> CIO (Array (sh:.Int) e)+cuFFT mode st arr =+ withScalarArrayPtr arr st $ \d_arr -> liftIO $ 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)+ -> CIO (Array (sh:.Int) e) -- this is really a packed array of (Vec2 e) type+a2c st arr | FloatingDict <- floatingDict (floatingType :: FloatingType e) = do+ let+ sh :. sz = shape arr+ n = size sh * sz+ --+ cs <- allocateArray (sh :. 2*sz)+ withComplexArrayPtrs arr st $ \d_re d_im -> do+ withScalarArrayPtr cs st $ \d_cs -> liftIO $ 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+ -> CIO (Array (sh:.Int) (Complex e))+c2a st cs | FloatingDict <- floatingDict (floatingType :: FloatingType e) = do+ let+ sh :. sz2 = shape cs+ sz = sz2 `quot` 2+ n = size sh * sz+ --+ arr <- allocateArray (sh :. sz)+ withComplexArrayPtrs arr st $ \d_re d_im -> do+ withScalarArrayPtr cs st $ \d_cs -> liftIO $ 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 -> CIO a)+ -> CIO a+withComplexArrayPtrs arr st k+ = case floatingType :: FloatingType e of+ TypeFloat{} -> withDevicePtrs arr (Just st) $ \(((),p1),p2) -> k p1 p2+ TypeDouble{} -> withDevicePtrs arr (Just st) $ \(((),p1),p2) -> k p1 p2+ TypeCDouble{} -> withDevicePtrs arr (Just st) $ \(((),p1),p2) -> k p1 p2+ TypeCFloat{} -> withDevicePtrs arr (Just st) $ \(((),p1),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 -> CIO a)+ -> CIO a+withScalarArrayPtr arr st k+ = case floatingType :: FloatingType e of+ TypeFloat{} -> withDevicePtrs arr (Just st) $ \p -> k p+ TypeDouble{} -> withDevicePtrs arr (Just st) $ \p -> k p+ TypeCDouble{} -> withDevicePtrs arr (Just st) $ \p -> k p+ TypeCFloat{} -> withDevicePtrs arr (Just st) $ \p -> k p+++-- Cache the FFT planning step for faster repeat evaluations.+{-# NOINLINE fft_plans #-}+fft_plans :: MVar [((FFT.Type, [Int]), FFT.Handle)]+fft_plans = unsafePerformIO $ do+ mv <- newMVar []+ _ <- mkWeakMVar mv+ $ withMVar mv+ $ mapM_ (\(_,p) -> FFT.destroy p)+ return mv++-- Cache the functions which convert between SoA and AoS format.+{-# NOINLINE ptx_twine_modules #-}+ptx_twine_modules :: MVar [((Int, CUDA.Context), CUDA.Module)]+ptx_twine_modules = unsafePerformIO $ do+ mv <- newMVar []+ _ <- mkWeakMVar mv+ $ withMVar mv+ $ mapM_ (\((_,ctx),mdl) -> bracket_ (CUDA.push ctx) CUDA.pop (CUDA.unload mdl))+ return mv+
+ Data/Array/Accelerate/Math/FFT/LLVM/Native.hs view
@@ -0,0 +1,255 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TemplateHaskell #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+-- |+-- Module : Data.Array.Accelerate.Math.FFT.LLVM.Native+-- Copyright : [2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module Data.Array.Accelerate.Math.FFT.LLVM.Native (++ fft1D,+ fft2D,+ fft3D,++) where++import Data.Array.Accelerate.Math.FFT.Mode++import Data.Array.Accelerate as A+import Data.Array.Accelerate.Type as A+import Data.Array.Accelerate.Array.Sugar as S+import Data.Array.Accelerate.Error as A+import Data.Array.Accelerate.Array.Data as A+import Data.Array.Accelerate.Array.Unique as A+import Data.Array.Accelerate.Data.Complex as A++import Data.Array.Accelerate.LLVM.Native.Foreign++import Data.Ix ( Ix )+import Data.Array.CArray ( CArray )+import qualified Data.Array.CArray as C++import Math.FFT.Base ( FFTWReal, Sign(..), Flag, measure, destroyInput )+import qualified Math.FFT as FFT++import Foreign.Ptr+import Foreign.Storable+import Foreign.Storable.Complex ()++import Data.Bits+import Text.Printf+import Prelude as P+++fft1D :: forall e. (Elt e, IsFloating e)+ => Mode+ -> ForeignAcc (Vector (Complex e) -> Vector (Complex e))+fft1D mode+ = ForeignAcc (nameOf mode (undefined::DIM1))+ $ case floatingType :: FloatingType e of+ TypeFloat{} -> liftIO . liftAtoC go+ TypeDouble{} -> liftIO . liftAtoC go+ TypeCFloat{} -> liftIO . liftAtoC go+ TypeCDouble{} -> liftIO . liftAtoC go+ where+ go :: FFTWReal r => CArray Int (Complex r) -> CArray Int (Complex r)+ go = FFT.dftGU (signOf mode) flags [0]+++fft2D :: forall e. (Elt e, IsFloating e)+ => Mode+ -> ForeignAcc (Array DIM2 (Complex e) -> Array DIM2 (Complex e))+fft2D mode+ = ForeignAcc (nameOf mode (undefined::DIM2))+ $ case floatingType :: FloatingType e of+ TypeFloat{} -> liftIO . liftAtoC go+ TypeDouble{} -> liftIO . liftAtoC go+ TypeCFloat{} -> liftIO . liftAtoC go+ TypeCDouble{} -> liftIO . liftAtoC go+ where+ go :: FFTWReal r => CArray (Int,Int) (Complex r) -> CArray (Int,Int) (Complex r)+ go = FFT.dftGU (signOf mode) flags [0,1]+++fft3D :: forall e. (Elt e, IsFloating e)+ => Mode+ -> ForeignAcc (Array DIM3 (Complex e) -> Array DIM3 (Complex e))+fft3D mode+ = ForeignAcc (nameOf mode (undefined::DIM3))+ $ case floatingType :: FloatingType e of+ TypeFloat{} -> liftIO . liftAtoC go+ TypeDouble{} -> liftIO . liftAtoC go+ TypeCFloat{} -> liftIO . liftAtoC go+ TypeCDouble{} -> liftIO . liftAtoC go+ where+ go :: FFTWReal r => CArray (Int,Int,Int) (Complex r) -> CArray (Int,Int,Int) (Complex r)+ go = FFT.dftGU (signOf mode) flags [0,1,2]+++signOf :: Mode -> Sign+signOf Forward = DFTForward+signOf _ = DFTBackward++flags :: Flag+flags = measure .|. destroyInput++nameOf :: forall sh. Shape sh => Mode -> sh -> String+nameOf Forward _ = printf "FFTW.dft%dD" (rank (undefined::sh))+nameOf _ _ = printf "FFTW.idft%dD" (rank (undefined::sh))+++-- | Lift an operation on CArray into an operation on Accelerate arrays+--+liftAtoC+ :: (IxShapeRepr (EltRepr ix) ~ EltRepr sh, Shape sh, Ix ix, Elt ix, Elt e, IsFloating e, Storable e', ArrayPtrs e ~ Ptr e')+ => (CArray ix (Complex e') -> CArray ix (Complex e'))+ -> Array sh (Complex e)+ -> IO (Array sh (Complex e))+liftAtoC f a = c2a . f =<< a2c a+++-- | Convert a multidimensional Accelerate array of complex numbers into+-- a packed CArray of complex numbers suitable for use by FFTW.+--+a2c :: forall ix sh e e'. (IxShapeRepr (EltRepr ix) ~ EltRepr sh, Ix ix, Elt ix, Shape sh, IsFloating e, Storable e', ArrayPtrs e ~ Ptr e')+ => Array sh (Complex e)+ -> IO (CArray ix (Complex e'))+a2c arr+ | FloatingDict <- floatingDict (floatingType :: FloatingType e)+ = let+ (lo,hi) = shapeToRange (arrayShape arr)+ bnds = (fromIxShapeRepr lo, fromIxShapeRepr hi)+ n = S.size (arrayShape arr)+ in+ C.createCArray bnds $ \p_cs ->+ withComplexArrayPtrs arr $ \p_re p_im ->+ let+ -- TLM: Should we execute this in parallel using the worker threads of+ -- the current target? (Native)+ go !i | i P.>= n = return ()+ go !i = do+ re <- peekElemOff p_re i+ im <- peekElemOff p_im i+ pokeElemOff p_cs i (re :+ im)+ go (i+1)+ in+ go 0+++-- | Convert a packed CArray of complex numbers into an unzipped (SoA)+-- multidimensional Accelerate array of complex numbers.+--+c2a :: forall ix sh e e'. (IxShapeRepr (EltRepr ix) ~ EltRepr sh, Ix ix, Elt ix, Shape sh, Elt e, IsFloating e, Storable e', ArrayPtrs e ~ Ptr e')+ => CArray ix (Complex e')+ -> IO (Array sh (Complex e))+c2a carr+ | FloatingDict <- floatingDict (floatingType :: FloatingType e)+ = let+ (lo,hi) = C.bounds carr+ n = C.rangeSize (lo,hi)+ sh = rangeToShape (toIxShapeRepr lo, toIxShapeRepr hi)+ in do+ arr <- allocateArray sh+ C.withCArray carr $ \p_cs -> do+ withComplexArrayPtrs arr $ \p_re p_im -> do+ let+ -- TLM: Should we execute this in parallel using the worker threads+ -- of the current target? (Native)+ go !i | i P.>= n = return ()+ go !i = do+ re :+ im <- peekElemOff p_cs i+ pokeElemOff p_re i re+ pokeElemOff p_im i im+ go (i+1)+ --+ go 0+ return arr+++-- Converting between Accelerate multidimensional shapes/indices and those used+-- by the CArray package (Data.Ix)+--++type family IxShapeRepr e where+ IxShapeRepr () = ()+ IxShapeRepr Int = ((),Int)+ IxShapeRepr (t,h) = (IxShapeRepr t, h)++fromIxShapeRepr+ :: forall ix sh. (IxShapeRepr (EltRepr ix) ~ EltRepr sh, Shape sh, Elt ix)+ => sh+ -> ix+fromIxShapeRepr = liftToElt (go (eltType (undefined::ix)))+ where+ go :: forall ix'. TupleType ix' -> IxShapeRepr ix' -> ix'+ go UnitTuple () = ()+ go (PairTuple tt _) (t, h) = (go tt t, h)+ go (SingleTuple (NumScalarType (IntegralNumType TypeInt{}))) ((),h) = h+ go _ _+ = $internalError "fromIxShapeRepr" "expected Int dimensions"++toIxShapeRepr+ :: forall ix sh. (IxShapeRepr (EltRepr ix) ~ EltRepr sh, Shape sh, Elt ix)+ => ix+ -> sh+toIxShapeRepr = liftToElt (go (eltType (undefined::ix)))+ where+ go :: forall ix'. TupleType ix' -> ix' -> IxShapeRepr ix'+ go UnitTuple () = ()+ go (SingleTuple (NumScalarType (IntegralNumType TypeInt{}))) h = ((), h)+ go (PairTuple tt _) (t, h) = (go tt t, h)+ go _ _+ = error "toIxShapeRepr: not a valid Data.Ix index"+++-- Dig out the underlying pointers of the Accelerate SoA data structure+--++withComplexArrayPtrs+ :: forall sh e a. IsFloating e+ => Array sh (Complex e)+ -> (ArrayPtrs e -> ArrayPtrs e -> IO a)+ -> IO a+withComplexArrayPtrs (Array _ adata) k+ | AD_Pair (AD_Pair AD_Unit ad1) ad2 <- adata+ = case floatingType :: FloatingType e of+ TypeFloat{} -> withArrayData arrayElt ad1 $ \p1 -> withArrayData arrayElt ad2 $ \p2 -> k p1 p2+ TypeDouble{} -> withArrayData arrayElt ad1 $ \p1 -> withArrayData arrayElt ad2 $ \p2 -> k p1 p2+ TypeCFloat{} -> withArrayData arrayElt ad1 $ \p1 -> withArrayData arrayElt ad2 $ \p2 -> k p1 p2+ TypeCDouble{} -> withArrayData arrayElt ad1 $ \p1 -> withArrayData arrayElt ad2 $ \p2 -> k p1 p2++-- withScalarArrayPtrs+-- :: forall sh e a. IsFloating e+-- => Array sh e+-- -> (ArrayPtrs e -> IO a)+-- -> IO a+-- withScalarArrayPtrs (Array _ adata) =+-- case floatingType :: FloatingType e of+-- TypeFloat{} -> withArrayData arrayElt adata+-- TypeDouble{} -> withArrayData arrayElt adata+-- TypeCFloat{} -> withArrayData arrayElt adata+-- TypeCDouble{} -> withArrayData arrayElt adata++withArrayData+ :: (ArrayPtrs e ~ Ptr a)+ => ArrayEltR e+ -> ArrayData e+ -> (Ptr a -> IO b)+ -> IO b+withArrayData ArrayEltRfloat (AD_Float ua) = withUniqueArrayPtr ua+withArrayData ArrayEltRdouble (AD_Double ua) = withUniqueArrayPtr ua+withArrayData ArrayEltRcfloat (AD_CFloat ua) = withUniqueArrayPtr ua+withArrayData ArrayEltRcdouble (AD_CDouble ua) = withUniqueArrayPtr ua+withArrayData _ _ =+ $internalError "withArrayData" "expected array of [C]Float or [C]Double"+
+ Data/Array/Accelerate/Math/FFT/LLVM/PTX.hs view
@@ -0,0 +1,279 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE PatternGuards #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TupleSections #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE ViewPatterns #-}+-- |+-- Module : Data.Array.Accelerate.Math.FFT.LLVM.PTX+-- Copyright : [2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module Data.Array.Accelerate.Math.FFT.LLVM.PTX (++ fft1D,+ fft2D,+ fft3D,++) where++import Data.Array.Accelerate.Math.FFT.Mode+import Data.Array.Accelerate.Math.FFT.Twine++import Data.Array.Accelerate.Array.Data+import Data.Array.Accelerate.Lifetime+import Data.Array.Accelerate.Array.Sugar+import Data.Array.Accelerate.Data.Complex+import Data.Array.Accelerate.Type++import Data.Array.Accelerate.LLVM.PTX.Foreign++import Foreign.CUDA.Ptr ( DevicePtr )+import Foreign.Ptr+import Foreign.Storable+import Foreign.CUDA.Analysis+import qualified Foreign.CUDA.FFT as FFT+import qualified Foreign.CUDA.Driver as CUDA hiding ( device )+import qualified Foreign.CUDA.Driver.Context as CUDA ( device )++import Control.Concurrent.MVar+import Control.Exception+import Control.Monad+import Data.Maybe+import Data.Typeable+import System.IO.Unsafe+++fft1D :: IsFloating e+ => Mode+ -> ForeignAcc (Vector (Complex e) -> (Vector (Complex e)))+fft1D mode = ForeignAcc "fft1D" $ liftAtoC (cuFFT mode)++fft2D :: IsFloating e+ => Mode+ -> ForeignAcc (Array DIM2 (Complex e) -> (Array DIM2 (Complex e)))+fft2D mode = ForeignAcc "fft2D" $ liftAtoC (cuFFT mode)++fft3D :: IsFloating e+ => Mode+ -> ForeignAcc (Array DIM3 (Complex e) -> (Array DIM3 (Complex e)))+fft3D mode = ForeignAcc "fft3D" $ liftAtoC (cuFFT mode)+++liftAtoC+ :: forall sh e. (Shape sh, IsFloating e)+ => (Stream -> Array (sh:.Int) e -> LLVM PTX (Array (sh:.Int) e))+ -> Stream+ -> Array (sh:.Int) (Complex e)+ -> LLVM PTX (Array (sh:.Int) (Complex e))+liftAtoC f s =+ case floatingType :: FloatingType e of+ TypeFloat{} -> c2a s <=< f s <=< a2c s+ TypeDouble{} -> c2a s <=< f s <=< a2c s+ TypeCFloat{} -> c2a s <=< f s <=< a2c s+ TypeCDouble{} -> c2a s <=< f s <=< a2c s+++-- | Call the cuFFT library to execute the FFT (inplace)+--+cuFFT :: forall sh e. (Shape sh, IsFloating e)+ => Mode+ -> Stream+ -> Array (sh:.Int) e+ -> LLVM PTX (Array (sh:.Int) e)+cuFFT mode stream arr =+ withScalarArrayPtr arr stream $ \d_arr -> liftIO $+ withLifetime stream $ \st -> do+ let sh :. sz = shape arr+ p <- plan (sh :. sz `quot` 2) (undefined::e) -- recall this is an array of packed (Vec2 e)+ FFT.setStream p st+ case floatingType :: FloatingType e of+ TypeFloat{} -> FFT.execC2C p d_arr d_arr (signOfMode mode) >> return arr+ TypeDouble{} -> FFT.execZ2Z p d_arr d_arr (signOfMode mode) >> return arr+ TypeCFloat{} -> FFT.execC2C p d_arr d_arr (signOfMode mode) >> return arr+ TypeCDouble{} -> FFT.execZ2Z p d_arr d_arr (signOfMode mode) >> return arr+++-- | Convert an unzipped Accelerate array of complex numbers into a (new) packed+-- array suitable for use with CUFFT.+--+a2c :: forall sh e. (Shape sh, Elt e, IsFloating e, Storable (DevicePtrs e))+ => Stream+ -> Array (sh:.Int) (Complex e)+ -> LLVM PTX (Array (sh:.Int) e) -- this is really a packed array of (Vec2 e) type+a2c stream arr | FloatingDict <- floatingDict (floatingType :: FloatingType e) = do+ let+ sh :. sz = shape arr+ n = size sh * sz+ --+ cs <- allocateRemote (sh :. 2*sz)+ withComplexArrayPtrs arr stream $ \d_re d_im -> do+ withScalarArrayPtr cs stream $ \d_cs -> liftIO $ do+ withLifetime stream $ \st -> do+ mdl <- twine (sizeOf (undefined::e))+ pack <- CUDA.getFun mdl "interleave"+ dev <- CUDA.device+ prp <- CUDA.props dev+ regs <- CUDA.requires pack CUDA.NumRegs+ let+ blockSize = 256+ sharedMem = 0+ maxBlocks = maxResidentBlocks prp blockSize regs sharedMem+ numBlocks = maxBlocks `min` ((n + blockSize - 1) `div` blockSize)+ --+ CUDA.launchKernel pack (numBlocks,1,1) (blockSize,1,1) sharedMem (Just st)+ [ CUDA.VArg d_cs, CUDA.VArg d_re, CUDA.VArg d_im, CUDA.IArg (fromIntegral n) ]+ return cs++-- | Convert a packed array of complex numbers into a (new) unzipped Accelerate+-- array.+--+c2a :: forall sh e. (Shape sh, Elt e, IsFloating e, Storable (DevicePtrs e))+ => Stream+ -> Array (sh:.Int) e+ -> LLVM PTX (Array (sh:.Int) (Complex e))+c2a stream cs | FloatingDict <- floatingDict (floatingType :: FloatingType e) = do+ let+ sh :. sz2 = shape cs+ sz = sz2 `quot` 2+ n = size sh * sz+ --+ arr <- allocateRemote (sh :. sz)+ withComplexArrayPtrs arr stream $ \d_re d_im -> do+ withScalarArrayPtr cs stream $ \d_cs -> liftIO $ do+ withLifetime stream $ \st -> do+ mdl <- twine (sizeOf (undefined::e))+ unpack <- CUDA.getFun mdl "deinterleave"+ dev <- CUDA.device+ prp <- CUDA.props dev+ regs <- CUDA.requires unpack CUDA.NumRegs+ let+ blockSize = 256+ sharedMem = 0+ maxBlocks = maxResidentBlocks prp blockSize regs sharedMem+ numBlocks = maxBlocks `min` ((n + blockSize - 1) `div` blockSize)+ --+ CUDA.launchKernel unpack (numBlocks,1,1) (blockSize,1,1) sharedMem (Just st)+ [ CUDA.VArg d_re, CUDA.VArg d_im, CUDA.VArg d_cs, CUDA.IArg (fromIntegral n) ]+ return arr+++-- | Generate an execute plan for a given type and size of FFT. These plans are+-- cached so that subsequent invocations are quicker.+--+plan :: forall sh e. (Shape sh, IsFloating e) => sh -> e -> IO FFT.Handle+plan (shapeToList -> sh) _ =+ modifyMVar fft_plans $ \ps ->+ case lookup (ty, sh) ps of+ Just p -> return (ps, p)+ Nothing -> do+ p <- case sh of+ [w] -> FFT.plan1D w ty 1+ [w,h] -> FFT.plan2D h w ty+ [w,h,d] -> FFT.plan3D d h w ty+ _ -> error "cuFFT only supports 1D, 2D, and 3D transforms"+ return (((ty,sh),p) : ps, p)+ where+ ty = case floatingType :: FloatingType e of+ TypeFloat{} -> FFT.C2C+ TypeDouble{} -> FFT.Z2Z+ TypeCFloat{} -> FFT.C2C+ TypeCDouble{} -> FFT.Z2Z+++-- | Load the module to convert between SoA and AoS representation for the given+-- type. This is cached for subsequent reuse.+--+twine :: Int -> IO CUDA.Module+twine bitsize = do+ ctx <- fromMaybe (error "could not determine current CUDA context") `fmap` CUDA.get+ modifyMVar ptx_twine_modules $ \ms -> do+ case lookup (bitsize,ctx) ms of+ Just m -> return (ms, m)+ Nothing -> do+ m <- CUDA.loadData $ case bitsize of+ 4 -> ptx_twine_f32+ 8 -> ptx_twine_f64+ _ -> error "cuFFT only supports Float and Double"+ return (((bitsize,ctx), m) : ms, m)+++-- | Dig out the two device pointers for an unzipped array of complex numbers.+--+withComplexArrayPtrs+ :: forall sh e a. IsFloating e+ => Array sh (Complex e)+ -> Stream+ -> (DevicePtrs e -> DevicePtrs e -> LLVM PTX a)+ -> LLVM PTX a+withComplexArrayPtrs (Array _ adata) st k+ | AD_Pair (AD_Pair AD_Unit ad1) ad2 <- adata+ = case floatingType :: FloatingType e of+ TypeFloat{} -> withArrayData arrayElt ad1 st $ \p1 -> withArrayData arrayElt ad2 st $ \p2 -> k p1 p2+ TypeDouble{} -> withArrayData arrayElt ad1 st $ \p1 -> withArrayData arrayElt ad2 st $ \p2 -> k p1 p2+ TypeCDouble{} -> withArrayData arrayElt ad1 st $ \p1 -> withArrayData arrayElt ad2 st $ \p2 -> k p1 p2+ TypeCFloat{} -> withArrayData arrayElt ad1 st $ \p1 -> withArrayData arrayElt ad2 st $ \p2 -> k p1 p2++-- | Dig out the device pointer for a scalar array+--+withScalarArrayPtr+ :: forall sh e a. IsFloating e+ => Array sh e+ -> Stream+ -> (DevicePtrs e -> LLVM PTX a)+ -> LLVM PTX a+withScalarArrayPtr (Array _ ad) st k+ = case floatingType :: FloatingType e of+ TypeFloat{} -> withArrayData arrayElt ad st $ \p -> k p+ TypeDouble{} -> withArrayData arrayElt ad st $ \p -> k p+ TypeCDouble{} -> withArrayData arrayElt ad st $ \p -> k p+ TypeCFloat{} -> withArrayData arrayElt ad st $ \p -> k p++withArrayData+ :: (Typeable e, Typeable a, ArrayElt e, Storable a, ArrayPtrs e ~ Ptr a)+ => ArrayEltR e+ -> ArrayData e+ -> Stream+ -> (DevicePtr a -> LLVM PTX b)+ -> LLVM PTX b+withArrayData _ ad s k =+ withDevicePtr ad $ \p -> do+ r <- k p+ e <- checkpoint s+ return (Just e,r)++type family DevicePtrs e :: *++type instance DevicePtrs Float = DevicePtr Float+type instance DevicePtrs Double = DevicePtr Double+type instance DevicePtrs CFloat = DevicePtr Float+type instance DevicePtrs CDouble = DevicePtr Double+++-- Cache the FFT planning step for faster repeat evaluations.+{-# NOINLINE fft_plans #-}+fft_plans :: MVar [((FFT.Type, [Int]), FFT.Handle)]+fft_plans = unsafePerformIO $ do+ mv <- newMVar []+ _ <- mkWeakMVar mv+ $ withMVar mv+ $ mapM_ (\(_,p) -> FFT.destroy p)+ return mv++-- Cache the functions which convert between SoA and AoS format.+{-# NOINLINE ptx_twine_modules #-}+ptx_twine_modules :: MVar [((Int, CUDA.Context), CUDA.Module)]+ptx_twine_modules = unsafePerformIO $ do+ mv <- newMVar []+ _ <- mkWeakMVar mv+ $ withMVar mv+ $ mapM_ (\((_,ctx),mdl) -> bracket_ (CUDA.push ctx) CUDA.pop (CUDA.unload mdl))+ return mv+
+ Data/Array/Accelerate/Math/FFT/Mode.hs view
@@ -0,0 +1,28 @@+-- |+-- Module : Data.Array.Accelerate.Math.FFT.Mode+-- Copyright : [2012..2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- [2013..2017] Robert Clifton-Everest+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module Data.Array.Accelerate.Math.FFT.Mode+ where+++data Mode+ = Forward -- ^ Forward DFT+ | Reverse -- ^ Inverse DFT, un-normalised+ | Inverse -- ^ Inverse DFT, normalised+ deriving (Eq, Show)++signOfMode :: Num a => Mode -> a+signOfMode m+ = case m of+ Forward -> -1+ Reverse -> 1+ Inverse -> 1+
+ Data/Array/Accelerate/Math/FFT/Twine.hs view
@@ -0,0 +1,88 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TemplateHaskell #-}+-- |+-- Module : Data.Array.Accelerate.Math.FFT.Twine+-- Copyright : [2017] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--+--++module Data.Array.Accelerate.Math.FFT.Twine+ where++import Data.Array.Accelerate as A+import Data.Array.Accelerate.Data.Complex++#if defined(ACCELERATE_CUDA_BACKEND) || defined(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+ #-}+++#if defined(ACCELERATE_CUDA_BACKEND) || defined(ACCELERATE_LLVM_PTX_BACKEND)++-- Embedded PTX code for interleave and deinterleave for 32- and 64-bit floating+-- point numbers respectively. These can be loaded and executed by the CUDA+-- driver at runtime as required.+--+-- The PTX code was compiled for SM-2.0 and 64-bit address space (the default+-- settings of nvcc-7.5), but the code is simple enough that the CUDA device+-- driver should be able to compile it for the actual target architecture+-- without issue. This has been confirmed with respect to SM, but I don't have+-- a 32-bit machine available to test that aspect with.+--++ptx_twine_f32 :: ByteString+ptx_twine_f32 = $(makeRelativeToProject "cubits/twine_f32.ptx" >>= embedFile)++ptx_twine_f64 :: ByteString+ptx_twine_f64 = $(makeRelativeToProject "cubits/twine_f64.ptx" >>= embedFile)++#endif+
accelerate-fft.cabal view
@@ -1,5 +1,5 @@ Name: accelerate-fft-Version: 0.15.1.0+Version: 1.0.0.0 Cabal-version: >= 1.6 Tested-with: GHC >= 7.8 Build-type: Simple@@ -7,7 +7,9 @@ Synopsis: FFT using the Accelerate library Description: Rank-polymorphic discrete Fourier transform (DFT), computed with a fast- Fourier transform (FFT) algorithm using the Accelerate library+ Fourier transform (FFT) algorithm using the Accelerate library. Note that+ optimised implementations are available via foreign libraries, but must be+ explicitly enabled. . Refer to the main /Accelerate/ package for more information: <http://hackage.haskell.org/package/accelerate>@@ -17,44 +19,101 @@ License-file: LICENSE Author: Manuel M T Chakravarty, Gabriele Keller,- Trevor L. McDonell-Maintainer: Manuel M T Chakravarty <chak@cse.unsw.edu.au>+ Trevor L. McDonell,+ Robert Clifton-Everest+Maintainer: Trevor L. McDonell <tmcdonell@cse.unsw.edu.au> Homepage: https://github.com/AccelerateHS/accelerate-fft Bug-reports: https://github.com/AccelerateHS/accelerate/issues Category: Compilers/Interpreters, Concurrency, Data, Parallelism Stability: Experimental +extra-source-files:+ cubits/twine_f32.ptx+ cubits/twine_f64.ptx+ cubits/twine_f32.cu+ cubits/twine_f64.cu+ Flag cuda- Description: Enable support for using CUFFT via the CUDA backend's FFI- Default: True+ Description: Use CUFFT-based implementation in the CUDA backend+ Default: False +Flag llvm-ptx+ Description: Use CUFFT-based implementation in the LLVM.PTX backend+ Default: False++Flag llvm-cpu+ Description: Use FFTW-based implementation in the LLVM.Native backend+ Default: False++ Library- Build-depends: accelerate == 0.15.*,- base >= 4.7 && < 4.9+ build-depends:+ base >= 4.7 && < 4.10+ , accelerate == 1.0.*+ , bytestring >= 0.9 - Exposed-modules: Data.Array.Accelerate.Math.FFT- Data.Array.Accelerate.Math.DFT- Data.Array.Accelerate.Math.DFT.Centre- Data.Array.Accelerate.Math.DFT.Roots+ exposed-modules:+ Data.Array.Accelerate.Math.FFT+ Data.Array.Accelerate.Math.DFT+ Data.Array.Accelerate.Math.DFT.Centre+ Data.Array.Accelerate.Math.DFT.Roots + other-modules:+ Data.Array.Accelerate.Math.FFT.Mode+ Data.Array.Accelerate.Math.FFT.Twine+ ghc-options: -O2 -Wall -funbox-strict-fields if flag(cuda)- CPP-options: -DACCELERATE_CUDA_BACKEND- Build-depends: accelerate-cuda == 0.15.*,- cuda >= 0.5,- cufft >= 0.1.2+ 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.*+ , accelerate-llvm-native == 1.0.*+ , carray >= 0.1.5+ , fft >= 0.1.8+ , storable-complex >= 0.2++ other-modules:+ Data.Array.Accelerate.Math.FFT.LLVM.Native++ if flag(llvm-ptx)+ cpp-options: -DACCELERATE_LLVM_PTX_BACKEND+ build-depends:+ accelerate-llvm == 1.0.*+ , accelerate-llvm-ptx == 1.0.*+ , cuda >= 0.5+ , cufft >= 0.1.2+ , file-embed >= 0.0.10++ other-modules:+ Data.Array.Accelerate.Math.FFT.LLVM.PTX+ -- Don't add the extensions list here. Instead, place individual LANGUAGE -- pragmas in the files that require a specific extension. This means the -- project loads in GHCi, and avoids extension clashes. -- -- Extensions: +Source-repository head+ Type: git+ Location: git://github.com/AccelerateHS/accelerate-fft.git+ Source-repository this- type: git- location: git://github.com/AccelerateHS/accelerate-fft.git- branch: release/0.15- tag: 0.15.1.0+ Type: git+ Tag: 1.0.0.0+ Location: git://github.com/AccelerateHS/accelerate-fft.git++-- vim: nospell
+ cubits/twine_f32.cu view
@@ -0,0 +1,63 @@+/*+ * Module : Twine+ * Copyright : [2016] Trevor L. McDonell+ * License : BSD3+ *+ * Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+ * Stability : experimental+ * Portability : non-portable (GHC extensions)+ *+ * Convert between Accelerate's Struct-of-Array representation of complex+ * numbers and the Array-of-Struct representation necessary for CUBLAS.+ *+ */++#include <cuda.h>+#include <cuComplex.h>++#ifdef __cplusplus+extern "C" {+#endif++__global__ void interleave+(+ cuFloatComplex * __restrict__ cplx,+ const float * __restrict__ real,+ const float * __restrict__ imag,+ const int size+)+{+ const int gridSize = blockDim.x * gridDim.x;+ int ix;++ for (ix = blockDim.x * blockIdx.x + threadIdx.x; ix < size; ix += gridSize) {+ const float re = real[ix];+ const float im = imag[ix];++ cplx[ix] = make_cuFloatComplex(re, im);+ }+}++__global__ void deinterleave+(+ float * __restrict__ real,+ float * __restrict__ imag,+ const cuFloatComplex * __restrict__ cplx,+ const int size+)+{+ const int gridSize = blockDim.x * gridDim.x;+ int ix;++ for (ix = blockDim.x * blockIdx.x + threadIdx.x; ix < size; ix += gridSize) {+ const cuFloatComplex c = cplx[ix];++ real[ix] = cuCrealf(c);+ imag[ix] = cuCimagf(c);+ }+}++#ifdef __cplusplus+}+#endif+
+ cubits/twine_f32.ptx view
@@ -0,0 +1,108 @@+//+// Generated by NVIDIA NVVM Compiler+//+// Compiler Build ID: CL-20633761+// Cuda compilation tools, release 7.5, V7.5.26+// Based on LLVM 3.4svn+//++.version 4.3+.target sm_20+.address_size 64++ // .globl interleave++.visible .entry interleave(+ .param .u64 interleave_param_0,+ .param .u64 interleave_param_1,+ .param .u64 interleave_param_2,+ .param .u32 interleave_param_3+)+{+ .reg .pred %p<3>;+ .reg .f32 %f<3>;+ .reg .b32 %r<11>;+ .reg .b64 %rd<12>;+++ ld.param.u64 %rd4, [interleave_param_0];+ ld.param.u64 %rd5, [interleave_param_1];+ ld.param.u64 %rd6, [interleave_param_2];+ ld.param.u32 %r5, [interleave_param_3];+ cvta.to.global.u64 %rd1, %rd4;+ cvta.to.global.u64 %rd2, %rd6;+ cvta.to.global.u64 %rd3, %rd5;+ mov.u32 %r6, %nctaid.x;+ mov.u32 %r7, %ntid.x;+ mul.lo.s32 %r1, %r6, %r7;+ mov.u32 %r8, %ctaid.x;+ mov.u32 %r9, %tid.x;+ mad.lo.s32 %r10, %r8, %r7, %r9;+ setp.ge.s32 %p1, %r10, %r5;+ @%p1 bra BB0_2;++BB0_1:+ mul.wide.s32 %rd7, %r10, 4;+ add.s64 %rd8, %rd3, %rd7;+ add.s64 %rd9, %rd2, %rd7;+ mul.wide.s32 %rd10, %r10, 8;+ add.s64 %rd11, %rd1, %rd10;+ ld.global.f32 %f1, [%rd9];+ ld.global.f32 %f2, [%rd8];+ st.global.v2.f32 [%rd11], {%f2, %f1};+ add.s32 %r10, %r10, %r1;+ setp.lt.s32 %p2, %r10, %r5;+ @%p2 bra BB0_1;++BB0_2:+ ret;+}++ // .globl deinterleave+.visible .entry deinterleave(+ .param .u64 deinterleave_param_0,+ .param .u64 deinterleave_param_1,+ .param .u64 deinterleave_param_2,+ .param .u32 deinterleave_param_3+)+{+ .reg .pred %p<3>;+ .reg .f32 %f<5>;+ .reg .b32 %r<11>;+ .reg .b64 %rd<12>;+++ ld.param.u64 %rd4, [deinterleave_param_0];+ ld.param.u64 %rd5, [deinterleave_param_1];+ ld.param.u64 %rd6, [deinterleave_param_2];+ ld.param.u32 %r5, [deinterleave_param_3];+ cvta.to.global.u64 %rd1, %rd5;+ cvta.to.global.u64 %rd2, %rd4;+ cvta.to.global.u64 %rd3, %rd6;+ mov.u32 %r6, %nctaid.x;+ mov.u32 %r7, %ntid.x;+ mul.lo.s32 %r1, %r6, %r7;+ mov.u32 %r8, %ctaid.x;+ mov.u32 %r9, %tid.x;+ mad.lo.s32 %r10, %r8, %r7, %r9;+ setp.ge.s32 %p1, %r10, %r5;+ @%p1 bra BB1_2;++BB1_1:+ mul.wide.s32 %rd7, %r10, 8;+ add.s64 %rd8, %rd3, %rd7;+ ld.global.v2.f32 {%f1, %f2}, [%rd8];+ mul.wide.s32 %rd9, %r10, 4;+ add.s64 %rd10, %rd2, %rd9;+ st.global.f32 [%rd10], %f1;+ add.s64 %rd11, %rd1, %rd9;+ st.global.f32 [%rd11], %f2;+ add.s32 %r10, %r10, %r1;+ setp.lt.s32 %p2, %r10, %r5;+ @%p2 bra BB1_1;++BB1_2:+ ret;+}++
+ cubits/twine_f64.cu view
@@ -0,0 +1,63 @@+/*+ * Module : Twine+ * Copyright : [2016] Trevor L. McDonell+ * License : BSD3+ *+ * Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+ * Stability : experimental+ * Portability : non-portable (GHC extensions)+ *+ * Convert between Accelerate's Struct-of-Array representation of complex+ * numbers and the Array-of-Struct representation necessary for CUBLAS.+ *+ */++#include <cuda.h>+#include <cuComplex.h>++#ifdef __cplusplus+extern "C" {+#endif++__global__ void interleave+(+ cuDoubleComplex * __restrict__ cplx,+ const double * __restrict__ real,+ const double * __restrict__ imag,+ const int size+)+{+ const int gridSize = blockDim.x * gridDim.x;+ int ix;++ for (ix = blockDim.x * blockIdx.x + threadIdx.x; ix < size; ix += gridSize) {+ const double re = real[ix];+ const double im = imag[ix];++ cplx[ix] = make_cuDoubleComplex(re, im);+ }+}++__global__ void deinterleave+(+ double * __restrict__ real,+ double * __restrict__ imag,+ const cuDoubleComplex * __restrict__ cplx,+ const int size+)+{+ const int gridSize = blockDim.x * gridDim.x;+ int ix;++ for (ix = blockDim.x * blockIdx.x + threadIdx.x; ix < size; ix += gridSize) {+ const cuDoubleComplex c = cplx[ix];++ real[ix] = cuCreal(c);+ imag[ix] = cuCimag(c);+ }+}++#ifdef __cplusplus+}+#endif+
+ cubits/twine_f64.ptx view
@@ -0,0 +1,108 @@+//+// 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;+}++