diff --git a/Data/Array/Accelerate/Math/DFT.hs b/Data/Array/Accelerate/Math/DFT.hs
--- a/Data/Array/Accelerate/Math/DFT.hs
+++ b/Data/Array/Accelerate/Math/DFT.hs
@@ -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'
 
diff --git a/Data/Array/Accelerate/Math/DFT/Centre.hs b/Data/Array/Accelerate/Math/DFT/Centre.hs
--- a/Data/Array/Accelerate/Math/DFT/Centre.hs
+++ b/Data/Array/Accelerate/Math/DFT/Centre.hs
@@ -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)
 
diff --git a/Data/Array/Accelerate/Math/DFT/Roots.hs b/Data/Array/Accelerate/Math/DFT/Roots.hs
--- a/Data/Array/Accelerate/Math/DFT/Roots.hs
+++ b/Data/Array/Accelerate/Math/DFT/Roots.hs
@@ -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 =
diff --git a/Data/Array/Accelerate/Math/FFT.hs b/Data/Array/Accelerate/Math/FFT.hs
--- a/Data/Array/Accelerate/Math/FFT.hs
+++ b/Data/Array/Accelerate/Math/FFT.hs
@@ -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
diff --git a/Data/Array/Accelerate/Math/FFT/CUDA.hs b/Data/Array/Accelerate/Math/FFT/CUDA.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Accelerate/Math/FFT/CUDA.hs
@@ -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
+
diff --git a/Data/Array/Accelerate/Math/FFT/LLVM/Native.hs b/Data/Array/Accelerate/Math/FFT/LLVM/Native.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Accelerate/Math/FFT/LLVM/Native.hs
@@ -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"
+
diff --git a/Data/Array/Accelerate/Math/FFT/LLVM/PTX.hs b/Data/Array/Accelerate/Math/FFT/LLVM/PTX.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Accelerate/Math/FFT/LLVM/PTX.hs
@@ -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
+
diff --git a/Data/Array/Accelerate/Math/FFT/Mode.hs b/Data/Array/Accelerate/Math/FFT/Mode.hs
new file mode 100644
--- /dev/null
+++ b/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 <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
+
diff --git a/Data/Array/Accelerate/Math/FFT/Twine.hs b/Data/Array/Accelerate/Math/FFT/Twine.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Accelerate/Math/FFT/Twine.hs
@@ -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
+
diff --git a/accelerate-fft.cabal b/accelerate-fft.cabal
--- a/accelerate-fft.cabal
+++ b/accelerate-fft.cabal
@@ -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
 
diff --git a/cubits/twine_f32.cu b/cubits/twine_f32.cu
new file mode 100644
--- /dev/null
+++ b/cubits/twine_f32.cu
@@ -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
+
diff --git a/cubits/twine_f32.ptx b/cubits/twine_f32.ptx
new file mode 100644
--- /dev/null
+++ b/cubits/twine_f32.ptx
@@ -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;
+}
+
+
diff --git a/cubits/twine_f64.cu b/cubits/twine_f64.cu
new file mode 100644
--- /dev/null
+++ b/cubits/twine_f64.cu
@@ -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
+
diff --git a/cubits/twine_f64.ptx b/cubits/twine_f64.ptx
new file mode 100644
--- /dev/null
+++ b/cubits/twine_f64.ptx
@@ -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;
+}
+
+