diff --git a/Data/Array/Accelerate/CUDA.hs b/Data/Array/Accelerate/CUDA.hs
--- a/Data/Array/Accelerate/CUDA.hs
+++ b/Data/Array/Accelerate/CUDA.hs
@@ -136,12 +136,50 @@
 --
 --  * @-ddump-cc:@ Print information about the CUDA kernels as they are compiled
 --    and run. Using this option will indicate whether your program is
---    generating the number of kernels that you were expecting.
+--    generating the number of kernels that you were expecting. Note that
+--    compiled kernels are cached in your home directory, and the generated code
+--    will only be displayed if it was not located in this persistent cache. To
+--    clear the cache and always print the generated code, use @-fflush-cache@
+--    as well.
 --
 --  * @-ddump-exec:@ Print each kernel as it is being executed, with timing
 --    information.
 --
 -- See the @accelerate-cuda.cabal@ file for the full list of options.
+--
+--
+-- [/Automatic Graphics Switching on Mac OS X:/]
+--
+-- Some Apple computers contain two graphics processors: a low-power integrated
+-- graphics chipset, as well as a higher-performance NVIDIA GPU. The latter is
+-- of course the one we want to use. Usually Mac OS X detects whenever a program
+-- attempts to run a CUDA function and switches to the NVIDIA GPU automatically.
+--
+-- However, sometimes this does not work correctly and the problem can manifest
+-- in several ways:
+--
+--  * The program may report an error such as \"No CUDA-capable device is
+--    available\" or \"invalid context handle\".
+--
+--  * For programs that also use OpenGL, the graphics switching might occur and
+--    the Accelerate computation complete as expected, but no OpenGL updates
+--    appear on screen.
+--
+-- There are several solutions:
+--
+--  * Use a tool such as /gfxCardStatus/ to manually select either the
+--    integrated or discrete GPU: <http://gfx.io>
+--
+--  * Disable automatic graphics switching in the Energy Saver pane of System
+--    Preferences. Since this disables use of the low-power integrated GPU, this
+--    can decrease battery life.
+--
+--  * When executing the program, disable the RTS clock by appending @+RTS -V0@
+--    to the command line arguments. This disables the RTS clock and all timers
+--    that depend on it: the context switch timer and the heap profiling timer.
+--    Context switches still happen, but deterministically and at a rate much
+--    faster than normal. Automatic graphics switching will work correctly, but
+--    this method has the disadvantage of reducing performance of the program.
 --
 
 module Data.Array.Accelerate.CUDA (
diff --git a/Data/Array/Accelerate/CUDA/Analysis/Launch.hs b/Data/Array/Accelerate/CUDA/Analysis/Launch.hs
--- a/Data/Array/Accelerate/CUDA/Analysis/Launch.hs
+++ b/Data/Array/Accelerate/CUDA/Analysis/Launch.hs
@@ -137,37 +137,38 @@
 --
 sharedMem :: CUDA.DeviceProperties -> PreOpenAcc DelayedOpenAcc aenv a -> Int -> Int
 -- non-computation forms
-sharedMem _ (Alet _ _)      _   = INTERNAL_ERROR(error) "sharedMem" "Let"
-sharedMem _ (Avar _)        _   = INTERNAL_ERROR(error) "sharedMem" "Avar"
-sharedMem _ (Apply _ _)     _   = INTERNAL_ERROR(error) "sharedMem" "Apply"
-sharedMem _ (Acond _ _ _)   _   = INTERNAL_ERROR(error) "sharedMem" "Acond"
-sharedMem _ (Atuple _)      _   = INTERNAL_ERROR(error) "sharedMem" "Atuple"
-sharedMem _ (Aprj _ _)      _   = INTERNAL_ERROR(error) "sharedMem" "Aprj"
-sharedMem _ (Use _)         _   = INTERNAL_ERROR(error) "sharedMem" "Use"
-sharedMem _ (Unit _)        _   = INTERNAL_ERROR(error) "sharedMem" "Unit"
-sharedMem _ (Reshape _ _)   _   = INTERNAL_ERROR(error) "sharedMem" "Reshape"
-sharedMem _ (Aforeign _ _ _) _  = INTERNAL_ERROR(error) "sharedMem" "Aforeign"
+sharedMem _ Alet{}     _ = INTERNAL_ERROR(error) "sharedMem" "Let"
+sharedMem _ Avar{}     _ = INTERNAL_ERROR(error) "sharedMem" "Avar"
+sharedMem _ Apply{}    _ = INTERNAL_ERROR(error) "sharedMem" "Apply"
+sharedMem _ Acond{}    _ = INTERNAL_ERROR(error) "sharedMem" "Acond"
+sharedMem _ Awhile{}   _ = INTERNAL_ERROR(error) "sharedMem" "Awhile"
+sharedMem _ Atuple{}   _ = INTERNAL_ERROR(error) "sharedMem" "Atuple"
+sharedMem _ Aprj{}     _ = INTERNAL_ERROR(error) "sharedMem" "Aprj"
+sharedMem _ Use{}      _ = INTERNAL_ERROR(error) "sharedMem" "Use"
+sharedMem _ Unit{}     _ = INTERNAL_ERROR(error) "sharedMem" "Unit"
+sharedMem _ Reshape{}  _ = INTERNAL_ERROR(error) "sharedMem" "Reshape"
+sharedMem _ Aforeign{} _ = INTERNAL_ERROR(error) "sharedMem" "Aforeign"
 
 -- skeleton nodes
-sharedMem _ (Generate _ _)       _        = 0
-sharedMem _ (Transform _ _ _ _)  _        = 0
-sharedMem _ (Replicate _ _ _)    _        = 0
-sharedMem _ (Slice _ _ _)        _        = 0
-sharedMem _ (Map _ _)            _        = 0
-sharedMem _ (ZipWith _ _ _)      _        = 0
-sharedMem _ (Permute _ _ _ _)    _        = 0
-sharedMem _ (Backpermute _ _ _)  _        = 0
-sharedMem _ (Stencil _ _ _)      _        = 0
-sharedMem _ (Stencil2 _ _ _ _ _) _        = 0
-sharedMem _ (Fold  _ x _)        blockDim = sizeOf (delayedExpType x) * blockDim
-sharedMem _ (Scanl _ x _)        blockDim = sizeOf (delayedExpType x) * blockDim
-sharedMem _ (Scanr _ x _)        blockDim = sizeOf (delayedExpType x) * blockDim
-sharedMem _ (Scanl' _ x _)       blockDim = sizeOf (delayedExpType x) * blockDim
-sharedMem _ (Scanr' _ x _)       blockDim = sizeOf (delayedExpType x) * blockDim
-sharedMem _ (Fold1 _ a)          blockDim = sizeOf (delayedAccType a) * blockDim
-sharedMem _ (Scanl1 _ a)         blockDim = sizeOf (delayedAccType a) * blockDim
-sharedMem _ (Scanr1 _ a)         blockDim = sizeOf (delayedAccType a) * blockDim
-sharedMem p (FoldSeg _ x _ _)    blockDim =
+sharedMem _ Generate{}          _        = 0
+sharedMem _ Transform{}         _        = 0
+sharedMem _ Replicate{}         _        = 0
+sharedMem _ Slice{}             _        = 0
+sharedMem _ Map{}               _        = 0
+sharedMem _ ZipWith{}           _        = 0
+sharedMem _ Permute{}           _        = 0
+sharedMem _ Backpermute{}       _        = 0
+sharedMem _ Stencil{}           _        = 0
+sharedMem _ Stencil2{}          _        = 0
+sharedMem _ (Fold  _ x _)       blockDim = sizeOf (delayedExpType x) * blockDim
+sharedMem _ (Scanl _ x _)       blockDim = sizeOf (delayedExpType x) * blockDim
+sharedMem _ (Scanr _ x _)       blockDim = sizeOf (delayedExpType x) * blockDim
+sharedMem _ (Scanl' _ x _)      blockDim = sizeOf (delayedExpType x) * blockDim
+sharedMem _ (Scanr' _ x _)      blockDim = sizeOf (delayedExpType x) * blockDim
+sharedMem _ (Fold1 _ a)         blockDim = sizeOf (delayedAccType a) * blockDim
+sharedMem _ (Scanl1 _ a)        blockDim = sizeOf (delayedAccType a) * blockDim
+sharedMem _ (Scanr1 _ a)        blockDim = sizeOf (delayedAccType a) * blockDim
+sharedMem p (FoldSeg _ x _ _)   blockDim =
   (blockDim `div` CUDA.warpSize p) * 8 + blockDim * sizeOf (delayedExpType x)  -- TLM: why 8? I can't remember...
 sharedMem p (Fold1Seg _ a _) blockDim =
   (blockDim `div` CUDA.warpSize p) * 8 + blockDim * sizeOf (delayedAccType a)
diff --git a/Data/Array/Accelerate/CUDA/Array/Data.hs b/Data/Array/Accelerate/CUDA/Array/Data.hs
--- a/Data/Array/Accelerate/CUDA/Array/Data.hs
+++ b/Data/Array/Accelerate/CUDA/Array/Data.hs
@@ -7,6 +7,7 @@
 -- Module      : Data.Array.Accelerate.CUDA.Array.Data
 -- Copyright   : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee
 --               [2009..2013] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell
+--               [2013]       Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell, Robert Clifton-Everest
 -- License     : BSD3
 --
 -- Maintainer  : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>
@@ -19,11 +20,13 @@
   -- * Array operations and representations
   mallocArray, indexArray,
   useArray,  useArrayAsync,
+  useDevicePtrs,
   copyArray, copyArrayPeer, copyArrayPeerAsync,
   peekArray, peekArrayAsync,
   pokeArray, pokeArrayAsync,
   marshalArrayData, marshalTextureData, marshalDevicePtrs,
   devicePtrsOfArrayData, advancePtrsOfArrayData,
+  devicePtrsFromList, devicePtrsToWordPtrs,
 
   -- * Garbage collection
   cleanupArrayData
@@ -32,15 +35,17 @@
 
 -- libraries
 import Prelude                                          hiding ( fst, snd )
+import qualified Prelude                                as P
 import Control.Applicative
 import Control.Monad.Reader                             ( asks )
 import Control.Monad.State                              ( gets )
 import Control.Monad.Trans                              ( liftIO )
 import Foreign.C.Types
+import Foreign.Ptr
 
 -- friends
 import Data.Array.Accelerate.Array.Data
-import Data.Array.Accelerate.Array.Sugar                ( Array(..), Shape, Elt, toElt )
+import Data.Array.Accelerate.Array.Sugar                ( Array(..), Shape, Elt, toElt, EltRepr )
 import Data.Array.Accelerate.Array.Representation       ( size )
 import Data.Array.Accelerate.CUDA.State
 import Data.Array.Accelerate.CUDA.Array.Table
@@ -159,6 +164,37 @@
         mkPrimDispatch(usePrim,Prim.useArrayAsync)
 
 
+useDevicePtrs :: (Shape sh, Elt e) => EltRepr sh -> Prim.DevicePtrs (EltRepr e) -> CIO (Array sh e)
+useDevicePtrs sh ptrs = run doUse
+  where
+    !n             = size sh
+    doUse !ctx !mt = Array sh <$> useD arrayElt ptrs
+      where
+        useD :: ArrayEltR e -> Prim.DevicePtrs e -> IO (ArrayData e)
+        useD ArrayEltRunit             _  = return AD_Unit
+        useD (ArrayEltRpair aeR1 aeR2) ps = AD_Pair <$> useD aeR1 (P.fst ps)
+                                                    <*> useD aeR2 (P.snd ps)
+        useD aer                       ps = usePrim aer ctx mt ps n
+        --
+        usePrim :: ArrayEltR e -> Context -> MemoryTable -> Prim.DevicePtrs e -> Int -> IO (ArrayData e)
+        mkPrimDispatch(usePrim,Prim.useDevicePtrs)
+
+devicePtrsFromList :: ArrayEltR e -> [WordPtr] -> Prim.DevicePtrs e
+devicePtrsFromList aeR = P.fst . (devP aeR)
+  where
+    devP :: ArrayEltR e -> [WordPtr] -> (Prim.DevicePtrs e, [WordPtr])
+    devP ArrayEltRunit             ps     = ((),ps)
+    devP (ArrayEltRpair aeR1 aeR2) ps     = let
+        (d1, ps')  = devP aeR1 ps
+        (d2, ps'') = devP aeR2 ps'
+      in ((d1,d2), ps'')
+    devP aer                       (p:ps) = (devPrim aer p, ps)
+    devP _                         []     = error "devicePtrsFromList: incorrect number of device pointers for element type"
+    --
+    devPrim :: ArrayEltR e -> WordPtr -> Prim.DevicePtrs e
+    mkPrimDispatch(devPrim,CUDA.wordPtrToDevPtr)
+
+
 -- |Read a single element from an array at the given row-major index. This is a
 -- synchronous operation.
 --
@@ -326,21 +362,25 @@
         mkPrimDispatch(pokePrim,Prim.pokeArrayAsync)
 
 
--- |Wrap device pointers into arguments that can be passed to a kernel
--- invocation
+-- |Convert the device pointers into a list of word pointers
 --
-marshalDevicePtrs :: ArrayElt e => ArrayData e -> Prim.DevicePtrs e -> [CUDA.FunParam]
-marshalDevicePtrs !adata = marshalR arrayElt adata
+devicePtrsToWordPtrs :: ArrayElt e => ArrayData e -> Prim.DevicePtrs e -> [WordPtr]
+devicePtrsToWordPtrs !adata = wordR arrayElt adata
   where
-    marshalR :: ArrayEltR e -> ArrayData e -> Prim.DevicePtrs e -> [CUDA.FunParam]
-    marshalR ArrayEltRunit             _  _       = []
-    marshalR (ArrayEltRpair aeR1 aeR2) ad (p1,p2) = marshalR aeR1 (fst ad) p1 ++
-                                                    marshalR aeR2 (snd ad) p2
-    marshalR aer                       ad ptr     = [marshalPrim aer ad ptr]
+    wordR :: ArrayEltR e -> ArrayData e -> Prim.DevicePtrs e -> [WordPtr]
+    wordR ArrayEltRunit             _  _       = []
+    wordR (ArrayEltRpair aeR1 aeR2) ad (p1,p2) = wordR aeR1 (fst ad) p1 ++
+                                                 wordR aeR2 (snd ad) p2
+    wordR aer                       ad ptr     = [wordPrim aer ad ptr]
     --
-    marshalPrim :: ArrayEltR e -> ArrayData e -> Prim.DevicePtrs e -> CUDA.FunParam
-    mkPrimDispatch(marshalPrim,Prim.marshalDevicePtrs)
+    wordPrim :: ArrayEltR e -> ArrayData e -> Prim.DevicePtrs e -> WordPtr
+    mkPrimDispatch(wordPrim,const CUDA.devPtrToWordPtr)
 
+-- |Wrap device pointers into arguments that can be passed to a kernel
+-- invocation
+--
+marshalDevicePtrs :: ArrayElt e => ArrayData e -> Prim.DevicePtrs e -> [CUDA.FunParam]
+marshalDevicePtrs !adata ptrs = map (CUDA.VArg . CUDA.wordPtrToDevPtr) $ devicePtrsToWordPtrs adata ptrs
 
 -- |Wrap the device pointers corresponding to a host-side array into arguments
 -- that can be passed to a kernel upon invocation.
@@ -413,4 +453,3 @@
     --
     advancePrim :: ArrayEltR e -> ArrayData e -> Prim.DevicePtrs e -> Prim.DevicePtrs e
     mkPrimDispatch(advancePrim,Prim.advancePtrsOfArrayData n)
-
diff --git a/Data/Array/Accelerate/CUDA/Array/Nursery.hs b/Data/Array/Accelerate/CUDA/Array/Nursery.hs
--- a/Data/Array/Accelerate/CUDA/Array/Nursery.hs
+++ b/Data/Array/Accelerate/CUDA/Array/Nursery.hs
@@ -14,7 +14,7 @@
 
 module Data.Array.Accelerate.CUDA.Array.Nursery (
 
-  Nursery(..), NRS, new, lookup, insert, flush,
+  Nursery(..), NRS, new, malloc, stash, flush,
 
 ) where
 
@@ -24,10 +24,10 @@
 import qualified Data.Array.Accelerate.CUDA.Debug               as D
 
 -- libraries
-import Prelude                                                  hiding ( lookup )
-import Data.IORef
+import Prelude
 import Data.Hashable
 import Control.Exception                                        ( bracket_ )
+import Control.Concurrent.MVar                                  ( MVar, newMVar, withMVar, mkWeakMVar )
 import System.Mem.Weak                                          ( Weak )
 import Foreign.Ptr                                              ( ptrToIntPtr )
 import Foreign.CUDA.Ptr                                         ( DevicePtr )
@@ -49,13 +49,14 @@
 --
 type HashTable key val  = HT.BasicHashTable key val
 
-type NRS                = IORef ( HashTable (CUDA.Context, Int) (FullList () (DevicePtr ())) )
+type NRS                = MVar ( HashTable (CUDA.Context, Int) (FullList () (DevicePtr ())) )
 data Nursery            = Nursery {-# UNPACK #-} !NRS
                                   {-# UNPACK #-} !(Weak NRS)
 
 instance Hashable CUDA.Context where
   {-# INLINE hashWithSalt #-}
-  hashWithSalt s (CUDA.Context ctx) = hashWithSalt s (fromIntegral (ptrToIntPtr ctx) :: Int)
+  hashWithSalt salt (CUDA.Context ctx)
+    = salt `hashWithSalt` (fromIntegral (ptrToIntPtr ctx) :: Int)
 
 
 -- Generate a fresh nursery
@@ -63,8 +64,8 @@
 new :: IO Nursery
 new = do
   tbl    <- HT.new
-  ref    <- newIORef tbl
-  weak   <- mkWeakIORef ref (flush tbl)
+  ref    <- newMVar tbl
+  weak   <- mkWeakMVar ref (flush tbl)
   return $! Nursery ref weak
 
 
@@ -72,12 +73,11 @@
 -- larger). If found, it is removed from the nursery and a pointer to it
 -- returned.
 --
-{-# INLINE lookup #-}
-lookup :: Int -> CUDA.Context -> Nursery -> IO (Maybe (DevicePtr ()))
-lookup !n !ctx (Nursery !ref _) = do
+{-# INLINE malloc #-}
+malloc :: Int -> CUDA.Context -> Nursery -> IO (Maybe (DevicePtr ()))
+malloc !n !ctx (Nursery !ref _) = withMVar ref $ \tbl -> do
   let !key = (ctx,n)
   --
-  tbl <- readIORef ref
   mp  <- HT.lookup tbl key
   case mp of
     Nothing               -> return Nothing
@@ -89,12 +89,11 @@
 
 -- Add a device pointer to the nursery.
 --
-{-# INLINE insert #-}
-insert :: Int -> CUDA.Context -> NRS -> DevicePtr a -> IO ()
-insert !n !ctx !ref (CUDA.castDevPtr -> !ptr) = do
+{-# INLINE stash #-}
+stash :: Int -> CUDA.Context -> NRS -> DevicePtr a -> IO ()
+stash !n !ctx !ref (CUDA.castDevPtr -> !ptr) = withMVar ref $ \tbl -> do
   let !key = (ctx, n)
   --
-  tbl <- readIORef ref
   mp  <- HT.lookup tbl key
   case mp of
     Nothing     -> HT.insert tbl key (FL.singleton () ptr)
diff --git a/Data/Array/Accelerate/CUDA/Array/Prim.hs b/Data/Array/Accelerate/CUDA/Array/Prim.hs
--- a/Data/Array/Accelerate/CUDA/Array/Prim.hs
+++ b/Data/Array/Accelerate/CUDA/Array/Prim.hs
@@ -3,6 +3,7 @@
 {-# LANGUAGE GADTs               #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE TypeFamilies        #-}
+{-# LANGUAGE TupleSections       #-}
 -- |
 -- Module      : Data.Array.Accelerate.CUDA.Array.Prim
 -- Copyright   : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee
@@ -20,6 +21,7 @@
 
   mallocArray, indexArray,
   useArray,  useArrayAsync,
+  useDevicePtrs,
   copyArray, copyArrayPeer, copyArrayPeerAsync,
   peekArray, peekArrayAsync,
   pokeArray, pokeArrayAsync,
@@ -29,11 +31,7 @@
 ) where
 
 -- libraries
-#if MIN_VERSION_base(4,6,0)
 import Prelude                                          hiding ( lookup )
-#else
-import Prelude                                          hiding ( catch, lookup )
-#endif
 import Data.Int
 import Data.Word
 import Data.Maybe
@@ -189,9 +187,8 @@
   in do
     exists <- isJust <$> (lookup ctx mt ad :: IO (Maybe (CUDA.DevicePtr a)))
     unless exists $ do
-      message $ "useArray/malloc: " ++ showBytes bytes
       dst <- malloc ctx mt ad n
-      CUDA.pokeArray n src dst
+      transfer "useArray/malloc" bytes $ CUDA.pokeArray n src dst
 
 
 useArrayAsync
@@ -209,11 +206,27 @@
   in do
     exists <- isJust <$> (lookup ctx mt ad :: IO (Maybe (CUDA.DevicePtr a)))
     unless exists $ do
-      message $ "useArrayAsync/malloc: " ++ showBytes bytes
       dst <- malloc ctx mt ad n
-      CUDA.pokeArrayAsync n src dst ms
+      transfer "useArrayAsync/malloc" bytes $ CUDA.pokeArrayAsync n src dst ms
 
 
+useDevicePtrs
+    :: forall e a. (ArrayElt e, ArrayPtrs e ~ Ptr a, DevicePtrs e ~ CUDA.DevicePtr a, Typeable e, Typeable a, Storable a)
+    => Context
+    -> MemoryTable
+    -> DevicePtrs e
+    -> Int
+    -> IO (ArrayData e)
+useDevicePtrs !ctx !mt !ptr !n0 =
+  let !n         = 1 `max` n0
+      !bytes     = n * sizeOf (undefined :: a)
+      (adata, _) = runArrayData $ (,undefined) `fmap` newArrayData n
+  in do
+    message $ "useDevicePtrs: " ++ showBytes bytes
+    insertRemote ctx mt adata ptr
+    return adata
+
+
 -- Read a single element from an array at the given row-major index
 --
 indexArray
@@ -243,10 +256,10 @@
     -> Int                      -- number of array elements
     -> IO ()
 copyArray !ctx !mt !from !to !n = do
-  message $ "copyArrayAsync: " ++ showBytes (n * sizeOf (undefined :: b))
   src <- devicePtrsOfArrayData ctx mt from
   dst <- devicePtrsOfArrayData ctx mt to
-  CUDA.copyArrayAsync n src dst
+  transfer "copyArrayAsync" (n * sizeOf (undefined :: b)) $
+    CUDA.copyArrayAsync n src dst
 
 
 -- Copy data between two device arrays that exist in different contexts and/or
@@ -260,10 +273,10 @@
     -> Int                      -- number of array elements
     -> IO ()
 copyArrayPeer !mt !from !ctxSrc !to !ctxDst !n = do
-  message $ "copyArrayPeer: " ++ showBytes (n * sizeOf (undefined :: b))
   src <- devicePtrsOfArrayData ctxSrc mt from
   dst <- devicePtrsOfArrayData ctxDst mt to
-  CUDA.copyArrayPeer n src (deviceContext ctxSrc) dst (deviceContext ctxDst)
+  transfer "copyArrayPeer" (n * sizeOf (undefined :: b)) $
+    CUDA.copyArrayPeer n src (deviceContext ctxSrc) dst (deviceContext ctxDst)
 
 copyArrayPeerAsync
     :: forall e a b. (ArrayElt e, ArrayPtrs e ~ Ptr a, DevicePtrs e ~ CUDA.DevicePtr b, Typeable a, Typeable b, Typeable e, Storable b)
@@ -274,10 +287,10 @@
     -> Maybe CUDA.Stream
     -> IO ()
 copyArrayPeerAsync !mt !from !ctxSrc !to !ctxDst !n !st = do
-  message $ "copyArrayPeerAsync: " ++ showBytes (n * sizeOf (undefined :: b))
   src <- devicePtrsOfArrayData ctxSrc mt from
   dst <- devicePtrsOfArrayData ctxDst mt to
-  CUDA.copyArrayPeerAsync n src (deviceContext ctxSrc) dst (deviceContext ctxDst) st
+  transfer "copyArrayPeerAsync" (n * sizeOf (undefined :: b)) $
+    CUDA.copyArrayPeerAsync n src (deviceContext ctxSrc) dst (deviceContext ctxDst) st
 
 
 -- Copy data from the device into the associated Accelerate host-side array
@@ -290,9 +303,9 @@
     -> Int
     -> IO ()
 peekArray !ctx !mt !ad !n =
-  devicePtrsOfArrayData ctx mt ad >>= \src -> do
-    message $ "peekArray: " ++ showBytes (n * sizeOf (undefined :: a))
-    CUDA.peekArray n src (ptrsOfArrayData ad)
+  devicePtrsOfArrayData ctx mt ad >>= \src ->
+    transfer "peekArray" (n * sizeOf (undefined :: a)) $
+      CUDA.peekArray n src (ptrsOfArrayData ad)
 
 peekArrayAsync
     :: forall e a. (ArrayElt e, ArrayPtrs e ~ Ptr a, DevicePtrs e ~ CUDA.DevicePtr a, Typeable a, Typeable e, Storable a)
@@ -303,9 +316,9 @@
     -> Maybe CUDA.Stream
     -> IO ()
 peekArrayAsync !ctx !mt !ad !n !st =
-  devicePtrsOfArrayData ctx mt ad >>= \src -> do
-    message $ "peekArrayAsync: " ++ showBytes (n * sizeOf (undefined :: a))
-    CUDA.peekArrayAsync n src (CUDA.HostPtr $ ptrsOfArrayData ad) st
+  devicePtrsOfArrayData ctx mt ad >>= \src ->
+    transfer "peekArrayAsync" (n * sizeOf (undefined :: a)) $
+      CUDA.peekArrayAsync n src (CUDA.HostPtr $ ptrsOfArrayData ad) st
 
 
 -- Copy data from an Accelerate array into the associated device array
@@ -318,9 +331,9 @@
     -> Int
     -> IO ()
 pokeArray !ctx !mt !ad !n =
-  devicePtrsOfArrayData ctx mt ad >>= \dst -> do
-    message $ "pokeArray: " ++ showBytes (n * sizeOf (undefined :: a))
-    CUDA.pokeArray n (ptrsOfArrayData ad) dst
+  devicePtrsOfArrayData ctx mt ad >>= \dst ->
+    transfer "pokeArray: " (n * sizeOf (undefined :: a)) $
+      CUDA.pokeArray n (ptrsOfArrayData ad) dst
 
 pokeArrayAsync
     :: forall e a. (ArrayElt e, ArrayPtrs e ~ Ptr a, DevicePtrs e ~ CUDA.DevicePtr a, Typeable a, Typeable e, Storable a)
@@ -331,9 +344,9 @@
     -> Maybe CUDA.Stream
     -> IO ()
 pokeArrayAsync !ctx !mt !ad !n !st =
-  devicePtrsOfArrayData ctx mt ad >>= \dst -> do
-    message $ "pokeArrayAsync: " ++ showBytes (n * sizeOf (undefined :: a))
-    CUDA.pokeArrayAsync n (CUDA.HostPtr $ ptrsOfArrayData ad) dst st
+  devicePtrsOfArrayData ctx mt ad >>= \dst ->
+    transfer "pokeArrayAsync: " (n * sizeOf (undefined :: a)) $
+      CUDA.pokeArrayAsync n (CUDA.HostPtr $ ptrsOfArrayData ad) dst st
 
 
 -- Marshal device pointers to arguments that can be passed to kernel invocation
@@ -417,4 +430,14 @@
 {-# INLINE message #-}
 message :: String -> IO ()
 message s = s `trace` return ()
+
+{-# INLINE transfer #-}
+transfer :: String -> Int -> IO () -> IO ()
+transfer name bytes action
+  = let showRate x t        = D.showFFloatSIBase (Just 3) 1024 (fromIntegral x / t) "B/s"
+        msg gpuTime cpuTime = "gc: " ++ name ++ ": "
+                                     ++ showBytes bytes ++ " @ " ++ showRate bytes gpuTime ++ ", "
+                                     ++ D.elapsed gpuTime cpuTime
+    in
+    D.timed D.dump_gc msg action
 
diff --git a/Data/Array/Accelerate/CUDA/Array/Sugar.hs b/Data/Array/Accelerate/CUDA/Array/Sugar.hs
--- a/Data/Array/Accelerate/CUDA/Array/Sugar.hs
+++ b/Data/Array/Accelerate/CUDA/Array/Sugar.hs
@@ -12,7 +12,7 @@
 module Data.Array.Accelerate.CUDA.Array.Sugar (
 
   module Data.Array.Accelerate.Array.Sugar,
-  newArray, allocateArray, useArray
+  newArray, allocateArray, useArray, useArrayAsync,
 
 ) where
 
diff --git a/Data/Array/Accelerate/CUDA/Array/Table.hs b/Data/Array/Accelerate/CUDA/Array/Table.hs
--- a/Data/Array/Accelerate/CUDA/Array/Table.hs
+++ b/Data/Array/Accelerate/CUDA/Array/Table.hs
@@ -17,21 +17,17 @@
 module Data.Array.Accelerate.CUDA.Array.Table (
 
   -- Tables for host/device memory associations
-  MemoryTable, new, lookup, malloc, insert, reclaim
+  MemoryTable, new, lookup, malloc, insert, insertRemote, reclaim
 
 ) where
 
-#if !MIN_VERSION_base(4,6,0)
-import Prelude                                                  hiding ( lookup, catch )
-#else
 import Prelude                                                  hiding ( lookup )
-#endif
-import Data.IORef                                               ( IORef, newIORef, readIORef, mkWeakIORef )
 import Data.Maybe                                               ( isJust )
 import Data.Hashable                                            ( Hashable(..) )
 import Data.Typeable                                            ( Typeable, gcast )
 import Control.Monad                                            ( unless )
 import Control.Concurrent                                       ( yield )
+import Control.Concurrent.MVar                                  ( MVar, newMVar, withMVar, mkWeakMVar )
 import Control.Exception                                        ( bracket_, catch, throwIO )
 import Control.Applicative                                      ( (<$>) )
 import System.Mem                                               ( performGC )
@@ -68,16 +64,18 @@
 -- semantics of weak pointers in the documentation).
 --
 type HashTable key val  = HT.BasicHashTable key val
-type MT                 = IORef ( HashTable HostArray DeviceArray )
+type MT                 = MVar ( HashTable HostArray DeviceArray )
 data MemoryTable        = MemoryTable {-# UNPACK #-} !MT
                                       {-# UNPACK #-} !(Weak MT)
                                       {-# UNPACK #-} !Nursery
 
 -- Arrays on the host and device
 --
+type ContextId = Int
+
 data HostArray where
   HostArray :: Typeable e
-            => {-# UNPACK #-} !Int      -- unique ID relating to the parent context
+            => {-# UNPACK #-} !ContextId        -- unique ID relating to the parent context
             -> {-# UNPACK #-} !(StableName (ArrayData e))
             -> HostArray
 
@@ -91,6 +89,7 @@
     = maybe False (== a2) (gcast a1)
 
 instance Hashable HostArray where
+  {-# INLINE hashWithSalt #-}
   hashWithSalt salt (HostArray cid sn)
     = salt `hashWithSalt` cid `hashWithSalt` sn
 
@@ -108,9 +107,9 @@
 new = do
   message "initialise memory table"
   tbl  <- HT.new
-  ref  <- newIORef tbl
+  ref  <- newMVar tbl
   nrs  <- N.new
-  weak <- mkWeakIORef ref (table_finalizer tbl)
+  weak <- mkWeakMVar ref (table_finalizer tbl)
   return $! MemoryTable ref weak nrs
 
 
@@ -119,7 +118,7 @@
 lookup :: (Typeable a, Typeable b) => Context -> MemoryTable -> ArrayData a -> IO (Maybe (DevicePtr b))
 lookup ctx (MemoryTable !ref _ _) !arr = do
   sa <- makeStableArray ctx arr
-  mw <- withIORef ref (`HT.lookup` sa)
+  mw <- withMVar ref (`HT.lookup` sa)
   case mw of
     Nothing              -> trace ("lookup/not found: " ++ show sa) $ return Nothing
     Just (DeviceArray w) -> do
@@ -127,6 +126,19 @@
       case mv of
         Just v | Just p <- gcast v -> trace ("lookup/found: " ++ show sa) $ return (Just p)
                | otherwise         -> INTERNAL_ERROR(error) "lookup" $ "type mismatch"
+
+        -- Note: [Weak pointer weirdness]
+        --
+        -- After the lookup is successful, there might conceivably be no further
+        -- references to 'arr'. If that is so, and a garbage collection
+        -- intervenes, the weak pointer might get tombstoned before 'deRefWeak'
+        -- gets to it. In that case we throw an error (below). However, because
+        -- we have used 'arr' in the continuation, this ensures that 'arr' is
+        -- reachable in the continuation of 'deRefWeak' and thus 'deRefWeak'
+        -- always succeeds. This sort of weirdness, typical of the world of weak
+        -- pointers, is why we can not reuse the stable name 'sa' computed
+        -- above in the error message.
+        --
         Nothing                    ->
           makeStableArray ctx arr >>= \x -> INTERNAL_ERROR(error) "lookup" $ "dead weak pair: " ++ show x
 
@@ -150,7 +162,7 @@
       !n'               = chunk * multiple (fromIntegral n) (fromIntegral chunk)
       !bytes            = n' * sizeOf (undefined :: b)
   --
-  mp  <- N.lookup bytes (deviceContext ctx) nursery
+  mp  <- N.malloc bytes (deviceContext ctx) nursery
   ptr <- case mp of
            Just p       -> trace "malloc/nursery" $ return (CUDA.castDevPtr p)
            Nothing      -> trace "malloc/new"     $
@@ -169,11 +181,22 @@
 insert !ctx (MemoryTable !ref !weak_ref (Nursery _ !weak_nrs)) !arr !ptr !bytes = do
   key  <- makeStableArray ctx arr
   dev  <- DeviceArray `fmap` mkWeak arr ptr (Just $ finalizer (weakContext ctx) weak_ref weak_nrs key ptr bytes)
-  tbl  <- readIORef ref
-  message $ "insert: " ++ show key
-  HT.insert tbl key dev
+  message      $ "insert: " ++ show key
+  withMVar ref $ \tbl -> HT.insert tbl key dev
 
 
+-- Record an association between a host-side array and a device memory area that was
+-- not allocated by accelerate. The device memory will NOT be freed when the host
+-- array is garbage collected.
+--
+insertRemote :: (Typeable a, Typeable b) => Context -> MemoryTable -> ArrayData a -> DevicePtr b -> IO ()
+insertRemote !ctx (MemoryTable !ref !weak_ref _) !arr !ptr = do
+  key  <- makeStableArray ctx arr
+  dev  <- DeviceArray `fmap` mkWeak arr ptr (Just $ remoteFinalizer weak_ref key)
+  message      $ "insert/remote: " ++ show key
+  withMVar ref $ \tbl -> HT.insert tbl key dev
+
+
 -- Removing entries
 -- ----------------
 
@@ -185,11 +208,11 @@
   (free, total) <- CUDA.getMemInfo
   performGC
   yield
-  withIORef nrs N.flush
+  withMVar nrs N.flush
   mr <- deRefWeak weak_ref
   case mr of
     Nothing  -> return ()
-    Just ref -> withIORef ref $ \tbl ->
+    Just ref -> withMVar ref $ \tbl ->
       flip HT.mapM_ tbl $ \(_,DeviceArray w) -> do
         alive <- isJust `fmap` deRefWeak w
         unless alive $ finalize w
@@ -214,7 +237,7 @@
   mr <- deRefWeak weak_ref
   case mr of
     Nothing  -> message ("finalise/dead table: " ++ show key)
-    Just ref -> withIORef ref (`HT.delete` key)
+    Just ref -> withMVar ref (`HT.delete` key)
   --
   mc <- deRefWeak weak_ctx
   case mc of
@@ -224,8 +247,14 @@
       mn <- deRefWeak weak_nrs
       case mn of
         Nothing  -> trace ("finalise/free: "     ++ show key) $ bracket_ (CUDA.push ctx) CUDA.pop (CUDA.free ptr)
-        Just nrs -> trace ("finalise/nursery: "  ++ show key) $ N.insert bytes ctx nrs ptr
+        Just nrs -> trace ("finalise/nursery: "  ++ show key) $ N.stash bytes ctx nrs ptr
 
+remoteFinalizer :: Weak MT -> HostArray -> IO ()
+remoteFinalizer !weak_ref !key = do
+  mr <- deRefWeak weak_ref
+  case mr of
+    Nothing  -> message ("finalise/dead table: " ++ show key)
+    Just ref -> trace   ("finalise: "            ++ show key) $ withMVar ref (`HT.delete` key)
 
 table_finalizer :: HashTable HostArray DeviceArray -> IO ()
 table_finalizer !tbl
@@ -243,10 +272,6 @@
       !cid              = fromIntegral (ptrToIntPtr p)
   in
   HostArray cid <$> makeStableName arr
-
-{-# INLINE withIORef #-}
-withIORef :: IORef a -> (a -> IO b) -> IO b
-withIORef ref f = readIORef ref >>= f
 
 
 -- Debug
diff --git a/Data/Array/Accelerate/CUDA/Async.hs b/Data/Array/Accelerate/CUDA/Async.hs
--- a/Data/Array/Accelerate/CUDA/Async.hs
+++ b/Data/Array/Accelerate/CUDA/Async.hs
@@ -12,9 +12,6 @@
 module Data.Array.Accelerate.CUDA.Async
   where
 
-#if !MIN_VERSION_base(4,6,0)
-import Prelude                                          hiding ( catch )
-#endif
 import Control.Exception
 import Control.Concurrent
 
diff --git a/Data/Array/Accelerate/CUDA/CodeGen.hs b/Data/Array/Accelerate/CUDA/CodeGen.hs
--- a/Data/Array/Accelerate/CUDA/CodeGen.hs
+++ b/Data/Array/Accelerate/CUDA/CodeGen.hs
@@ -4,6 +4,7 @@
 {-# LANGUAGE QuasiQuotes         #-}
 {-# LANGUAGE RecordWildCards     #-}
 {-# LANGUAGE ScopedTypeVariables #-}
+{-# OPTIONS -fno-warn-name-shadowing #-}
 -- |
 -- Module      : Data.Array.Accelerate.CUDA.CodeGen
 -- Copyright   : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee
@@ -22,8 +23,8 @@
 ) where
 
 -- libraries
-import Prelude                                                  hiding ( id, exp, replicate, iterate )
-import Control.Applicative                                      ( (<$>), (<*>), (<*) )
+import Prelude                                                  hiding ( id, exp, replicate )
+import Control.Applicative                                      ( (<$>), (<*>) )
 import Control.Monad.State.Strict
 import Data.Loc
 import Data.Char
@@ -45,7 +46,7 @@
 import qualified Data.Array.Accelerate.Analysis.Type            as Sugar
 
 import Data.Array.Accelerate.CUDA.AST                           hiding ( Val(..), prj )
-import Data.Array.Accelerate.CUDA.CodeGen.Base                  hiding ( shapeSize )
+import Data.Array.Accelerate.CUDA.CodeGen.Base
 import Data.Array.Accelerate.CUDA.CodeGen.Type
 import Data.Array.Accelerate.CUDA.CodeGen.Monad
 import Data.Array.Accelerate.CUDA.CodeGen.Mapping
@@ -53,7 +54,7 @@
 import Data.Array.Accelerate.CUDA.CodeGen.PrefixSum
 import Data.Array.Accelerate.CUDA.CodeGen.Reduction
 import Data.Array.Accelerate.CUDA.CodeGen.Stencil
-import Data.Array.Accelerate.CUDA.Foreign                       ( canExecuteExp )
+import Data.Array.Accelerate.CUDA.Foreign.Import                ( canExecuteExp )
 
 #include "accelerate.h"
 
@@ -113,20 +114,21 @@
       Stencil2 f b1 a1 b2 a2    -> mkStencil2 dev aenv  <$> travF2 f <*> travB a1 b1 <*> travB a2 b2
 
       -- Non-computation forms -> sadness
-      Alet _ _                  -> unexpectedError
-      Avar _                    -> unexpectedError
-      Apply _ _                 -> unexpectedError
-      Acond _ _ _               -> unexpectedError
-      Atuple _                  -> unexpectedError
-      Aprj _ _                  -> unexpectedError
-      Use _                     -> unexpectedError
-      Unit _                    -> unexpectedError
-      Aforeign _ _ _            -> unexpectedError
-      Reshape _ _               -> unexpectedError
+      Alet{}                    -> unexpectedError
+      Avar{}                    -> unexpectedError
+      Apply{}                   -> unexpectedError
+      Acond{}                   -> unexpectedError
+      Awhile{}                  -> unexpectedError
+      Atuple{}                  -> unexpectedError
+      Aprj{}                    -> unexpectedError
+      Use{}                     -> unexpectedError
+      Unit{}                    -> unexpectedError
+      Aforeign{}                -> unexpectedError
+      Reshape{}                 -> unexpectedError
 
-      Replicate _ _ _           -> fusionError
-      Slice _ _ _               -> fusionError
-      ZipWith _ _ _             -> fusionError
+      Replicate{}               -> fusionError
+      Slice{}                   -> fusionError
+      ZipWith{}                 -> fusionError
 
   where
     codegen :: CUDA [CUTranslSkel aenv a] -> [CUTranslSkel aenv a]
@@ -177,12 +179,16 @@
 -- Generate code for scalar function abstractions.
 --
 -- This is quite awkward: we have an outer monad to generate fresh variable
--- names, but since we know that even if the function in applied many times
--- (for example, collective operations such as 'fold' and 'scan'), the variables
--- will not shadow each other. Thus, we don't need fresh names at _every_
--- invocation site, so we hack this a bit to return a pure closure.
+-- names, but since we know that even if the function in applied many times (for
+-- example, collective operations such as 'fold' and 'scan'), the variables will
+-- not shadow each other. Thus, we don't need fresh names at _every_ invocation
+-- site, so we hack this a bit to return a pure closure.
 --
--- Still, there has got to be a cleaner way to do this...
+-- Note that the implementation of def-use analysis used for dead code
+-- elimination requires that we always generate code for closed functions.
+-- Additionally, we require two passes over the function: once when performing
+-- the analysis, and a second time when instantiating the function in the
+-- skeleton.
 --
 codegenFun1
     :: forall aenv a b. DeviceProperties
@@ -195,19 +201,22 @@
         go :: Rvalue x => [x] -> Gen ([C.BlockItem], [C.Exp])
         go x = do
           code  <- mapM use =<< codegenOpenExp dev aenv f (Empty `Push` map rvalue x)
-          env   <- getEnv
-          return (env, code)
+          env'  <- getEnv
+          return (env', code)
 
+        -- Initial code generation proceeds with dummy variable names. The real
+        -- names are substituted later when we instantiate the skeleton.
         (_,u,_) = locals "undefined_x" (undefined :: a)
     in do
       n                 <- get
-      ExpST _ used _    <- execCGM (go u)
+      ExpST _ used      <- execCGM (go u)
       return $ CUFun1 (mark used u)
              $ \xs -> evalState (evalCGM (go xs)) n
   --
   | otherwise
   = INTERNAL_ERROR(error) "codegenFun1" "expected unary function"
 
+
 codegenFun2
     :: forall aenv a b c. DeviceProperties
     -> Gamma aenv
@@ -219,14 +228,14 @@
         go :: (Rvalue x, Rvalue y) => [x] -> [y] -> Gen ([C.BlockItem], [C.Exp])
         go x y = do
           code  <- mapM use =<< codegenOpenExp dev aenv f (Empty `Push` map rvalue x `Push` map rvalue y)
-          env   <- getEnv
-          return (env, code)
+          env'  <- getEnv
+          return (env', code)
 
         (_,u,_)  = locals "undefined_x" (undefined :: a)
         (_,v,_)  = locals "undefined_y" (undefined :: b)
     in do
       n                 <- get
-      ExpST _ used _    <- execCGM (go u v)
+      ExpST _ used      <- execCGM (go u v)
       return $ CUFun2 (mark used u) (mark used v)
              $ \xs ys -> evalState (evalCGM (go xs ys)) n
   --
@@ -311,15 +320,14 @@
         Tuple t                 -> cvtT t env
         Prj i t                 -> prjT i t exp env
         Cond p t e              -> cond p t e env
-        Iterate n f x           -> iterate n f x env
---        While p f x             -> while p f x env
+        While p f x             -> while p f x env
 
         -- Shapes and indices
         IndexNil                -> return []
         IndexAny                -> return []
         IndexCons sh sz         -> (++) <$> cvtE sh env <*> cvtE sz env
-        IndexHead ix            -> return . last <$> cvtE ix env
-        IndexTail ix            ->          init <$> cvtE ix env
+        IndexHead ix            -> return . cindexHead <$> cvtE ix env
+        IndexTail ix            ->          cindexTail <$> cvtE ix env
         IndexSlice ix slix sh   -> indexSlice ix slix sh env
         IndexFull  ix slix sl   -> indexFull  ix slix sl env
         ToIndex sh ix           -> toIndex   sh ix env
@@ -349,9 +357,8 @@
     --
     elet :: DelayedOpenExp env aenv bnd -> DelayedOpenExp (env, bnd) aenv body -> Val env -> Gen [C.Exp]
     elet bnd body env = do
-      bnd'      <- cvtE bnd env
-      x         <- pushEnv bnd bnd'
-      body'     <- cvtE body (env `Push` x)
+      bnd'      <- cvtE bnd env >>= pushEnv bnd
+      body'     <- cvtE body (env `Push` bnd')
       return body'
 
     -- Convert an open expression into a sequence of C expressions. We retain
@@ -409,78 +416,62 @@
       ok        <- single "Cond" <$> pushEnv p p'
       zipWith (\a b -> [cexp| $exp:ok ? $exp:a : $exp:b |]) <$> cvtE t env <*> cvtE e env
 
-    -- Value recursion. Two flavours.
+    -- Value recursion
     --
-    iterate :: DelayedOpenExp env     aenv Int          -- fixed iteration depth
-            -> DelayedOpenExp (env,a) aenv a            -- loop body
-            -> DelayedOpenExp env     aenv a            -- initial value
-            -> Val env
-            -> Gen [C.Exp]
-    iterate n f x env
-      = do [n']         <- cvtE n env
-           x'           <- cvtE x env
-           var_x        <- mapM (\_ -> lift fresh) x'
-           var_n        <- lift fresh
-           let seed      = [cdecl| const int $id:var_n = $exp:n'; |]
-                         : zipWith3 (\t a v -> [cdecl| $ty:t $id:a = $exp:v; |]) (expType x) var_x x'
-               acc       = map cvar var_x
-
-           -- generate the loop in a clean environment, so that the previous
-           -- environment fragments are not included in the body
-           outer        <- gets bindings <* modify (\st -> st { bindings = [] })
-           body         <- cvtE f (env `Push` acc)
-           inner        <- getEnv
-           i            <- lift fresh
-
-           let go        = C.BlockStm
-                         $ [cstm| for (int $id:i = 0; $id:i < $id:var_n; ++ $id:i) {
-                                      $items:inner
-                                      $items:(acc .=. body)
-                                  } |]
-
-           -- restore the outer environment, plus the new loop
-           modify (\st -> st { bindings = go : map C.BlockDecl (reverse seed) ++ outer })
-           return acc
-
-{--
-    while :: DelayedOpenExp (env,a) aenv Bool           -- continue while predicate returns true
-          -> DelayedOpenExp (env,a) aenv a              -- loop body
-          -> DelayedOpenExp env     aenv a              -- initial value
+    while :: forall env a. Elt a
+          => DelayedOpenFun env aenv (a -> Bool)        -- continue while predicate returns true
+          -> DelayedOpenFun env aenv (a -> a)           -- loop body
+          -> DelayedOpenExp env aenv a                  -- initial value
           -> Val env
           -> Gen [C.Exp]
-    while p f x env
-      = do x'           <- cvtE x env
-
-           var_x        <- mapM (\_ -> lift fresh) x'
+    while test step x env
+      | Lam (Body p)    <- test
+      , Lam (Body f)    <- step
+      = do
+           -- Generate code for the initial value, then bind this to a fresh
+           -- (mutable) variable. We need build the declarations ourselves, and
+           -- twiddle the names a bit to avoid clobbering.
+           --
+           x'           <- cvtE x env
+           var_acc      <- lift fresh
            var_ok       <- lift fresh
 
-           let seed      = [cdecl| int $id:var_ok; |]
-                         : zipWith3 (\t a v -> [cdecl| $ty:t $id:a = $exp:v; |]) (expType x) var_x x'
-               acc       = map cvar var_x
-               ok        = cvar var_ok
+           let (tn_acc, acc, _)  = locals ('l':var_acc) (undefined :: a)
+               (tn_ok,  ok,  _)  = locals ('l':var_ok)  (undefined :: Bool)
 
-           -- generate the loop functions in a clean environment
-           outer        <- gets bindings <* modify (\st -> st { bindings = [] })
-           [done]       <- cvtE p (env `Push` acc)
-           envP         <- getEnv <* modify (\st -> st { bindings = [] })
-           body         <- cvtE f (env `Push` acc)
-           envF         <- getEnv
+           -- Generate code for the predicate and body expressions, with the new
+           -- names baked in directly. We can't use 'codegenFun1', because
+           -- def-use analysis won't be able to see into this new function.
+           --
+           -- However, we do need to generate the function with a clean set of
+           -- local bindings, and extract and new declarations afterwards.
+           --
+           let cvtF :: forall env t. Elt t => DelayedOpenExp env aenv t -> Val env -> Gen ([C.BlockItem], [C.Exp])
+               cvtF e env = do
+                 old  <- state (\s -> ( localBindings s, s { localBindings = []  } ))
+                 e'   <- cvtE e env
+                 env' <- state (\s -> ( localBindings s, s { localBindings = old } ))
+                 return (reverse env', e')
 
-           let go        =  envP
-                         ++ (ok .=. done)
-                         ++ [C.BlockStm
-                            [cstm| while ( $exp:ok ) {
-                                       $items:envF
-                                       $items:(acc .=. body)
-                                       $items:envP
-                                       $items:(ok .=. done)
-                                   } |]]
+           p'   <- cvtF p (env `Push` acc)
+           f'   <- cvtF f (env `Push` acc)
 
-           -- restore the outer environment, plus the new loop
-           modify (\st -> st { bindings = reverse (map C.BlockDecl seed ++ go) ++ outer })
+           -- Piece it all together. Note that declarations are added to the
+           -- localBindings in reverse order.
+           let loop = [citem| while ( $exp:(single "while" ok) ) {
+                                  $items:(acc .=. f')
+                                  $items:(ok  .=. p')
+                              } |]
+                    : (ok .=. p')
+                   ++ map     (\(t,n)   -> [citem| $ty:t $id:n ; |])      tn_ok
+                   ++ zipWith (\(t,n) v -> [citem| $ty:t $id:n = $v ; |]) tn_acc x'
+
+           modify (\s -> s { localBindings = loop ++ localBindings s })
            return acc
---}
 
+      | otherwise
+      = error "Would you say we'd be venturing into a zone of danger?"
+
     -- Restrict indices based on a slice specification. In the SliceAll case we
     -- elide the presence of IndexAny from the head of slx, as this is not
     -- represented in by any C term (Any ~ [])
@@ -520,29 +511,22 @@
       in
       replicate <$> cvtE slix env <*> cvtE sl env
 
-    -- Convert between linear and multidimensional indices. For the
-    -- multidimensional case, we've inlined the definition of 'fromIndex'
-    -- because we need to return an expression for each component.
+    -- Convert between linear and multidimensional indices
     --
     toIndex :: DelayedOpenExp env aenv sh -> DelayedOpenExp env aenv sh -> Val env -> Gen [C.Exp]
     toIndex sh ix env = do
-      sh'   <- cvtE sh env
-      ix'   <- cvtE ix env
-      return [ ccall "toIndex" [ ccall "shape" sh', ccall "shape" ix' ] ]
+      sh'   <- mapM use =<< cvtE sh env
+      ix'   <- mapM use =<< cvtE ix env
+      return [ ctoIndex sh' ix' ]
 
     fromIndex :: DelayedOpenExp env aenv sh -> DelayedOpenExp env aenv Int -> Val env -> Gen [C.Exp]
     fromIndex sh ix env = do
-      sh'   <- cvtE sh env
+      sh'   <- mapM use =<< cvtE sh env
       ix'   <- cvtE ix env
-      reverse <$> fromIndex' (reverse sh') (single "fromIndex" ix')
-      where
-        fromIndex' :: [C.Exp] -> C.Exp -> Gen [C.Exp]
-        fromIndex' []     _     = return []
-        fromIndex' [_]    i     = return [i]
-        fromIndex' (d:ds) i     = do
-          i'    <- bind [cty| int |] i
-          ds'   <- fromIndex' ds [cexp| $exp:i' / $exp:d |]
-          return $ [cexp| $exp:i' % $exp:d |] : ds'
+      tmp   <- lift fresh
+      let (ls, sz) = cfromIndex sh' (single "fromIndex" ix') tmp
+      modify (\st -> st { localBindings = reverse ls ++ localBindings st })
+      return sz
 
     -- Project out a single scalar element from an array. The array expression
     -- does not contain any free scalar variables (strictly flat data
@@ -566,8 +550,8 @@
       = let (sh, arr)   = namesOfAvar aenv idx
             ty          = accType acc
         in do
-        ix'     <- cvtE ix env
-        i       <- bind [cty| int |] $ ccall "toIndex" [ cvar sh, ccall "shape" ix' ]
+        ix'     <- mapM use =<< cvtE ix env
+        i       <- bind cint $ ctoIndex (cshape (expDim ix) sh) ix'
         return   $ zipWith (\t a -> indexArray dev t (cvar a) i) ty arr
       --
       | otherwise
@@ -584,7 +568,7 @@
       = let (_, arr)    = namesOfAvar aenv idx
             ty          = accType acc
         in do
-        ix'     <- cvtE ix env
+        ix'     <- mapM use =<< cvtE ix env
         i       <- bind [cty| int |] $ single "LinearIndex" ix'
         return   $ zipWith (\t a -> indexArray dev t (cvar a) i) ty arr
       --
@@ -600,7 +584,7 @@
     shape :: (Shape sh, Elt e) => DelayedOpenAcc aenv (Array sh e) -> Val env -> Gen [C.Exp]
     shape acc _env
       | Manifest (Avar idx) <- acc
-      = return $ cshape (delayedDim acc) (cvar $ fst (namesOfAvar aenv idx))
+      = return $ cshape (delayedDim acc) (fst (namesOfAvar aenv idx))
 
       | otherwise
       = INTERNAL_ERROR(error) "Shape" "expected array variable"
@@ -609,11 +593,7 @@
     -- definition is inlined, but we could also call the C function helpers.
     --
     shapeSize :: DelayedOpenExp env aenv sh -> Val env -> Gen [C.Exp]
-    shapeSize sh env =
-      let size [] = return $ [cexp| 1 |]
-          size ss = return $ foldl1 (\a b -> [cexp| $exp:a * $exp:b |]) ss
-      in
-      size <$> cvtE sh env
+    shapeSize sh env = return . csize <$> cvtE sh env
 
     -- Intersection of two shapes, taken as the minimum in each dimension.
     --
@@ -621,10 +601,8 @@
               => DelayedOpenExp env aenv sh
               -> DelayedOpenExp env aenv sh
               -> Val env -> Gen [C.Exp]
-    intersect sh1 sh2 env = let
-        sh1' = ccastTup (Sugar.eltType (undefined::sh)) <$> cvtE sh1 env
-        sh2' = ccastTup (Sugar.eltType (undefined::sh)) <$> cvtE sh2 env
-      in zipWith (\a b -> ccall "min" [a,b]) <$> sh1' <*> sh2'
+    intersect sh1 sh2 env =
+      zipWith (\a b -> ccall "min" [a,b]) <$> cvtE sh1 env <*> cvtE sh2 env
 
     -- Foreign scalar functions. We need to extract any header files that might
     -- be required so they can be added to the top level definitions.
@@ -638,14 +616,12 @@
              -> Val env
              -> Gen [C.Exp]
     foreignE ff x env = case canExecuteExp ff of
-      Nothing   -> INTERNAL_ERROR(error) "codegenOpenExp" "Non-CUDA foreign expression encountered"
-      Just f    -> do
-        unless (null hdr) . lift $ modify (\st -> st { headers = Set.insert hdr (headers st) })
+      Nothing      -> INTERNAL_ERROR(error) "codegenOpenExp" "Non-CUDA foreign expression encountered"
+      Just (hs, f) -> do
+        lift $ modify (\st -> st { headers = foldl (flip Set.insert) (headers st) hs })
         args    <- cvtE x env
-        return  $  [ccall name (ccastTup (Sugar.eltType (undefined::a)) args)]
-        where
-          (hdr, rest)   = break isSpace f
-          name          = if null rest then f else tail rest
+        mapM_ use args
+        return  $  [ccall f (ccastTup (Sugar.eltType (undefined::a)) args)]
 
     -- Some terms demand we extract only singly typed expressions
     --
@@ -672,10 +648,8 @@
 codegenPrim (PrimSig             ty) [a]   = codegenSig ty a
 codegenPrim (PrimQuot             _) [a,b] = [cexp|$exp:a / $exp:b|]
 codegenPrim (PrimRem              _) [a,b] = [cexp|$exp:a % $exp:b|]
-codegenPrim (PrimIDiv            ty) [a,b] = ccall "idiv" [ccast (NumScalarType $ IntegralNumType ty) a,
-                                                           ccast (NumScalarType $ IntegralNumType ty) b]
-codegenPrim (PrimMod             ty) [a,b] = ccall "mod"  [ccast (NumScalarType $ IntegralNumType ty) a,
-                                                           ccast (NumScalarType $ IntegralNumType ty) b]
+codegenPrim (PrimIDiv             _) [a,b] = ccall "idiv" [a,b]
+codegenPrim (PrimMod              _) [a,b] = ccall "mod"  [a,b]
 codegenPrim (PrimBAnd             _) [a,b] = [cexp|$exp:a & $exp:b|]
 codegenPrim (PrimBOr              _) [a,b] = [cexp|$exp:a | $exp:b|]
 codegenPrim (PrimBXor             _) [a,b] = [cexp|$exp:a ^ $exp:b|]
@@ -807,7 +781,10 @@
     one  | IntegralDict <- integralDict ty = codegenIntegralScalar ty 1
 
 codegenFloatingSig :: FloatingType a -> C.Exp -> C.Exp
-codegenFloatingSig ty x = [cexp|$exp:x == $exp:zero ? $exp:zero : $exp:(ccall (FloatingNumType ty `postfix` "copysign") [one,x]) |]
+codegenFloatingSig ty x =
+  [cexp|$exp:x == $exp:zero
+            ? $exp:zero
+            : $exp:(ccall (FloatingNumType ty `postfix` "copysign") [one,x]) |]
   where
     zero | FloatingDict <- floatingDict ty = codegenFloatingScalar ty 0
     one  | FloatingDict <- floatingDict ty = codegenFloatingScalar ty 1
diff --git a/Data/Array/Accelerate/CUDA/CodeGen/Base.hs b/Data/Array/Accelerate/CUDA/CodeGen/Base.hs
--- a/Data/Array/Accelerate/CUDA/CodeGen/Base.hs
+++ b/Data/Array/Accelerate/CUDA/CodeGen/Base.hs
@@ -1,11 +1,11 @@
 {-# LANGUAGE CPP                   #-}
 {-# LANGUAGE FlexibleInstances     #-}
 {-# LANGUAGE GADTs                 #-}
+{-# LANGUAGE ImpredicativeTypes    #-}
 {-# LANGUAGE MultiParamTypeClasses #-}
 {-# LANGUAGE OverlappingInstances  #-}
 {-# LANGUAGE PatternGuards         #-}
 {-# LANGUAGE QuasiQuotes           #-}
-{-# LANGUAGE RankNTypes            #-}
 {-# LANGUAGE ScopedTypeVariables   #-}
 -- |
 -- Module      : Data.Array.Accelerate.CUDA.CodeGen.Base
@@ -22,11 +22,13 @@
 
   -- Names and Types
   CUTranslSkel(..), CUDelayedAcc(..), CUExp(..), CUFun1(..), CUFun2(..),
+  Eliminate, Instantiate1, Instantiate2,
   Name, namesOfArray, namesOfAvar, groupOfInt,
 
   -- Declaration generation
-  cvar, ccall, cchar, cintegral, cbool, cdim, cshape, getters, setters, shared,
-  indexArray, indexHead, shapeSize, environment, arrayAsTex, arrayAsArg,
+  cint, cvar, ccall, cchar, cintegral, cbool, cshape, csize, cindexHead, cindexTail, ctoIndex, cfromIndex,
+  readArray, writeArray, shared,
+  indexArray, environment, arrayAsTex, arrayAsArg,
   umul24, gridSize, threadIdx,
 
   -- Mutable operations
@@ -34,12 +36,19 @@
 
 ) where
 
+-- library
+import Prelude                                          hiding ( zipWith, zipWith3 )
+import Data.List                                        ( mapAccumR )
 import Text.PrettyPrint.Mainland
 import Language.C.Quote.CUDA
 import qualified Language.C.Syntax                      as C
 import qualified Data.HashMap.Strict                    as Map
 
+-- cuda
 import Foreign.CUDA.Analysis.Device
+
+-- friends
+import Data.Array.Accelerate.Type
 import Data.Array.Accelerate.Array.Sugar                ( Array, Shape, Elt )
 import Data.Array.Accelerate.Analysis.Shape
 import Data.Array.Accelerate.CUDA.CodeGen.Type
@@ -59,7 +68,7 @@
     -> (Name, [Name])   -- shape and array field names
 namesOfArray grp _
   = let ty      = eltType (undefined :: e)
-        arr x   = "arr" ++ grp ++ "_a" ++ show x
+        arr x   = "arr" ++ grp ++ '_':show x
         n       = length ty
     in
     ( "sh" ++ grp, map arr [n-1, n-2 .. 0] )
@@ -98,17 +107,21 @@
 
 -- Scalar functions of particular arity, with local bindings.
 --
+type Eliminate a        = forall x. [x] -> [(Bool,x)]
+type Instantiate1 a b   = forall x. Rvalue x => [x] -> ([C.BlockItem], [C.Exp])
+type Instantiate2 a b c = forall x y. (Rvalue x, Rvalue y) => [x] -> [y] -> ([C.BlockItem], [C.Exp])
+
 data CUFun1 aenv f where
   CUFun1 :: (Elt a, Elt b)
-         => (forall x. [x] -> [(Bool,x)])
-         -> (forall x. Rvalue x => [x] -> ([C.BlockItem], [C.Exp]))
+         => Eliminate a
+         -> Instantiate1 a b
          -> CUFun1 aenv (a -> b)
 
 data CUFun2 aenv f where
   CUFun2 :: (Elt a, Elt b, Elt c)
-         => (forall x. [x] -> [(Bool,x)])
-         -> (forall y. [y] -> [(Bool,y)])
-         -> (forall x y. (Rvalue x, Rvalue y) => [x] -> [y] -> ([C.BlockItem], [C.Exp]))
+         => Eliminate a
+         -> Eliminate b
+         -> Instantiate2 a b c
          -> CUFun2 aenv (a -> b -> c)
 
 -- Delayed arrays
@@ -120,9 +133,12 @@
             -> CUDelayedAcc aenv sh e
 
 
--- Expression and declaration generation
--- -------------------------------------
+-- Common expression forms
+-- -----------------------
 
+cint :: C.Type
+cint = codegenScalarType (scalarType :: ScalarType Int)
+
 cvar :: Name -> C.Exp
 cvar x = [cexp|$id:x|]
 
@@ -138,41 +154,69 @@
 cbool :: Bool -> C.Exp
 cbool = cintegral . fromEnum
 
-cdim :: Name -> Int -> C.Definition
-cdim name n = [cedecl|typedef typename $id:("DIM" ++ show n) $id:name;|]
+-- Generate all the names of a shape given a base name and dimensionality
+cshape :: Int -> Name -> [C.Exp]
+cshape dim sh = [ cvar x | x <- cshape' dim sh ]
 
--- Disassemble a struct-shape into a list of expressions accessing the fields
-cshape :: Int -> C.Exp -> [C.Exp]
-cshape dim sh
-  | dim == 0  = []
-  | dim == 1  = [sh]
-  | otherwise = map (\i -> [cexp|$exp:sh . $id:('a':show i)|]) [dim-1, dim-2 .. 0]
+cshape' :: Int -> Name -> [Name]
+cshape' dim sh = [ (sh ++ '_':show i) | i <- [dim-1, dim-2 .. 0] ]
 
--- Calculate the size of a shape from its component dimensions
-shapeSize :: Rvalue r => [r] -> C.Exp
-shapeSize [] = [cexp| 1 |]
-shapeSize ss = foldl1 (\a b -> [cexp| $exp:a * $exp:b |]) (map rvalue ss)
+-- Get the innermost index of a shape/index
+cindexHead :: Rvalue r => [r] -> C.Exp
+cindexHead = rvalue . last
 
-indexHead :: Rvalue r => [r] -> C.Exp
-indexHead = rvalue . last
+-- Get the tail of a shape/index
+cindexTail :: Rvalue r => [r] -> [C.Exp]
+cindexTail = map rvalue . init
 
+-- generate code that calculates the product of the list of expressions
+csize :: Rvalue r => [r] -> C.Exp
+csize [] = [cexp| 1 |]
+csize ss = foldr1 (\a b -> [cexp| $exp:a * $exp:b |]) (map rvalue ss)
 
--- Thread blocks and indices
+-- Generate code to calculate a linear from a multi-dimensional index (given an array shape).
 --
+ctoIndex :: (Rvalue sh, Rvalue ix) => [sh] -> [ix] -> C.Exp
+ctoIndex extent index
+  = toIndex (reverse $ map rvalue extent) (reverse $ map rvalue index)  -- we use a row-major representation
+  where
+    toIndex []      []     = [cexp| $int:(0::Int) |]
+    toIndex [_]     [i]    = i
+    toIndex (sz:sh) (i:ix) = [cexp| $exp:(toIndex sh ix) * $exp:sz + $exp:i |]
+    toIndex _       _      = INTERNAL_ERROR(error) "toIndex" "argument mismatch"
+
+-- Generate code to calculate a multi-dimensional index from a linear index and a given array shape.
+-- This version creates temporary values that are reused in the computation.
+--
+cfromIndex :: (Rvalue sh, Rvalue ix) => [sh] -> ix -> Name -> ([C.BlockItem], [C.Exp])
+cfromIndex shName ixName tmpName = fromIndex (map rvalue shName) (rvalue ixName)
+  where
+    fromIndex [sh]   ix = ([], [[cexp| ({ assert( $exp:ix >= 0 && $exp:ix < $exp:sh ); $exp:ix; }) |]])
+    fromIndex extent ix = let ((env, _, _), sh) = mapAccumR go ([], ix, 0) extent
+                          in  (reverse env, sh)
+
+    go (tmps,ix,n) d
+      = let tmp         = tmpName ++ '_':show (n::Int)
+            ix'         = [citem| const $ty:cint $id:tmp = $exp:ix ; |]
+        in
+        ((ix':tmps, [cexp| $id:tmp / $exp:d |], n+1), [cexp| $id:tmp % $exp:d |])
+
+
+-- Thread blocks and indices
+-- -------------------------
+
 umul24 :: DeviceProperties -> C.Exp -> C.Exp -> C.Exp
 umul24 dev x y
   | computeCapability dev < Compute 2 0 = [cexp| __umul24($exp:x, $exp:y) |]
   | otherwise                           = [cexp| $exp:x * $exp:y |]
 
 gridSize :: DeviceProperties -> C.Exp
-gridSize dev
-  | computeCapability dev < Compute 2 0 = [cexp| __umul24(blockDim.x, gridDim.x) |]
-  | otherwise                           = [cexp| blockDim.x * gridDim.x |]
+gridSize dev = umul24 dev [cexp|blockDim.x|] [cexp|gridDim.x|]
 
 threadIdx :: DeviceProperties -> C.Exp
-threadIdx dev
-  | computeCapability dev < Compute 2 0 = [cexp| __umul24(blockDim.x, blockIdx.x) + threadIdx.x |]
-  | otherwise                           = [cexp| blockDim.x * blockIdx.x + threadIdx.x |]
+threadIdx dev =
+  let block = umul24 dev [cexp|blockDim.x|] [cexp|blockIdx.x|]
+  in  [cexp| $exp:block + threadIdx.x |]
 
 
 -- Generate an array indexing expression. Depending on the hardware class, this
@@ -186,31 +230,34 @@
     -> C.Exp
 indexArray dev elt arr ix
   -- use the L2 cache of newer devices
-  | computeCapability dev >= Compute 2 0                = [cexp| $exp:arr [ $exp:ix ] |]
+  | computeCapability dev >= Compute 2 0 = [cexp| $exp:arr [ $exp:ix ] |]
 
   -- use the texture cache of compute 1.x devices
-  | C.Type (C.DeclSpec _ _ (C.Tdouble _) _) _ _ <- elt  = ccall "indexDArray" [arr, ix]
-  | otherwise                                           = ccall "indexArray"  [arr, ix]
+  | [cty|double|] <- elt                 = ccall "indexDArray" [arr, ix]
+  | otherwise                            = ccall "indexArray"  [arr, ix]
 
 
+-- Kernel function parameters
+-- --------------------------
+
 -- Generate kernel parameters for an array valued argument, and a function to
 -- linearly index this array. Note that dimensional indexing results in error.
 --
-getters
+readArray
     :: forall aenv sh e. (Shape sh, Elt e)
     => Name                             -- group names
     -> Array sh e                       -- dummy to fix types
     -> ( [C.Param], CUDelayedAcc aenv sh e )
-getters grp dummy
+readArray grp dummy
   = let (sh, arrs)      = namesOfArray grp (undefined :: e)
         args            = arrayAsArg dummy grp
 
         dim             = expDim (undefined :: Exp aenv sh)
-        sh'             = cshape dim (cvar sh)
+        sh'             = cshape dim sh
         get ix          = ([], map (\a -> [cexp| $id:a [ $exp:ix ] |]) arrs)
         manifest        = CUDelayed (CUExp ([], sh'))
-                                    (INTERNAL_ERROR(error) "getters" "linear indexing only")
-                                    (CUFun1 (zip (repeat True)) (get . rvalue . head))
+                                    (INTERNAL_ERROR(error) "readArray" "linear indexing only")
+                                    (CUFun1 (zip (repeat True)) (\[i] -> get (rvalue i)))
     in ( args, manifest )
 
 
@@ -220,20 +267,24 @@
 -- name (say "Out") to be welded with a shape name "shOut" followed by the
 -- non-parametric array data "arrOut_aX".
 --
-setters
+writeArray
     :: forall sh e. (Shape sh, Elt e)
     => Name                             -- group names
     -> Array sh e                       -- dummy to fix types
-    -> ([C.Param], Name -> [C.Exp])
-setters grp _
-  = let (sh, arrs)      = namesOfArray grp (undefined :: e)
-        dim             = expDim (undefined :: Exp aenv sh)
-        sh'             = [cparam| const typename $id:("DIM" ++ show dim) $id:sh |]
-        arrs'           = zipWith (\t n -> [cparam| $ty:t * __restrict__ $id:n |]) (eltType (undefined :: e)) arrs
-    in
-    ( sh' : arrs'
-    , \ix -> map (\a -> [cexp| $id:a [ $id:ix ] |]) arrs
-    )
+    -> ( [C.Param]                      -- function parameters to marshal the output array
+       , [C.Exp]                        -- the shape of the output array
+       , Rvalue x => x -> [C.Exp] )     -- write an element at a given index
+writeArray grp _ =
+  let (sh, arrs)        = namesOfArray grp (undefined :: e)
+      dim               = expDim (undefined :: Exp aenv sh)
+      sh'               = cshape' dim sh
+      extent            = [ [cparam| const $ty:cint $id:i |] | i <- sh' ]
+      adata             = zipWith (\t n -> [cparam| $ty:t * __restrict__ $id:n |]) (eltType (undefined :: e)) arrs
+  in
+  ( extent ++ adata
+  , map cvar sh'
+  , \ix -> map (\a -> [cexp| $id:a [ $exp:(rvalue ix) ] |]) arrs
+  )
 
 
 -- All dynamically allocated __shared__ memory will begin at the same base
@@ -247,7 +298,8 @@
     -> Name                             -- group name
     -> C.Exp                            -- how much shared memory per type
     -> Maybe C.Exp                      -- (optional) initialise from this base address
-    -> ([C.InitGroup], Name -> [C.Exp]) -- shared memory declaration and indexing function
+    -> ( [C.InitGroup]                  -- shared memory declaration and...
+       , Rvalue x => x -> [C.Exp])      -- ...indexing function
 shared _ grp size mprev
   = let e:es                    = eltType (undefined :: e)
         x:xs                    = let k = length es in map (\n -> grp ++ show n) [k, k-1 .. 0]
@@ -257,8 +309,8 @@
           | Just p <- mprev     = [cdecl| volatile $ty:t * $id:v = ($ty:t *) $exp:p; |]
           | otherwise           = [cdecl| extern volatile __shared__ $ty:t $id:v [] ; |]
     in
-    ( sbase e x : zipWith3 sdata es xs (x:xs)
-    , \ix -> map (\v -> [cexp| $id:v [ $id:ix ] |]) (x:xs)
+    ( sbase e x : zipWith3 sdata es xs (init (x:xs))
+    , \ix -> map (\v -> [cexp| $id:v [ $exp:(rvalue ix) ] |]) (x:xs)
     )
 
 -- Array environment references. The method in which arrays are accessed depends
@@ -281,36 +333,36 @@
     -> ([C.Definition], [C.Param])
 environment dev gamma@(Gamma aenv)
   | computeCapability dev < Compute 2 0
-  = Map.foldrWithKey (\(Idx_ v) _ (ds,ps) -> let (d,p) = asTex v in (d++ds, p:ps)) ([],[]) aenv
+  = Map.foldrWithKey (\(Idx_ v) _ (ds,ps) -> let (d,p) = asTex v in (d++ds, p++ps)) ([],[]) aenv
 
   | otherwise
   = ([], Map.foldrWithKey (\(Idx_ v) _ vs -> asArg v ++ vs) [] aenv)
 
   where
-    asTex :: forall sh e. (Shape sh, Elt e) => Idx aenv (Array sh e) -> ([C.Definition], C.Param)
+    asTex :: forall sh e. (Shape sh, Elt e) => Idx aenv (Array sh e) -> ([C.Definition], [C.Param])
     asTex ix = arrayAsTex (undefined :: Array sh e) (groupOfAvar gamma ix)
 
     asArg :: forall sh e. (Shape sh, Elt e) => Idx aenv (Array sh e) -> [C.Param]
     asArg ix = arrayAsArg (undefined :: Array sh e) (groupOfAvar gamma ix)
 
 
-arrayAsTex :: forall sh e. (Shape sh, Elt e) => Array sh e -> Name -> ([C.Definition], C.Param)
+arrayAsTex :: forall sh e. (Shape sh, Elt e) => Array sh e -> Name -> ([C.Definition], [C.Param])
 arrayAsTex _ grp =
   let (sh, arrs)        = namesOfArray grp (undefined :: e)
       dim               = expDim (undefined :: Exp aenv sh)
-      sh'               = [cparam| const typename $id:("DIM" ++ show dim) $id:sh |]
-      arrs'             = zipWith (\t a -> [cedecl| static $ty:t $id:a; |]) (eltTypeTex (undefined :: e)) arrs
+      extent            = [ [cparam| const $ty:cint $id:i |] | i <- cshape' dim sh ]
+      adata             = zipWith (\t a -> [cedecl| static $ty:t $id:a; |]) (eltTypeTex (undefined :: e)) arrs
   in
-  (arrs', sh')
+  (adata, extent)
 
 arrayAsArg :: forall sh e. (Shape sh, Elt e) => Array sh e -> Name -> [C.Param]
 arrayAsArg _ grp =
   let (sh, arrs)        = namesOfArray grp (undefined :: e)
       dim               = expDim (undefined :: Exp aenv sh)
-      sh'               = [cparam| const typename $id:("DIM" ++ show dim) $id:sh |]
-      arrs'             = zipWith (\t n -> [cparam| const $ty:t * __restrict__ $id:n |]) (eltType (undefined :: e)) arrs
+      extent            = [ [cparam| const $ty:cint $id:i |] | i <- cshape' dim sh ]
+      adata             = zipWith (\t n -> [cparam| const $ty:t * __restrict__ $id:n |]) (eltType (undefined :: e)) arrs
   in
-  sh' : arrs'
+  extent ++ adata
 
 
 -- Mutable operations
@@ -337,10 +389,10 @@
   lvalue :: a -> C.Exp -> C.BlockItem
 
 instance Lvalue C.Exp where
-  lvalue x y = C.BlockStm  [cstm| $exp:x = $exp:y; |]
+  lvalue x y = [citem| $exp:x = $exp:y; |]
 
 instance Lvalue (C.Type, Name) where
-  lvalue (t,x) y = C.BlockDecl [cdecl| const $ty:t $id:x = $exp:y; |]
+  lvalue (t,x) y = [citem| const $ty:t $id:x = $exp:y; |]
 
 
 class Rvalue a where
@@ -349,8 +401,11 @@
 instance Rvalue C.Exp where
   rvalue = id
 
+instance Rvalue Name where
+  rvalue = cvar
+
 instance Rvalue (C.Type, Name) where
-  rvalue (_,x) = cvar x
+  rvalue (_,x) = rvalue x
 
 
 infixr 0 .=.
@@ -361,7 +416,7 @@
   assign :: l -> r -> [C.BlockItem]
 
 instance (Lvalue l, Rvalue r) => Assign l r where
-  assign lhs rhs = return $ lvalue lhs (rvalue rhs)
+  assign lhs rhs = [ lvalue lhs (rvalue rhs) ]
 
 instance Assign l r => Assign (Bool,l) r where
   assign (used,lhs) rhs
@@ -375,4 +430,20 @@
 
 instance Assign l r => Assign l ([C.BlockItem], r) where
   assign lhs (env, rhs) = env ++ assign lhs rhs
+
+
+-- Prelude'
+-- --------
+
+-- A version of zipWith that requires the lists to be equal length
+--
+zipWith :: (a -> b -> c) -> [a] -> [b] -> [c]
+zipWith f (x:xs) (y:ys) = f x y : zipWith f xs ys
+zipWith _ []     []     = []
+zipWith _ _      _      = INTERNAL_ERROR(error) "zipWith" "argument mismatch"
+
+zipWith3 :: (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
+zipWith3 f (x:xs) (y:ys) (z:zs) = f x y z : zipWith3 f xs ys zs
+zipWith3 _ []     []     []     = []
+zipWith3 _ _      _      _      = INTERNAL_ERROR(error) "zipWith3" "argument mismatch"
 
diff --git a/Data/Array/Accelerate/CUDA/CodeGen/IndexSpace.hs b/Data/Array/Accelerate/CUDA/CodeGen/IndexSpace.hs
--- a/Data/Array/Accelerate/CUDA/CodeGen/IndexSpace.hs
+++ b/Data/Array/Accelerate/CUDA/CodeGen/IndexSpace.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE CPP                 #-}
 {-# LANGUAGE GADTs               #-}
 {-# LANGUAGE PatternGuards       #-}
 {-# LANGUAGE QuasiQuotes         #-}
@@ -28,13 +29,14 @@
 import Foreign.CUDA.Analysis.Device
 import qualified Language.C.Syntax                      as C
 
-import Data.Array.Accelerate.Array.Sugar                ( Array, Shape, Elt )
-import Data.Array.Accelerate.Analysis.Shape
-import Data.Array.Accelerate.CUDA.AST                   ( Gamma, Exp )
+import Data.Array.Accelerate.Array.Sugar                ( Array, Shape, Elt, ignore, shapeToList )
+import Data.Array.Accelerate.CUDA.AST                   ( Gamma )
 import Data.Array.Accelerate.CUDA.CodeGen.Type
 import Data.Array.Accelerate.CUDA.CodeGen.Base
 
+#include "accelerate.h"
 
+
 -- Construct a new array by applying a function to each index. Each thread
 -- processes multiple elements, striding the array by the grid size.
 --
@@ -49,12 +51,11 @@
     -> Gamma aenv
     -> CUFun1 aenv (sh -> e)
     -> [CUTranslSkel aenv (Array sh e)]
-mkGenerate dev aenv (CUFun1 _ f)
+mkGenerate dev aenv (CUFun1 dce f)
   = return
   $ CUTranslSkel "generate" [cunit|
 
-    $esc:("#include <accelerate_cuda_extras.h>")
-    $edecl:(cdim "DimOut" dim)
+    $esc:("#include <accelerate_cuda.h>")
     $edecls:texIn
 
     extern "C" __global__ void
@@ -64,7 +65,7 @@
         $params:argOut
     )
     {
-        const int shapeSize     = size(shOut);
+        const int shapeSize     = $exp:(csize shOut);
         const int gridSize      = $exp:(gridSize dev);
               int ix;
 
@@ -72,17 +73,15 @@
             ; ix <  shapeSize
             ; ix += gridSize )
         {
-            const typename DimOut sh = fromIndex(shOut, ix);
-
+            $items:(dce sh      .=. cfromIndex shOut "ix" "tmp")
             $items:(setOut "ix" .=. f sh)
         }
     }
   |]
   where
-    dim                 = expDim (undefined :: Exp aenv sh)
-    sh                  = cshape dim (cvar "sh")
-    (texIn, argIn)      = environment dev aenv
-    (argOut, setOut)    = setters "Out" (undefined :: Array sh e)
+    (sh, _, _)                  = locals "sh" (undefined :: sh)
+    (texIn, argIn)              = environment dev aenv
+    (argOut, shOut, setOut)     = writeArray "Out" (undefined :: Array sh e)
 
 
 -- A combination map/backpermute, where the index and value transformations have
@@ -104,15 +103,13 @@
     -> CUDelayedAcc aenv sh a
     -> [CUTranslSkel aenv (Array sh' b)]
 mkTransform dev aenv perm fun arr
-  | CUFun1 _ p                   <- perm
-  , CUFun1 dce f                 <- fun
-  , CUDelayed _ (CUFun1 _ get) _ <- arr
+  | CUFun1 dce_p p                      <- perm
+  , CUFun1 dce_f f                      <- fun
+  , CUDelayed _ (CUFun1 dce_g get) _    <- arr
   = return
   $ CUTranslSkel "transform" [cunit|
 
-    $esc:("#include <accelerate_cuda_extras.h>")
-    $edecl:(cdim "DimOut" dimOut)
-    $edecl:(cdim "DimIn"  dimIn)
+    $esc:("#include <accelerate_cuda.h>")
     $edecls:texIn
 
     extern "C" __global__ void
@@ -122,7 +119,7 @@
         $params:argOut
     )
     {
-        const int shapeSize     = size(shOut);
+        const int shapeSize     = $exp:(csize shOut);
         const int gridSize      = $exp:(gridSize dev);
               int ix;
 
@@ -130,21 +127,19 @@
             ; ix <  shapeSize
             ; ix += gridSize )
         {
-            const typename DimOut sh_ = fromIndex(shOut, ix);
-            $items:(sh          .=. p sh_)
-            $items:(dce x0      .=. get sh)
+            $items:(dce_p sh'   .=. cfromIndex shOut "ix" "tmp")
+            $items:(dce_g sh    .=. p sh')
+            $items:(dce_f x0    .=. get sh)
             $items:(setOut "ix" .=. f x0)
         }
     }
   |]
   where
-    dimIn               = expDim (undefined :: Exp aenv sh)
-    dimOut              = expDim (undefined :: Exp aenv sh')
-    sh_                 = cshape dimOut (cvar "sh_")
-    (texIn, argIn)      = environment dev aenv
-    (argOut, setOut)    = setters "Out" (undefined :: Array sh' b)
-    (x0, _, _)          = locals "x"  (undefined :: a)
-    (sh, _, _)          = locals "sh" (undefined :: sh)
+    (texIn, argIn)              = environment dev aenv
+    (argOut, shOut, setOut)     = writeArray "Out" (undefined :: Array sh' b)
+    (x0, _, _)                  = locals "x"   (undefined :: a)
+    (sh, _, _)                  = locals "sh"  (undefined :: sh)
+    (sh', _, _)                 = locals "sh_" (undefined :: sh')
 
 
 -- Forward permutation specified by an index mapping that determines for each
@@ -171,14 +166,12 @@
     -> CUFun1 aenv (sh -> sh')
     -> CUDelayedAcc aenv sh e
     -> [CUTranslSkel aenv (Array sh' e)]
-mkPermute dev aenv (CUFun2 dcex dcey combine) (CUFun1 _ prj) arr
+mkPermute dev aenv (CUFun2 dce_x dce_y combine) (CUFun1 dce_p prj) arr
   | CUDelayed (CUExp shIn) _ (CUFun1 _ get) <- arr
   = return
   $ CUTranslSkel "permute" [cunit|
 
-    $esc:("#include <accelerate_cuda_extras.h>")
-    $edecl:(cdim "DimOut" dimOut)
-    $edecl:(cdim "DimIn"  dimIn)
+    $esc:("#include <accelerate_cuda.h>")
     $edecls:texIn
 
     extern "C" __global__ void
@@ -188,9 +181,12 @@
         $params:argOut
     )
     {
+        /*
+         * The input shape might be a complex expression. Evaluate it first to reuse the result.
+         */
         $items:(sh .=. shIn)
-        const typename DimIn shIn       = $exp:(ccall "shape" (map rvalue sh));
-        const int shapeSize             = $exp:(shapeSize sh);
+
+        const int shapeSize             = $exp:(csize sh);
         const int gridSize              = $exp:(gridSize dev);
               int ix;
 
@@ -198,19 +194,17 @@
             ; ix <  shapeSize
             ; ix += gridSize )
         {
-            typename DimOut dst;
-
-            const typename DimIn src = fromIndex( shIn, ix );
-            $items:(dst .=. prj src)
+            $items:(dce_p src   .=. cfromIndex sh "ix" "srcTmp")
+            $items:(dst         .=. prj src)
 
-            if ( !ignore(dst) )
+            if ( ! $exp:(cignore dst) )
             {
                 $decls:decly
                 $decls:decly'
 
-                const int jx = toIndex(shOut, dst);
-                $items:(dcex x .=. get ix)
-                $items:(dcey y .=. arrOut "jx")
+                $items:(jx      .=. ctoIndex shOut dst)
+                $items:(dce_x x .=. get ix)
+                $items:(dce_y y .=. arrOut jx)
 
                 $items:write
             }
@@ -218,20 +212,27 @@
     }
   |]
   where
-    dimIn               = expDim (undefined :: Exp aenv sh)
-    dimOut              = expDim (undefined :: Exp aenv sh')
-    sizeof              = eltSizeOf (undefined::e)
-    (texIn, argIn)      = environment dev aenv
-    (argOut, arrOut)    = setters "Out" (undefined :: Array sh' e)
-    (sh, _, _)          = locals "sh" (undefined :: sh)
-    (x, _, _)           = locals "x"  (undefined :: e)
-    (_, y,  decly)      = locals "y"  (undefined :: e)
-    (_, y', decly')     = locals "_y" (undefined :: e)
-    ix                  = [cvar "ix"]
-    src                 = cshape dimIn  (cvar "src")
-    dst                 = cshape dimOut (cvar "dst")
-    sm                  = computeCapability dev
+    sizeof                      = eltSizeOf (undefined::e)
+    (texIn, argIn)              = environment dev aenv
+    (argOut, shOut, arrOut)     = writeArray "Out" (undefined :: Array sh' e)
+    (x, _, _)                   = locals "x" (undefined :: e)
+    (_, y,  decly)              = locals "y" (undefined :: e)
+    (_, y', decly')             = locals "_y" (undefined :: e)
+    (sh, _, _)                  = locals "shIn" (undefined :: sh)
+    (src, _, _)                 = locals "sh" (undefined :: sh)
+    (dst, _, _)                 = locals "sh_" (undefined :: sh')
+    ([jx], _, _)                = locals "jx" (undefined :: Int)
+    ix                          = [cvar "ix"]
+    sm                          = computeCapability dev
 
+    -- If the destination index resolves to the magic index "ignore", the result
+    -- is dropped from the output array.
+    cignore :: Rvalue x => [x] -> C.Exp
+    cignore []  = INTERNAL_ERROR(error) "permute" "singleton arrays not supported"
+    cignore xs  = foldl1 (\a b -> [cexp| $exp:a && $exp:b |])
+                $ zipWith (\a b -> [cexp| $exp:(rvalue a) == $int:b |]) xs
+                $ shapeToList (ignore :: sh')
+
     -- Apply the combining function between old and new values. If multiple
     -- threads attempt to write to the same location, the hardware
     -- write-combining mechanism will accept one transaction and all other
@@ -245,21 +246,20 @@
     -- Each element of a tuple is necessarily written individually, so the tuple
     -- as a whole is not stored atomically.
     --
-    write               = env ++ zipWith6 apply sizeof (arrOut "jx") fun x (dcey y) y'
-    (env, fun)          = combine x y
+    write       = env ++ zipWith6 apply sizeof (arrOut jx) fun x (dce_y y) y'
+    (env, fun)  = combine x y
 
     apply size out f x1 (used,y1) y1'
       | used
       , Just atomicCAS <- reinterpret size
-      = C.BlockStm
-        [cstm| do {
-                      $exp:y1' = $exp:y1;
-                      $exp:y1  = $exp:atomicCAS ( & $exp:out, $exp:y1', $exp:f );
+      = [citem| do {
+                       $exp:y1' = $exp:y1;
+                       $exp:y1  = $exp:atomicCAS ( & $exp:out, $exp:y1', $exp:f );
 
-                  } while ( $exp:y1 != $exp:y1' ); |]
+                   } while ( $exp:y1 != $exp:y1' ); |]
 
       | otherwise
-      = C.BlockStm [cstm| $exp:out = $exp:(rvalue x1); |]
+      = [citem| $exp:out = $exp:(rvalue x1); |]
 
     --
     reinterpret :: Int -> Maybe C.Exp
diff --git a/Data/Array/Accelerate/CUDA/CodeGen/Mapping.hs b/Data/Array/Accelerate/CUDA/CodeGen/Mapping.hs
--- a/Data/Array/Accelerate/CUDA/CodeGen/Mapping.hs
+++ b/Data/Array/Accelerate/CUDA/CodeGen/Mapping.hs
@@ -47,7 +47,7 @@
   = return
   $ CUTranslSkel "map" [cunit|
 
-    $esc:("#include <accelerate_cuda_extras.h>")
+    $esc:("#include <accelerate_cuda.h>")
     $edecls:texIn
 
     extern "C" __global__ void
@@ -57,7 +57,7 @@
         $params:argOut
     )
     {
-        const int shapeSize     = size(shOut);
+        const int shapeSize     = $exp:(csize shOut);
         const int gridSize      = $exp:(gridSize dev);
               int ix;
 
@@ -71,8 +71,8 @@
     }
   |]
   where
-    (texIn, argIn)      = environment dev aenv
-    (argOut, setOut)    = setters "Out" (undefined :: Array sh b)
-    (x, _, _)           = locals "x" (undefined :: a)
-    ix                  = [cvar "ix"]
+    (texIn, argIn)              = environment dev aenv
+    (argOut, shOut, setOut)     = writeArray "Out" (undefined :: Array sh b)
+    (x, _, _)                   = locals "x" (undefined :: a)
+    ix                          = [cvar "ix"]
 
diff --git a/Data/Array/Accelerate/CUDA/CodeGen/Monad.hs b/Data/Array/Accelerate/CUDA/CodeGen/Monad.hs
--- a/Data/Array/Accelerate/CUDA/CodeGen/Monad.hs
+++ b/Data/Array/Accelerate/CUDA/CodeGen/Monad.hs
@@ -20,17 +20,13 @@
 
 import Prelude                                          hiding ( exp )
 import Data.HashSet                                     ( HashSet )
-import Data.HashMap.Strict                              ( HashMap )
 import Data.Hashable
 import Control.Monad
 import Control.Monad.State.Strict
-import Control.Applicative
 import Language.C.Quote.CUDA
 import qualified Language.C                             as C
 import qualified Data.HashSet                           as Set
-import qualified Data.HashMap.Strict                    as Map
 
-import Data.Array.Accelerate.AST
 import Data.Array.Accelerate.Trafo
 import Data.Array.Accelerate.CUDA.CodeGen.Type
 
@@ -40,8 +36,9 @@
 
 -- The state of the code generator monad. The outer monad is used to generate
 -- fresh variable names and collect any headers required for foreign functions.
--- The inner is used to collect local environment bindings when generating code
--- for each individual scalar expression.
+-- The inner is used to collect local bindings (for let-bound expressions) and
+-- to record when previously bound variables are later used for dead-code
+-- analysis.
 --
 -- This separation is required so that names are unique across all generated
 -- code fragments of a skeleton.
@@ -55,10 +52,17 @@
   }
 
 data ExpST = ExpST
-  { bindings    :: [C.BlockItem]
-  , terms       :: !(HashSet C.Exp)
-  , letterms    :: !(HashMap C.Exp C.Exp)
-    -- TODO: this should be a set of reverse dependencies: HashMap C.Exp [C.Exp]
+  -- A stack of (typically) declarations that must be evaluated before
+  -- computation of the main expression. These are typically introduced by
+  -- let-bindings. The list is kept in reverse order, with the last to be
+  -- evaluated (newest addition) at the front of the list.
+  --
+  { localBindings       :: [C.BlockItem]
+
+  -- A set of the Var's that we know have been used in the expression. With this
+  -- we can do def-use analysis for simple dead-code elimination.
+  --
+  , usedTerms           :: !(HashSet C.Exp)
   }
 
 
@@ -69,7 +73,7 @@
 runCUDA a = runState a (AccST 0 Set.empty)
 
 runCGM :: Gen a -> CUDA (a, ExpST)
-runCGM a = runStateT a (ExpST [] Set.empty Map.empty)
+runCGM a = runStateT a (ExpST [] Set.empty)
 
 evalCGM :: Gen a -> CUDA a
 evalCGM = fmap fst . runCGM
@@ -81,53 +85,51 @@
 -- Create new binding points for the C expressions associated with the given AST
 -- term, unless the term is itself a variable.
 --
--- Additionally, add these new terms to a map from the variable name to original
--- binding expression. This will be used as a reverse lookup when marking terms
--- as used.
---
 pushEnv :: DelayedOpenExp env aenv t -> [C.Exp] -> Gen [C.Exp]
 pushEnv exp cs =
-  case exp of
-    Var _       -> return cs
-    Prj _ _     -> return cs
-    _           -> do
-      vs <- zipWithM bind (expType exp) cs
-      modify (\st -> st { letterms = Map.union (Map.fromList (zip vs cs)) (letterms st) })
-      return vs
+  let tys               = expType exp
+      zipWithM' xs ys f = zipWithM f xs ys
+  in
+  zipWithM' tys cs $ \ty c ->
+    case c of
+      C.Var{}   -> return c
+      C.Const{} -> return c
+      _         -> bind ty c
 
+
 -- Return the local environment code, consisting of a list of initialisation
 -- declarations and statements. During construction, these are introduced to the
 -- front of the list, so reverse to get in execution order.
 --
 getEnv :: Gen [C.BlockItem]
-getEnv = reverse <$> gets bindings
+getEnv = reverse `fmap` gets localBindings
 
 -- Generate a fresh variable name
 --
 fresh :: CUDA String
-fresh = do
-  n     <- gets counter <* modify (\st -> st { counter = counter st + 1 })
-  return $ 'v' : show n
+fresh = state $ \s -> let n = counter s in ('v':show n, s { counter = n+1 })
 
 -- Add an expression of given type to the environment and return the (new,
 -- unique) binding name that can be used in place of the thing just bound.
 --
 bind :: C.Type -> C.Exp -> Gen C.Exp
-bind t e = do
-  name <- lift fresh
-  modify (\st -> st { bindings = C.BlockDecl [cdecl| const $ty:t $id:name = $exp:e;|] : bindings st })
-  return [cexp| $id:name |]
+bind t e =
+  case e of
+    C.Var{}     -> return e
+    _           -> do
+      name <- lift fresh
+      modify (\st -> st { localBindings = [citem| const $ty:t $id:name = $exp:e;|] : localBindings st })
+      return [cexp| $id:name |]
 
+
 -- Add an expression to the set marking that it will be used to generate the
 -- output value(s). If the term exists in the reverse let-map, add that binding
 -- instead.
 --
 use :: C.Exp -> Gen C.Exp
-use e = do
-  m <- gets letterms
-  case Map.lookup e m of
-    Nothing     -> modify (\st -> st { terms = Set.insert e (terms st) })
-    Just x      -> modify (\st -> st { terms = Set.insert x (terms st) })
-  --
+use e@(C.Var{}) = do
+  modify (\st -> st { usedTerms = Set.insert e (usedTerms st) })
   return e
+--
+use e           = return e
 
diff --git a/Data/Array/Accelerate/CUDA/CodeGen/PrefixSum.hs b/Data/Array/Accelerate/CUDA/CodeGen/PrefixSum.hs
--- a/Data/Array/Accelerate/CUDA/CodeGen/PrefixSum.hs
+++ b/Data/Array/Accelerate/CUDA/CodeGen/PrefixSum.hs
@@ -170,7 +170,7 @@
 mkScan dir dev aenv fun@(CUFun2 _ _ combine) mseed (CUDelayed (CUExp shIn) _ (CUFun1 _ get)) =
   CUTranslSkel scan [cunit|
 
-    $esc:("#include <accelerate_cuda_extras.h>")
+    $esc:("#include <accelerate_cuda.h>")
     $edecls:texIn
 
     extern "C" __global__ void
@@ -188,7 +188,7 @@
         $decls:declz
         $items:(sh .=. shIn)
 
-        const int shapeSize     = $exp:(shapeSize sh);
+        const int shapeSize     = $exp:(csize sh);
         const int intervalSize  = (shapeSize + gridDim.x - 1) / gridDim.x;
 
         /*
@@ -211,7 +211,7 @@
             ; seg < numElements
             ; seg += blockDim.x )
         {
-            const int ix = $id:(if dir == L then "start + seg" else "end - seg - 1") ;
+            const int ix = $exp:firstIndex;
 
             /*
              * Generate the next set of values
@@ -232,7 +232,7 @@
             $items:(sdata "threadIdx.x" .=. x)
             __syncthreads();
 
-            $stms:(scanBlock dev fun x y sdata Nothing)
+            $items:(scanBlock dev fun x y sdata Nothing)
 
             /*
              * Exclusive scans write the result of the previous thread to global
@@ -258,36 +258,49 @@
                 const int last = min(numElements - seg, blockDim.x) - 1;
                 $items:(z .=. sdata "last")
             }
-            $id:( if isNothing mseed then "carryIn = 1" else [] ) ;
+            $items:setCarry
         }
 
         /*
          * Finally, exclusive scans set the overall scan result (reduction value)
          */
-        if ( $exp:(cbool (isJust mseed)) && threadIdx.x == 0 && blockIdx.x == $id:lastBlock ) {
-            $items:(setSum .=. z)
-        }
+        $items:setFinal
     }
   |]
   where
-    scan                = "scan" ++ show dir ++ maybe "1" (const []) mseed
-    (texIn, argIn)      = environment dev aenv
-    (argOut, setOut)    = setters "Out" (undefined :: Vector e)
-    (argSum, totalSum)  = setters "Sum" (undefined :: Vector e)
-    (argBlk, blkSum)    = setters "Blk" (undefined :: Vector e)
-    (_, x, declx)       = locals "x" (undefined :: e)
-    (_, y, decly)       = locals "y" (undefined :: e)
-    (_, z, declz)       = locals "z" (undefined :: e)
-    (sh, _, _)          = locals "sh" (undefined :: DIM1)
-    (smem, sdata)       = shared (undefined :: e) "sdata" [cexp| blockDim.x |] Nothing
-    ix                  = [cvar "ix"]
-    setSum              = totalSum "0"
+    scan                        = "scan" ++ show dir ++ maybe "1" (const []) mseed
+    (texIn, argIn)              = environment dev aenv
+    (argOut, _, setOut)         = writeArray "Out" (undefined :: Vector e)
+    (argSum, _, totalSum)       = writeArray "Sum" (undefined :: Vector e)
+    (argBlk, _, blkSum)         = writeArray "Blk" (undefined :: Vector e)
+    (_, x, declx)               = locals "x" (undefined :: e)
+    (_, y, decly)               = locals "y" (undefined :: e)
+    (_, z, declz)               = locals "z" (undefined :: e)
+    (sh, _, _)                  = locals "sh" (undefined :: DIM1)
+    (smem, sdata)               = shared (undefined :: e) "sdata" [cexp| blockDim.x |] Nothing
+    ix                          = [cvar "ix"]
+    setSum                      = totalSum "0"
 
+    -- depending on whether we are inclusive/exclusive scans
+    setCarry
+      | isNothing mseed         = [[citem| carryIn = 1; |]]
+      | otherwise               = []
+
+    setFinal
+      | isNothing mseed         = []
+      | otherwise               = [[citem| if ( threadIdx.x == 0 && blockIdx.x == $id:lastBlock ) {
+                                               $items:(setSum .=. z)
+                                           } |]]
+
     -- accessing neighbouring blocks
     firstBlock          = if dir == L then "0" else "gridDim.x - 1"
     lastBlock           = if dir == R then "0" else "gridDim.x - 1"
     prevBlock           = if dir == L then "blockIdx.x - 1" else "blockIdx.x + 1"
 
+    firstIndex
+      | dir == L        = [cexp| start + seg |]
+      | otherwise       = [cexp| end - seg - 1 |]
+
     carryIn
       | isJust mseed    = [cexp| threadIdx.x == 0 |]
       | otherwise       = [cexp| threadIdx.x == 0 && carryIn |]
@@ -327,7 +340,7 @@
 mkScanUp1 dir dev aenv fun@(CUFun2 _ _ combine) (CUDelayed (CUExp shIn) _ (CUFun1 _ get)) =
   CUTranslSkel scan [cunit|
 
-    $esc:("#include <accelerate_cuda_extras.h>")
+    $esc:("#include <accelerate_cuda.h>")
     $edecls:texIn
 
     extern "C" __global__ void
@@ -342,7 +355,7 @@
         $decls:decly
         $items:(sh .=. shIn)
 
-        const int shapeSize     = $exp:(shapeSize sh);
+        const int shapeSize     = $exp:(csize sh);
         const int intervalSize  = (shapeSize + gridDim.x - 1) / gridDim.x;
 
         const int start         = blockIdx.x * intervalSize;
@@ -355,7 +368,7 @@
             ; seg < numElements
             ; seg += blockDim.x )
         {
-            const int ix = $id:(if dir == L then "start + seg" else "end - seg - 1") ;
+            const int ix = $exp:firstIndex ;
 
             /*
              * Read in new values, combine with carry-in
@@ -372,7 +385,7 @@
             $items:(sdata "threadIdx.x" .=. x)
             __syncthreads();
 
-            $stms:(scanBlock dev fun x y sdata Nothing)
+            $items:(scanBlock dev fun x y sdata Nothing)
 
             /*
              * Store the final result of the block to be carried in
@@ -393,16 +406,20 @@
     }
   |]
   where
-    scan                = "scan" ++ show dir ++ "Up"
-    (texIn, argIn)      = environment dev aenv
-    (argOut, setOut)    = setters "Out" (undefined :: Vector e)
-    (_, x, declx)       = locals "x" (undefined :: e)
-    (_, y, decly)       = locals "y" (undefined :: e)
-    (sh, _, _)          = locals "sh" (undefined :: DIM1)
-    (smem, sdata)       = shared (undefined :: e) "sdata" [cexp| blockDim.x |] Nothing
-    ix                  = [cvar "ix"]
+    scan                        = "scan" ++ show dir ++ "Up"
+    (texIn, argIn)              = environment dev aenv
+    (argOut, _, setOut)         = writeArray "Out" (undefined :: Vector e)
+    (_, x, declx)               = locals "x" (undefined :: e)
+    (_, y, decly)               = locals "y" (undefined :: e)
+    (sh, _, _)                  = locals "sh" (undefined :: DIM1)
+    (smem, sdata)               = shared (undefined :: e) "sdata" [cexp| blockDim.x |] Nothing
+    ix                          = [cvar "ix"]
 
+    firstIndex
+      | dir == L                = [cexp| start + seg |]
+      | otherwise               = [cexp| end - seg - 1 |]
 
+
 -- Second step of the upsweep phase: scan the interval sums to produce carry-in
 -- values for each block of the final downsweep step
 --
@@ -415,7 +432,7 @@
     -> Maybe (CUExp aenv e)
     -> CUTranslSkel aenv (Vector e)
 mkScanUp2 dir dev aenv f z
-  = let (_, get) = getters "Blk" (undefined :: Vector e)
+  = let (_, get) = readArray "Blk" (undefined :: Vector e)
     in  mkScan dir dev aenv f z get
 
 
@@ -429,7 +446,7 @@
     -> [C.Exp] -> [C.Exp]
     -> (Name -> [C.Exp])
     -> Maybe C.Exp
-    -> [C.Stm]
+    -> [C.BlockItem]
 scanBlock dev f x0 x1 sdata mlim
   | shflOK dev (undefined :: e) = error "shfl-scan"
   | otherwise                   = scanBlockTree dev f x0 x1 sdata mlim
@@ -447,7 +464,7 @@
     -> [C.Exp] -> [C.Exp]               -- temporary variables x0 and x1
     -> (Name -> [C.Exp])                -- index elements from shared memory
     -> Maybe C.Exp                      -- partially full block bounds check?
-    -> [C.Stm]
+    -> [C.BlockItem]
 scanBlockTree dev (CUFun2 _ _ f) x0 x1 sdata mlim = map (scan . pow2) [ 0 .. maxThreads ]
   where
     pow2 :: Int -> Int
@@ -458,7 +475,7 @@
       | Just m <- mlim  = [cexp| threadIdx.x >= $int:n && threadIdx.x < $exp:m |]
       | otherwise       = [cexp| threadIdx.x >= $int:n |]
 
-    scan n = [cstm|
+    scan n = [citem|
       if ( blockDim.x > $int:n ) {
           if ( $exp:(inrange n) ) {
               $items:(x1 .=. sdata ("threadIdx.x - " ++ show n))
diff --git a/Data/Array/Accelerate/CUDA/CodeGen/Reduction.hs b/Data/Array/Accelerate/CUDA/CodeGen/Reduction.hs
--- a/Data/Array/Accelerate/CUDA/CodeGen/Reduction.hs
+++ b/Data/Array/Accelerate/CUDA/CodeGen/Reduction.hs
@@ -93,7 +93,7 @@
     -> CUDelayedAcc aenv (sh :. Int) e
     -> [ CUTranslSkel aenv (Array sh e) ]
 mkFoldAll dev aenv f z a
-  = let (_, rec) = getters "Rec" (undefined :: Array (sh:.Int) e)
+  = let (_, rec) = readArray "Rec" (undefined :: Array (sh:.Int) e)
     in
     [ mkFoldAll' False dev aenv f z a
     , mkFoldAll' True  dev aenv f z rec ]
@@ -108,10 +108,10 @@
     -> Maybe (CUExp aenv e)
     -> CUDelayedAcc aenv (sh :. Int) e
     -> CUTranslSkel aenv (Array sh e)
-mkFoldAll' recursive dev aenv fun@(CUFun2 _ _ combine) mseed (CUDelayed (CUExp sh) _ (CUFun1 _ get))
+mkFoldAll' recursive dev aenv fun@(CUFun2 _ _ combine) mseed (CUDelayed (CUExp shIn) _ (CUFun1 _ get))
   = CUTranslSkel foldAll [cunit|
 
-    $esc:("#include <accelerate_cuda_extras.h>")
+    $esc:("#include <accelerate_cuda.h>")
     $edecls:texIn
 
     extern "C" __global__ void
@@ -126,8 +126,9 @@
         $decls:declx
         $decls:decly
 
-        $items:(shIn .=. sh)
-        const int shapeSize     = $exp:(shapeSize shIn);
+        $items:(sh .=. shIn)
+
+        const int shapeSize     = $exp:(csize sh);
         const int gridSize      = $exp:(gridSize dev);
               int ix            = $exp:(threadIdx dev);
 
@@ -169,7 +170,7 @@
         $items:(sdata "threadIdx.x" .=. y)
 
         ix = min(shapeSize - blockIdx.x * blockDim.x, blockDim.x);
-        $stms:(reduceBlock dev fun x y sdata (cvar "ix"))
+        $items:(reduceBlock dev fun x y sdata (cvar "ix"))
 
         /*
          * Write the results of this block back to global memory. If we are the last
@@ -182,21 +183,21 @@
     }
   |]
   where
-    foldAll             = maybe "fold1All" (const "foldAll") mseed
-    (texIn, argIn)      = environment dev aenv
-    (argOut, setOut)    = setters "Out" (undefined :: Array sh e)
+    foldAll                     = maybe "fold1All" (const "foldAll") mseed
+    (texIn, argIn)              = environment dev aenv
+    (argOut, _, setOut)         = writeArray "Out" (undefined :: Array (sh :. Int) e)
     (argRec, _)
-      | recursive       = getters "Rec" (undefined :: Array (sh:.Int) e)
-      | otherwise       = ([], undefined)
+      | recursive               = readArray "Rec" (undefined :: Array (sh :. Int) e)
+      | otherwise               = ([], undefined)
 
-    (_, x, declx)       = locals "x" (undefined :: e)
-    (_, y, decly)       = locals "y" (undefined :: e)
-    (shIn, _, _)        = locals "sh" (undefined :: sh :. Int)
-    ix                  = [cvar "ix"]
-    (smem, sdata)       = shared (undefined :: e) "sdata" [cexp| blockDim.x |] Nothing
+    (_, x, declx)               = locals "x" (undefined :: e)
+    (_, y, decly)               = locals "y" (undefined :: e)
+    (sh, _, _)                  = locals "sh" (undefined :: sh :. Int)
+    ix                          = [cvar "ix"]
+    (smem, sdata)               = shared (undefined :: e) "sdata" [cexp| blockDim.x |] Nothing
     --
     inclusive_finish                    = setOut "blockIdx.x" .=. y
-    exclusive_finish (CUExp seed)       = C.BlockStm [cstm|
+    exclusive_finish (CUExp seed)       = [[citem|
       if ( shapeSize > 0 ) {
           if ( gridDim.x == 1 ) {
               $items:(x .=. seed)
@@ -207,7 +208,7 @@
       else {
           $items:(setOut "blockIdx.x" .=. seed)
       }
-    |] : []
+    |]]
 
 
 -- Reduction of the innermost dimension of an array of arbitrary rank. Each
@@ -221,11 +222,11 @@
     -> Maybe (CUExp aenv e)
     -> CUDelayedAcc aenv (sh :. Int) e
     -> [ CUTranslSkel aenv (Array sh e) ]
-mkFoldDim dev aenv fun@(CUFun2 _ _ combine) mseed (CUDelayed (CUExp sh) _ (CUFun1 _ get))
+mkFoldDim dev aenv fun@(CUFun2 _ _ combine) mseed (CUDelayed (CUExp shIn) _ (CUFun1 _ get))
   = return
   $ CUTranslSkel fold [cunit|
 
-    $esc:("#include <accelerate_cuda_extras.h>")
+    $esc:("#include <accelerate_cuda.h>")
     $edecls:texIn
 
     extern "C" __global__ void
@@ -239,10 +240,10 @@
         $decls:declx
         $decls:decly
 
-        $items:(shIn .=. sh)
+        $items:(sh .=. shIn)
 
-        const int numIntervals  = size(shOut);
-        const int intervalSize  = $exp:(indexHead shIn);
+        const int numIntervals  = $exp:(csize (cindexTail sh));
+        const int intervalSize  = $exp:(cindexHead sh);
               int ix;
               int seg;
 
@@ -250,7 +251,7 @@
          * If the intervals of an exclusive fold are empty, use all threads to
          * map the seed value to the output array and exit.
          */
-        $stms:(maybe [] mapseed mseed)
+        $items:(maybe [] mapseed mseed)
 
         /*
          * Threads in a block cooperatively reduce all elements in an interval.
@@ -316,7 +317,7 @@
              * cooperatively reduces this to a single value.
              */
             $items:(sdata "threadIdx.x" .=. y)
-            $stms:(reduceBlock dev fun x y sdata (cvar "n"))
+            $items:(reduceBlock dev fun x y sdata (cvar "n"))
 
             /*
              * Finally, the first thread writes the result for this segment. For
@@ -330,26 +331,27 @@
     }
   |]
   where
-    fold                = maybe "fold1" (const "fold") mseed
-    (texIn, argIn)      = environment dev aenv
-    (argOut, setOut)    = setters "Out" (undefined :: Array sh e)
-    (_, x, declx)       = locals "x" (undefined :: e)
-    (_, y, decly)       = locals "y" (undefined :: e)
-    (shIn, _, _)        = locals "sh" (undefined :: sh :. Int)
-    ix                  = [cvar "ix"]
-    (smem, sdata)       = shared (undefined :: e) "sdata" [cexp| blockDim.x |] Nothing
+    fold                        = maybe "fold1" (const "fold") mseed
+    (texIn, argIn)              = environment dev aenv
+    (argOut, shOut, setOut)     = writeArray "Out" (undefined :: Array sh e)
+    (_, x, declx)               = locals "x" (undefined :: e)
+    (_, y, decly)               = locals "y" (undefined :: e)
+    (sh, _, _)                  = locals "sh" (undefined :: sh :. Int)
+    ix                          = [cvar "ix"]
+    (smem, sdata)               = shared (undefined :: e) "sdata" [cexp| blockDim.x |] Nothing
     --
     mapseed (CUExp seed)
-      = [cstm|  if ( intervalSize == 0 ) {
-                    const int gridSize  = $exp:(gridSize dev);
+      = [citem|  if ( intervalSize == 0 || numIntervals == 0 ) {
+                     const int gridSize  = $exp:(gridSize dev);
 
-                    for ( ix = $exp:(threadIdx dev)
-                        ; ix < numIntervals
-                        ; ix += gridSize )
-                    {
-                        $items:(setOut "ix" .=. seed)
-                    }
-                } |] :[]
+                     for ( ix = $exp:(threadIdx dev)
+                         ; ix < $exp:(csize shOut)
+                         ; ix += gridSize )
+                     {
+                         $items:(setOut "ix" .=. seed)
+                     }
+                     return;
+                 } |] :[]
     --
     exclusive_finish (CUExp seed)
       = concat [ x .=. seed
@@ -421,7 +423,7 @@
   (CUDelayed _            _ (CUFun1 _ offset))
   = CUTranslSkel foldSeg [cunit|
 
-    $esc:("#include <accelerate_cuda_extras.h>")
+    $esc:("#include <accelerate_cuda.h>")
     $edecls:texIn
 
     extern "C"
@@ -439,8 +441,8 @@
         const int thread_lane           = threadIdx.x & (warpSize - 1);
         const int vector_lane           = threadIdx.x / warpSize;
 
-        const int num_segments          = indexHead(shOut);
-        const int total_segments        = size(shOut);
+        const int num_segments          = $exp:(cindexHead shOut);
+        const int total_segments        = $exp:(csize shOut);
               int seg;
               int ix;
 
@@ -459,7 +461,7 @@
             ; seg += num_vectors )
         {
             const int s    =  seg % num_segments;
-            const int base = (seg / num_segments) * $exp:(indexHead sh);
+            const int base = (seg / num_segments) * $exp:(cindexHead sh);
 
             /*
              * Use two threads to fetch the indices of the start and end of this
@@ -523,7 +525,7 @@
              */
             ix = min(num_elements, warpSize);
             $items:(sdata "threadIdx.x" .=. y)
-            $stms:(reduceWarp dev fun x y sdata (cvar "ix") (cvar "thread_lane"))
+            $items:(reduceWarp dev fun x y sdata (cvar "ix") (cvar "thread_lane"))
 
             /*
              * Finally, the first thread writes the result for this segment
@@ -537,25 +539,25 @@
     }
   |]
   where
-    foldSeg             = maybe "fold1Seg" (const "foldSeg") mseed
-    (texIn, argIn)      = environment dev aenv
-    (argOut, setOut)    = setters "Out" (undefined :: Array (sh :. Int) e)
-    (_, x, declx)       = locals "x" (undefined :: e)
-    (_, y, decly)       = locals "y" (undefined :: e)
-    (sh, _, _)          = locals "sh" (undefined :: sh :. Int)
-    (smem, sdata)       = shared (undefined :: e) "sdata" [cexp| blockDim.x |] (Just $ [cexp| &s_ptrs[vectors_per_block][2] |])
+    foldSeg                     = maybe "fold1Seg" (const "foldSeg") mseed
+    (texIn, argIn)              = environment dev aenv
+    (argOut, shOut, setOut)     = writeArray "Out" (undefined :: Array (sh :. Int) e)
+    (_, x, declx)               = locals "x" (undefined :: e)
+    (_, y, decly)               = locals "y" (undefined :: e)
+    (sh, _, _)                  = locals "sh" (undefined :: sh :. Int)
+    (smem, sdata)               = shared (undefined :: e) "sdata" [cexp| blockDim.x |] (Just $ [cexp| &s_ptrs[vectors_per_block][2] |])
     --
-    ix                  = [cvar "ix"]
-    vectors_per_block   = cvar "vectors_per_block"
-    gridDim             = cvar "gridDim.x"
+    ix                          = [cvar "ix"]
+    vectors_per_block           = cvar "vectors_per_block"
+    gridDim                     = cvar "gridDim.x"
     --
     exclusive_finish (CUExp seed)
-      = C.BlockStm [cstm| if ( num_elements > 0 ) {
-                              $items:(x .=. seed)
-                              $items:(y .=. combine x y)
-                          } else {
-                              $items:(y .=. seed)
-                          } |] :[]
+      = [[citem| if ( num_elements > 0 ) {
+                     $items:(x .=. seed)
+                     $items:(y .=. combine x y)
+                 } else {
+                     $items:(y .=. seed)
+                 } |]]
 
 
 -- Reducers
@@ -574,7 +576,7 @@
     -> (Name -> [C.Exp])                -- index elements from shared memory
     -> C.Exp                            -- number of elements
     -> C.Exp                            -- thread identifier: usually lane or thread ID
-    -> [C.Stm]
+    -> [C.BlockItem]
 reduceWarp dev fun x0 x1 sdata n tid
   | shflOK dev (undefined :: e) = return
                                 $ reduceWarpShfl dev fun x0 x1       n tid
@@ -588,7 +590,7 @@
     -> [C.Exp] -> [C.Exp]               -- temporary variables x0 and x1
     -> (Name -> [C.Exp])                -- index elements from shared memory
     -> C.Exp                            -- number of elements
-    -> [C.Stm]
+    -> [C.BlockItem]
 reduceBlock dev fun x0 x1 sdata n
   | shflOK dev (undefined :: e) = reduceBlockShfl dev fun x0 x1 sdata n
   | otherwise                   = reduceBlockTree dev fun x0 x1 sdata n
@@ -605,7 +607,7 @@
     -> (Name -> [C.Exp])                -- index elements from shared memory
     -> C.Exp                            -- number of elements
     -> C.Exp                            -- thread identifier: usually lane or thread ID
-    -> [C.Stm]
+    -> [C.BlockItem]
 reduceWarpTree dev (CUFun2 _ _ f) x0 x1 sdata n tid
   = map (reduce . pow2) [v, v-1 .. 0]
   where
@@ -614,18 +616,18 @@
     pow2 :: Int -> Int
     pow2 x = 2 ^ x
 
-    reduce :: Int -> C.Stm
+    reduce :: Int -> C.BlockItem
     reduce 0
-      = [cstm| if ( $exp:tid < $exp:n ) {
-                   $items:(x0 .=. sdata "threadIdx.x + 1")
-                   $items:(x1 .=. f x1 x0)
-               } |]
+      = [citem| if ( $exp:tid < $exp:n ) {
+                    $items:(x0 .=. sdata "threadIdx.x + 1")
+                    $items:(x1 .=. f x1 x0)
+                } |]
     reduce i
-      = [cstm| if ( $exp:tid + $int:i < $exp:n ) {
-                   $items:(x0 .=. sdata ("threadIdx.x + " ++ show i))
-                   $items:(x1 .=. f x1 x0)
-                   $items:(sdata "threadIdx.x" .=. x1)
-               } |]
+      = [citem| if ( $exp:tid + $int:i < $exp:n ) {
+                    $items:(x0 .=. sdata ("threadIdx.x + " ++ show i))
+                    $items:(x1 .=. f x1 x0)
+                    $items:(sdata "threadIdx.x" .=. x1)
+                } |]
 
 reduceBlockTree
     :: Elt e
@@ -634,7 +636,7 @@
     -> [C.Exp] -> [C.Exp]               -- temporary variables x0 and x1
     -> (Name -> [C.Exp])                -- index elements from shared memory
     -> C.Exp                            -- number of elements
-    -> [C.Stm]
+    -> [C.BlockItem]
 reduceBlockTree dev fun@(CUFun2 _ _ f) x0 x1 sdata n
   = flip (foldr1 (.)) []
   $ map (reduce . pow2) [u-1, u-2 .. v]
@@ -646,32 +648,37 @@
     pow2 :: Int -> Int
     pow2 x = 2 ^ x
 
-    reduce :: Int -> [C.Stm] -> [C.Stm]
+    reduce :: Int -> [C.BlockItem] -> [C.BlockItem]
     reduce i rest
       -- Ensure that threads synchronise before reading from or writing to
       -- shared memory. Synchronising after each reduction step is not enough,
       -- because one warp could update the partial results before a different
       -- warp has read in their data for this step.
       --
+      -- Additionally, note that all threads of a warp must participate in the
+      -- synchronisation. Thus, this must go outside of the test against the
+      -- bounds of the array. We do a bit of extra work here, with all threads
+      -- writing into shared memory whether they updated their value or not.
+      --
       | i > warpSize dev
-      = [cstm| if ( threadIdx.x + $int:i < $exp:n ) {
-                   __syncthreads();
-                   $items:(x0 .=. sdata ("threadIdx.x + " ++ show i))
-                   $items:(x1 .=. f x1 x0)
-                   __syncthreads();
-                   $items:(sdata "threadIdx.x" .=. x1)
-               } |]
-      : rest
+      = [citem| __syncthreads(); |]
+      : [citem| if ( threadIdx.x + $int:i < $exp:n ) {
+                    $items:(x0 .=. sdata ("threadIdx.x + " ++ show i))
+                    $items:(x1 .=. f x1 x0)
+                } |]
+      : [citem| __syncthreads(); |]
+      : (sdata "threadIdx.x" .=. x1)
+      ++ rest
 
       -- The threads of a warp execute in lockstep, so it is only necessary to
       -- synchronise at the top, to ensure all threads have written their
       -- results into shared memory.
       --
       | otherwise
-      = [cstm| __syncthreads(); |]
-      : [cstm| if ( threadIdx.x < $int:(warpSize dev) ) {
-                   $stms:(reduceWarpTree dev fun x0 x1 sdata n (cvar "threadIdx.x"))
-               } |]
+      = [citem| __syncthreads(); |]
+      : [citem| if ( threadIdx.x < $int:(warpSize dev) ) {
+                    $items:(reduceWarpTree dev fun x0 x1 sdata n (cvar "threadIdx.x"))
+                } |]
       : rest
 
 
@@ -696,15 +703,15 @@
     -> [C.Exp] -> [C.Exp]
     -> C.Exp
     -> C.Exp
-    -> C.Stm
+    -> C.BlockItem
 reduceWarpShfl _dev (CUFun2 _ _ f) x0 x1 n tid
-  = [cstm| for ( int z = warpSize/2; z >= 1; z /= 2 ) {
-               $items:(x0 .=. shfl_xor x1)
+  = [citem| for ( int z = warpSize/2; z >= 1; z /= 2 ) {
+                $items:(x0 .=. shfl_xor x1)
 
-               if ( $exp:tid + z < $exp:n ) {
-                   $items:(x1 .=. f x1 x0)
-               }
-           } |]
+                if ( $exp:tid + z < $exp:n ) {
+                    $items:(x1 .=. f x1 x0)
+                }
+            } |]
   where
     sizeof      = eltSizeOf (undefined :: e)
     shfl_xor    = zipWith (\s x -> ccall (shfl s) [ x, cvar "z" ]) sizeof
@@ -726,17 +733,17 @@
     -> [C.Exp] -> [C.Exp]
     -> (Name -> [C.Exp])
     -> C.Exp
-    -> [C.Stm]
+    -> [C.BlockItem]
 reduceBlockShfl dev fun x0 x1 sdata n
   = reduceWarpShfl dev fun x0 x1 n (cvar "threadIdx.x")
-  : [cstm|  if ( (threadIdx.x & warpSize - 1) == 0 ) {
-                $items:(sdata "threadIdx.x / warpSize" .=. x1)
-            } |]
-  : [cstm|  __syncthreads(); |]
-  : [cstm|  if ( threadIdx.x < warpSize ) {
-                $items:(x1 .=. sdata "threadIdx.x")
-                $exp:n = ($exp:n + warpSize - 1) / warpSize;
-                $stm:(reduceWarpShfl dev fun x0 x1 n (cvar "threadIdx.x"))
-            } |]
+  : [citem|  if ( (threadIdx.x & warpSize - 1) == 0 ) {
+                 $items:(sdata "threadIdx.x / warpSize" .=. x1)
+             } |]
+  : [citem|  __syncthreads(); |]
+  : [citem|  if ( threadIdx.x < warpSize ) {
+                 $items:(x1 .=. sdata "threadIdx.x")
+                 $exp:n = ($exp:n + warpSize - 1) / warpSize;
+                 $item:(reduceWarpShfl dev fun x0 x1 n (cvar "threadIdx.x"))
+             } |]
   : []
 
diff --git a/Data/Array/Accelerate/CUDA/CodeGen/Stencil.hs b/Data/Array/Accelerate/CUDA/CodeGen/Stencil.hs
--- a/Data/Array/Accelerate/CUDA/CodeGen/Stencil.hs
+++ b/Data/Array/Accelerate/CUDA/CodeGen/Stencil.hs
@@ -19,19 +19,15 @@
 
 ) where
 
-import Control.Applicative
-import Control.Monad.State.Strict
 import Foreign.CUDA.Analysis
 import Language.C.Quote.CUDA
-import qualified Language.C.Syntax                      as C
 
 import Data.Array.Accelerate.Type                       ( Boundary(..) )
-import Data.Array.Accelerate.Array.Sugar                ( Array, Shape, Elt, shapeToList )
-import Data.Array.Accelerate.Analysis.Shape
+import Data.Array.Accelerate.Array.Sugar                ( Array, Elt )
 import Data.Array.Accelerate.Analysis.Stencil
 import Data.Array.Accelerate.CUDA.AST                   hiding ( stencil, stencilAccess )
 import Data.Array.Accelerate.CUDA.CodeGen.Base
-import Data.Array.Accelerate.CUDA.CodeGen.Type
+import Data.Array.Accelerate.CUDA.CodeGen.Stencil.Extra
 
 
 -- Map a stencil over an array.  In contrast to 'map', the domain of a stencil
@@ -66,8 +62,7 @@
   = return
   $ CUTranslSkel "stencil" [cunit|
 
-    $esc:("#include <accelerate_cuda_extras.h>")
-    $edecl:(cdim "Shape" dim)
+    $esc:("#include <accelerate_cuda.h>")
     $edecls:texIn
     $edecls:texStencil
 
@@ -79,7 +74,7 @@
         $params:argStencil
     )
     {
-        const int shapeSize     = size(shOut);
+        const int shapeSize     = $exp:(csize shOut);
         const int gridSize      = $exp:(gridSize dev);
               int ix;
 
@@ -87,24 +82,38 @@
             ; ix <  shapeSize
             ; ix += gridSize )
         {
-            const typename Shape sh = fromIndex( shOut, ix );
-            $items:(dce xs      .=. stencil sh)
-            $items:(setOut "ix" .=. f xs)
+            $items:(sh .=. cfromIndex shOut "ix" "tmp")
+            $items:stencilBody
         }
     }
   |]
   where
-    dim                 = expDim (undefined :: Exp aenv sh)
-    (texIn,  argIn)     = environment dev aenv
-    (argOut, setOut)    = setters "Out" (undefined :: Array sh b)
-    ix                  = cvar "ix"
-    sh                  = cvar "sh"
-    (xs,_,_)            = locals "x" (undefined :: stencil)
-    dx                  = offsets (undefined :: Fun aenv (stencil -> b)) (undefined :: OpenAcc aenv (Array sh a))
+    (texIn,  argIn)             = environment dev aenv
+    (argOut, shOut, setOut)     = writeArray "Out" (undefined :: Array sh b)
+    (sh, _, _)                  = locals "sh" (undefined :: sh)
+    (xs,_,_)                    = locals "x" (undefined :: stencil)
 
-    (texStencil, argStencil, stencil) = stencilAccess True "Stencil" "w" dev dx ix boundary dce
+    dx  = offsets (undefined :: Fun aenv (stencil -> b))
+                  (undefined :: OpenAcc aenv (Array sh a))
 
+    (texStencil, argStencil, safeIndex)   = stencilAccess dev True True  "Stencil" "Stencil" "w" "ix" dx boundary dce
+    (_,          _,          unsafeIndex) = stencilAccess dev True False "Stencil" "Stencil" "w" "ix" dx boundary dce
 
+    stencilBody
+      | computeCapability dev < Compute 1 2     = with safeIndex
+      | otherwise                               =
+          [[citem| if ( __all( $exp:(insideRegion shOut (borderRegion dx) (map rvalue sh)) ) ) {
+                       $items:(with unsafeIndex)
+                   } else {
+                       $items:(with safeIndex)
+                   } |]]
+
+      where
+        with stencil = (dce xs      .=. stencil sh) ++
+                       (setOut "ix" .=. f xs)
+
+
+
 -- Map a binary stencil of an array.  The extent of the resulting array is the
 -- intersection of the extents of the two source arrays.
 --
@@ -127,12 +136,25 @@
     -> Boundary (CUExp aenv a)
     -> Boundary (CUExp aenv b)
     -> [CUTranslSkel aenv (Array sh c)]
-mkStencil2 dev aenv (CUFun2 dce1 dce2 f) boundary1 boundary2
-  = return
-  $ CUTranslSkel "stencil2" [cunit|
+mkStencil2 dev aenv stencil boundary1 boundary2
+  = [ mkStencil2' dev False aenv stencil boundary1 boundary2
+    , mkStencil2' dev True  aenv stencil boundary1 boundary2
+    ]
 
-    $esc:("#include <accelerate_cuda_extras.h>")
-    $edecl:(cdim "Shape" dim)
+mkStencil2'
+    :: forall aenv sh stencil1 stencil2 a b c.
+       (Stencil sh a stencil1, Stencil sh b stencil2, Elt c)
+    => DeviceProperties
+    -> Bool                                     -- are the source arrays the same extent?
+    -> Gamma aenv
+    -> CUFun2 aenv (stencil1 -> stencil2 -> c)
+    -> Boundary (CUExp aenv a)
+    -> Boundary (CUExp aenv b)
+    -> CUTranslSkel aenv (Array sh c)
+mkStencil2' dev sameExtent aenv (CUFun2 dce1 dce2 f) boundary1 boundary2
+  = CUTranslSkel "stencil2" [cunit|
+
+    $esc:("#include <accelerate_cuda.h>")
     $edecls:texIn
     $edecls:texS1
     $edecls:texS2
@@ -146,7 +168,7 @@
         $params:argS2
     )
     {
-        const int shapeSize     = size(shOut);
+        const int shapeSize     = $exp:(csize shOut);
         const int gridSize      = $exp:(gridSize dev);
               int ix;
 
@@ -154,135 +176,53 @@
             ; ix <  shapeSize
             ; ix += gridSize )
         {
-            const typename Shape sh = fromIndex( shOut, ix );
-
-            $items:(dce1 xs     .=. stencil1 sh)
-            $items:(dce2 ys     .=. stencil2 sh)
-            $items:(setOut "ix" .=. f xs ys)
+            $items:(sh        .=. cfromIndex shOut "ix" "tmp")
+            $items:stencilBody
         }
     }
   |]
   where
-    dim                 = expDim (undefined :: Exp aenv sh)
-    (texIn,  argIn)     = environment dev aenv
-    (argOut, setOut)    = setters "Out" (undefined :: Array sh c)
-    ix                  = cvar "ix"
-    sh                  = cvar "sh"
-    (xs,_,_)            = locals "x" (undefined :: stencil1)
-    (ys,_,_)            = locals "y" (undefined :: stencil2)
-
-    (dx1, dx2)          = offsets2 (undefined :: Fun aenv (stencil1 -> stencil2 -> c))
-                                   (undefined :: OpenAcc aenv (Array sh a))
-                                   (undefined :: OpenAcc aenv (Array sh b))
-
-    (texS1, argS1, stencil1) = stencilAccess False "Stencil1" "w" dev dx1 ix boundary1 dce1
-    (texS2, argS2, stencil2) = stencilAccess False "Stencil2" "z" dev dx2 ix boundary2 dce2
-
-
--- Generate declarations for reading in a stencil pattern surrounding a given
--- focal point. The first parameter determines whether it is safe to use linear
--- indexing at the centroid position. This is true for:
---
---  * stencil1
---  * stencil2 if both input stencil have the same dimensionality
---
-stencilAccess
-    :: forall aenv sh e. (Shape sh, Elt e)
-    => Bool                                     -- linear indexing at centroid?
-    -> Name                                     -- array group name
-    -> Name                                     -- secondary group name, for fresh variables
-    -> DeviceProperties                         -- properties of currently executing device
-    -> [sh]                                     -- list of offset indices
-    -> C.Exp                                    -- linear index of the centroid
-    -> Boundary (CUExp aenv e)                  -- stencil boundary condition
-    -> ([C.Exp] -> [(Bool,C.Exp)])              -- dead code elimination flags for this var
-    -> ( [C.Definition]                         -- input arrays as texture references; or
-       , [C.Param]                              -- function arguments
-       , (C.Exp -> ([C.BlockItem], [C.Exp])) )  -- read data at a given shape centroid
-stencilAccess linear grp grp' dev shx centroid boundary dce
-  = (texStencil, argStencil, stencil)
-  where
-    stencil ix = flip evalState 0 $ do
-      (envs, xs) <- mapAndUnzipM (access ix . shapeToList) shx
-
-      let (envs', xs') = unzip
-                       $ eliminate
-                       $ zip envs
-                       $ unconcat (map length xs)
-                       $ dce (concat xs)
-
-      return ( concat envs', concat xs' )
-
-    -- Filter unused components of the stencil. Environment bindings are shared
-    -- between tuple components of each cursor position, so filter these out
-    -- only if all elements of that position are unused.
-    --
-    unconcat :: [Int] -> [a] -> [[a]]
-    unconcat []     _  = []
-    unconcat (n:ns) xs = let (h,t) = splitAt n xs in h : unconcat ns t
-
-    eliminate :: [ ([C.BlockItem], [(Bool, C.Exp)]) ] -> [ ([C.BlockItem], [C.Exp]) ]
-    eliminate []         = []
-    eliminate ((e,v):xs) = (e', x) : eliminate xs
-      where
-        (flags, x)      = unzip v
-        e' | or flags   = e
-           | otherwise  = []
+    (texIn,  argIn)             = environment dev aenv
+    (argOut, shOut, setOut)     = writeArray "Out" (undefined :: Array sh c)
+    (sh, _, _)                  = locals "sh" (undefined :: sh)
+    (xs,_,_)                    = locals "x" (undefined :: stencil1)
+    (ys,_,_)                    = locals "y" (undefined :: stencil2)
+    grp1                        = "Stencil1"
+    grp2                        = "Stencil2"
 
-    -- Generate the entire stencil, including any local environment bindings
+    -- If the source arrays are the same extent, twiddle the names a bit so that
+    -- code generation refers to the same set of shape variables. Then, if there
+    -- are any duplicate calculations, hope that the CUDA compiler is smart
+    -- enough and spots this.
     --
-    access :: C.Exp -> [Int] -> State Int ([C.BlockItem], [C.Exp])
-    access ix dx = case boundary of
-      Clamp                     -> bounded "clamp"
-      Mirror                    -> bounded "mirror"
-      Wrap                      -> bounded "wrap"
-      Constant (CUExp (_,c))    -> inrange c            -- constant value: no environment possible
+    sh1                         = grp1
+    sh2 | sameExtent            = sh1
+        | otherwise             = grp2
 
-      where
-        focus                   = all (==0) dx
-        dim                     = expDim (undefined :: Exp aenv sh)
-        cursor
-          | all (==0) dx        = ix
-          | otherwise           = ccall "shape"
-                                $ zipWith (\a b -> [cexp| $exp:a + $int:b |]) (cshape dim ix) (reverse dx)
+    (dx1, dx2)  = offsets2 (undefined :: Fun aenv (stencil1 -> stencil2 -> c))
+                           (undefined :: OpenAcc aenv (Array sh a))
+                           (undefined :: OpenAcc aenv (Array sh b))
 
-        bounded f
-          | focus && linear     = return $ ( [], getStencil centroid )
-          | otherwise           = do
-              j <- fresh
-              return ( if focus then [C.BlockDecl [cdecl| const int $id:j = toIndex( $id:shIn, $exp:ix ); |]]
-                                else [C.BlockDecl [cdecl| const int $id:j = toIndex( $id:shIn, $exp:(ccall f [cvar shIn, cursor]) ); |]]
-                     , getStencil (cvar j) )
+    border      = zipWith max (borderRegion dx1) (borderRegion dx2)
 
-        inrange cs
-          | focus && linear     = return ( [], getStencil centroid )
-          | focus               = do
-              j <- fresh
-              return ( [C.BlockDecl [cdecl| const int $id:j = toIndex( $id:shIn, $exp:ix ); |]]
-                     , getStencil (cvar j) )
+    (texS1, argS1, safeIndex1)   = stencilAccess dev sameExtent True  grp1 sh1 "w" "ix" dx1 boundary1 dce1
+    (_, _,         unsafeIndex1) = stencilAccess dev sameExtent False grp1 sh1 "w" "ix" dx1 boundary1 dce1
 
-          | otherwise           = do
-              j     <- fresh
-              i     <- fresh
-              p     <- fresh
-              return $ ( [ C.BlockDecl [cdecl| const typename Shape $id:j = $exp:cursor; |]
-                         , C.BlockDecl [cdecl| const typename bool  $id:p = inRange( $id:shIn, $id:j ); |]
-                         , C.BlockDecl [cdecl| const int            $id:i = toIndex( $id:shIn, $id:j ); |] ]
-                       , zipWith (\a c -> [cexp| $id:p ? $exp:a : $exp:c |]) (getStencil (cvar i)) cs )
+    (texS2, argS2, safeIndex2)   = stencilAccess dev sameExtent True  grp2 sh2 "z" "ix" dx2 boundary2 dce2
+    (_, _,         unsafeIndex2) = stencilAccess dev sameExtent False grp2 sh2 "z" "ix" dx2 boundary2 dce2
 
-    -- Extra parameters for accessing the stencil data. We are doing things a
-    -- little out of the ordinary, so don't get this "for free". sadface.
-    --
-    getStencil ix       = zipWith (\t a -> indexArray dev t a ix) (eltType (undefined :: e)) (map cvar stencilIn)
-    (shIn, stencilIn)   = namesOfArray grp (undefined :: e)
-    (texStencil, argStencil)
-      | computeCapability dev < Compute 2 0 = let (d,p) = arrayAsTex (undefined :: Array sh e) grp in (d,[p])
-      | otherwise                           = ([], arrayAsArg (undefined :: Array sh e) grp)
+    stencilBody
+      | computeCapability dev < Compute 1 2     = with safeIndex1 safeIndex2
+      | otherwise                               =
+          [[citem| if ( __all( $exp:(insideRegion shOut border (map rvalue sh)) ) ) {
+                       $items:(with unsafeIndex1 unsafeIndex2)
+                   } else {
+                       $items:(with safeIndex1 safeIndex2)
+                   } |]]
 
-    -- Generate a fresh variable name
-    --
-    fresh :: State Int Name
-    fresh = do
-      n <- get <* modify (+1)
-      return $ grp' ++ show n
+      where
+        with stencil1 stencil2 =
+          (dce1 xs     .=. stencil1 sh) ++
+          (dce2 ys     .=. stencil2 sh) ++
+          (setOut "ix" .=. f xs ys)
 
diff --git a/Data/Array/Accelerate/CUDA/CodeGen/Stencil/Extra.hs b/Data/Array/Accelerate/CUDA/CodeGen/Stencil/Extra.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Accelerate/CUDA/CodeGen/Stencil/Extra.hs
@@ -0,0 +1,280 @@
+{-# LANGUAGE CPP                 #-}
+{-# LANGUAGE ImpredicativeTypes  #-}
+{-# LANGUAGE QuasiQuotes         #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# LANGUAGE ViewPatterns        #-}
+-- |
+-- Module      : Data.Array.Accelerate.CUDA.CodeGen.Stencil.Extra
+-- Copyright   : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee
+--               [2009..2013] 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.CUDA.CodeGen.Stencil.Extra (
+
+  stencilAccess,
+
+  cinRange, cclamp, cmirror, cwrap, insideRegion, borderRegion,
+
+) where
+
+-- standard library
+import Prelude                                          hiding ( and, zipWith, zipWith3 )
+import Data.List                                        ( transpose )
+import Control.Monad
+import Control.Monad.State.Strict
+
+-- language-c
+import Language.C.Quote.CUDA
+import qualified Language.C.Syntax                      as C
+
+-- friends
+import Data.Array.Accelerate.Type                       ( Boundary(..) )
+import Data.Array.Accelerate.Array.Sugar                ( Array, Shape, Elt, shapeToList )
+import Data.Array.Accelerate.Analysis.Shape
+
+import Foreign.CUDA.Analysis
+import Data.Array.Accelerate.CUDA.AST                   hiding ( stencil, stencilAccess )
+import Data.Array.Accelerate.CUDA.CodeGen.Base
+import Data.Array.Accelerate.CUDA.CodeGen.Type
+
+#include "accelerate.h"
+
+
+-- Stencil Access
+-- --------------
+
+-- Generate declarations for reading in a stencil pattern surrounding a given
+-- focal point.
+--
+stencilAccess
+    :: forall aenv sh e. (Shape sh, Elt e)
+    => DeviceProperties
+    -> Bool                             -- can we use linear indexing?
+    -> Bool                             -- do we need to do bounds checking?
+    -> Name                             -- array group name
+    -> Name                             -- group name for array shape (hax!)
+    -> Name                             -- seed name for temporary variables
+    -> Name                             -- linear index at the focus
+    -> [sh]                             -- list of offset indices
+    -> Boundary (CUExp aenv e)          -- stencil boundary condition
+    -> Eliminate e                      -- dead code elimination flags
+    -> ( [C.Definition]                 -- kernel texture reference definitions
+       , [C.Param]                      -- kernel function arguments
+       , Instantiate1 sh e )            -- access stencil at given multidimensional index
+stencilAccess dev doLinearIndexing doBoundsChecks grp sh tmp centroid positions boundary dce =
+  ( decls, params, stencil )
+  where
+    (decls, params, _, getIn)   = readStencil dev grp (undefined :: Array sh e)
+    (_, _, shIn, _)             = readStencil dev sh  (undefined :: Array sh e)
+
+    getInAt ix                  = fresh >>= \j ->
+                                  return ( [[citem| const $ty:cint $id:j = $exp:ix; |]], getIn j )
+
+    -- Generate the entire stencil, reading elements from those positions of the
+    -- pattern that are used and eliminating reads from those that are not.
+    --
+    stencil ix = withNameGen tmp $ do
+      (envs, xs) <- mapAndUnzipM (access ix . shapeToList) positions
+
+      let (envs', xs') = unzip
+                       $ eliminate
+                       $ zipWith (,) envs               -- our version of zipwith that checks lengths
+                       $ unconcat (map length xs)
+                       $ dce (concat xs)
+
+      return ( concat envs', concat xs' )
+
+    -- Read the stencil component at the given offset (second argument). This
+    -- may generate additional environment terms, such as for the index
+    -- calculations.
+    --
+    access :: Rvalue x => [x] -> [Int] -> Gen ([C.BlockItem], [C.Exp])
+    access (map rvalue -> ix) dx
+      | doBoundsChecks          = safeAccess
+      | otherwise               = unsafeAccess
+      where
+        focus                   = all (==0) dx
+
+        -- The current stencil position into the array, as a multidimensional index
+        cursor | focus          = ix
+               | otherwise      = zipWith (\i d -> [cexp| $exp:i + $int:d |]) ix (reverse dx)
+
+        -- Read the array position without any bounds checks
+        unsafeAccess
+          | doLinearIndexing && focus   = return $ ([], getIn centroid)
+          | otherwise                   = getInAt (ctoIndex shIn cursor)
+
+        -- Read the array, applying appropriate bounds checks
+        safeAccess = case boundary of
+          Clamp                  -> bounded cclamp
+          Mirror                 -> bounded cmirror
+          Wrap                   -> bounded cwrap
+          Constant (CUExp (_,c)) -> inrange c
+
+        bounded f
+          | focus               = unsafeAccess
+          | otherwise           = getInAt (ctoIndex shIn (f shIn cursor))
+
+        inrange cs
+          | focus               = unsafeAccess
+          | otherwise           = do
+              (env, as) <- unsafeAccess
+              p         <- fresh
+              return ( [citem| const int $id:p = $exp:(cinRange shIn cursor); |] : env
+                     , zipWith (\a c -> [cexp| $id:p ? $exp:a : $exp:c |]) as cs )
+
+
+-- Filter unused components of the stencil. Environment bindings are shared
+-- between tuple components of each cursor position, so filter these out only if
+-- all elements of that position are unused.
+--
+eliminate :: [ ([a], [(Bool,b)]) ] -> [ ([a],[b]) ]
+eliminate []         = []
+eliminate ((e,v):xs) = (e', x) : eliminate xs
+  where
+    (flags, x)          = unzip v
+    e' | or flags       = e
+       | otherwise      = []
+
+
+-- A simple fresh name supply
+--
+type Gen = State (Name,Int)
+
+withNameGen :: Name -> Gen a -> a
+withNameGen base f = evalState f (base,0)
+
+fresh :: Gen Name
+fresh = state $ \(base,n) -> (base ++ show n, (base,n+1))
+
+
+-- Boundary conditions
+-- -------------------
+
+-- Test whether the given multidimensional index lies in the inside region of
+-- the stencil.
+--
+insideRegion
+    :: [C.Exp]                  -- The shape of the array
+    -> [Int]                    -- The width of the stencil in each direction
+    -> [C.Exp]                  -- The index in question
+    -> C.Exp
+insideRegion shape border index = foldl1 and (zipWith3 inside shape border index)
+  where
+    inside sz dx i      = [cexp| $exp:i >= $int:dx && $exp:i < $exp:sz - $int:dx |]
+    and x y             = [cexp| $exp:x && $exp:y |]
+
+
+-- Given a list of stencil offset positions, calculate the size of the border
+-- region along each dimension.
+--
+-- Note that this does not consider any positions of the stencil that are not
+-- actually used. We assume the user is sensible and uses the minimally sized
+-- stencil for their application, but this can still be problematic for
+-- non-symmetric stencils. For example, a large stencil that uses elements from
+-- only one quadrant.
+--
+borderRegion :: Shape sh => [sh] -> [Int]
+borderRegion
+  = reverse
+  . map maximum
+  . transpose
+  . map shapeToList
+
+
+-- Test whether an index lies within the boundaries of a shape (first argument)
+--
+cinRange :: [C.Exp] -> [C.Exp] -> C.Exp
+cinRange []    []    = INTERNAL_ERROR(error) "inRange" "singleton index"
+cinRange shape index = foldl1 and (zipWith inside shape index)
+  where
+    inside sz i = [cexp| ({ const $ty:cint _i = $exp:i; _i >= 0 && _i < $exp:sz; }) |]
+    and x y     = [cexp| $exp:x && $exp:y |]
+
+-- Clamp an index to the boundary of the shape (first argument)
+--
+cclamp :: [C.Exp] -> [C.Exp] -> [C.Exp]
+cclamp = zipWith f
+  where
+    f sz i = [cexp| max(($ty:cint) 0, min( $exp:i, $exp:sz - 1 )) |]
+
+-- Indices out of bounds of the shape are mirrored back in range. Assumes that
+-- the array is at least as large as the stencil.
+--
+cmirror :: [C.Exp] -> [C.Exp] -> [C.Exp]
+cmirror = zipWith f
+  where
+    f sz i = [cexp| ({ const $ty:cint _i  = $exp:i;
+                       const $ty:cint _sz = $exp:sz;
+                      _i < 0    ? -_i
+                    : _i >= _sz ?  _sz - (_i - _sz + 2)
+                    : _i; }) |]
+
+-- Indices out of bounds are wrapped to the opposite edge of the shape
+--
+cwrap :: [C.Exp] -> [C.Exp] -> [C.Exp]
+cwrap = zipWith f
+  where
+    f sz i = [cexp| ({ const $ty:cint _i  = $exp:i;
+                       const $ty:cint _sz = $exp:sz;
+                    _i < 0    ? _sz + _i
+                  : _i >= _sz ? _i  - _sz
+                  : _i; }) |]
+
+
+-- Kernel parameters
+-- -----------------
+
+-- Generate kernel parameters for input arrays. This is similar to 'readArray',
+-- but we force compute 1.x devices to read through the texture cache as well.
+--
+readStencil
+    :: forall sh e. (Shape sh, Elt e)
+    => DeviceProperties
+    -> Name                                     -- group names
+    -> Array sh e                               -- dummy to fix the types
+    -> ( [C.Definition]                         -- global definitions for stencils read via texture references (compute < 2.0)
+       , [C.Param]                              -- function arguments for stencils read as arrays (compute >= 2.0)
+       , [C.Exp]                                -- shape of the array
+       , forall x. Rvalue x => x -> [C.Exp]     -- read elements from a linear index
+       )
+readStencil dev grp dummy
+  = let (sh, arrs)      = namesOfArray grp (undefined :: e)
+        (decl, args)
+          | computeCapability dev < Compute 2 0 = arrayAsTex dummy grp
+          | otherwise                           = ([], arrayAsArg dummy grp)
+
+        dim             = expDim (undefined :: Exp aenv sh)
+        sh'             = cshape dim sh
+        fetch ix        = zipWith (\t a -> indexArray dev t (cvar a) (rvalue ix)) (eltType (undefined :: e)) arrs
+    in
+    ( decl, args, sh', fetch )
+
+
+-- Prelude'
+-- --------
+
+-- A version of 'zipWith' that requires the lists to be equal length
+--
+zipWith :: (a -> b -> c) -> [a] -> [b] -> [c]
+zipWith f (x:xs) (y:ys) = f x y : zipWith f xs ys
+zipWith _ []     []     = []
+zipWith _ _      _      = INTERNAL_ERROR(error) "zipWith" "argument mismatch"
+
+zipWith3 :: (a -> b -> c -> d) -> [a] -> [b] -> [c] -> [d]
+zipWith3 f (x:xs) (y:ys) (z:zs) = f x y z : zipWith3 f xs ys zs
+zipWith3 _ []     []     []     = []
+zipWith3 _ _      _      _      = INTERNAL_ERROR(error) "zipWith3" "argument mismatch"
+
+
+-- Split a list into segments of given length
+--
+unconcat :: [Int] -> [a] -> [[a]]
+unconcat []     _  = []
+unconcat (n:ns) xs = let (h,t) = splitAt n xs in h : unconcat ns t
+
diff --git a/Data/Array/Accelerate/CUDA/Compile.hs b/Data/Array/Accelerate/CUDA/Compile.hs
--- a/Data/Array/Accelerate/CUDA/Compile.hs
+++ b/Data/Array/Accelerate/CUDA/Compile.hs
@@ -33,7 +33,7 @@
 import Data.Array.Accelerate.CUDA.CodeGen
 import Data.Array.Accelerate.CUDA.Array.Sugar
 import Data.Array.Accelerate.CUDA.Analysis.Launch
-import Data.Array.Accelerate.CUDA.Foreign                       ( canExecute, canExecuteExp )
+import Data.Array.Accelerate.CUDA.Foreign.Import                ( canExecuteAcc, canExecuteExp )
 import Data.Array.Accelerate.CUDA.Persistent                    as KT
 import qualified Data.Array.Accelerate.CUDA.FullList            as FL
 import qualified Data.Array.Accelerate.CUDA.Debug               as D
@@ -119,9 +119,10 @@
       case pacc of
         -- Environment and control flow
         Avar ix                 -> node $ pure (Avar ix)
-        Alet a b                -> node . pure =<< Alet         <$> traverseAcc a     <*> traverseAcc b
-        Apply f a               -> node . pure =<< Apply        <$> compileOpenAfun f <*> traverseAcc a
-        Acond p t e             -> node =<< liftA3 Acond        <$> travE p <*> travA t <*> travA e
+        Alet a b                -> node . pure =<< Alet         <$> traverseAcc a <*> traverseAcc b
+        Apply f a               -> node =<< liftA2 Apply        <$> travAF f <*> travA a
+        Awhile p f a            -> node =<< liftA3 Awhile       <$> travAF p <*> travAF f <*> travA a
+        Acond p t e             -> node =<< liftA3 Acond        <$> travE  p <*> travA  t <*> travA e
         Atuple tup              -> node =<< liftA Atuple        <$> travAtup tup
         Aprj ix tup             -> node =<< liftA (Aprj ix)     <$> travA    tup
 
@@ -163,7 +164,7 @@
       where
         use :: ArraysR a -> a -> CIO ()
         use ArraysRunit         ()       = return ()
-        use ArraysRarray        arr      = useArray arr
+        use ArraysRarray        arr      = useArrayAsync arr Nothing
         use (ArraysRpair r1 r2) (a1, a2) = use r1 a1 >> use r2 a2
 
         exec :: (Free aenv, PreOpenAcc ExecOpenAcc aenv arrs) -> CIO (ExecOpenAcc aenv arrs)
@@ -183,6 +184,9 @@
           Manifest{}    -> pure                    <$> traverseAcc acc
           Delayed{..}   -> liftA2 (const EmbedAcc) <$> travF indexD <*> travE extentD
 
+        travAF :: DelayedOpenAfun aenv f -> CIO (Free aenv, PreOpenAfun ExecOpenAcc aenv f)
+        travAF afun = pure <$> compileOpenAfun afun
+
         travAtup :: Atuple (DelayedOpenAcc aenv) a -> CIO (Free aenv, Atuple (ExecOpenAcc aenv) a)
         travAtup NilAtup        = return (pure NilAtup)
         travAtup (SnocAtup t a) = liftA2 SnocAtup <$> travAtup t <*> travA a
@@ -207,7 +211,7 @@
                  -> DelayedAfun (a -> r)
                  -> DelayedOpenAcc aenv a
                  -> CIO (Free aenv, PreOpenAcc ExecOpenAcc aenv r)
-        foreignA ff afun a = case canExecute ff of
+        foreignA ff afun a = case canExecuteAcc ff of
           Nothing       -> liftA2 (Aforeign ff)          <$> pure <$> compileAfun afun <*> travA a
           Just _        -> liftA  (Aforeign ff err)      <$> travA a
             where
@@ -237,8 +241,7 @@
         Tuple t                 -> liftA  Tuple                 <$> travT t
         Prj ix e                -> liftA  (Prj ix)              <$> travE e
         Cond p t e              -> liftA3 Cond                  <$> travE p <*> travE t <*> travE e
-        Iterate n f x           -> liftA3 Iterate               <$> travE n <*> travE f <*> travE x
---        While p f x             -> liftA3 While                 <$> travE p <*> travE f <*> travE x
+        While p f x             -> liftA3 While                 <$> travF p <*> travF f <*> travE x
         PrimApp f e             -> liftA  (PrimApp f)           <$> travE e
         Index a e               -> liftA2 Index                 <$> travA a <*> travE e
         LinearIndex a e         -> liftA2 LinearIndex           <$> travA a <*> travE e
@@ -453,15 +456,21 @@
     [ "-I", ddir </> "cubits"
     , "-arch=sm_" ++ show m ++ show n
     , "-cubin"
+--    , "--restrict"    -- requires nvcc >= 5.0
+--    , "--maxrregcount", "32"
     , "-o", cufile `replaceExtension` "cubin"
-    , if D.mode D.dump_cc  then ""   else "--disable-warnings"
-    , if D.mode D.debug_cc then "-G" else "-O3"
+    , if warnings then ""   else "--disable-warnings"
+    , if debug    then ""   else "-DNDEBUG"
+    , if debug    then "-G" else "-O3"
     , machine
     , cufile ]
   where
+    warnings    = D.mode D.dump_cc && D.mode D.verbose
+    debug       = D.mode D.debug_cc
     machine     = case (SIZEOF_HTYPE_INT :: Int) of
                     4   -> "-m32"
                     8   -> "-m64"
+                    _   -> INTERNAL_ERROR(error) "compileFlags" "unknown 'Int' size"
 
 
 -- Open a unique file in the temporary directory used for compilation
@@ -483,13 +492,13 @@
 -- Worker pool
 -- -----------
 
-{-# NOINLINE pool #-}
-pool :: Q.MSem Int
-pool = unsafePerformIO $ Q.new =<< getNumProcessors
+{-# NOINLINE workers #-}
+workers :: Q.MSem Int
+workers = unsafePerformIO $ Q.new =<< getNumProcessors
 
 -- Queue a system process to be executed and return an MVar flag that will be
--- filled once the process completes. The task will only be launched once there
--- is a worker available from the pool. This ensures we don't run out of process
+-- filled once the process completes. The task will only begin once there is a
+-- worker available from the pool. This ensures we don't run out of process
 -- handles or flood the IO bus, degrading performance.
 --
 enqueueProcess :: FilePath -> [String] -> IO (MVar ())
@@ -497,35 +506,27 @@
   mvar  <- newEmptyMVar
   _     <- forkIO $ do
 
-    -- wait for a worker to become available
-    (_, queueT) <- time $ Q.wait pool
-    ccBegin     <- getTime
-    (_,_,_,pid) <- createProcess (proc nvcc flags)
+    -- Wait for a worker to become available
+    (ccT, queueT) <- time $ Q.with workers $ do
 
-    -- asynchronously notify the queue when the compiler has completed
-    _           <- forkIO $ do
+      -- Initiate the external process...
+      ccBegin           <- getTime
+      (_,_,_,pid)       <- createProcess (proc nvcc flags)
 
-      -- Wait for the process to complete
-      --
+      -- ... and wait for it to complete
       waitFor pid
-      ccEnd     <- getTime
-
-      let ccT   = diffTime ccBegin ccEnd
-          msg2  = nvcc ++ " " ++ unwords flags
-          msg1  = "queue: " ++ D.showFFloatSIBase (Just 3) 1000 queueT "s, "
-             ++ "execute: " ++ D.showFFloatSIBase (Just 3) 1000 ccT    "s"
+      ccEnd             <- getTime
 
-      message $ intercalate "\n     ... " [msg1, msg2]
+      return (diffTime ccBegin ccEnd)
+    --
+    let msg2  = nvcc ++ " " ++ unwords flags
+        msg1  = "queue: " ++ D.showFFloatSIBase (Just 3) 1000 queueT "s, "
+           ++ "execute: " ++ D.showFFloatSIBase (Just 3) 1000 ccT    "s"
 
-      -- If there was an error (compilation failed) then this spot in the queue
-      -- is never released and the MVar never signalled that compilation is
-      -- done. This means you'll get a "blocked indefinitely on MVar" error, but
-      -- since the compiler failed in the first places that is somewhat moot.
-      --
-      Q.signal pool
-      putMVar mvar ()
+    message $ intercalate "\n     ... " [msg1, msg2]
 
-    return ()
+    -- Signal to the host thread that the compiled result is available
+    putMVar mvar ()
   --
   return mvar
 
diff --git a/Data/Array/Accelerate/CUDA/Context.hs b/Data/Array/Accelerate/CUDA/Context.hs
--- a/Data/Array/Accelerate/CUDA/Context.hs
+++ b/Data/Array/Accelerate/CUDA/Context.hs
@@ -18,7 +18,7 @@
 
   -- An execution context
   Context(..), create, push, pop, destroy,
-  keepAlive,
+  keepAlive, fromDeviceContext
 
 ) where
 
@@ -55,12 +55,8 @@
 --
 create :: CUDA.Device -> [CUDA.ContextFlag] -> IO Context
 create dev flags = do
-  ctx           <- CUDA.create dev flags >> CUDA.pop >>= keepAlive
-  prp           <- CUDA.props dev
-  weak          <- mkWeakContext ctx $ do
-    message dump_gc $ "gc: finalise context #" ++ show (CUDA.useContext ctx)
-    CUDA.destroy ctx
-  message dump_gc $ "gc: initialise context #" ++ show (CUDA.useContext ctx)
+  ctx                    <- CUDA.create dev flags >> CUDA.pop >>= keepAlive
+  actx@(Context prp _ _) <- fromDeviceContext dev ctx
 
   -- Generated code does not take particular advantage of shared memory, so
   -- for devices that support it use those banks as an L1 cache instead.
@@ -72,8 +68,18 @@
      $ bracket_ (CUDA.push ctx) CUDA.pop (CUDA.setCacheConfig CUDA.PreferL1)
 
   message verbose (deviceInfo dev prp)
-  return $! Context prp ctx weak
+  return actx
 
+-- |Given a device context, construct a new context around it.
+fromDeviceContext :: CUDA.Device -> CUDA.Context -> IO Context
+fromDeviceContext dev ctx = do
+  prp           <- CUDA.props dev
+  weak          <- mkWeakContext ctx $ do
+    message dump_gc $ "gc: finalise context #" ++ show (CUDA.useContext ctx)
+    CUDA.destroy ctx
+  message dump_gc $ "gc: initialise context #" ++ show (CUDA.useContext ctx)
+
+  return $! Context prp ctx weak
 
 -- | Destroy the specified context. This will fail if the context is more than
 -- single attachment.
diff --git a/Data/Array/Accelerate/CUDA/Debug.hs b/Data/Array/Accelerate/CUDA/Debug.hs
--- a/Data/Array/Accelerate/CUDA/Debug.hs
+++ b/Data/Array/Accelerate/CUDA/Debug.hs
@@ -1,8 +1,9 @@
 {-# LANGUAGE CPP             #-}
 {-# LANGUAGE TemplateHaskell #-}
 {-# LANGUAGE TypeOperators   #-}
-{-# OPTIONS -fno-warn-unused-imports #-}
-{-# OPTIONS -fno-warn-unused-binds   #-}
+{-# OPTIONS -fno-warn-incomplete-patterns #-}
+{-# OPTIONS -fno-warn-unused-binds        #-}
+{-# OPTIONS -fno-warn-unused-imports      #-}
 -- |
 -- Module      : Data.Array.Accelerate.CUDA.Debug
 -- Copyright   : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee
@@ -21,7 +22,7 @@
 
   showFFloatSIBase,
 
-  message, trace, event, when, unless, mode,
+  message, trace, event, when, unless, mode, timed, elapsed,
   verbose, flush_cache,
   dump_gc, dump_cc, debug_cc, dump_exec,
 
@@ -31,23 +32,17 @@
 import Data.List
 import Data.Label
 import Data.IORef
-import Control.Monad.IO.Class
+import Debug.Trace                                      ( traceIO, traceEventIO )
+import Control.Monad                                    ( void )
+import Control.Monad.IO.Class                           ( liftIO, MonadIO )
+import Control.Concurrent                               ( forkIO )
 import System.CPUTime
 import System.IO.Unsafe
 import System.Environment
 import System.Console.GetOpt
-
-#if MIN_VERSION_base(4,5,0)
-import Debug.Trace                              ( traceIO, traceEventIO )
-#else
-import Debug.Trace                              ( putTraceMsg )
-
-traceIO :: String -> IO ()
-traceIO = putTraceMsg
+import qualified Foreign.CUDA.Driver.Event              as Event
 
-traceEventIO :: String -> IO ()
-traceEventIO = traceIO
-#endif
+import GHC.Float
 
 
 -- -----------------------------------------------------------------------------
@@ -56,12 +51,20 @@
 showFFloatSIBase :: RealFloat a => Maybe Int -> a -> a -> ShowS
 showFFloatSIBase p b n
   = showString
-  . nubBy (\x y -> x == ' ' && y == ' ')
-  $ showFFloat p n' [ ' ', si_unit ]
+  $ showFFloat p n' (' ':si_unit)
   where
-    n'          = n / (b ^^ (pow-4))
-    pow         = max 0 . min 8 . (+) 4 . floor $ logBase b n
-    si_unit     = "pnµm kMGT" !! pow
+    n'          = n / (b ^^ pow)
+    pow         = (-4) `max` floor (logBase b n) `min` 4        :: Int
+    si_unit     = case pow of
+                       -4 -> "p"
+                       -3 -> "n"
+                       -2 -> "µ"
+                       -1 -> "m"
+                       0  -> ""
+                       1  -> "k"
+                       2  -> "M"
+                       3  -> "G"
+                       4  -> "T"
 
 
 -- -----------------------------------------------------------------------------
@@ -169,4 +172,50 @@
 #else
 unless _ action = action
 #endif
+
+{-# INLINE timed #-}
+timed
+    :: MonadIO m
+    => (Flags :-> Bool)
+    -> (Double -> Double -> String)
+    -> m ()
+    -> m ()
+timed _f _str action
+#ifdef ACCELERATE_DEBUG
+  | mode _f
+  = do
+      gpuBegin  <- liftIO $ Event.create []
+      gpuEnd    <- liftIO $ Event.create []
+      cpuBegin  <- liftIO getCPUTime
+      liftIO $ Event.record gpuBegin Nothing
+      action
+      liftIO $ Event.record gpuEnd Nothing
+      cpuEnd    <- liftIO getCPUTime
+
+      -- Wait for the GPU to finish executing then display the timing execution
+      -- message. Do this in a separate thread so that the remaining kernels can
+      -- be queued asynchronously.
+      --
+      _         <- liftIO . forkIO $ do
+        Event.block gpuEnd
+        diff    <- Event.elapsedTime gpuBegin gpuEnd
+        let gpuTime = float2Double $ diff * 1E-3                         -- milliseconds
+            cpuTime = fromIntegral (cpuEnd - cpuBegin) * 1E-12 :: Double -- picoseconds
+
+        Event.destroy gpuBegin
+        Event.destroy gpuEnd
+        --
+        message _f (_str gpuTime cpuTime)
+      --
+      return ()
+
+  | otherwise
+#endif
+  = action
+
+{-# INLINE elapsed #-}
+elapsed :: Double -> Double -> String
+elapsed gpuTime cpuTime
+  = "gpu: " ++ showFFloatSIBase (Just 3) 1000 gpuTime "s, " ++
+    "cpu: " ++ showFFloatSIBase (Just 3) 1000 cpuTime "s"
 
diff --git a/Data/Array/Accelerate/CUDA/Execute.hs b/Data/Array/Accelerate/CUDA/Execute.hs
--- a/Data/Array/Accelerate/CUDA/Execute.hs
+++ b/Data/Array/Accelerate/CUDA/Execute.hs
@@ -1,14 +1,14 @@
-{-# LANGUAGE BangPatterns         #-}
-{-# LANGUAGE CPP                  #-}
-{-# LANGUAGE FlexibleInstances    #-}
-{-# LANGUAGE GADTs                #-}
-{-# LANGUAGE IncoherentInstances  #-}
-{-# LANGUAGE PatternGuards        #-}
-{-# LANGUAGE ScopedTypeVariables  #-}
-{-# LANGUAGE TypeOperators        #-}
-{-# LANGUAGE TypeSynonymInstances #-}
-{-# LANGUAGE UndecidableInstances #-}
-{-# OPTIONS_GHC -fno-warn-orphans #-}
+{-# LANGUAGE BangPatterns               #-}
+{-# LANGUAGE CPP                        #-}
+{-# LANGUAGE FlexibleInstances          #-}
+{-# LANGUAGE GADTs                      #-}
+{-# LANGUAGE IncoherentInstances        #-}
+{-# LANGUAGE NoForeignFunctionInterface #-}
+{-# LANGUAGE PatternGuards              #-}
+{-# LANGUAGE ScopedTypeVariables        #-}
+{-# LANGUAGE TypeOperators              #-}
+{-# LANGUAGE TypeSynonymInstances       #-}
+{-# LANGUAGE UndecidableInstances       #-}
 -- |
 -- Module      : Data.Array.Accelerate.CUDA.Execute
 -- Copyright   : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee
@@ -33,12 +33,14 @@
 import Data.Array.Accelerate.CUDA.FullList                      ( FullList(..), List(..) )
 import Data.Array.Accelerate.CUDA.Array.Data
 import Data.Array.Accelerate.CUDA.Array.Sugar
-import Data.Array.Accelerate.CUDA.Foreign                       ( canExecute )
+import Data.Array.Accelerate.CUDA.Foreign.Import                ( canExecuteAcc )
 import Data.Array.Accelerate.CUDA.CodeGen.Base                  ( Name, namesOfArray, groupOfInt )
+import Data.Array.Accelerate.CUDA.Execute.Event                 ( Event )
+import Data.Array.Accelerate.CUDA.Execute.Stream                ( Stream )
 import qualified Data.Array.Accelerate.CUDA.Array.Prim          as Prim
-#ifdef ACCELERATE_DEBUG
 import qualified Data.Array.Accelerate.CUDA.Debug               as D
-#endif
+import qualified Data.Array.Accelerate.CUDA.Execute.Event       as Event
+import qualified Data.Array.Accelerate.CUDA.Execute.Stream      as Stream
 
 import Data.Array.Accelerate.Tuple
 import Data.Array.Accelerate.Interpreter                        ( evalPrim, evalPrimConst, evalPrj )
@@ -52,29 +54,67 @@
 import Control.Applicative                                      hiding ( Const )
 import Control.Monad                                            ( join, when, liftM )
 import Control.Monad.Reader                                     ( asks )
+import Control.Monad.State                                      ( gets )
 import Control.Monad.Trans                                      ( MonadIO, liftIO )
 import System.IO.Unsafe                                         ( unsafeInterleaveIO )
 import Data.Int
 import Data.Word
 import Data.Maybe
 
-import Foreign.Ptr                                              ( Ptr, castPtr )
-import Foreign.Storable                                         ( Storable(..) )
 import Foreign.CUDA.Analysis.Device                             ( DeviceProperties, computeCapability, Compute(..) )
 import qualified Foreign.CUDA.Driver                            as CUDA
-import qualified Foreign.Marshal.Array                          as F
 import qualified Data.HashMap.Strict                            as Map
 
-#ifdef ACCELERATE_DEBUG
-import Control.Monad                                            ( void )
-import Control.Concurrent                                       ( forkIO )
-import System.CPUTime
-import qualified Foreign.CUDA.Driver.Event                      as Event
-#endif
-
 #include "accelerate.h"
 
 
+-- Asynchronous kernel execution
+-- -----------------------------
+
+-- Arrays with an associated CUDA Event that will be signalled once the
+-- computation has completed.
+--
+data Async a = Async {-# UNPACK #-} !Event !a
+
+-- Valuation for an environment of asynchronous array computations
+--
+data Aval env where
+  Aempty :: Aval ()
+  Apush  :: Aval env -> Async t -> Aval (env, t)
+
+
+-- Projection of a value from a valuation using a de Bruijn index.
+--
+aprj :: Idx env t -> Aval env -> Async t
+aprj ZeroIdx       (Apush _   x) = x
+aprj (SuccIdx idx) (Apush val _) = aprj idx val
+aprj _             _             = INTERNAL_ERROR(error) "aprj" "inconsistent valuation"
+
+
+-- All work submitted to the given stream will occur after the asynchronous
+-- event for the given array has been fulfilled. Synchronisation is performed
+-- efficiently on the device. This function returns immediately.
+--
+after :: MonadIO m => Stream -> Async a -> m a
+after stream (Async event arr) = liftIO $ Event.after event stream >> return arr
+
+
+-- Block the calling thread until the event for the given array computation
+-- is recorded.
+--
+wait :: MonadIO m => Async a -> m a
+wait (Async e x) = liftIO $ Event.block e >> return x
+
+
+-- Execute the given computation in a unique execution stream.
+--
+streaming :: (Stream -> CIO a) -> CIO (Async a)
+streaming action = do
+  context   <- asks activeContext
+  reservoir <- gets streamReservoir
+  uncurry Async <$> Stream.streaming context reservoir action
+
+
 -- Array expression evaluation
 -- ---------------------------
 
@@ -92,20 +132,20 @@
 --    skeleton are invoked
 --
 executeAcc :: Arrays a => ExecAcc a -> CIO a
-executeAcc !acc = executeOpenAcc acc Empty
+executeAcc !acc = wait =<< streaming (executeOpenAcc acc Aempty)
 
 executeAfun1 :: (Arrays a, Arrays b) => ExecAfun (a -> b) -> a -> CIO b
 executeAfun1 !afun !arrs = do
-  useArrays (arrays arrs) (fromArr arrs)
-  executeOpenAfun1 afun Empty arrs
+  Async event () <- streaming $ \st -> useArrays (arrays arrs) (fromArr arrs) st
+  executeOpenAfun1 afun Aempty (Async event arrs)
   where
-    useArrays :: ArraysR arrs -> arrs -> CIO ()
-    useArrays ArraysRunit         ()       = return ()
-    useArrays (ArraysRpair r1 r0) (a1, a0) = useArrays r1 a1 >> useArrays r0 a0
-    useArrays ArraysRarray        arr      = useArray arr
+    useArrays :: ArraysR arrs -> arrs -> Stream -> CIO ()
+    useArrays ArraysRunit         ()       _  = return ()
+    useArrays (ArraysRpair r1 r0) (a1, a0) st = useArrays r1 a1 st >> useArrays r0 a0 st
+    useArrays ArraysRarray        arr      st = useArrayAsync arr (Just st)
 
-executeOpenAfun1 :: PreOpenAfun ExecOpenAcc aenv (a -> b) -> Val aenv -> a -> CIO b
-executeOpenAfun1 (Alam (Abody f)) aenv x = executeOpenAcc f (aenv `Push` x)
+executeOpenAfun1 :: PreOpenAfun ExecOpenAcc aenv (a -> b) -> Aval aenv -> Async a -> CIO b
+executeOpenAfun1 (Alam (Abody f)) aenv x = wait =<< streaming (executeOpenAcc f (aenv `Apush` x))
 executeOpenAfun1 _                _    _ = error "the sword comes out after you swallow it, right?"
 
 
@@ -114,11 +154,12 @@
 executeOpenAcc
     :: forall aenv arrs.
        ExecOpenAcc aenv arrs
-    -> Val aenv
+    -> Aval aenv
+    -> Stream
     -> CIO arrs
-executeOpenAcc EmbedAcc{} _
+executeOpenAcc EmbedAcc{} _ _
   = INTERNAL_ERROR(error) "execute" "unexpected delayed array"
-executeOpenAcc (ExecAcc (FL () kernel more) !gamma !pacc) !aenv
+executeOpenAcc (ExecAcc (FL () kernel more) !gamma !pacc) !aenv !stream
   = case pacc of
 
       -- Array introduction
@@ -126,15 +167,16 @@
       Unit x                    -> newArray Z . const =<< travE x
 
       -- Environment manipulation
-      Avar ix                   -> return (prj ix aenv)
-      Alet bnd body             -> executeOpenAcc body . (aenv `Push`) =<< travA bnd
+      Avar ix                   -> after stream (aprj ix aenv)
+      Alet bnd body             -> streamA bnd >>= \x -> executeOpenAcc body (aenv `Apush` x) stream
       Atuple tup                -> toTuple <$> travT tup
       Aprj ix tup               -> evalPrj ix . fromTuple <$> travA tup
-      Apply f a                 -> executeOpenAfun1 f aenv =<< travA a
+      Apply f a                 -> executeOpenAfun1 f aenv =<< streamA a
       Acond p t e               -> travE p >>= \x -> if x then travA t else travA e
+      Awhile p f a              -> awhile p f =<< travA a
 
       -- Foreign
-      Aforeign ff afun a        -> fromMaybe (executeAfun1 afun) (canExecute ff) =<< travA a
+      Aforeign ff afun a        -> fromMaybe (executeAfun1 afun) (canExecuteAcc ff) =<< travA a
 
       -- Producers
       Map _ a                   -> executeOp =<< extent a
@@ -168,15 +210,26 @@
 
     -- term traversals
     travA :: ExecOpenAcc aenv a -> CIO a
-    travA !acc = executeOpenAcc acc aenv
+    travA !acc = executeOpenAcc acc aenv stream
 
+    streamA :: ExecOpenAcc aenv a -> CIO (Async a)
+    streamA !acc = streaming $ \st -> executeOpenAcc acc aenv st
+
     travE :: ExecExp aenv t -> CIO t
-    travE !exp = executeExp exp aenv
+    travE !exp = executeExp exp aenv stream
 
     travT :: Atuple (ExecOpenAcc aenv) t -> CIO t
     travT NilAtup          = return ()
     travT (SnocAtup !t !a) = (,) <$> travT t <*> travA a
 
+    awhile :: PreOpenAfun ExecOpenAcc aenv (a -> Scalar Bool) -> PreOpenAfun ExecOpenAcc aenv (a -> a) -> a -> CIO a
+    awhile p f a = do
+      nop <- liftIO Event.create                -- record event never call, so this is a functional no-op
+      r   <- executeOpenAfun1 p aenv (Async nop a)
+      ok  <- indexArray r 0                     -- TLM TODO: memory manager should remember what is already on the host
+      if ok then awhile p f =<< executeOpenAfun1 f aenv (Async nop a)
+            else return a
+
     -- get the extent of an embedded array
     extent :: Shape sh => ExecOpenAcc aenv (Array sh e) -> CIO sh
     extent ExecAcc{}     = INTERNAL_ERROR(error) "executeOpenAcc" "expected delayed array"
@@ -191,7 +244,7 @@
     executeOp :: (Shape sh, Elt e) => sh -> CIO (Array sh e)
     executeOp !sh = do
       out       <- allocateArray sh
-      execute kernel gamma aenv (size sh) out
+      execute kernel gamma aenv (size sh) out stream
       return out
 
     -- Change the shape of an array without altering its contents. This does not
@@ -209,17 +262,21 @@
     fold1Op :: (Shape sh, Elt e) => (sh :. Int) -> CIO (Array sh e)
     fold1Op !sh@(_ :. sz)
       = BOUNDS_CHECK(check) "fold1" "empty array" (sz > 0)
-      $ foldOp sh
+      $ foldCore sh
 
     foldOp :: (Shape sh, Elt e) => (sh :. Int) -> CIO (Array sh e)
     foldOp !(!sh :. sz)
+      = foldCore ((listToShape . map (max 1) . shapeToList $ sh) :. sz)
+
+    foldCore :: (Shape sh, Elt e) => (sh :. Int) -> CIO (Array sh e)
+    foldCore !(!sh :. sz)
       | dim sh > 0              = executeOp sh
       | otherwise
       = let !numElements        = size sh * sz
             (_,!numBlocks,_)    = configure kernel numElements
         in do
           out   <- allocateArray (sh :. numBlocks)
-          execute kernel gamma aenv numElements out
+          execute kernel gamma aenv numElements out stream
           foldRec out
 
     -- Recursive step(s) of a multi-block reduction
@@ -234,7 +291,7 @@
               then return $ Array (fromElt sh) adata
               else do
                 out     <- allocateArray (sh :. numBlocks)
-                execute rec gamma aenv numElements (out, arr)
+                execute rec gamma aenv numElements (out, arr) stream
                 foldRec out
 
       | otherwise
@@ -297,13 +354,13 @@
           --          the final scan result from each interval.
           --
           when (numIntervals > 1) $ do
-            execute upsweep1 gamma aenv numElements blk
-            execute upsweep2 gamma aenv numIntervals (blk, blk, d_sum)
+            execute upsweep1 gamma aenv numElements blk stream
+            execute upsweep2 gamma aenv numIntervals (blk, blk, d_sum) stream
 
           -- Phase 2: Re-scan the input using the carry-in value from each
           --          interval sum calculated in phase 1.
           --
-          execute kernel gamma aenv numElements (Z :. numElements, d_body, blk, d_sum)
+          execute kernel gamma aenv numElements (numElements, d_body, blk, d_sum) stream
 
       | otherwise
       = INTERNAL_ERROR(error) "scanOp" "missing multi-block kernel module(s)"
@@ -314,7 +371,7 @@
     permuteOp !sh !dfs = do
       out <- allocateArray (shape dfs)
       copyArray dfs out
-      execute kernel gamma aenv (size sh) out
+      execute kernel gamma aenv (size sh) out stream
       return out
 
     -- Stencil operations. NOTE: the arguments to 'namesOfArray' must be the
@@ -328,50 +385,56 @@
 
       if computeCapability dev < Compute 2 0
          then marshalAccTex (namesOfArray "Stencil" (undefined :: a)) kernel arr >>
-              execute kernel gamma aenv (size sh) (out, sh)
-         else execute kernel gamma aenv (size sh) (out, arr)
+              execute kernel gamma aenv (size sh) (out, sh) stream
+         else execute kernel gamma aenv (size sh) (out, arr) stream
       --
       return out
 
     stencil2Op :: forall sh a b c. (Shape sh, Elt a, Elt b, Elt c)
                => Array sh a -> Array sh b -> CIO (Array sh c)
-    stencil2Op !arr1 !arr2 = do
-      let sh1   =  shape arr1
-          sh2   =  shape arr2
-          sh    =  sh1 `intersect` sh2
-      out       <- allocateArray sh
-      dev       <- asks deviceProperties
+    stencil2Op !arr1 !arr2
+      | Cons _ spec _ <- more
+      = let sh1         =  shape arr1
+            sh2         =  shape arr2
+            (sh, op)
+              | fromElt sh1 == fromElt sh2      = (sh1,                 spec)
+              | otherwise                       = (sh1 `intersect` sh2, kernel)
+        in do
+          out   <- allocateArray sh
+          dev   <- asks deviceProperties
 
-      if computeCapability dev < Compute 2 0
-         then marshalAccTex (namesOfArray "Stencil1" (undefined :: a)) kernel arr1 >>
-              marshalAccTex (namesOfArray "Stencil2" (undefined :: b)) kernel arr2 >>
-              execute kernel gamma aenv (size sh) (out, sh1,  sh2)
-         else execute kernel gamma aenv (size sh) (out, arr1, arr2)
-      --
-      return out
+          if computeCapability dev < Compute 2 0
+             then marshalAccTex (namesOfArray "Stencil1" (undefined :: a)) op arr1 >>
+                  marshalAccTex (namesOfArray "Stencil2" (undefined :: b)) op arr2 >>
+                  execute op gamma aenv (size sh) (out, sh1,  sh2) stream
+             else execute op gamma aenv (size sh) (out, arr1, arr2) stream
+          --
+          return out
 
+      | otherwise
+      = INTERNAL_ERROR(error) "stencil2Op" "missing stencil specialisation kernel"
 
+
 -- Scalar expression evaluation
 -- ----------------------------
 
-executeExp :: ExecExp aenv t -> Val aenv -> CIO t
-executeExp !exp !aenv = executeOpenExp exp Empty aenv
+executeExp :: ExecExp aenv t -> Aval aenv -> Stream -> CIO t
+executeExp !exp !aenv !stream = executeOpenExp exp Empty aenv stream
 
-executeOpenExp :: forall env aenv exp. ExecOpenExp env aenv exp -> Val env -> Val aenv -> CIO exp
-executeOpenExp !rootExp !env !aenv = travE rootExp
+executeOpenExp :: forall env aenv exp. ExecOpenExp env aenv exp -> Val env -> Aval aenv -> Stream -> CIO exp
+executeOpenExp !rootExp !env !aenv !stream = travE rootExp
   where
     travE :: ExecOpenExp env aenv t -> CIO t
     travE exp = case exp of
       Var ix                    -> return (prj ix env)
-      Let bnd body              -> travE bnd >>= \x -> executeOpenExp body (env `Push` x) aenv
+      Let bnd body              -> travE bnd >>= \x -> executeOpenExp body (env `Push` x) aenv stream
       Const c                   -> return (toElt c)
       PrimConst c               -> return (evalPrimConst c)
       PrimApp f x               -> evalPrim f <$> travE x
       Tuple t                   -> toTuple <$> travT t
       Prj ix e                  -> evalPrj ix . fromTuple <$> travE e
       Cond p t e                -> travE p >>= \x -> if x then travE t else travE e
-      Iterate n f x             -> join $ iterate f <$> travE n <*> travE x
---      While p f x               -> while p f =<< travE x
+      While p f x               -> while p f =<< travE x
       IndexAny                  -> return Any
       IndexNil                  -> return Z
       IndexCons sh sz           -> (:.) <$> travE sh <*> travE sz
@@ -386,7 +449,7 @@
       Shape acc                 -> shape <$> travA acc
       Index acc ix              -> join $ index      <$> travA acc <*> travE ix
       LinearIndex acc ix        -> join $ indexArray <$> travA acc <*> travE ix
-      Foreign _ f x             -> travF1 f x
+      Foreign _ f x             -> foreign f x
 
     -- Helpers
     -- -------
@@ -397,30 +460,21 @@
       SnocTup !t !e     -> (,) <$> travT t <*> travE e
 
     travA :: ExecOpenAcc aenv a -> CIO a
-    travA !acc = executeOpenAcc acc aenv
+    travA !acc = executeOpenAcc acc aenv stream
 
-    travF1 :: ExecFun () (a -> b) -> ExecOpenExp env aenv a -> CIO b
-    travF1 (Lam (Body f)) x = travE x >>= \a -> executeOpenExp f (Empty `Push` a) Empty
-    travF1 _              _ = error "I bless the rains down in Africa"
+    foreign :: ExecFun () (a -> b) -> ExecOpenExp env aenv a -> CIO b
+    foreign (Lam (Body f)) x = travE x >>= \e -> executeOpenExp f (Empty `Push` e) Aempty stream
+    foreign _              _ = error "I bless the rains down in Africa"
 
-    iterate :: ExecOpenExp (env,a) aenv a -> Int -> a -> CIO a
-    iterate !f !limit !x
-      = let go !i !acc
-              | i >= limit      = return acc
-              | otherwise       = go (i+1) =<< executeOpenExp f (env `Push` acc) aenv
-        in
-        go 0 x
+    travF1 :: ExecOpenFun env aenv (a -> b) -> a -> CIO b
+    travF1 (Lam (Body f)) x = executeOpenExp f (env `Push` x) aenv stream
+    travF1 _              _ = error "Gonna take some time to do the things we never have"
 
-{--
-    while :: ExecOpenExp (env,a) aenv Bool -> ExecOpenExp (env,a) aenv a -> a -> CIO a
-    while !p !f !x
-      = let go !acc = do
-              done <- executeOpenExp p (env `Push` acc) aenv
-              if done then return x
-                      else go =<< executeOpenExp f (env `Push` acc) aenv
-        in
-        go x
---}
+    while :: ExecOpenFun env aenv (a -> Bool) -> ExecOpenFun env aenv (a -> a) -> a -> CIO a
+    while !p !f !x = do
+      ok <- travF1 p x
+      if ok then while p f =<< travF1 f x
+            else return x
 
     indexSlice :: (Elt slix, Elt sh, Elt sl)
                => SliceIndex (EltRepr slix) (EltRepr sl) co (EltRepr sh)
@@ -468,7 +522,7 @@
   marshal !ad = marshalArrayData ad
 
 instance Shape sh => Marshalable sh where
-  marshal !sh = return [CUDA.VArg sh]
+  marshal !sh = marshal (reverse (shapeToList sh))
 
 instance Marshalable a => Marshalable [a] where
   marshal = concatMapM marshal
@@ -505,26 +559,9 @@
 primMarshalable(Word64)
 primMarshalable(Float)
 primMarshalable(Double)
-primMarshalable(Ptr a)
 primMarshalable(CUDA.DevicePtr a)
 
-instance Shape sh => Storable sh where  -- undecidable, incoherent
-  sizeOf sh     = sizeOf    (undefined :: Int32) * (dim sh)
-  alignment _   = alignment (undefined :: Int32)
-  poke !p !sh   = F.pokeArray (castPtr p) (convertShape (shapeToList sh))
 
-
--- Convert shapes into 32-bit integers for marshalling onto the device
---
-convertShape :: [Int] -> [Int32]
-convertShape [] = [1]
-convertShape sh = reverse (map convertIx sh)
-
-convertIx :: Int -> Int32
-convertIx !ix = INTERNAL_ASSERT "convertIx" (ix <= fromIntegral (maxBound :: Int32))
-              $ fromIntegral ix
-
-
 -- Note [Array references in scalar code]
 -- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 --
@@ -542,11 +579,11 @@
 -- as to use the larger available caches.
 --
 
-marshalAccEnvTex :: AccKernel a -> Val aenv -> Gamma aenv -> CIO [CUDA.FunParam]
-marshalAccEnvTex !kernel !aenv (Gamma !gamma)
+marshalAccEnvTex :: AccKernel a -> Aval aenv -> Gamma aenv -> Stream -> CIO [CUDA.FunParam]
+marshalAccEnvTex !kernel !aenv (Gamma !gamma) !stream
   = flip concatMapM (Map.toList gamma)
   $ \(Idx_ !(idx :: Idx aenv (Array sh e)), i) ->
-        do let arr = prj idx aenv
+        do arr <- after stream (aprj idx aenv)
            marshalAccTex (namesOfArray (groupOfInt i) (undefined :: e)) kernel arr
            marshal (shape arr)
 
@@ -554,9 +591,9 @@
 marshalAccTex (_, !arrIn) (AccKernel _ _ !mdl _ _ _ _) (Array !sh !adata)
   = marshalTextureData adata (R.size sh) =<< liftIO (sequence' $ map (CUDA.getTex mdl) (reverse arrIn))
 
-marshalAccEnvArg :: Val aenv -> Gamma aenv -> CIO [CUDA.FunParam]
-marshalAccEnvArg !aenv (Gamma !gamma)
-  = concatMapM (\(Idx_ !idx) -> marshal (prj idx aenv)) (Map.keys gamma)
+marshalAccEnvArg :: Aval aenv -> Gamma aenv -> Stream -> CIO [CUDA.FunParam]
+marshalAccEnvArg !aenv (Gamma !gamma) !stream
+  = concatMapM (\(Idx_ !idx) -> marshal =<< after stream (aprj idx aenv)) (Map.keys gamma)
 
 
 -- A lazier version of 'Control.Monad.sequence'
@@ -587,16 +624,17 @@
 --
 arguments :: Marshalable args
           => AccKernel a
-          -> Val aenv
+          -> Aval aenv
           -> Gamma aenv
           -> args
+          -> Stream
           -> CIO [CUDA.FunParam]
-arguments !kernel !aenv !gamma !a = do
+arguments !kernel !aenv !gamma !a !stream = do
   dev <- asks deviceProperties
   let marshaller | computeCapability dev < Compute 2 0   = marshalAccEnvTex kernel
                  | otherwise                             = marshalAccEnvArg
   --
-  (++) <$> marshaller aenv gamma <*> marshal a
+  (++) <$> marshaller aenv gamma stream <*> marshal a
 
 
 -- Link the binary object implementing the computation, configure the kernel
@@ -606,63 +644,25 @@
 execute :: Marshalable args
         => AccKernel a                  -- The binary module implementing this kernel
         -> Gamma aenv                   -- variables of arrays embedded in scalar expressions
-        -> Val aenv                     -- the environment
+        -> Aval aenv                    -- the environment
         -> Int                          -- a "size" parameter, typically number of elements in the output
         -> args                         -- arguments to marshal to the kernel function
+        -> Stream                       -- Compute stream to execute in
         -> CIO ()
-execute !kernel !gamma !aenv !n !a = do
-  args <- arguments kernel aenv gamma a
-  launch kernel (configure kernel n) args
+execute !kernel !gamma !aenv !n !a !stream = do
+  args <- arguments kernel aenv gamma a stream
+  launch kernel (configure kernel n) args stream
 
 
 -- Execute a device function, with the given thread configuration and function
 -- parameters. The tuple contains (threads per block, grid size, shared memory)
 --
-launch :: AccKernel a -> (Int,Int,Int) -> [CUDA.FunParam] -> CIO ()
-launch (AccKernel _entry !fn _ _ _ _ _) !(cta, grid, smem) !args
-#ifdef ACCELERATE_DEBUG
-  | D.mode D.dump_exec
-  = liftIO $ do
-      gpuBegin  <- Event.create []
-      gpuEnd    <- Event.create []
-      cpuBegin  <- getCPUTime
-      Event.record gpuBegin Nothing
-      CUDA.launchKernel fn (grid,1,1) (cta,1,1) smem Nothing args
-      Event.record gpuEnd Nothing
-      cpuEnd    <- getCPUTime
-
-      -- Wait for the GPU to finish executing then display the timing execution
-      -- message. Do this in a separate thread so that the remaining kernels can
-      -- be queued asynchronously.
-      --
-      void . forkIO $ do
-        Event.block gpuEnd
-        diff    <- Event.elapsedTime gpuBegin gpuEnd
-        let gpuTime = diff * 1E-3                                        -- milliseconds
-            cpuTime = fromIntegral (cpuEnd - cpuBegin) * 1E-12 :: Double -- picoseconds
-
-        Event.destroy gpuBegin
-        Event.destroy gpuEnd
-        --
-        message $
-          _entry ++ "<<< " ++ shows grid ", " ++ shows cta ", " ++ shows smem " >>> "
-                 ++ "gpu: " ++ D.showFFloatSIBase (Just 3) 1000 gpuTime "s, "
-                 ++ "cpu: " ++ D.showFFloatSIBase (Just 3) 1000 cpuTime "s"
-#endif
-  | otherwise
-  = liftIO $ CUDA.launchKernel fn (grid,1,1) (cta,1,1) smem Nothing args
-
-
--- Debugging
--- ---------
-
-#ifdef ACCELERATE_DEBUG
-{-# INLINE trace #-}
-trace :: MonadIO m => String -> m a -> m a
-trace msg next = D.message D.dump_exec ("exec: " ++ msg) >> next
-
-{-# INLINE message #-}
-message :: MonadIO m => String -> m ()
-message s = s `trace` return ()
-#endif
+launch :: AccKernel a -> (Int,Int,Int) -> [CUDA.FunParam] -> Stream -> CIO ()
+launch (AccKernel entry !fn _ _ _ _ _) !(cta, grid, smem) !args !stream
+  = D.timed D.dump_exec msg
+  $ liftIO $ CUDA.launchKernel fn (grid,1,1) (cta,1,1) smem (Just stream) args
+  where
+    msg gpuTime cpuTime
+      = "exec: " ++ entry ++ "<<< " ++ shows grid ", " ++ shows cta ", " ++ shows smem " >>> "
+                 ++ D.elapsed gpuTime cpuTime
 
diff --git a/Data/Array/Accelerate/CUDA/Execute/Event.hs b/Data/Array/Accelerate/CUDA/Execute/Event.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Accelerate/CUDA/Execute/Event.hs
@@ -0,0 +1,86 @@
+{-# LANGUAGE BangPatterns  #-}
+{-# LANGUAGE MagicHash     #-}
+{-# LANGUAGE UnboxedTuples #-}
+-- |
+-- Module      : Data.Array.Accelerate.CUDA.Execute.Event
+-- Copyright   : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee
+--               [2009..2013] 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.CUDA.Execute.Event (
+
+  Event, create, waypoint, after, block,
+
+) where
+
+-- friends
+import qualified Data.Array.Accelerate.CUDA.Debug               as D
+
+-- libraries
+import Foreign.CUDA.Driver.Event                                ( Event(..) )
+import Foreign.CUDA.Driver.Stream                               ( Stream )
+import qualified Foreign.CUDA.Driver.Event                      as Event
+
+import GHC.Base
+import GHC.Ptr
+
+
+-- Create a new event that will be automatically garbage collected. The event is
+-- not suitable for timing purposes.
+--
+{-# INLINE create #-}
+create :: IO Event
+create = do
+  event <- Event.create [Event.DisableTiming]
+  addEventFinalizer event $ trace ("destroy " ++ show event) (Event.destroy event)
+  return event
+
+-- Create a new event marker that will be filled once execution in the specified
+-- stream has completed all previously submitted work.
+--
+{-# INLINE waypoint #-}
+waypoint :: Stream -> IO Event
+waypoint stream = do
+  event <- create
+  Event.record event (Just stream)
+--  message $ "waypoint " ++ show event ++ " in " ++ show stream
+  return event
+
+-- Make all future work submitted to the given stream wait until the event
+-- reports completion before beginning execution.
+--
+{-# INLINE after #-}
+after :: Event -> Stream -> IO ()
+after event stream = do
+--  message $ "after " ++ show event ++ " in " ++ show stream
+  Event.wait event (Just stream) []
+
+-- Block the calling thread until the event is recorded
+--
+{-# INLINE block #-}
+block :: Event -> IO ()
+block = Event.block
+
+
+-- Add a finaliser to an event token
+--
+addEventFinalizer :: Event -> IO () -> IO ()
+addEventFinalizer e@(Event (Ptr e#)) f = IO $ \s ->
+  case mkWeak# e# e f s of (# s', _w #) -> (# s', () #)
+
+
+-- Debug
+-- -----
+
+{-# INLINE trace #-}
+trace :: String -> IO a -> IO a
+trace msg next = D.message D.dump_exec ("event: " ++ msg) >> next
+
+-- {-# INLINE message #-}
+-- message :: String -> IO ()
+-- message s = s `trace` return ()
+
diff --git a/Data/Array/Accelerate/CUDA/Execute/Stream.hs b/Data/Array/Accelerate/CUDA/Execute/Stream.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Accelerate/CUDA/Execute/Stream.hs
@@ -0,0 +1,162 @@
+{-# LANGUAGE BangPatterns #-}
+-- |
+-- Module      : Data.Array.Accelerate.CUDA.Execute.Stream
+-- Copyright   : [2013] 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.CUDA.Execute.Stream (
+
+  Stream, Reservoir, new, streaming,
+
+) where
+
+-- friends
+import Data.Array.Accelerate.CUDA.Array.Nursery                 ( ) -- hashable CUDA.Context instance
+import Data.Array.Accelerate.CUDA.Context                       ( Context(..) )
+import Data.Array.Accelerate.CUDA.FullList                      ( FullList(..) )
+import Data.Array.Accelerate.CUDA.Execute.Event                 ( Event )
+import qualified Data.Array.Accelerate.CUDA.Execute.Event       as Event
+import qualified Data.Array.Accelerate.CUDA.FullList            as FL
+import qualified Data.Array.Accelerate.CUDA.Debug               as D
+
+-- libraries
+import Control.Monad                                            ( void )
+import Control.Monad.Trans                                      ( MonadIO, liftIO )
+import Control.Exception                                        ( bracket_ )
+import Control.Concurrent                                       ( forkIO )
+import Control.Concurrent.MVar                                  ( MVar, newMVar, withMVar, mkWeakMVar )
+import System.Mem.Weak                                          ( Weak, deRefWeak )
+import Foreign.CUDA.Driver.Stream                               ( Stream )
+import qualified Foreign.CUDA.Driver                            as CUDA
+import qualified Foreign.CUDA.Driver.Stream                     as Stream
+
+import qualified Data.HashTable.IO                              as HT
+
+
+-- Representation
+-- --------------
+
+-- The Reservoir is a place to store CUDA execution streams that are currently
+-- inactive. When a new stream is requested one is provided from the reservoir
+-- if available, otherwise a fresh execution stream is created.
+--
+type HashTable key val  = HT.BasicHashTable key val
+
+type RSV                = MVar ( HashTable CUDA.Context (FullList () Stream) )
+data Reservoir          = Reservoir {-# UNPACK #-} !RSV
+                                    {-# UNPACK #-} !(Weak RSV)
+
+
+-- Executing operations in streams
+-- -------------------------------
+
+-- Execute an operation in a unique execution stream.
+--
+{-# INLINE streaming #-}
+streaming :: MonadIO m => Context -> Reservoir -> (Stream -> m a) -> m (Event, a)
+streaming !ctx !rsv@(Reservoir !_ !weak_rsv) !action = do
+  stream <- liftIO $ create ctx rsv
+  result <- action stream
+  end    <- liftIO $ Event.waypoint stream
+  liftIO $ merge (weakContext ctx) weak_rsv stream
+  return $ (end, result)
+
+
+-- Primitive operations
+-- --------------------
+
+-- Generate a new empty reservoir. It is not necessary to pre-populate it with
+-- any streams because stream creation does not cause a device synchronisation.
+--
+new :: IO Reservoir
+new = do
+  tbl    <- HT.new
+  ref    <- newMVar tbl
+  weak   <- mkWeakMVar ref (flush tbl)
+  return $! Reservoir ref weak
+
+
+-- Create a CUDA execution stream. If an inactive stream is available for use,
+-- that is returned, else a fresh stream is created.
+--
+{-# INLINE create #-}
+create :: Context -> Reservoir -> IO Stream
+create !ctx (Reservoir !ref _) = withMVar ref $ \tbl -> do
+  let key = deviceContext ctx
+  --
+  ms    <- HT.lookup tbl key
+  case ms of
+    Nothing -> do
+        stream <- Stream.create []
+        message ("new " ++ showStream stream)
+        return stream
+
+    Just (FL () stream rest) -> do
+      case rest of
+        FL.Nil          -> HT.delete tbl key
+        FL.Cons () s ss -> HT.insert tbl key (FL () s ss)
+      --
+      return stream
+
+
+-- Merge a stream back into the reservoir. This is done asynchronously, once all
+-- pending operations in the stream have completed.
+--
+{-# INLINE merge #-}
+merge :: Weak CUDA.Context -> Weak RSV -> Stream -> IO ()
+merge !weak_ctx !weak_rsv !stream
+  = void . forkIO
+  $ do  -- wait for all operations to complete
+        -- Stream.block stream
+
+        -- Now check whether the context and reservoir are still active. Return
+        -- the stream back to the reservoir for later reuse if we can, otherwise
+        -- destroy it.
+        --
+        mc <- deRefWeak weak_ctx
+        case mc of
+          Nothing       -> message ("finalise/dead context " ++ showStream stream)
+          Just ctx      -> do
+            --
+            mr <- deRefWeak weak_rsv
+            case mr of
+              Nothing   -> trace ("merge/free " ++ showStream stream) $ Stream.destroy stream
+              Just ref  -> trace ("merge/save " ++ showStream stream) $ withMVar ref $ \tbl -> do
+                --
+                ms <- HT.lookup tbl ctx
+                case ms of
+                  Nothing       -> HT.insert tbl ctx (FL.singleton () stream)
+                  Just ss       -> HT.insert tbl ctx (FL.cons () stream ss)
+
+
+-- Destroy all streams in the reservoir.
+--
+flush :: HashTable CUDA.Context (FullList () Stream) -> IO ()
+flush !tbl =
+  let clean (!ctx,!ss) = do
+        bracket_ (CUDA.push ctx) CUDA.pop (FL.mapM_ (const Stream.destroy) ss)
+        HT.delete tbl ctx
+  in
+  message "flush reservoir" >> HT.mapM_ clean tbl
+
+
+-- Debug
+-- -----
+
+{-# INLINE trace #-}
+trace :: String -> IO a -> IO a
+trace msg next = D.message D.dump_exec ("stream: " ++ msg) >> next
+
+{-# INLINE message #-}
+message :: String -> IO ()
+message s = s `trace` return ()
+
+{-# INLINE showStream #-}
+showStream :: Stream -> String
+showStream (Stream.Stream s) = show s
+
diff --git a/Data/Array/Accelerate/CUDA/Foreign.hs b/Data/Array/Accelerate/CUDA/Foreign.hs
--- a/Data/Array/Accelerate/CUDA/Foreign.hs
+++ b/Data/Array/Accelerate/CUDA/Foreign.hs
@@ -1,13 +1,3 @@
-{-# LANGUAGE CPP                  #-}
-{-# LANGUAGE ConstraintKinds      #-}
-{-# LANGUAGE DeriveDataTypeable   #-}
-{-# LANGUAGE FlexibleContexts     #-}
-{-# LANGUAGE FlexibleInstances    #-}
-{-# LANGUAGE RankNTypes           #-}
-{-# LANGUAGE ScopedTypeVariables  #-}
-{-# LANGUAGE TypeFamilies         #-}
-{-# LANGUAGE UndecidableInstances #-}
-{-# OPTIONS_GHC -fno-warn-orphans #-}
 -- |
 -- Module      : Data.Array.Accelerate.CUDA.Foreign
 -- Copyright   : [2013] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell, Robert Clifton-Everest
@@ -17,140 +7,14 @@
 -- Stability   : experimental
 -- Portability : non-portable (GHC extensions)
 --
--- This module provides the CUDA backend's implementation of Accelerate's
--- foreign function interface. Also provided are a series of utility functions
--- for transferring arrays from the device to the host (and vice-versa),
--- allocating new arrays, getting the CUDA device pointers of a given array, and
--- executing IO actions within a CUDA context.
---
--- [/NOTE:/]
---
--- When arrays are passed to the foreign function there is no guarantee that the
--- host side data matches the device side data. If the data is needed host side
--- 'peekArray' or 'peekArrayAsync' must be called.
---
--- Arrays of tuples are represented as tuples of arrays so for example an array
--- of type @Array DIM1 (Float, Float)@ would have two device pointers associated
--- with it.
---
 
 module Data.Array.Accelerate.CUDA.Foreign (
 
-  -- * Backend representation
-  cudaAcc, canExecute, CuForeignAcc, CuForeignExp, CIO,
-  liftIO, canExecuteExp, cudaExp,
-
-  -- * Manipulating arrays
-  DevicePtrs,
-  devicePtrsOfArray,
-  indexArray, copyArray,
-  useArray,  useArrayAsync,
-  peekArray, peekArrayAsync,
-  pokeArray, pokeArrayAsync,
-  allocateArray, newArray,
-
-  -- * Running IO actions in a CUDA context
-  inContext, inDefaultContext
+  module Data.Array.Accelerate.CUDA.Foreign.Export,
+  module Data.Array.Accelerate.CUDA.Foreign.Import,
 
 ) where
 
-import Data.Array.Accelerate.CUDA.State
-import Data.Array.Accelerate.CUDA.Context
-import Data.Array.Accelerate.CUDA.Array.Sugar
-import Data.Array.Accelerate.CUDA.Array.Data
-import Data.Array.Accelerate.CUDA.Array.Prim            ( DevicePtrs )
-
-import Data.Dynamic
-import Control.Applicative
-import Control.Exception                                ( bracket_ )
-import Control.Monad.Trans                              ( liftIO )
-import System.IO.Unsafe                                 ( unsafePerformIO )
-import System.Mem.StableName
-
-
--- CUDA backend representation of foreign functions
--- ------------------------------------------------
-
--- CUDA foreign Acc functions are just CIO functions.
---
-newtype CuForeignAcc args results = CuForeignAcc (args -> CIO results)
-  deriving (Typeable)
-
-instance Foreign CuForeignAcc where
-  -- Using the hash of the StableName in order to uniquely identify the function
-  -- when it is pretty printed.
-  --
-  strForeign ff =
-    let sn = unsafePerformIO $ makeStableName ff
-    in
-    "cudaAcc<" ++ (show (hashStableName sn)) ++ ">"
-
--- |Gives the executable form of a foreign function if it can be executed by the
--- CUDA backend.
---
-canExecute :: forall ff args results. (Foreign ff, Typeable args, Typeable results)
-           => ff args results
-           -> Maybe (args -> CIO results)
-canExecute ff =
-  let
-    df = toDyn ff
-    fd = fromDynamic :: Dynamic -> Maybe (CuForeignAcc args results)
-  in (\(CuForeignAcc ff') -> ff') <$> fd df
-
--- CUDA foreign Exp functions are just strings with the header filename and the name of the
--- function separated by a space.
---
-newtype CuForeignExp args results = CuForeignExp String
-  deriving (Typeable)
-
-instance Foreign CuForeignExp where
-  strForeign (CuForeignExp n) = "cudaExp<" ++ n ++ ">"
-
--- |Gives the foreign function name as a string if it is a foreign Exp function
--- for the CUDA backend.
---
-canExecuteExp :: forall ff args results. (Foreign ff, Typeable results, Typeable args)
-              => ff args results
-              -> Maybe String
-canExecuteExp ff =
-  let
-    df = toDyn ff
-    fd = fromDynamic :: Dynamic -> Maybe (CuForeignExp args results)
-  in (\(CuForeignExp ff') -> ff') <$> fd df
-
-
--- User facing utility functions
--- -----------------------------
-
--- |Create a CUDA foreign function
---
-cudaAcc :: (Arrays args, Arrays results)
-        => (args -> CIO results)
-        -> CuForeignAcc args results
-cudaAcc = CuForeignAcc
-
--- |Create a CUDA foreign scalar function. The string needs to be formatted in
--- the same way as for the Haskell FFI. That is, the header file name and the
--- name of the function separated by a space. i.e cudaExp "stdlib.h min".
---
-cudaExp :: (Elt args, Elt results)
-        => String
-        -> CuForeignExp args results
-cudaExp = CuForeignExp
-
--- |Get the raw CUDA device pointers associated with an array
---
-devicePtrsOfArray :: Array sh e -> CIO (DevicePtrs (EltRepr e))
-devicePtrsOfArray (Array _ adata) = devicePtrsOfArrayData adata
-
--- |Run an IO action within the given CUDA context
---
-inContext :: Context -> IO a -> IO a
-inContext ctx action =
-  bracket_ (push ctx) pop action
-
--- |Run an IO action in the default CUDA context
---
-inDefaultContext :: IO a -> IO a
-inDefaultContext = inContext defaultContext
+import Data.Array.Accelerate.CUDA.Foreign.Export
+import Data.Array.Accelerate.CUDA.Foreign.Import
 
diff --git a/Data/Array/Accelerate/CUDA/Foreign/Export.hs b/Data/Array/Accelerate/CUDA/Foreign/Export.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Accelerate/CUDA/Foreign/Export.hs
@@ -0,0 +1,276 @@
+{-# LANGUAGE RankNTypes               #-}
+{-# LANGUAGE ScopedTypeVariables      #-}
+{-# LANGUAGE BangPatterns             #-}
+{-# LANGUAGE GADTs                    #-}
+{-# LANGUAGE CPP                      #-}
+{-# LANGUAGE ForeignFunctionInterface #-}
+{-# LANGUAGE TemplateHaskell          #-}
+{-# LANGUAGE QuasiQuotes              #-}
+{-# LANGUAGE TypeFamilies             #-}
+{-# LANGUAGE FlexibleInstances        #-}
+{-# LANGUAGE ImpredicativeTypes       #-}
+{-# OPTIONS_GHC -fno-warn-name-shadowing -fno-warn-orphans #-}
+-- |
+-- Module      : Data.Array.Accelerate.CUDA.Foreign.Export
+-- Copyright   : [2013] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell, Robert Clifton-Everest
+-- License     : BSD3
+--
+-- Maintainer  : Robert Clifton-Everest <robertce@cse.unsw.edu.au>
+-- Stability   : experimental
+-- Portability : non-portable (GHC extensions)
+--
+-- This module is intended to allow the calling of Accelerate functions from
+-- within CUDA C/C++code. See the nbody/visualizer example in the accelerate-examples
+-- package to see how it is used.
+--
+
+module Data.Array.Accelerate.CUDA.Foreign.Export (
+
+  -- ** Functions callable from foreign code
+  -- In order to call these from from C, see the corresponding C function signature.
+  accelerateCreate, accelerateDestroy, freeOutput, freeProgram,
+
+  -- ** Exporting
+  exportAfun, buildExported,
+
+  -- ** Types
+  InputArray, OutputArray, ShapeBuffer, DevicePtrBuffer,
+
+) where
+
+import Prelude                                          as P
+import Data.Functor
+import Control.Applicative
+import Foreign.StablePtr
+import Foreign.C.Types
+import Foreign.Ptr
+import Foreign.Storable                                 ( Storable(..) )
+import Foreign.Marshal.Array                            ( peekArray, pokeArray, mallocArray )
+import Foreign.Marshal.Alloc                            ( free )
+import Control.Monad.State                              ( liftIO )
+import qualified Foreign.CUDA.Driver                    as CUDA
+import Language.Haskell.TH                              hiding ( ppr )
+
+-- friends
+import Data.Array.Accelerate.Smart                      ( Acc )
+import Data.Array.Accelerate.Type
+import Data.Array.Accelerate.Array.Data
+import Data.Array.Accelerate.CUDA                       ( run1In )
+import Data.Array.Accelerate.CUDA.Array.Sugar           hiding ( shape, size )
+import Data.Array.Accelerate.CUDA.Array.Data            hiding ( pokeArray, peekArray, mallocArray )
+import Data.Array.Accelerate.CUDA.State
+import Data.Array.Accelerate.CUDA.Context
+
+-- |A handle foreign code can use to call accelerate functions.
+type AccHandle = StablePtr Context
+
+-- |A foreign buffer that represents a shape as an array of ints.
+type ShapeBuffer = Ptr CInt
+
+-- |A buffer of device pointers
+type DevicePtrBuffer = Ptr WordPtr
+
+-- |The input required from foreign code.
+type InputArray = (ShapeBuffer, DevicePtrBuffer)
+
+-- |A result array from an accelerate program.
+type OutputArray = (ShapeBuffer, DevicePtrBuffer, StablePtr EArray)
+
+-- |Foreign exportable representation of a CUDA device
+type Device = Int32
+
+-- |Foreign representation of a CUDA context.
+type ForeignContext = Ptr ()
+
+-- We need to capture the Arrays constraint
+data Afun where
+  Afun :: (Arrays a, Arrays b)
+       => (a -> b)
+       -> a {- dummy -}
+       -> b {- dummy -}
+       -> Afun
+
+data EArray where
+  EArray :: (Shape sh, Elt e) => Array sh e -> EArray
+
+-- We need to export these
+foreign export ccall accelerateCreate  :: Device -> ForeignContext -> IO AccHandle
+foreign export ccall accelerateDestroy :: AccHandle -> IO ()
+foreign export ccall runProgram        :: AccHandle -> StablePtr Afun -> Ptr InputArray -> Ptr OutputArray -> IO ()
+foreign export ccall freeOutput        :: Ptr OutputArray -> IO ()
+foreign export ccall freeProgram       :: StablePtr a -> IO ()
+
+instance Storable InputArray where
+  sizeOf (sh, ptrs) = sizeOf sh + sizeOf ptrs
+
+  alignment _ = 0
+
+  peek ptr = do
+    let p_sh   = castPtr ptr :: Ptr ShapeBuffer
+    sh         <- peek p_sh
+    let p_ptrs = (castPtr p_sh :: Ptr DevicePtrBuffer) `plusPtr` sizeOf sh
+    ptrs       <- peek p_ptrs
+    return (sh, ptrs)
+
+  poke ptr (sh, ptrs) = do
+    let p_sh   = castPtr ptr :: Ptr ShapeBuffer
+        p_ptrs = (castPtr p_sh :: Ptr DevicePtrBuffer) `plusPtr` sizeOf sh
+    poke p_sh sh
+    poke p_ptrs ptrs
+
+instance Storable OutputArray where
+  sizeOf (sh, ptrs, sa) = sizeOf sh + sizeOf ptrs + sizeOf sa
+  alignment _ = 0
+  peek ptr = do
+    let p_sh   = castPtr ptr :: Ptr ShapeBuffer
+    sh         <- peek p_sh
+    let p_ptrs = (castPtr p_sh :: Ptr DevicePtrBuffer) `plusPtr` sizeOf sh
+    ptrs       <- peek p_ptrs
+    let p_sa   = (castPtr p_ptrs :: Ptr (StablePtr a)) `plusPtr` sizeOf ptrs
+    sa         <- peek p_sa
+    return (sh, ptrs, sa)
+  poke ptr (sh, ptrs, sa) = do
+    let p_sh   = castPtr ptr :: Ptr ShapeBuffer
+        p_ptrs = (castPtr p_sh :: Ptr DevicePtrBuffer) `plusPtr` sizeOf sh
+        p_sa   = (castPtr p_ptrs :: Ptr (StablePtr a)) `plusPtr` sizeOf ptrs
+    poke p_sh sh
+    poke p_ptrs ptrs
+    poke p_sa sa
+
+-- |Create an Accelerate handle given a device and a cuda context.
+--
+-- @AccHandle accelerateCreate(int device, CUcontext ctx);@
+accelerateCreate :: Device -> ForeignContext -> IO AccHandle
+accelerateCreate dev ctx = fromDeviceContext (CUDA.Device $ CInt dev) (CUDA.Context ctx) >>= newStablePtr
+
+-- |Releases all resources used by the accelerate library.
+--
+-- @void accelerateDestroy(AccHandle hndl);@
+accelerateDestroy :: AccHandle -> IO ()
+accelerateDestroy = freeStablePtr
+
+-- |Function callable from foreign code to 'free' a OutputArray returned after executing
+-- an Accelerate computation.
+--
+-- Once freed, the device pointers associated with an array are no longer valid.
+--
+-- @void freeOutput(OutputArray arr);@
+freeOutput :: Ptr OutputArray -> IO ()
+freeOutput o = do
+  (sh, dptrs, sa) <- peek o
+  free sh
+  free dptrs
+  freeStablePtr sa
+
+-- |Free a compiled accelerate program.
+--
+-- @void freeProgram(Program prg);@
+freeProgram :: StablePtr a -> IO ()
+freeProgram = freeStablePtr
+
+-- |Execute the given accelerate program with @is@ as the input and @os@ as the output.
+--
+-- @void runProgram(AccHandle hndl, AccProgram p, InputArray* is, OutputArray* os);@
+runProgram :: AccHandle -> StablePtr Afun -> Ptr InputArray -> Ptr OutputArray -> IO ()
+runProgram hndl fun input output = do
+  ctx <- deRefStablePtr hndl
+  af <- deRefStablePtr fun
+  run ctx af
+  where
+    run :: Context -> Afun -> IO ()
+    run ctx (Afun f (_ :: a) (_ :: b)) = do
+      _ <- evalCUDA ctx $ do
+        (a, _) <- marshalIn (arrays (undefined :: a)) input
+        let !b = f (toArr a)
+        marshalOut (arrays (undefined :: b)) (fromArr b) output
+      return ()
+
+    marshalIn :: ArraysR a -> Ptr InputArray -> CIO (a, Ptr InputArray)
+    marshalIn ArraysRunit  ptr = return ((), ptr)
+    marshalIn ArraysRarray ptr = do
+      (sh, ptrs) <- liftIO (peek ptr)
+      a <- arrayFromForeignData ptrs sh
+      let ptr' = plusPtr ptr (sizeOf (sh, ptrs))
+      return (a, ptr')
+    marshalIn (ArraysRpair aR1 aR2) ptr = do
+      (x, ptr')  <- marshalIn aR1 ptr
+      (y, ptr'') <- marshalIn aR2 ptr'
+      return ((x,y), ptr'')
+
+    marshalOut :: ArraysR b -> b -> Ptr OutputArray -> CIO (Ptr OutputArray)
+    marshalOut ArraysRunit  () ptr = return ptr
+    marshalOut ArraysRarray a  ptr = do
+      oarr <- mkOutput a
+      liftIO $ poke ptr oarr
+      return (plusPtr ptr (sizeOf oarr))
+      where
+        mkOutput :: forall sh e. Shape sh => Array sh e -> CIO OutputArray
+        mkOutput (Array sh adata) = do
+          let sh' = shapeToList (toElt sh :: sh)
+          shbuf <- liftIO $ mallocArray (P.length sh')
+          liftIO $ pokeArray shbuf (map fromIntegral sh')
+
+          dptrs <- devicePtrsToWordPtrs adata <$> devicePtrsOfArrayData adata
+          pbuf  <- liftIO $ mallocArray (P.length dptrs)
+          liftIO $ pokeArray pbuf dptrs
+
+          sa <- liftIO $ newStablePtr (EArray a)
+          return (shbuf, pbuf, sa)
+
+    marshalOut (ArraysRpair aR1 aR2) (x,y) ptr = do
+      ptr' <- marshalOut aR1 x ptr
+      marshalOut aR2 y ptr'
+
+
+-- |Given the 'Name' of an Accelerate function (a function of type ''Acc a -> Acc b'') generate a
+-- a function callable from foreign code with the second argument specifying it's name.
+exportAfun :: Name -> String -> Q [Dec]
+exportAfun fname ename = do
+  (VarI n ty _ _) <- reify fname
+
+  -- Generate initialisation function
+  genCompileFun n ename ty
+
+genCompileFun :: Name -> String -> Type -> Q [Dec]
+genCompileFun fname ename ty@(AppT (AppT ArrowT (AppT _ a)) (AppT _ b))
+  = sequence [sig, dec, expt]
+  where
+    initName = mkName ename
+
+    body = [| \hndl -> buildExported hndl $(varE fname) |]
+    dec  = FunD initName . (:[]) <$> cls
+    cls  = Clause [] <$> (NormalB <$> body) <*> return []
+    sig  = SigD initName <$> ety
+    expt = ForeignD <$> (ExportF cCall (nameBase initName) initName <$> ety)
+
+    ety   = [t| AccHandle -> IO (StablePtr Afun) |]
+genCompileFun _     _     _
+  = error "Invalid accelerate function"
+
+-- |Given a handle and an Accelerate function, generate an exportable version.
+buildExported :: forall a b. (Arrays a, Arrays b) => AccHandle -> (Acc a -> Acc b) -> IO (StablePtr Afun)
+buildExported hndl f = ef
+  where
+    ef :: IO (StablePtr Afun)
+    ef = do
+      ctx <- deRefStablePtr hndl
+      newStablePtr (Afun (run1In ctx f) (undefined :: a) (undefined :: b))
+
+-- Utility functions
+-- ------------------
+
+arrayFromForeignData :: forall sh e. (Shape sh, Elt e) => DevicePtrBuffer -> ShapeBuffer -> CIO (Array sh e)
+arrayFromForeignData ptrs shape = do
+   let d  = dim (ignore :: sh) -- Using ignore as using dim requires a non-dummy argument
+   let sz = eltSize (eltType (undefined :: e))
+   lst <- liftIO (peekArray d shape)
+   let sh = listToShape (map fromIntegral lst) :: sh
+   plst <- liftIO $ peekArray sz ptrs
+   let ptrs' = devicePtrsFromList (arrayElt :: ArrayEltR (EltRepr e)) plst
+   useDevicePtrs (fromElt sh) ptrs'
+
+eltSize :: TupleType e -> Int
+eltSize UnitTuple         = 0
+eltSize (SingleTuple _  ) = 1
+eltSize (PairTuple   a b) = eltSize a + eltSize b
+
diff --git a/Data/Array/Accelerate/CUDA/Foreign/Import.hs b/Data/Array/Accelerate/CUDA/Foreign/Import.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Accelerate/CUDA/Foreign/Import.hs
@@ -0,0 +1,146 @@
+{-# LANGUAGE CPP                  #-}
+{-# LANGUAGE ConstraintKinds      #-}
+{-# LANGUAGE DeriveDataTypeable   #-}
+{-# LANGUAGE FlexibleContexts     #-}
+{-# LANGUAGE FlexibleInstances    #-}
+{-# LANGUAGE GADTs                #-}
+{-# LANGUAGE PatternGuards        #-}
+{-# LANGUAGE RankNTypes           #-}
+{-# LANGUAGE ScopedTypeVariables  #-}
+{-# LANGUAGE StandaloneDeriving   #-}
+{-# LANGUAGE TypeFamilies         #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+-- |
+-- Module      : Data.Array.Accelerate.CUDA.Foreign.Import
+-- Copyright   : [2013] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell, Robert Clifton-Everest
+-- License     : BSD3
+--
+-- Maintainer  : Robert Clifton-Everest <robertce@cse.unsw.edu.au>
+-- Stability   : experimental
+-- Portability : non-portable (GHC extensions)
+--
+-- This module provides the CUDA backend's implementation of Accelerate's
+-- foreign function interface. Also provided are a series of utility functions
+-- for transferring arrays from the device to the host (and vice-versa),
+-- allocating new arrays, getting the CUDA device pointers of a given array, and
+-- executing IO actions within a CUDA context.
+--
+-- [/NOTE:/]
+--
+-- When arrays are passed to the foreign function there is no guarantee that the
+-- host side data matches the device side data. If the data is needed host side
+-- 'peekArray' or 'peekArrayAsync' must be called.
+--
+-- Arrays of tuples are represented as tuples of arrays so for example an array
+-- of type @Array DIM1 (Float, Float)@ would have two device pointers associated
+-- with it.
+--
+
+module Data.Array.Accelerate.CUDA.Foreign.Import (
+
+  -- * Backend representation
+  CUDAForeignAcc(..), canExecuteAcc,
+  CUDAForeignExp(..), canExecuteExp,
+
+  -- * Manipulating arrays
+  DevicePtrs,
+  devicePtrsOfArray,
+  indexArray, copyArray,
+  useArray,  useArrayAsync,
+  peekArray, peekArrayAsync,
+  pokeArray, pokeArrayAsync,
+  allocateArray, newArray,
+
+  -- * Running IO actions in an Accelerate context
+  CIO, liftIO, inContext, inDefaultContext
+
+) where
+
+import Data.Array.Accelerate.Type
+import Data.Array.Accelerate.CUDA.State
+import Data.Array.Accelerate.CUDA.Context
+import Data.Array.Accelerate.CUDA.Array.Sugar
+import Data.Array.Accelerate.CUDA.Array.Data
+import Data.Array.Accelerate.CUDA.Array.Prim            ( DevicePtrs )
+
+import Data.Typeable
+import Control.Exception                                ( bracket_ )
+import Control.Monad.Trans                              ( liftIO )
+
+
+-- CUDA backend representation of foreign functions
+-- ------------------------------------------------
+
+-- |CUDA foreign Acc functions are just CIO functions.
+--
+data CUDAForeignAcc as bs where
+  CUDAForeignAcc :: String                      -- name of the function
+                 -> (as -> CIO bs)              -- operation to execute
+                 -> CUDAForeignAcc as bs
+
+deriving instance Typeable2 CUDAForeignAcc
+
+instance Foreign CUDAForeignAcc where
+  strForeign (CUDAForeignAcc n _) = n
+
+-- |Gives the executable form of a foreign function if it can be executed by the
+-- CUDA backend.
+--
+canExecuteAcc
+    :: (Foreign f, Typeable as, Typeable bs)
+    => f as bs
+    -> Maybe (as -> CIO bs)
+canExecuteAcc ff
+  | Just (CUDAForeignAcc _ fun) <- cast ff
+  = Just fun
+
+  | otherwise
+  = Nothing
+
+-- |CUDA foreign Exp functions consist of a list of C header files necessary to call the function
+-- and the name of the function to call.
+--
+data CUDAForeignExp x y where
+  CUDAForeignExp :: IsScalar y
+                 => [String]                    -- header files to be imported
+                 -> String                      -- name of the foreign function
+                 -> CUDAForeignExp x y
+
+deriving instance Typeable2 CUDAForeignExp
+
+instance Foreign CUDAForeignExp where
+  strForeign (CUDAForeignExp _ n) = n
+
+-- |Gives the foreign function name as a string if it is a foreign Exp function
+-- for the CUDA backend.
+--
+canExecuteExp
+    :: forall f x y. (Foreign f, Typeable y, Typeable x)
+    => f x y
+    -> Maybe ([String], String)
+canExecuteExp ff
+  | Just (CUDAForeignExp hdr fun) <- cast ff    :: Maybe (CUDAForeignExp x y)
+  = Just (hdr, fun)
+
+  | otherwise
+  = Nothing
+
+
+-- User facing utility functions
+-- -----------------------------
+
+-- |Get the raw CUDA device pointers associated with an array
+--
+devicePtrsOfArray :: Array sh e -> CIO (DevicePtrs (EltRepr e))
+devicePtrsOfArray (Array _ adata) = devicePtrsOfArrayData adata
+
+-- |Run an IO action within the given Acclerate context
+--
+inContext :: Context -> IO a -> IO a
+inContext ctx = bracket_ (push ctx) pop
+
+-- |Run an IO action in the default Acclerate context
+--
+inDefaultContext :: IO a -> IO a
+inDefaultContext = inContext defaultContext
+
diff --git a/Data/Array/Accelerate/CUDA/Persistent.hs b/Data/Array/Accelerate/CUDA/Persistent.hs
--- a/Data/Array/Accelerate/CUDA/Persistent.hs
+++ b/Data/Array/Accelerate/CUDA/Persistent.hs
@@ -67,8 +67,8 @@
 -- Interface -------------------------------------------------------------------
 -- ---------                                                                  --
 
-data KernelTable = KT {-# UNPACK #-} !ProgramCache      -- first level cache
-                      {-# UNPACK #-} !PersistentCache   -- second level cache
+data KernelTable = KT {-# UNPACK #-} !ProgramCache      -- first level, in-memory cache
+                      {-# UNPACK #-} !PersistentCache   -- second level, on-disk cache
 
 new :: IO KernelTable
 new = do
@@ -131,11 +131,11 @@
 -- Unload a kernel module from the specified context
 --
 module_finalizer :: Weak CUDA.Context -> KernelKey -> CUDA.Module -> IO ()
-module_finalizer weak_ctx (_,key) mdl = do
+module_finalizer weak_ctx key mdl = do
   mc <- deRefWeak weak_ctx
   case mc of
-    Nothing     -> D.message D.dump_gc ("gc: finalise module/dead context: "       ++ show key)
-    Just ctx    -> D.message D.dump_gc ("gc: finalise module: " ++ show ctx ++ "/" ++ show key)
+    Nothing     -> D.message D.dump_gc ("gc: finalise module/dead context: " ++ cacheFilePath key)
+    Just ctx    -> D.message D.dump_gc ("gc: finalise module: "              ++ cacheFilePath key)
                 >> bracket_ (CUDA.push ctx) CUDA.pop (CUDA.unload mdl)
 
 
diff --git a/Data/Array/Accelerate/CUDA/State.hs b/Data/Array/Accelerate/CUDA/State.hs
--- a/Data/Array/Accelerate/CUDA/State.hs
+++ b/Data/Array/Accelerate/CUDA/State.hs
@@ -23,21 +23,19 @@
   CIO, Context, evalCUDA,
 
   -- Querying execution state
-  defaultContext, deviceProperties, activeContext, kernelTable, memoryTable
+  defaultContext, deviceProperties, activeContext, kernelTable, memoryTable, streamReservoir,
 
 ) where
 
 -- friends
 import Data.Array.Accelerate.CUDA.Context
 import Data.Array.Accelerate.CUDA.Debug                 ( message, dump_gc )
-import Data.Array.Accelerate.CUDA.Persistent            as KT
-import Data.Array.Accelerate.CUDA.Array.Table           as MT
+import Data.Array.Accelerate.CUDA.Persistent            as KT ( KernelTable, new )
+import Data.Array.Accelerate.CUDA.Array.Table           as MT ( MemoryTable, new )
+import Data.Array.Accelerate.CUDA.Execute.Stream        as ST ( Reservoir, new )
 import Data.Array.Accelerate.CUDA.Analysis.Device
 
 -- library
-#if !MIN_VERSION_base(4,6,0)
-import Prelude                                          hiding ( catch )
-#endif
 import Control.Applicative                              ( Applicative )
 import Control.Concurrent                               ( runInBoundThread )
 import Control.Exception                                ( catch, bracket_ )
@@ -61,7 +59,8 @@
 --
 data State = State {
     memoryTable         :: {-# UNPACK #-} !MemoryTable,                 -- host/device memory associations
-    kernelTable         :: {-# UNPACK #-} !KernelTable                  -- compiled kernel object code
+    kernelTable         :: {-# UNPACK #-} !KernelTable,                 -- compiled kernel object code
+    streamReservoir     :: {-# UNPACK #-} !Reservoir                    -- kernel execution streams
   }
 
 newtype CIO a = CIO {
@@ -105,7 +104,8 @@
   $ do  message dump_gc "gc: initialise CUDA state"
         mtb     <- keepAlive =<< MT.new
         ktb     <- keepAlive =<< KT.new
-        return  $! State mtb ktb
+        rsv     <- keepAlive =<< ST.new
+        return  $! State mtb ktb rsv
 
 
 -- Select and initialise a default CUDA device, and create a new execution
diff --git a/accelerate-cuda.cabal b/accelerate-cuda.cabal
--- a/accelerate-cuda.cabal
+++ b/accelerate-cuda.cabal
@@ -1,7 +1,7 @@
 Name:                   accelerate-cuda
-Version:                0.13.0.4
+Version:                0.14.0.0
 Cabal-version:          >= 1.6
-Tested-with:            GHC >= 7.4
+Tested-with:            GHC == 7.6.*
 Build-type:             Custom
 
 Synopsis:               Accelerate backend for NVIDIA GPUs
@@ -40,27 +40,20 @@
 Category:               Compilers/Interpreters, Concurrency, Data, Parallelism
 Stability:              Experimental
 
--- We require 'accelerate_cuda_shape.h' to be in this list so that it is copied
--- as part of the installation, although 'cabal sdist' does not grok that it is
--- generated by the configure script.
---
-Data-files:             cubits/accelerate_cuda_extras.h
+Data-files:             cubits/accelerate_cuda.h
+                        cubits/accelerate_cuda_assert.h
                         cubits/accelerate_cuda_function.h
-                        cubits/accelerate_cuda_shape.h
-                        cubits/accelerate_cuda_stencil.h
                         cubits/accelerate_cuda_texture.h
                         cubits/accelerate_cuda_type.h
-                        cubits/accelerate_cuda_util.h
+                        include/AccFFI.h
 
 Extra-tmp-files:        config.status
                         config.log
                         autom4te.cache
                         accelerate-cuda.buildinfo               -- generated by configure
-                        cubits/accelerate_cuda_shape.h          -- generated by configure
 
 Extra-source-files:     configure
                         accelerate-cuda.buildinfo.in
-                        cubits/accelerate_cuda_shape.h.in
                         include/accelerate.h
 
 Flag debug
@@ -97,31 +90,33 @@
   Default:              False
 
 Library
-  include-dirs:         include
+  Include-Dirs:         include
 
-  build-depends:        accelerate              == 0.13.*,
-                        base                    == 4.*,
-                        array                   >= 0.3      && < 0.5,
-                        binary                  >= 0.5      && < 0.7,
-                        bytestring              >= 0.9      && < 0.11,
-                        cryptohash              >= 0.7      && < 0.10,
-                        cuda                    >= 0.5.0.2  && < 0.6,
-                        directory               >= 1.0      && < 1.3,
-                        fclabels                >= 1.0      && < 1.2,
-                        filepath                >= 1.0      && < 1.4,
-                        hashable                >= 1.1      && < 1.3,
-                        hashtables              >= 1.0      && < 1.2,
-                        language-c-quote        >= 0.4.4    && < 0.8,
-                        mainland-pretty         >= 0.2      && < 0.3,
-                        mtl                     >= 2.0      && < 2.2,
-                        old-time                >= 1.0      && < 1.2,
-                        pretty                  >= 1.0      && < 1.2,
-                        process                 >= 1.0      && < 1.2,
-                        SafeSemaphore           >= 0.9      && < 0.10,
-                        srcloc                  >= 0.2      && < 0.5,
-                        text                    >= 0.11     && < 0.12,
-                        transformers            >= 0.2      && < 0.4,
-                        unordered-containers    >= 0.1.4    && < 0.3
+  Build-depends:        accelerate              == 0.14.*,
+                        array                   >= 0.3,
+                        base                    == 4.6.*,
+                        binary                  >= 0.5          && < 0.7,
+                            -- binary 0.7 untested
+                        bytestring              >= 0.9          && < 0.11,
+                        cryptohash              >= 0.7          && < 0.12,
+                        cuda                    >= 0.5.1.1      && < 0.6,
+                        directory               >= 1.0,
+                        fclabels                >= 2.0          && < 2.1,
+                        filepath                >= 1.0,
+                        hashable                >= 1.1          && < 1.3,
+                        hashtables              >= 1.0.1        && < 1.2,
+                        language-c-quote        >= 0.4.4        && < 0.8,
+                        mainland-pretty         >= 0.2          && < 0.3,
+                        mtl                     >= 2.0          && < 2.2,
+                        old-time                >= 1.0,
+                        pretty                  >= 1.0,
+                        process                 >= 1.0,
+                        SafeSemaphore           >= 0.9          && < 0.10,
+                        srcloc                  >= 0.2          && < 0.5,
+                        text                    >= 0.11         && < 0.12,
+                        template-haskell        >= 2.2,
+                        transformers            >= 0.2          && < 0.4,
+                        unordered-containers    >= 0.1.4        && < 0.3
 
   if os(windows)
     cpp-options:        -DWIN32
@@ -150,11 +145,16 @@
                         Data.Array.Accelerate.CUDA.CodeGen.PrefixSum
                         Data.Array.Accelerate.CUDA.CodeGen.Reduction
                         Data.Array.Accelerate.CUDA.CodeGen.Stencil
+                        Data.Array.Accelerate.CUDA.CodeGen.Stencil.Extra
                         Data.Array.Accelerate.CUDA.CodeGen.Type
                         Data.Array.Accelerate.CUDA.Compile
                         Data.Array.Accelerate.CUDA.Context
                         Data.Array.Accelerate.CUDA.Debug
                         Data.Array.Accelerate.CUDA.Execute
+                        Data.Array.Accelerate.CUDA.Execute.Event
+                        Data.Array.Accelerate.CUDA.Execute.Stream
+                        Data.Array.Accelerate.CUDA.Foreign.Export
+                        Data.Array.Accelerate.CUDA.Foreign.Import
                         Data.Array.Accelerate.CUDA.FullList
                         Data.Array.Accelerate.CUDA.Persistent
                         Data.Array.Accelerate.CUDA.State
@@ -176,6 +176,8 @@
   ghc-options:          -O2
                         -Wall
                         -fwarn-tabs
+
+  ghc-prof-options:     -auto-all
 
   -- Don't add the extensions list here. Instead, place individual LANGUAGE
   -- pragmas in the files that require a specific extension. This means the
diff --git a/configure b/configure
--- a/configure
+++ b/configure
@@ -589,7 +589,6 @@
 LIBOBJS
 cpp_flags
 ghc_flags
-type_hs_int
 OBJEXT
 EXEEXT
 ac_ct_CXX
@@ -1743,7 +1742,7 @@
 
 
 
-ac_config_files="$ac_config_files accelerate-cuda.buildinfo cubits/accelerate_cuda_shape.h"
+ac_config_files="$ac_config_files accelerate-cuda.buildinfo"
 
 
 
@@ -2409,7 +2408,6 @@
     ghc_flags="$ghc_flags -optP$def"
 
 
-         # array index type in accelerate_cuda_shape.h.in
 
 
 cat >confcache <<\_ACEOF
@@ -3118,7 +3116,6 @@
 do
   case $ac_config_target in
     "accelerate-cuda.buildinfo") CONFIG_FILES="$CONFIG_FILES accelerate-cuda.buildinfo" ;;
-    "cubits/accelerate_cuda_shape.h") CONFIG_FILES="$CONFIG_FILES cubits/accelerate_cuda_shape.h" ;;
 
   *) as_fn_error $? "invalid argument: \`$ac_config_target'" "$LINENO" 5;;
   esac
diff --git a/cubits/accelerate_cuda.h b/cubits/accelerate_cuda.h
new file mode 100644
--- /dev/null
+++ b/cubits/accelerate_cuda.h
@@ -0,0 +1,22 @@
+/* -----------------------------------------------------------------------------
+ *
+ * Module      : Accelerate-CUDA
+ * Copyright   : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee
+ *               [2009..2012] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell
+ * License     : BSD3
+ *
+ * Maintainer  : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>
+ * Stability   : experimental
+ *
+ * ---------------------------------------------------------------------------*/
+
+#ifndef __ACCELERATE_CUDA_H__
+#define __ACCELERATE_CUDA_H__
+
+#include "accelerate_cuda_assert.h"
+#include "accelerate_cuda_function.h"
+#include "accelerate_cuda_texture.h"
+#include "accelerate_cuda_type.h"
+
+#endif
+
diff --git a/cubits/accelerate_cuda_assert.h b/cubits/accelerate_cuda_assert.h
new file mode 100644
--- /dev/null
+++ b/cubits/accelerate_cuda_assert.h
@@ -0,0 +1,28 @@
+/* -----------------------------------------------------------------------------
+ *
+ * Module      : Assert
+ * Copyright   : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee
+ *               [2009..2012] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell
+ * License     : BSD3
+ *
+ * Maintainer  : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>
+ * Stability   : experimental
+ *
+ * ---------------------------------------------------------------------------*/
+
+#ifndef __ACCELERATE_CUDA_ASSERT_H__
+#define __ACCELERATE_CUDA_ASSERT_H__
+
+#include <assert.h>
+#include <cuda_runtime.h>
+
+/*
+ * assert() is only supported for devices of compute capability 2.0 and higher
+ */
+#if defined(__CUDA_ARCH__) && (__CUDA_ARCH__ < 200)
+#undef  assert
+#define assert(arg)
+#endif
+
+#endif  // __ACCELERATE_CUDA_ASSERT_H__
+
diff --git a/cubits/accelerate_cuda_extras.h b/cubits/accelerate_cuda_extras.h
deleted file mode 100644
--- a/cubits/accelerate_cuda_extras.h
+++ /dev/null
@@ -1,23 +0,0 @@
-/* -----------------------------------------------------------------------------
- *
- * Module      : Extras
- * Copyright   : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee
- *               [2009..2012] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell
- * License     : BSD3
- *
- * Maintainer  : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>
- * Stability   : experimental
- *
- * ---------------------------------------------------------------------------*/
-
-#ifndef __ACCELERATE_CUDA_EXTRAS_H__
-#define __ACCELERATE_CUDA_EXTRAS_H__
-
-#include "accelerate_cuda_function.h"
-#include "accelerate_cuda_shape.h"
-#include "accelerate_cuda_stencil.h"
-#include "accelerate_cuda_texture.h"
-#include "accelerate_cuda_type.h"
-
-#endif
-
diff --git a/cubits/accelerate_cuda_shape.h b/cubits/accelerate_cuda_shape.h
deleted file mode 100644
--- a/cubits/accelerate_cuda_shape.h
+++ /dev/null
@@ -1,338 +0,0 @@
-/* -----------------------------------------------------------------------------
- *
- * Module      : Shape
- * Copyright   : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee
- *               [2009..2012] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell
- * License     : BSD3
- *
- * Maintainer  : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>
- * Stability   : experimental
- *
- * ---------------------------------------------------------------------------*/
-
-#ifndef __ACCELERATE_CUDA_SHAPE_H__
-#define __ACCELERATE_CUDA_SHAPE_H__
-
-#include <cuda_runtime.h>
-#include "accelerate_cuda_type.h"
-
-/*
- * For the time being, the D.A.A.CUDA.Execute converts shape components from
- * Haskell Ints (which may be 32- or 64-bits) into CUDA integers which are
- * 32-bit. This is to avoid undue register pressure. Lift this restriction if
- * future hardware gains better 64-bit support and/or we need to access very
- * large arrays.
- *
- * typedef Int64                             Ix;
- */
-typedef Int32                                     Ix;
-typedef void*                                     DIM0;
-typedef Ix                                        DIM1;
-typedef struct { Ix a1,a0; }                      DIM2;
-typedef struct { Ix a2,a1,a0; }                   DIM3;
-typedef struct { Ix a3,a2,a1,a0; }                DIM4;
-typedef struct { Ix a4,a3,a2,a1,a0; }             DIM5;
-typedef struct { Ix a5,a4,a3,a2,a1,a0; }          DIM6;
-typedef struct { Ix a6,a5,a4,a3,a2,a1,a0; }       DIM7;
-typedef struct { Ix a7,a6,a5,a4,a3,a2,a1,a0; }    DIM8;
-typedef struct { Ix a8,a7,a6,a5,a4,a3,a2,a1,a0; } DIM9;
-
-#ifdef __cplusplus
-
-/* -----------------------------------------------------------------------------
- * Shape construction and destruction
- */
-
-/*
- * Convert the individual dimensions of a linear array into a shape
- */
-static __inline__ __device__ DIM0 shape()
-{
-    return NULL;
-}
-
-static __inline__ __device__ DIM1 shape(const Ix a)
-{
-    return a;
-}
-
-static __inline__ __device__ DIM2 shape(const Ix b, const Ix a)
-{
-    DIM2 sh = { b, a };
-    return sh;
-}
-
-static __inline__ __device__ DIM3 shape(const Ix c, const Ix b, const Ix a)
-{
-    DIM3 sh = { c, b, a };
-    return sh;
-}
-
-static __inline__ __device__ DIM4 shape(const Ix d, const Ix c, const Ix b, const Ix a)
-{
-    DIM4 sh = { d, c, b, a };
-    return sh;
-}
-
-static __inline__ __device__ DIM5 shape(const Ix e, const Ix d, const Ix c, const Ix b, const Ix a)
-{
-    DIM5 sh = { e, d, c, b, a };
-    return sh;
-}
-
-static __inline__ __device__ DIM6 shape(const Ix f, const Ix e, const Ix d, const Ix c, const Ix b, const Ix a)
-{
-    DIM6 sh = { f, e, d, c, b, a };
-    return sh;
-}
-
-static __inline__ __device__ DIM7 shape(const Ix g, const Ix f, const Ix e, const Ix d, const Ix c, const Ix b, const Ix a)
-{
-    DIM7 sh = { g, f, e, d, c, b, a };
-    return sh;
-}
-
-static __inline__ __device__ DIM8 shape(const Ix h, const Ix g, const Ix f, const Ix e, const Ix d, const Ix c, const Ix b, const Ix a)
-{
-    DIM8 sh = { h, g, f, e, d, c, b, a };
-    return sh;
-}
-
-static __inline__ __device__ DIM9 shape(const Ix i, const Ix h, const Ix g, const Ix f, const Ix e, const Ix d, const Ix c, const Ix b, const Ix a)
-{
-    DIM9 sh = { i, h, g, f, e, d, c, b, a };
-    return sh;
-}
-
-/*
- * Yield the inner-most dimension of a shape only
- */
-template <typename Shape>
-static __inline__ __device__ Ix indexHead(const Shape ix)
-{
-    return ix.a0;
-}
-
-template <>
-static __inline__ __device__ Ix indexHead(const DIM0 ix)
-{
-    return 0;
-}
-
-template <>
-static __inline__ __device__ Ix indexHead(const DIM1 ix)
-{
-    return ix;
-}
-
-
-/*
- * Yield all but the inner-most dimension of a shape
- */
-static __inline__ __device__ DIM0 indexTail(const DIM1 ix)
-{
-    return 0;
-}
-
-static __inline__ __device__ DIM1 indexTail(const DIM2 ix)
-{
-    return ix.a1;
-}
-
-static __inline__ __device__ DIM2 indexTail(const DIM3 ix)
-{
-    return shape(ix.a2, ix.a1);
-}
-
-static __inline__ __device__ DIM3 indexTail(const DIM4 ix)
-{
-    return shape(ix.a3, ix.a2, ix.a1);
-}
-
-static __inline__ __device__ DIM4 indexTail(const DIM5 ix)
-{
-    return shape(ix.a4, ix.a3, ix.a2, ix.a1);
-}
-
-static __inline__ __device__ DIM5 indexTail(const DIM6 ix)
-{
-    return shape(ix.a5, ix.a4, ix.a3, ix.a2, ix.a1);
-}
-
-static __inline__ __device__ DIM6 indexTail(const DIM7 ix)
-{
-    return shape(ix.a6, ix.a5, ix.a4, ix.a3, ix.a2, ix.a1);
-}
-
-static __inline__ __device__ DIM7 indexTail(const DIM8 ix)
-{
-    return shape(ix.a7, ix.a6, ix.a5, ix.a4, ix.a3, ix.a2, ix.a1);
-}
-
-static __inline__ __device__ DIM8 indexTail(const DIM9 ix)
-{
-    return shape(ix.a8, ix.a7, ix.a6, ix.a5, ix.a4, ix.a3, ix.a2, ix.a1);
-}
-
-
-/* -----------------------------------------------------------------------------
- * Shape methods
- */
-
-/*
- * Number of dimensions of a shape
- */
-template <typename Shape>
-static __inline__ __device__ int dim(const Shape sh)
-{
-    return dim(indexTail(sh)) + 1;
-}
-
-template <>
-static __inline__ __device__ int dim(const DIM0 sh)
-{
-    return 0;
-}
-
-
-/*
- * Yield the total number of elements in a shape
- */
-template <typename Shape>
-static __inline__ __device__ int size(const Shape sh)
-{
-    return size(indexTail(sh)) * indexHead(sh);
-}
-
-template <>
-static __inline__ __device__ int size(const DIM0 sh)
-{
-    return 1;
-}
-
-
-/*
- * Add an index to the head of a shape
- */
-static __inline__ __device__ DIM1 indexCons(const DIM0 sh, const Ix ix)
-{
-    return shape(ix);
-}
-
-static __inline__ __device__ DIM2 indexCons(const DIM1 sh, const Ix ix)
-{
-    return shape(sh, ix);
-}
-
-static __inline__ __device__ DIM3 indexCons(const DIM2 sh, const Ix ix)
-{
-    return shape(sh.a1, sh.a0, ix);
-}
-
-static __inline__ __device__ DIM4 indexCons(const DIM3 sh, const Ix ix)
-{
-    return shape(sh.a2, sh.a1, sh.a0, ix);
-}
-
-static __inline__ __device__ DIM5 indexCons(const DIM4 sh, const Ix ix)
-{
-    return shape(sh.a3, sh.a2, sh.a1, sh.a0, ix);
-}
-
-static __inline__ __device__ DIM6 indexCons(const DIM5 sh, const Ix ix)
-{
-    return shape(sh.a4, sh.a3, sh.a2, sh.a1, sh.a0, ix);
-}
-
-static __inline__ __device__ DIM7 indexCons(const DIM6 sh, const Ix ix)
-{
-    return shape(sh.a5, sh.a4, sh.a3, sh.a2, sh.a1, sh.a0, ix);
-}
-
-static __inline__ __device__ DIM8 indexCons(const DIM7 sh, const Ix ix)
-{
-    return shape(sh.a6, sh.a5, sh.a4, sh.a3, sh.a2, sh.a1, sh.a0, ix);
-}
-
-static __inline__ __device__ DIM9 indexCons(const DIM8 sh, const Ix ix)
-{
-    return shape(sh.a7, sh.a6, sh.a5, sh.a4, sh.a3, sh.a2, sh.a1, sh.a0, ix);
-}
-
-
-/*
- * Yield the index position in a linear, row-major representation of the array.
- * First argument is the shape of the array, the second the index
- */
-template <typename Shape>
-static __inline__ __device__ Ix toIndex(const Shape sh, const Shape ix)
-{
-    return toIndex(indexTail(sh), indexTail(ix)) * indexHead(sh) + indexHead(ix);
-}
-
-template <>
-static __inline__ __device__ Ix toIndex(const DIM0 sh, const DIM0 ix)
-{
-    return 0;
-}
-
-template <>
-static __inline__ __device__ Ix toIndex(const DIM1 sh, const DIM1 ix)
-{
-    return ix;
-}
-
-
-/*
- * Inverse of 'toIndex'
- */
-template <typename Shape>
-static __inline__ __device__ Shape fromIndex(const Shape sh, const Ix ix)
-{
-    const Ix d = indexHead(sh);
-    return indexCons(fromIndex(indexTail(sh), ix / d), ix % d);
-}
-
-template <>
-static __inline__ __device__ DIM0 fromIndex(const DIM0 sh, const Ix ix)
-{
-    return 0;
-}
-
-template <>
-static __inline__ __device__ DIM1 fromIndex(const DIM1 sh, const Ix ix)
-{
-    return ix;
-}
-
-
-/*
- * Test for the magic index `ignore`
- */
-template <typename Shape>
-static __inline__ __device__ int ignore(const Shape ix)
-{
-    return indexHead(ix) == -1 && ignore(indexTail(ix));
-}
-
-template <>
-static __inline__ __device__ int ignore(const DIM0 ix)
-{
-    return 1;
-}
-
-
-#else
-
-static __inline__ __device__ int dim(const Ix sh);
-static __inline__ __device__ int size(const Ix sh);
-static __inline__ __device__ int shape(const Ix sh);
-static __inline__ __device__ int toIndex(const Ix sh, const Ix ix);
-static __inline__ __device__ int fromIndex(const Ix sh, const Ix ix);
-static __inline__ __device__ int indexHead(const Ix ix);
-static __inline__ __device__ int indexTail(const Ix ix);
-static __inline__ __device__ int indexCons(const Ix sh, const Ix ix);
-
-#endif  // __cplusplus
-#endif  // __ACCELERATE_CUDA_SHAPE_H__
-
diff --git a/cubits/accelerate_cuda_shape.h.in b/cubits/accelerate_cuda_shape.h.in
deleted file mode 100644
--- a/cubits/accelerate_cuda_shape.h.in
+++ /dev/null
@@ -1,338 +0,0 @@
-/* -----------------------------------------------------------------------------
- *
- * Module      : Shape
- * Copyright   : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee
- *               [2009..2012] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell
- * License     : BSD3
- *
- * Maintainer  : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>
- * Stability   : experimental
- *
- * ---------------------------------------------------------------------------*/
-
-#ifndef __ACCELERATE_CUDA_SHAPE_H__
-#define __ACCELERATE_CUDA_SHAPE_H__
-
-#include <cuda_runtime.h>
-#include "accelerate_cuda_type.h"
-
-/*
- * For the time being, the D.A.A.CUDA.Execute converts shape components from
- * Haskell Ints (which may be 32- or 64-bits) into CUDA integers which are
- * 32-bit. This is to avoid undue register pressure. Lift this restriction if
- * future hardware gains better 64-bit support and/or we need to access very
- * large arrays.
- *
- * typedef @type_hs_int@                             Ix;
- */
-typedef Int32                                     Ix;
-typedef void*                                     DIM0;
-typedef Ix                                        DIM1;
-typedef struct { Ix a1,a0; }                      DIM2;
-typedef struct { Ix a2,a1,a0; }                   DIM3;
-typedef struct { Ix a3,a2,a1,a0; }                DIM4;
-typedef struct { Ix a4,a3,a2,a1,a0; }             DIM5;
-typedef struct { Ix a5,a4,a3,a2,a1,a0; }          DIM6;
-typedef struct { Ix a6,a5,a4,a3,a2,a1,a0; }       DIM7;
-typedef struct { Ix a7,a6,a5,a4,a3,a2,a1,a0; }    DIM8;
-typedef struct { Ix a8,a7,a6,a5,a4,a3,a2,a1,a0; } DIM9;
-
-#ifdef __cplusplus
-
-/* -----------------------------------------------------------------------------
- * Shape construction and destruction
- */
-
-/*
- * Convert the individual dimensions of a linear array into a shape
- */
-static __inline__ __device__ DIM0 shape()
-{
-    return NULL;
-}
-
-static __inline__ __device__ DIM1 shape(const Ix a)
-{
-    return a;
-}
-
-static __inline__ __device__ DIM2 shape(const Ix b, const Ix a)
-{
-    DIM2 sh = { b, a };
-    return sh;
-}
-
-static __inline__ __device__ DIM3 shape(const Ix c, const Ix b, const Ix a)
-{
-    DIM3 sh = { c, b, a };
-    return sh;
-}
-
-static __inline__ __device__ DIM4 shape(const Ix d, const Ix c, const Ix b, const Ix a)
-{
-    DIM4 sh = { d, c, b, a };
-    return sh;
-}
-
-static __inline__ __device__ DIM5 shape(const Ix e, const Ix d, const Ix c, const Ix b, const Ix a)
-{
-    DIM5 sh = { e, d, c, b, a };
-    return sh;
-}
-
-static __inline__ __device__ DIM6 shape(const Ix f, const Ix e, const Ix d, const Ix c, const Ix b, const Ix a)
-{
-    DIM6 sh = { f, e, d, c, b, a };
-    return sh;
-}
-
-static __inline__ __device__ DIM7 shape(const Ix g, const Ix f, const Ix e, const Ix d, const Ix c, const Ix b, const Ix a)
-{
-    DIM7 sh = { g, f, e, d, c, b, a };
-    return sh;
-}
-
-static __inline__ __device__ DIM8 shape(const Ix h, const Ix g, const Ix f, const Ix e, const Ix d, const Ix c, const Ix b, const Ix a)
-{
-    DIM8 sh = { h, g, f, e, d, c, b, a };
-    return sh;
-}
-
-static __inline__ __device__ DIM9 shape(const Ix i, const Ix h, const Ix g, const Ix f, const Ix e, const Ix d, const Ix c, const Ix b, const Ix a)
-{
-    DIM9 sh = { i, h, g, f, e, d, c, b, a };
-    return sh;
-}
-
-/*
- * Yield the inner-most dimension of a shape only
- */
-template <typename Shape>
-static __inline__ __device__ Ix indexHead(const Shape ix)
-{
-    return ix.a0;
-}
-
-template <>
-static __inline__ __device__ Ix indexHead(const DIM0 ix)
-{
-    return 0;
-}
-
-template <>
-static __inline__ __device__ Ix indexHead(const DIM1 ix)
-{
-    return ix;
-}
-
-
-/*
- * Yield all but the inner-most dimension of a shape
- */
-static __inline__ __device__ DIM0 indexTail(const DIM1 ix)
-{
-    return 0;
-}
-
-static __inline__ __device__ DIM1 indexTail(const DIM2 ix)
-{
-    return ix.a1;
-}
-
-static __inline__ __device__ DIM2 indexTail(const DIM3 ix)
-{
-    return shape(ix.a2, ix.a1);
-}
-
-static __inline__ __device__ DIM3 indexTail(const DIM4 ix)
-{
-    return shape(ix.a3, ix.a2, ix.a1);
-}
-
-static __inline__ __device__ DIM4 indexTail(const DIM5 ix)
-{
-    return shape(ix.a4, ix.a3, ix.a2, ix.a1);
-}
-
-static __inline__ __device__ DIM5 indexTail(const DIM6 ix)
-{
-    return shape(ix.a5, ix.a4, ix.a3, ix.a2, ix.a1);
-}
-
-static __inline__ __device__ DIM6 indexTail(const DIM7 ix)
-{
-    return shape(ix.a6, ix.a5, ix.a4, ix.a3, ix.a2, ix.a1);
-}
-
-static __inline__ __device__ DIM7 indexTail(const DIM8 ix)
-{
-    return shape(ix.a7, ix.a6, ix.a5, ix.a4, ix.a3, ix.a2, ix.a1);
-}
-
-static __inline__ __device__ DIM8 indexTail(const DIM9 ix)
-{
-    return shape(ix.a8, ix.a7, ix.a6, ix.a5, ix.a4, ix.a3, ix.a2, ix.a1);
-}
-
-
-/* -----------------------------------------------------------------------------
- * Shape methods
- */
-
-/*
- * Number of dimensions of a shape
- */
-template <typename Shape>
-static __inline__ __device__ int dim(const Shape sh)
-{
-    return dim(indexTail(sh)) + 1;
-}
-
-template <>
-static __inline__ __device__ int dim(const DIM0 sh)
-{
-    return 0;
-}
-
-
-/*
- * Yield the total number of elements in a shape
- */
-template <typename Shape>
-static __inline__ __device__ int size(const Shape sh)
-{
-    return size(indexTail(sh)) * indexHead(sh);
-}
-
-template <>
-static __inline__ __device__ int size(const DIM0 sh)
-{
-    return 1;
-}
-
-
-/*
- * Add an index to the head of a shape
- */
-static __inline__ __device__ DIM1 indexCons(const DIM0 sh, const Ix ix)
-{
-    return shape(ix);
-}
-
-static __inline__ __device__ DIM2 indexCons(const DIM1 sh, const Ix ix)
-{
-    return shape(sh, ix);
-}
-
-static __inline__ __device__ DIM3 indexCons(const DIM2 sh, const Ix ix)
-{
-    return shape(sh.a1, sh.a0, ix);
-}
-
-static __inline__ __device__ DIM4 indexCons(const DIM3 sh, const Ix ix)
-{
-    return shape(sh.a2, sh.a1, sh.a0, ix);
-}
-
-static __inline__ __device__ DIM5 indexCons(const DIM4 sh, const Ix ix)
-{
-    return shape(sh.a3, sh.a2, sh.a1, sh.a0, ix);
-}
-
-static __inline__ __device__ DIM6 indexCons(const DIM5 sh, const Ix ix)
-{
-    return shape(sh.a4, sh.a3, sh.a2, sh.a1, sh.a0, ix);
-}
-
-static __inline__ __device__ DIM7 indexCons(const DIM6 sh, const Ix ix)
-{
-    return shape(sh.a5, sh.a4, sh.a3, sh.a2, sh.a1, sh.a0, ix);
-}
-
-static __inline__ __device__ DIM8 indexCons(const DIM7 sh, const Ix ix)
-{
-    return shape(sh.a6, sh.a5, sh.a4, sh.a3, sh.a2, sh.a1, sh.a0, ix);
-}
-
-static __inline__ __device__ DIM9 indexCons(const DIM8 sh, const Ix ix)
-{
-    return shape(sh.a7, sh.a6, sh.a5, sh.a4, sh.a3, sh.a2, sh.a1, sh.a0, ix);
-}
-
-
-/*
- * Yield the index position in a linear, row-major representation of the array.
- * First argument is the shape of the array, the second the index
- */
-template <typename Shape>
-static __inline__ __device__ Ix toIndex(const Shape sh, const Shape ix)
-{
-    return toIndex(indexTail(sh), indexTail(ix)) * indexHead(sh) + indexHead(ix);
-}
-
-template <>
-static __inline__ __device__ Ix toIndex(const DIM0 sh, const DIM0 ix)
-{
-    return 0;
-}
-
-template <>
-static __inline__ __device__ Ix toIndex(const DIM1 sh, const DIM1 ix)
-{
-    return ix;
-}
-
-
-/*
- * Inverse of 'toIndex'
- */
-template <typename Shape>
-static __inline__ __device__ Shape fromIndex(const Shape sh, const Ix ix)
-{
-    const Ix d = indexHead(sh);
-    return indexCons(fromIndex(indexTail(sh), ix / d), ix % d);
-}
-
-template <>
-static __inline__ __device__ DIM0 fromIndex(const DIM0 sh, const Ix ix)
-{
-    return 0;
-}
-
-template <>
-static __inline__ __device__ DIM1 fromIndex(const DIM1 sh, const Ix ix)
-{
-    return ix;
-}
-
-
-/*
- * Test for the magic index `ignore`
- */
-template <typename Shape>
-static __inline__ __device__ int ignore(const Shape ix)
-{
-    return indexHead(ix) == -1 && ignore(indexTail(ix));
-}
-
-template <>
-static __inline__ __device__ int ignore(const DIM0 ix)
-{
-    return 1;
-}
-
-
-#else
-
-static __inline__ __device__ int dim(const Ix sh);
-static __inline__ __device__ int size(const Ix sh);
-static __inline__ __device__ int shape(const Ix sh);
-static __inline__ __device__ int toIndex(const Ix sh, const Ix ix);
-static __inline__ __device__ int fromIndex(const Ix sh, const Ix ix);
-static __inline__ __device__ int indexHead(const Ix ix);
-static __inline__ __device__ int indexTail(const Ix ix);
-static __inline__ __device__ int indexCons(const Ix sh, const Ix ix);
-
-#endif  // __cplusplus
-#endif  // __ACCELERATE_CUDA_SHAPE_H__
-
diff --git a/cubits/accelerate_cuda_stencil.h b/cubits/accelerate_cuda_stencil.h
deleted file mode 100644
--- a/cubits/accelerate_cuda_stencil.h
+++ /dev/null
@@ -1,93 +0,0 @@
-/* -----------------------------------------------------------------------------
- *
- * Module      : Stencil
- * Copyright   : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee
- *               [2009..2012] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell
- * License     : BSD3
- *
- * Maintainer  : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>
- * Stability   : experimental
- *
- * ---------------------------------------------------------------------------*/
-
-#ifndef __ACCELERATE_CUDA_STENCIL_H__
-#define __ACCELERATE_CUDA_STENCIL_H__
-
-#include <accelerate_cuda_shape.h>
-
-#ifdef __cplusplus
-
-/*
- * Test if an index lies within the boundaries of a shape
- */
-template <typename Shape>
-static __inline__ __device__ int inRange(const Shape sh, const Shape ix)
-{
-    return inRange(indexHead(sh), indexHead(ix)) && inRange(indexTail(sh), indexTail(ix));
-}
-
-template <>
-static __inline__ __device__ int inRange(const DIM1 sz, const DIM1 i)
-{
-    return i >= 0 && i < sz;
-}
-
-template <>
-static __inline__ __device__ int inRange(const DIM0 sz, const DIM0 i)
-{
-    return i == 0;
-}
-
-
-/*
- * Boundary condition handlers
- */
-template <typename Shape>
-static __inline__ __device__ Shape clamp(const Shape sh, const Shape ix)
-{
-    return indexCons( clamp(indexTail(sh), indexTail(ix))
-                    , clamp(indexHead(sh), indexHead(ix)) );
-}
-
-template <>
-static __inline__ __device__ DIM1 clamp(const DIM1 sz, const DIM1 i)
-{
-    // CUDA-4.0 does not include 64-bit min/max functions for Fermi (?)
-    return max(0, min(i, sz-1));
-}
-
-
-template <typename Shape>
-static __inline__ __device__ Shape mirror(const Shape sh, const Shape ix)
-{
-    return indexCons( mirror(indexTail(sh), indexTail(ix))
-                    , mirror(indexHead(sh), indexHead(ix)) );
-}
-
-template <>
-static __inline__ __device__ DIM1 mirror(const DIM1 sz, const DIM1 i)
-{
-    if      (i <  0)  return -i;
-    else if (i >= sz) return sz - (i-sz+2);
-    else              return i;
-}
-
-
-template <typename Shape>
-static __inline__ __device__ Shape wrap(const Shape sh, const Shape ix)
-{
-    return indexCons( wrap(indexTail(sh), indexTail(ix))
-                    , wrap(indexHead(sh), indexHead(ix)) );
-}
-
-template <>
-static __inline__ __device__ DIM1 wrap(const DIM1 sz, const DIM1 i)
-{
-    if      (i <  0)  return sz+i;
-    else if (i >= sz) return i-sz;
-    else              return i;
-}
-
-#endif  // __cplusplus
-#endif  // __ACCELERATE_CUDA_STENCIL_H__
-
diff --git a/cubits/accelerate_cuda_util.h b/cubits/accelerate_cuda_util.h
deleted file mode 100644
--- a/cubits/accelerate_cuda_util.h
+++ /dev/null
@@ -1,70 +0,0 @@
-/* -----------------------------------------------------------------------------
- *
- * Module      : Util
- * Copyright   : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee
- *               [2009..2012] Manuel M T Chakravarty, Gabriele Keller, Trevor L. McDonell
- * License     : BSD3
- *
- * Maintainer  : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>
- * Stability   : experimental
- *
- * ---------------------------------------------------------------------------*/
-
-#ifndef __ACCELERATE_CUDA_UTIL_H__
-#define __ACCELERATE_CUDA_UTIL_H__
-
-#include <math.h>
-#include <cuda_runtime.h>
-
-/*
- * Core assert function. Don't let this escape...
- */
-#if defined(__CUDACC__) || !defined(__DEVICE_EMULATION__)
-#define __assert(e, file, line) ((void)0)
-#else
-#define __assert(e, file, line) \
-    ((void) fprintf (stderr, "%s:%u: failed assertion `%s'\n", file, line, e), abort())
-#endif
-
-/*
- * Test the given expression, and abort the program if it evaluates to false.
- * Only available in debug mode.
- */
-#ifndef _DEBUG
-#define assert(e)               ((void)0)
-#else
-#define assert(e)  \
-    ((void) ((e) ? (void(0)) : __assert (#e, __FILE__, __LINE__)))
-#endif
-
-/*
- * Macro to insert __syncthreads() in device emulation mode
- */
-#ifdef __DEVICE_EMULATION__
-#define __EMUSYNC               __syncthreads()
-#else
-#define __EMUSYNC
-#endif
-
-/*
- * Check the return status of CUDA API calls, and abort with an appropriate
- * error string on failure.
- */
-#define CUDA_SAFE_CALL_NO_SYNC(call)                                           \
-    do {                                                                       \
-        cudaError err = call;                                                  \
-        if(cudaSuccess != err) {                                               \
-            const char *str = cudaGetErrorString(err);                         \
-            __assert(str, __FILE__, __LINE__);                                 \
-        }                                                                      \
-    } while (0)
-
-#define CUDA_SAFE_CALL(call)                                                   \
-    do {                                                                       \
-        CUDA_SAFE_CALL_NO_SYNC(call);                                          \
-        CUDA_SAFE_CALL_NO_SYNC(cudaThreadSynchronize());                       \
-    } while (0)
-
-#undef __assert
-#endif  // __ACCELERATE_CUDA_UTIL_H__
-
diff --git a/include/AccFFI.h b/include/AccFFI.h
new file mode 100644
--- /dev/null
+++ b/include/AccFFI.h
@@ -0,0 +1,34 @@
+#include "HsFFI.h"
+#include <cuda.h>
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef HsStablePtr AccHandle;
+typedef HsStablePtr AccProgram;
+
+typedef struct {
+  int*    shape;
+  void**  adata;
+} InputArray;
+
+typedef struct {
+  int*    shape;
+  void**  adata;
+  HsStablePtr stable_ptr;
+} OutputArray;
+
+extern HsStablePtr accelerateCreate(HsInt32 device, HsPtr ctx);
+extern void accelerateDestroy(AccHandle hndl);
+extern void runProgram(AccHandle hndl, AccProgram p, InputArray* in, OutputArray* out);
+extern void freeOutput(OutputArray* out);
+extern void freeProgram(AccProgram a1);
+
+AccHandle accelerateInit(int argc, char** argv, int device, CUcontext ctx) {
+  hs_init(&argc, &argv);
+  return accelerateCreate(device, ctx);
+}
+
+#ifdef __cplusplus
+}
+#endif