diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,30 @@
+Copyright (c)2011, Trevor L. McDonell
+
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+      copyright notice, this list of conditions and the following
+      disclaimer in the documentation and/or other materials provided
+      with the distribution.
+
+    * Neither the name of Trevor L. McDonell nor the names of other
+      contributors may be used to endorse or promote products derived
+      from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/accelerate-examples.cabal b/accelerate-examples.cabal
new file mode 100644
--- /dev/null
+++ b/accelerate-examples.cabal
@@ -0,0 +1,92 @@
+Name:                accelerate-examples
+Version:             0.2.0.0
+Synopsis:            Examples using the Accelerate library
+Description:         Examples using the Accelerate library
+License:             BSD3
+License-file:        LICENSE
+Author:              The Accelerate Team
+Maintainer:          Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>
+Category:            Compilers/Interpreters,
+Build-type:          Simple
+Cabal-version:       >=1.6
+
+
+Flag cuda
+  Description:          Enable the CUDA parallel backend for NVIDIA GPUs
+  Default:              True
+
+Flag io
+  Description:          Provide access to the block copy I/O functionality
+  Default:              False
+
+
+Executable accelerate-examples
+  Main-is:              Main.hs
+  Other-modules:        Benchmark
+                        Test
+                        Validate
+                        Config
+                        Random
+                        Util
+                        Canny
+                        IntegralImage
+                        PGM
+                        BlockCopy
+                        Vector
+                        Backpermute
+                        Map
+                        ScanSeg
+                        Stencil2
+                        ZipWith
+                        Fold
+                        Permute
+                        Stencil
+                        Zip
+                        BlackScholes
+                        SASUM
+                        SharingRecovery
+                        DotP
+                        SAXPY
+                        SliceExamples
+                        Filter
+                        Radix
+                        SMVM
+                        SMVM.Matrix
+                        SMVM.MatrixMarket
+  hs-source-dirs:       src
+                        tests/primitives
+                        tests/simple
+                        tests/image-processing
+                        tests/io
+
+  c-sources:            tests/io/fill_with_values.cpp
+  extra-libraries:      stdc++
+
+  ghc-options:          -Wall -O2
+  if impl(ghc >= 7.0)
+    ghc-options:        -rtsopts
+
+  if flag(cuda)
+    CPP-options:        -DACCELERATE_CUDA_BACKEND
+
+  if flag(io)
+    CPP-options:        -DACCELERATE_IO
+
+  build-depends:        accelerate              == 0.10.*,
+                        array                   >= 0.3 && < 0.5,
+                        attoparsec              == 0.8.*,
+                        base                    == 4.*,
+                        bytestring              == 0.9.*,
+                        bytestring-lexing       == 0.2.*,
+                        cmdargs                 == 0.6.*,
+                        criterion               == 0.5.*,
+                        deepseq                 >= 1.1 && < 1.4,
+                        directory               >= 1.0 && < 1.2,
+                        filepath                >= 1.0 && < 1.4,
+                        mtl                     >= 1.1 && < 3.0,
+                        mwc-random              == 0.8.*,
+                        pgm                     == 0.1.*,
+                        pretty                  >= 1.0 && < 1.2,
+                        vector                  == 0.7.*,
+                        vector-algorithms       == 0.4.*
+
diff --git a/src/Benchmark.hs b/src/Benchmark.hs
new file mode 100644
--- /dev/null
+++ b/src/Benchmark.hs
@@ -0,0 +1,58 @@
+{-# LANGUAGE FlexibleContexts #-}
+{-# OPTIONS_GHC -fno-warn-orphans #-}
+
+module Benchmark where
+
+import Data.Array.IArray
+import Data.Array.Unboxed                       (UArray)
+
+import Data.List
+import Data.Monoid
+import Criterion
+import Criterion.IO
+import Criterion.Config
+import Criterion.Main
+import Criterion.Monad
+import Criterion.Types
+import Criterion.Environment
+import Control.Monad
+import Control.Monad.Trans                      (liftIO)
+import Control.DeepSeq
+import System.IO
+import System.Directory
+import System.Environment
+
+
+instance (Ix dim, IArray UArray e) => NFData (UArray dim e) where
+  rnf a = a ! head (indices a) `seq` ()
+
+
+-- Much like defaultMain, but we ignore any non-flag command line arguments,
+-- which we take as the inputs to the program itself (returned via getArg')
+--
+runBenchmark :: [Benchmark] -> IO ()
+runBenchmark = runBenchmarkWith defaultConfig (return ())
+
+runBenchmarkWith :: Config -> Criterion () -> [Benchmark] -> IO ()
+runBenchmarkWith defCfg prep bs = do
+  (cfg, _) <- parseArgs defCfg defaultOptions =<< getArgs
+  withConfig cfg $
+    if cfgPrintExit cfg == List
+    then do
+      _ <- note "Benchmarks:\n"
+      mapM_ (note "  %s\n") (sort $ concatMap benchNames bs)
+    else do
+      case getLast $ cfgSummaryFile cfg of
+        Just fn -> liftIO $ writeFileOnce fn "Name,Mean,MeanLB,MeanUB,Stddev,StddevLB,StddevUB\n"
+        Nothing -> return ()
+      env <- measureEnvironment
+      let shouldRun = const True
+      prep
+      runAndAnalyse shouldRun env $ BenchGroup "" bs
+
+writeFileOnce :: FilePath -> String -> IO ()
+writeFileOnce fn line = do
+  exists <- doesFileExist fn
+  size   <- withFile fn ReadWriteMode hFileSize
+  unless (exists && size > 0) $ writeFile fn line
+
diff --git a/src/Config.hs b/src/Config.hs
new file mode 100644
--- /dev/null
+++ b/src/Config.hs
@@ -0,0 +1,135 @@
+{-# LANGUAGE CPP, DeriveDataTypeable #-}
+
+module Config where
+
+import Data.Version
+import Text.PrettyPrint
+import System.Console.CmdArgs
+import Paths_accelerate_examples
+
+-- The Accelerate backends available to test, which should be no larger than the
+-- build configuration for the Accelerate library itself.
+--
+data Backend
+  = Interpreter
+#ifdef ACCELERATE_CUDA_BACKEND
+  | CUDA
+#endif
+  deriving (Show, Data, Typeable)
+
+
+-- Program configuration options
+--
+data Config = Config
+  {
+    -- common options
+    cfgBackend     :: Backend
+  , cfgVerify      :: Bool
+  , cfgElements    :: Int
+  , cfgImage       :: Maybe FilePath
+  , cfgMatrix      :: Maybe FilePath
+
+    -- criterion hooks
+  , cfgPerformGC   :: Bool
+  , cfgConfidence  :: Maybe Double
+  , cfgResamples   :: Maybe Int
+  , cfgSummaryFile :: Maybe FilePath
+
+    -- names of tests to run (all non-option arguments)
+  , cfgArgs        :: [String]
+  }
+  deriving (Show, Data, Typeable)
+
+
+-- With list of (name,description) pairs for the available tests
+--
+defaultConfig :: [(String,String)] -> Config
+defaultConfig testPrograms = Config
+  {
+    cfgBackend = enum
+    [ Interpreter
+        &= help "Reference implementation (sequential)"
+#ifdef ACCELERATE_CUDA_BACKEND
+    , CUDA
+        &= explicit
+        &= name "cuda"
+        &= help "Implementation for NVIDIA GPUs (parallel)"
+#endif
+    ]
+
+  , cfgVerify = def
+      &= explicit
+      &= name "k"
+      &= name "verify"
+      &= help "Only verify examples, do not run timing tests"
+
+  , cfgElements = 1000000
+      &= explicit
+      &= name "n"
+      &= name "size"
+      &= help "Canonical test data size (1000000)"
+
+  , cfgImage = def
+      &= explicit
+      &= name "i"
+      &= name "image"
+      &= help "PGM image file to use for image-processing tests"
+      &= typFile
+
+  , cfgMatrix = def
+      &= name "m"
+      &= name "matrix"
+      &= explicit
+      &= help "MatrixMarket file to use for SMVM test"
+      &= typFile
+
+  , cfgPerformGC = enum
+    [ False
+        &= name "G"
+        &= name "no-gc"
+        &= explicit
+        &= help "Do not collect garbage between iterations"
+    , True
+        &= name "g"
+        &= name "gc"
+        &= explicit
+        &= help "Collect garbage between iterations"
+    ]
+
+  , cfgConfidence = def
+      &= explicit
+      &= name "I"
+      &= name "ci"
+      &= help "Bootstrap confidence interval"
+      &= typ  "CI"
+
+  , cfgResamples = def
+      &= explicit
+      &= name "s"
+      &= name "resamples"
+      &= help "Number of bootstrap resamples to perform"
+
+  , cfgSummaryFile = def
+      &= name "u"
+      &= name "summary"
+      &= explicit
+      &= help "Produce a summary CSV file of all results"
+      &= typFile
+
+  , cfgArgs = def
+      &= args
+      &= typ  "TESTS"
+  }
+  &= program "accelerate-examples"
+  &= summary "accelerate-examples (c) 2011 The Accelerate Team"
+  &= versionArg [summary $ "accelerate-examples-" ++ showVersion version]
+  &= verbosityArgs [help "Print more output"] [help "Print less output"]
+  &= details (
+      [ "Available tests, by prefix match:"
+      , "  <default>             run all tests"
+      ]
+      ++
+      map (\(n,d) -> render . nest 2 $ text n $$ nest 22 (text d)) testPrograms)
+      --
+      -- magic number to make the second columns of the help text align
+
diff --git a/src/Main.hs b/src/Main.hs
new file mode 100644
--- /dev/null
+++ b/src/Main.hs
@@ -0,0 +1,74 @@
+{-# LANGUAGE TupleSections #-}
+
+module Main where
+
+import Test
+import Config
+import Benchmark
+
+import Prelude                  hiding (catch)
+import Data.List
+import Data.Maybe
+import Control.Arrow
+import Control.Monad
+import System.Environment
+import System.Console.CmdArgs   (cmdArgs, getVerbosity, Verbosity(..))
+
+
+-- Process command line options, prepare selected test programs for benchmarking
+-- or verification
+--
+processArgs :: IO (Config, [Test])
+processArgs = do
+  testInfo <- map (title &&& description) `fmap` allTests undefined
+  config   <- cmdArgs $ defaultConfig testInfo
+  tests    <- filter (selected config) `fmap` allTests config
+  --
+  return (config, tests)
+  where
+    selected a = case cfgArgs a of
+                   [] -> const True
+                   ps -> \x -> any (\p -> p `isPrefixOf` title x) ps
+
+
+-- Verify results with the chosen backend, turning exceptions into failures.
+-- Pass back the tests which succeeded.
+--
+runVerify :: Config -> [Test] -> IO [Test]
+runVerify cfg tests = do
+  results <- forM tests $ \t -> (t,) `fmap` verifyTest cfg t
+  return . map fst
+         $ filter (\(_,r) -> r `elem` [Ok, Skipped]) results
+
+
+-- Run criterion timing tests in the chosen backend
+--
+runTiming :: Config -> [Test] -> IO ()
+runTiming cfg tests = do
+  verbose <- getVerbosity
+  unless (verbose == Quiet) $ putStrLn ""
+  let args = [ maybe "" (\ci -> "--ci=" ++ show ci)       (cfgConfidence cfg)
+             , maybe "" (\r  -> "--resamples=" ++ show r) (cfgResamples cfg)
+             , maybe "" (\f  -> "--summary=" ++ f)        (cfgSummaryFile cfg)
+             , if cfgPerformGC cfg then "-g" else "-G"
+             , case verbose of
+                 Loud   -> "--verbose"
+                 Quiet  -> "--quiet"
+                 Normal -> ""
+             ]
+  --
+  withArgs args
+    . runBenchmark
+    . catMaybes
+    $ map (benchmarkTest cfg) tests
+
+
+-- Main
+-- ====
+
+main :: IO ()
+main = do
+  (config, tests) <- processArgs
+  valid           <- runVerify config tests
+  --
+  unless (null valid || cfgVerify config) $ runTiming config valid
diff --git a/src/Random.hs b/src/Random.hs
new file mode 100644
--- /dev/null
+++ b/src/Random.hs
@@ -0,0 +1,72 @@
+{-# LANGUAGE BangPatterns, FlexibleContexts #-}
+
+module Random where
+
+import System.Random.MWC
+import Data.Array.IArray
+import Data.Array.Unboxed			(UArray)
+import Data.Array.IO                    	(MArray, IOUArray)
+import Control.Exception                	(evaluate)
+import Data.Array.Accelerate            	(Z(..),(:.)(..))
+import qualified Data.Vector.Generic		as G
+import qualified Data.Vector.Generic.Mutable	as GM
+import qualified Data.Array.MArray		as M
+import qualified Data.Array.Accelerate  	as Acc
+
+
+-- Convert an Unboxed Data.Array to an Accelerate Array
+--
+convertUArray :: (IArray UArray e, Acc.Elt e) => UArray Int e -> IO (Acc.Vector e)
+convertUArray v =
+  let arr = Acc.fromIArray v
+  in  evaluate (arr `Acc.indexArray` (Z:.0)) >> return arr
+
+
+-- Convert a Data.Vector to an Accelerate Array
+--
+convertVector
+  :: (IArray UArray a, MArray IOUArray a IO, G.Vector v a, Acc.Elt a)
+  => v a
+  -> IO (Acc.Vector a)
+
+convertVector vec = do
+  arr <- Acc.fromIArray `fmap` toIArray vec
+  evaluate (arr `Acc.indexArray` (Z:.0)) >> return arr
+  where
+    toIArray :: (MArray IOUArray a IO, IArray UArray a, G.Vector v a) => v a -> IO (UArray Int a)
+    toIArray v = do
+      let n = G.length v
+      mu <- M.newArray_ (0,n-1) :: MArray IOUArray a IO => IO (IOUArray Int a)
+      let go !i | i < n     = M.writeArray mu i (G.unsafeIndex v i) >> go (i+1)
+                | otherwise = M.unsafeFreeze mu
+      go 0
+
+
+-- Generate a random, uniformly distributed vector of specified size over the
+-- range. For integral types the range is inclusive, for floating point numbers
+-- the range (a,b] is used, if one ignores rounding errors.
+--
+randomUArrayR
+  :: (Variate a, MArray IOUArray a IO, IArray UArray a)
+  => (a,a)
+  -> GenIO
+  -> Int
+  -> IO (UArray Int a)
+
+randomUArrayR lim gen n = do
+  mu  <- M.newArray_ (0,n-1) :: MArray IOUArray e IO => IO (IOUArray Int e)
+  let go !i | i < n     = uniformR lim gen >>= M.writeArray mu i >> go (i+1)
+            | otherwise = M.unsafeFreeze mu
+  go 0
+
+
+-- Generate a uniformly distributed Data.Vector of specified range and size
+--
+randomVectorR :: (G.Vector v a, Variate a) => (a,a) -> GenIO -> Int -> IO (v a)
+randomVectorR lim gen n = do
+  mu <- GM.unsafeNew n
+  let go !i | i < n     = uniformR lim gen >>= GM.unsafeWrite mu i >> go (i+1)
+            | otherwise = G.unsafeFreeze mu
+  go 0
+
+
diff --git a/src/Test.hs b/src/Test.hs
new file mode 100644
--- /dev/null
+++ b/src/Test.hs
@@ -0,0 +1,234 @@
+{-# LANGUAGE CPP, ExistentialQuantification #-}
+
+module Test (
+
+  Title, Description, Test(..), Status(..),
+  allTests, verifyTest, benchmarkTest
+
+) where
+
+-- individual test implementations
+import qualified Map
+import qualified Zip
+import qualified ZipWith
+import qualified Fold
+import qualified ScanSeg
+import qualified Stencil
+import qualified Stencil2
+import qualified Permute
+import qualified Backpermute
+
+import qualified SASUM
+import qualified SAXPY
+import qualified DotP
+import qualified Filter
+import qualified SMVM
+import qualified BlackScholes
+import qualified Radix
+import qualified SliceExamples
+
+#ifdef ACCELERATE_IO
+import qualified BlockCopy
+#endif
+
+import qualified Canny
+import qualified IntegralImage
+import qualified SharingRecovery
+
+-- friends
+import Util
+import Config
+import Validate
+
+-- libraries
+import Prelude                                          hiding (catch)
+import Criterion                                        (Benchmark, bench, whnf)
+import Data.Maybe
+import Data.Array.IArray
+import Control.Monad
+import Control.Exception
+import System.IO
+import System.IO.Unsafe
+import System.Console.CmdArgs                           (getVerbosity, Verbosity(..))
+
+import Data.Array.Accelerate                            (Acc)
+import qualified Data.Array.Accelerate                  as Acc
+import qualified Data.Array.Accelerate.Interpreter      as Interpreter
+
+#ifdef ACCELERATE_CUDA_BACKEND
+import qualified Data.Array.Accelerate.CUDA             as CUDA
+#endif
+
+
+data Status
+  = Ok
+  | Skipped
+  | Failed String
+
+instance Eq Status where
+  Ok      == Ok      = True
+  Skipped == Skipped = True
+  _       == _       = False
+
+instance Show Status where
+  show Ok         = "Ok"
+  show Skipped    = "Skipped"
+  show (Failed s) = "Failed: " ++ s
+
+
+type Title       = String
+type Description = String
+
+data Test
+  -- A cannonical test program, where we have a reference implementation that
+  -- the Accelerate program must match. The 'convert' field is slightly magic:
+  -- we need to carry it around as a proof that ELtRepr sh ~ EltRepr ix.
+  --
+  = forall array ix sh e. (Similar e, Acc.Elt e, Acc.Shape sh, Show ix, Show e, IArray array e, Ix ix) => Test
+  { title       :: Title
+  , description :: Description
+  , reference   :: () -> array ix e
+  , accelerate  :: () -> Acc (Acc.Array sh e)
+  , convert     :: Acc.Array sh e -> array ix e
+  }
+
+  -- No reference implementation, so the result can not be validated, but we can
+  -- check that no exceptions are thrown, and benchmark the operation.
+  --
+  | forall sh e. (Acc.Elt e, Acc.Shape sh) => TestNoRef
+  { title       :: Title
+  , description :: Description
+  , accelerate  :: () -> Acc (Acc.Array sh e)
+  }
+
+  -- An IO action. Run once to verify that no exceptions are thrown, do not
+  -- benchmark.
+  --
+  | forall a. TestIO
+  { title       :: Title
+  , description :: Description
+  , action      :: IO a
+  }
+
+
+allTests :: Config -> IO [Test]
+allTests cfg = sequence'
+  [
+
+    -- primitive functions
+    mkTest "map-abs"               "absolute value of each element"             $ Map.run "abs" n
+  , mkTest "map-plus"              "add a constant to each element"             $ Map.run "plus" n
+  , mkTest "map-square"            "square of each element"                     $ Map.run "square" n
+  , mkTest "zip"                   "vector zip"                                 $ Zip.run n
+  , mkTest "zipWith-plus"          "element-wise addition"                      $ ZipWith.run "plus" n
+  , mkTest "fold-sum"              "vector reduction: fold (+) 0"               $ Fold.run "sum" n
+  , mkTest "fold-product"          "vector product: fold (*) 1"                 $ Fold.run "product" n
+  , mkTest "fold-maximum"          "maximum of a vector: fold1 max"             $ Fold.run "maximum" n
+  , mkTest "fold-minimum"          "minimum of a vector: fold1 min"             $ Fold.run "minimum" n
+  , mkTest "fold-2d-sum"           "reduction along innermost matrix dimension" $ Fold.run2d "sum-2d" n
+  , mkTest "fold-2d-product"       "product along innermost matrix dimension"   $ Fold.run2d "product-2d" n
+  , mkTest "scanseg-sum"           "segmented reduction"                        $ ScanSeg.run "sum" n
+  , mkTest "stencil-1D"            "3-element vector"                           $ Stencil.run "1D" n
+  , mkTest "stencil-2D"            "3x3 pattern"                                $ Stencil.run2D "2D" n
+  , mkTest "stencil-3D"            "3x3x3 pattern"                              $ Stencil.run3D "3D" n
+  , mkTest "stencil-3x3-cross"     "3x3 cross pattern"                          $ Stencil.run2D "3x3-cross" n
+  , mkTest "stencil-3x3-pair"      "3x3 non-symmetric pattern with pairs"       $ Stencil.run2D "3x3-pair" n
+  , mkTest "stencil2-2D"           "3x3 pattern"                                $ Stencil2.run2D "2D" n
+  , mkTest "permute-hist"          "histogram"                                  $ Permute.run "histogram" n
+  , mkTest "backpermute-reverse"   "reverse a vector"                           $ Backpermute.run "reverse" n
+  , mkTest "backpermute-transpose" "transpose a matrix"                         $ Backpermute.run2d "transpose" n
+
+    -- simple examples
+  , mkTest "sasum"                 "sum of absolute values"                     $ SASUM.run n
+  , mkTest "saxpy"                 "scalar alpha*x + y"                         $ SAXPY.run n
+  , mkTest "dotp"                  "vector dot-product"                         $ DotP.run n
+  , mkTest "filter"                "return elements that satisfy a predicate"   $ Filter.run n
+  , mkTest "smvm"                  "sparse-matrix vector multiplication"        $ SMVM.run (cfgMatrix cfg)
+  , mkTest "black-scholes"         "Black-Scholes option pricing"               $ BlackScholes.run n
+  , mkTest "radixsort"             "radix sort"                                 $ Radix.run n
+
+#ifdef ACCELERATE_IO
+    -- Array IO
+  , mkIO   "io"                    "array IO test"                              $ BlockCopy.run
+#endif
+
+  --  image processing
+  , mkNoRef "canny"          "canny edge detection"                       $ Canny.run img
+  , mkNoRef "integral-image" "image integral (2D scan)"                   $ IntegralImage.run img
+  -- slices
+  , mkTest "slices"  "replicate (Z:.2:.All:.All)" $ SliceExamples.run1
+  , mkTest "slices"  "replicate (Z:.All:.2:.All)" $ SliceExamples.run2
+  , mkTest "slices"  "replicate (Z:.All:.All:.2)" $ SliceExamples.run3
+  , mkTest "slices"  "replicate (Any:.2)"         $ SliceExamples.run4  
+  , mkTest "slices"  "replicate (Z:.2:.2:.2)"     $ SliceExamples.run5
+  --
+  , mkIO "sharing-recovery" "simple"    $ return (show SharingRecovery.simple)
+  , mkIO "sharing-recovery" "orderFail" $ return (show SharingRecovery.orderFail)
+  , mkIO "sharing-recovery" "testSort"  $ return (show SharingRecovery.testSort)
+  , mkIO "sharing-recovery" "muchSharing" $ return (show $ SharingRecovery.muchSharing 20)
+  , mkIO "sharing-recovery" "bfsFail"   $ return (show SharingRecovery.bfsFail)
+  , mkIO "sharing-recovery" "twoLetsSameLevel"  $ return (show SharingRecovery.twoLetsSameLevel)
+  , mkIO "sharing-recovery" "twoLetsSameLevel2" $ return (show SharingRecovery.twoLetsSameLevel2)
+  , mkIO "sharing-recovery" "noLetAtTop"   $ return (show SharingRecovery.noLetAtTop)
+  , mkIO "sharing-recovery" "noLetAtTop2"   $ return (show SharingRecovery.noLetAtTop2)
+  , mkIO "sharing-recovery" "pipe" $ return (show SharingRecovery.pipe)
+  ]
+  where
+    n   = cfgElements cfg
+    img = fromMaybe (error "no image file specified") (cfgImage cfg)
+    --
+    mkTest name desc builder = do
+      ~(ref,acc) <- unsafeInterleaveIO builder  -- must be super lazy
+      return $ Test name desc ref acc Acc.toIArray
+
+    mkNoRef name desc builder = do
+      acc <- unsafeInterleaveIO builder
+      return $ TestNoRef name desc acc
+
+    mkIO name desc act = return $ TestIO name desc act
+
+
+-- How to evaluate Accelerate programs with the chosen backend?
+--
+backend :: Acc.Arrays a => Config -> Acc a -> a
+backend cfg =
+  case cfgBackend cfg of
+    Interpreter -> Interpreter.run
+#ifdef ACCELERATE_CUDA_BACKEND
+    CUDA        -> CUDA.run
+#endif
+
+
+-- Verify that the Accelerate and reference implementations yield the same
+-- result in the chosen backend
+--
+verifyTest :: Config -> Test -> IO Status
+verifyTest cfg test = do
+  quiet <- (==Quiet) `fmap` getVerbosity
+  verify quiet `catch` \e -> let r = Failed (show (e :: SomeException))
+                             in  putStrLn (show r) >> return r
+  where
+    run acc      = backend cfg $ acc ()
+    verify quiet = do
+      unless quiet $ putStr (title test ++ ": ") >> hFlush stdout
+      result <- case test of
+        Test _ _ ref acc cvt ->
+          return $ case validate (ref ()) (cvt $ run acc) of
+                     []   -> Ok
+                     errs -> Failed . unlines . ("":)
+                                    $ map (\(i,v) -> ">>> " ++ shows i " : " ++ show v) errs
+
+        TestNoRef _ _ acc -> return $ run acc `seq` Ok
+        TestIO _ _ act    -> act >>= \v -> v `seq` return Ok
+      --
+      unless quiet $ putStrLn (show result)
+      return result
+
+
+-- Benchmark a test with Criterion
+--
+benchmarkTest :: Config -> Test -> Maybe Benchmark
+benchmarkTest cfg (Test name _ _ acc _)  = Just . bench name $ whnf (backend cfg . acc) ()
+benchmarkTest cfg (TestNoRef name _ acc) = Just . bench name $ whnf (backend cfg . acc) ()
+benchmarkTest _   (TestIO _ _ _)         = Nothing
+
diff --git a/src/Util.hs b/src/Util.hs
new file mode 100644
--- /dev/null
+++ b/src/Util.hs
@@ -0,0 +1,17 @@
+
+module Util where
+
+import System.IO.Unsafe
+
+-- Lazier version of 'Control.Monad.sequence'
+--
+sequence' :: [IO a] -> IO [a]
+sequence' = foldr k (return [])
+  where k m ms = do { x <- m; xs <- unsafeInterleaveIO ms; return (x:xs) }
+
+mapM' :: (a -> IO b) -> [a] -> IO [b]
+mapM' f xs = sequence' $ map f xs
+
+forM' :: [a] -> (a -> IO b) -> IO [b]
+forM' = flip mapM'
+
diff --git a/src/Validate.hs b/src/Validate.hs
new file mode 100644
--- /dev/null
+++ b/src/Validate.hs
@@ -0,0 +1,102 @@
+{-# LANGUAGE FlexibleContexts, ParallelListComp #-}
+{-# OPTIONS_GHC -fno-warn-unused-binds #-}
+
+module Validate (Similar(..), validate, validate') where
+
+import Data.Int
+import Data.Word
+import Data.Array.IArray
+import Foreign.C.Types
+import Foreign.Storable
+import Control.Exception                (assert)
+import Unsafe.Coerce
+
+class Similar a where
+  sim :: a -> a -> Bool
+
+instance Similar Int     where sim = (==)
+instance Similar Int8    where sim = (==)
+instance Similar Int16   where sim = (==)
+instance Similar Int32   where sim = (==)
+instance Similar Int64   where sim = (==)
+instance Similar Word    where sim = (==)
+instance Similar Word8   where sim = (==)
+instance Similar Word16  where sim = (==)
+instance Similar Word32  where sim = (==)
+instance Similar Word64  where sim = (==)
+instance Similar CShort  where sim = (==)
+instance Similar CUShort where sim = (==)
+instance Similar CInt    where sim = (==)
+instance Similar CUInt   where sim = (==)
+instance Similar CLong   where sim = (==)
+instance Similar CULong  where sim = (==)
+instance Similar CLLong  where sim = (==)
+instance Similar CULLong where sim = (==)
+
+instance Similar Bool    where sim = (==)
+instance Similar Char    where sim = (==)
+instance Similar CChar   where sim = (==)
+instance Similar CSChar  where sim = (==)
+instance Similar CUChar  where sim = (==)
+
+instance Similar Float   where sim = absoluteOrRelative
+instance Similar CFloat  where sim = absoluteOrRelative
+instance Similar Double  where sim = absoluteOrRelative
+instance Similar CDouble where sim = absoluteOrRelative
+
+instance (Similar a, Similar b) => Similar (a,b) where
+  (x,y) `sim` (u,v) = x `sim` u && y `sim` v
+
+--
+-- http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
+--
+
+absoluteOrRelative :: (Fractional a, Ord a) => a -> a -> Bool
+absoluteOrRelative u v
+  | abs (u-v) < epsilonAbs = True
+  | abs u > abs v          = abs ((u-v) / u) < epsilonRel
+  | otherwise              = abs ((v-u) / v) < epsilonRel
+  where
+    epsilonRel = 0.001
+    epsilonAbs = 0.00001
+
+
+-- Comparisons using lexicographically ordered floating-point numbers
+-- reinterpreted as twos-complement integers.
+--
+lexicographic32 :: (Num a, Storable a) => Int -> a -> a -> Bool
+lexicographic32 maxUlps a b
+  = assert (sizeOf a == 4 && maxUlps > 0 && maxUlps < 4 * 1024 * 1024)
+  $ intDiff < fromIntegral maxUlps
+  where
+    intDiff = abs (toInt a - toInt b)
+    toInt x | x' < 0    = 0x80000000 - x'
+            | otherwise = x'
+            where x'    = unsafeCoerce x :: Int32
+
+
+lexicographic64 :: (Num a, Storable a) => Int -> a -> a -> Bool
+lexicographic64 maxUlps a b
+  = assert (sizeOf a == 8 && maxUlps > 0 && maxUlps < 8 * 1024 * 1024)
+  $ intDiff < fromIntegral maxUlps
+  where
+    intDiff = abs (toInt a - toInt b)
+    toInt x | x' < 0    = 0x8000000000000000 - x'
+            | otherwise = x'
+            where x'    = unsafeCoerce x :: Int64
+
+
+-- Compare two vectors element-wise for equality, for a given measure of
+-- similarity. The index and values are returned for pairs that fail.
+--
+validate
+  :: (IArray array e, Ix ix, Similar e)
+  => array ix e
+  -> array ix e
+  -> [(ix,(e,e))]
+validate ref arr = validate' (assocs ref) (elems arr)
+
+validate' :: (Ix ix, Similar e) => [(ix,e)] -> [e] -> [(ix,(e,e))]
+validate' ref arr =
+  filter (not . uncurry sim . snd) [ (i,(x,y)) | (i,x) <- ref | y <- arr ]
+
diff --git a/tests/image-processing/Canny.hs b/tests/image-processing/Canny.hs
new file mode 100644
--- /dev/null
+++ b/tests/image-processing/Canny.hs
@@ -0,0 +1,117 @@
+{-# LANGUAGE TypeOperators #-}
+
+module Canny where
+
+import PGM
+
+import Data.Array.Accelerate                hiding (zipWith, unindex2)
+import qualified Data.Array.Accelerate      as Acc
+
+
+type Image a      = Array DIM2 a
+type Stencil7x1 a = (Stencil3 a, Stencil7 a, Stencil3 a)
+type Stencil1x7 a = (Stencil3 a, Stencil3 a, Stencil3 a, Stencil3 a, Stencil3 a, Stencil3 a, Stencil3 a)
+
+
+convolve7x1 :: (Elt a, IsNum a) => [Exp a] -> Stencil7x1 a -> Exp a
+convolve7x1 kernel (_, (a,b,c,d,e,f,g), _) =
+  sum $ zipWith (*) kernel [a,b,c,d,e,f,g]
+
+convolve1x7 :: (Elt a, IsNum a) => [Exp a] -> Stencil1x7 a -> Exp a
+convolve1x7 kernel ((_,a,_), (_,b,_), (_,c,_), (_,d,_), (_,e,_), (_,f,_), (_,g,_)) =
+  sum $ zipWith (*) kernel [a,b,c,d,e,f,g]
+
+
+-- Gaussian smoothing
+--
+gaussian :: (Elt a, IsFloating a) => [Exp a]
+gaussian = [ 0.00442012927963
+           , 0.05384819825462
+           , 0.24133088157513
+           , 0.39788735772974
+           , 0.24133088157513
+           , 0.05384819825462
+           , 0.00442012927963 ]
+
+gaussianX :: (Elt a, IsFloating a) => Acc (Image a) -> Acc (Image a)
+gaussianX = stencil (convolve7x1 gaussian) (Constant 0)
+
+gaussianY :: (Elt a, IsFloating a) => Acc (Image a) -> Acc (Image a)
+gaussianY = stencil (convolve1x7 gaussian) (Constant 0)
+
+
+-- Gaussian derivative and gradient quantisation
+--
+gaussian' :: (Elt a, IsFloating a) => [Exp a]
+gaussian' = [ 0.02121662054222
+            , 0.17231423441479
+            , 0.38612941052022
+            , 0.0
+            ,-0.38612941052022
+            ,-0.17231423441479
+            ,-0.02121662054222 ]
+
+gradientX :: (Elt a, IsFloating a) => Acc (Image a) -> Acc (Image a)
+gradientX = stencil (convolve7x1 gaussian') (Constant 0)
+
+gradientY :: (Elt a, IsFloating a) => Acc (Image a) -> Acc (Image a)
+gradientY = stencil (convolve1x7 gaussian') (Constant 0)
+
+gradientMagnitude :: (Elt a, IsFloating a) => Acc (Image a) -> Acc (Image a) -> Acc (Image a)
+gradientMagnitude = Acc.zipWith magdir
+  where
+    magdir dx dy = let mag = sqrt (dx*dx + dy*dy)
+                    -- dir = atan2 dy dx
+                   in  mag -- lift (mag, dir)
+
+
+-- Non-maximum suppression
+--
+nonMaximumSuppression
+  :: (Elt a, IsFloating a)
+  => Exp a
+  -> Acc (Image a)
+  -> Acc (Image a)
+  -> Acc (Image a)
+  -> Acc (Image a)
+nonMaximumSuppression threshold gradX gradY gradM =
+  generate (shape gradX) $ \ix ->
+    let dx          = gradX ! ix
+        dy          = gradY ! ix
+        mag         = gradM ! ix
+        alpha       = 1.3065629648763766  -- 0.5 / sin (pi / 8.0)
+        offsetx     = Acc.round (alpha * dx / mag)
+        offsety     = Acc.round (alpha * dy / mag)
+        --
+        (m,n)       = unindex2 (shape gradX)
+        (x,y)       = unindex2 ix
+        fwd         = gradM ! lift (clamp (x+offsetx, y+offsety))
+        rev         = gradM ! lift (clamp (x-offsetx, y-offsety))
+        --
+        unindex2 uv = let Z:.u:.v = unlift uv in (u,v)
+        clamp (u,v) = lift (Z:. 0 `Acc.max` u `Acc.min` (m-1) :. 0 `Acc.max` v `Acc.min` (n-1))
+    in
+    (mag <* threshold ||* fwd >* mag ||* rev >* mag) ? (0, 1)
+
+
+-- Canny edge detection
+--
+canny :: (Elt a, IsFloating a) => Image a -> Acc (Image a)
+canny img =
+  let smooth  = gaussianX . gaussianY $ use img
+      gradX   = gradientX smooth
+      gradY   = gradientY smooth
+      gradMag = gradientMagnitude gradX gradY
+  in
+  nonMaximumSuppression 0.1 gradX gradY gradMag
+
+
+-- Main
+-- ----
+
+-- TLM: should compare to a pre-saved reference image
+run :: FilePath -> IO (() -> Acc (Array DIM2 Float))
+run file = do
+  pgm <- readPGM file
+  return (\() -> canny pgm)
+
diff --git a/tests/image-processing/IntegralImage.hs b/tests/image-processing/IntegralImage.hs
new file mode 100644
--- /dev/null
+++ b/tests/image-processing/IntegralImage.hs
@@ -0,0 +1,49 @@
+{-# LANGUAGE ScopedTypeVariables, TypeOperators #-}
+
+module IntegralImage where
+
+import PGM
+
+import Data.Array.Accelerate as Acc
+
+
+-- |The value of each element in an integral image is the sum of all input elements
+-- above and to the left, inclusive. It is calculated by performing an inclusive/post
+-- scan from left-to-right then top-to-bottom.
+--
+integralImage :: (Elt a, IsNum a) => Array DIM2 a -> Acc (Array DIM2 a)
+integralImage img = sumTable
+  where
+    -- scan rows
+    rowArr  = reshape (lift $ Z:.(w * h)) arr
+    rowSegs = Acc.replicate (lift $ Z:.h) $ unit w
+    rowSum  = reshape (lift (Z:.w:.h)) $ Acc.scanl1Seg (+) rowArr rowSegs
+
+    -- scan cols
+    colArr  = reshape (lift $ Z:.(h * w)) $ transpose2D rowSum
+    colSegs = Acc.replicate (lift $ Z:.w) $ unit h
+    colSum  = reshape (lift (Z:.h:.w)) $ Acc.scanl1Seg (+) colArr colSegs
+
+    -- transpose back
+    sumTable = transpose2D colSum
+
+    --
+    arr     = use img
+    Z:.w:.h = unlift $ shape arr
+
+
+-- |Simple 2D matrix transpose.
+--
+transpose2D :: Elt a => Acc (Array DIM2 a) -> Acc (Array DIM2 a)
+transpose2D arr = backpermute (swap $ shape arr) swap arr
+  where
+    swap = lift1 $ \(Z:.x:.y) -> Z:.y:.x :: Z :. Exp Int :. Exp Int
+
+
+-- Run integralImage over the input PGM
+--
+run :: FilePath -> IO (() -> Acc (Array DIM2 Float))
+run file = do
+  pgm <- readPGM file
+  return (\() -> integralImage pgm)
+
diff --git a/tests/image-processing/PGM.hs b/tests/image-processing/PGM.hs
new file mode 100644
--- /dev/null
+++ b/tests/image-processing/PGM.hs
@@ -0,0 +1,32 @@
+--
+-- Load a PGM file. MacOS X users might find the following quicklook plugin
+-- useful for viewing PGM files:
+--
+-- http://code.google.com/p/quicklook-pfm/
+--
+
+module PGM where
+
+import Control.Applicative
+import Graphics.Pgm
+
+import Prelude                   as P
+import Data.Array.Accelerate     as Acc
+import Data.Array.Unboxed        hiding (Array)
+import qualified Data.ByteString as B
+
+
+-- Read an 8-bit PGM file, and marshal to an Accelerate array as floating-point
+-- data in the range [0,1].
+--
+readPGM :: FilePath -> IO (Array DIM2 Float)
+readPGM fp = do
+  img <- either (error . show) head . pgmsToArrays <$> B.readFile fp :: IO (UArray (Int,Int) Word8)
+  return . fromIArray $ amap (\x -> P.fromIntegral x / 255) img
+
+
+writePGM :: FilePath -> Array DIM2 Float -> IO ()
+writePGM fp img =
+  let arr = toIArray img :: UArray (Int,Int) Float
+  in  arrayToFile fp $ amap (\x -> P.round (255 * x)) arr
+
diff --git a/tests/io/BlockCopy.hs b/tests/io/BlockCopy.hs
new file mode 100644
--- /dev/null
+++ b/tests/io/BlockCopy.hs
@@ -0,0 +1,90 @@
+{-# LANGUAGE ForeignFunctionInterface, ScopedTypeVariables, TypeOperators #-}
+
+module BlockCopy where
+
+-- standard libraries
+import Prelude as P
+import Foreign.Ptr
+import Control.Monad
+import Control.Exception
+
+-- friends
+import Data.Array.Accelerate
+import Data.Array.Accelerate.IO
+
+assertEqual :: (Eq a, Show a) => String -> a -> a -> IO ()
+assertEqual preface expected actual =
+  unless (actual == expected) (throw $ AssertionFailed msg)
+  where
+    msg = (if null preface then "" else preface ++ "\n")  ++
+          "expected: " ++ show expected ++ "\n but got: " ++ show actual
+
+run :: IO ()
+run =
+  mapM_ (\(msg,act) -> putStrLn ("test: " ++ msg) >> act)
+    [ ("fromPtr Int",          testBlockCopyPrim)
+    , ("fromPtr (Int,Double)", testBlockCopyTuples)
+    , ("toPtr Int16",          testBlockCopyFromArrayInt16)
+    , ("toPtr Int32",          testBlockCopyFromArrayInt32)
+    , ("toPtr Int64",          testBlockCopyFromArrayInt64)
+    , ("fromArray Int",        testBlockCopyFromArrayWithFunctions) ]
+
+
+testBlockCopyPrim :: IO ()
+testBlockCopyPrim = do
+  ptr <- oneToTen
+  (arr :: Vector Int32) <- fromPtr (Z :. 10) ((), ptr)
+  assertEqual "Not equal" [1..10] (toList arr)
+
+testBlockCopyTuples :: IO ()
+testBlockCopyTuples = do
+  intPtr    <- oneToTen
+  doublePtr <- tenToOne
+  (arr :: Vector (Int32, Double)) <- fromPtr (Z :. 10) (((), intPtr), doublePtr)
+  assertEqual "Not equal" [ (x, P.fromIntegral (11 - x)) | x <- [1..10]] (toList arr)
+
+testBlockCopyFromArrayWithFunctions :: IO ()
+testBlockCopyFromArrayWithFunctions = do
+  let n = 5^(3::Int)
+  let (arr :: Array (Z:.Int:.Int:.Int) Int32) = fromList (Z:.5:.5:.5) [2*x | x <- [0..n-1]]
+  ohi <- nInt32s (P.fromIntegral n)
+  fromArray arr ((), memcpy ohi)
+  b   <- isFilledWithEvens32 ohi (P.fromIntegral n)
+  assertEqual "Not equal" 1 b
+
+testBlockCopyFromArrayInt16 :: IO ()
+testBlockCopyFromArrayInt16 = do
+  let n = 50
+  let (arr :: Vector Int16) = fromList (Z:.n) [2 * P.fromIntegral x | x <- [0..n-1]]
+  ohi <- nInt16s (P.fromIntegral n)
+  toPtr arr ((), ohi)
+  b   <- isFilledWithEvens16 ohi (P.fromIntegral n)
+  assertEqual "Not equal" 1 b
+
+testBlockCopyFromArrayInt32 :: IO ()
+testBlockCopyFromArrayInt32 = do
+  let (arr :: Array (Z:.Int:.Int) Int32) = fromList (Z:.10:.10) [2*x | x <- [0..99]]
+  ohi <- nInt32s 100
+  toPtr arr ((), ohi)
+  b   <- isFilledWithEvens32 ohi 100
+  assertEqual "Not equal" 1 b
+
+testBlockCopyFromArrayInt64 :: IO ()
+testBlockCopyFromArrayInt64 = do
+  let n = 73
+  let (arr :: Vector Int64) = fromList (Z:.n) [2 * P.fromIntegral x | x <- [0..n-1]]
+  ohi <- nInt64s (P.fromIntegral n)
+  toPtr arr ((), ohi)
+  b   <- isFilledWithEvens64 ohi (P.fromIntegral n)
+  assertEqual "Not equal" 1 b
+
+foreign import ccall "one_to_ten" oneToTen :: IO (Ptr Int32)
+foreign import ccall "ten_to_one" tenToOne :: IO (Ptr Double)
+foreign import ccall "n_int_16s" nInt16s :: CInt -> IO (Ptr Int16)
+foreign import ccall "n_int_32s" nInt32s :: CInt -> IO (Ptr Int32)
+foreign import ccall "n_int_64s" nInt64s :: CInt -> IO (Ptr Int64)
+foreign import ccall "is_filled_with_evens_16" isFilledWithEvens16 :: Ptr Int16 -> CInt -> IO CInt
+foreign import ccall "is_filled_with_evens_32" isFilledWithEvens32 :: Ptr Int32 -> CInt -> IO CInt
+foreign import ccall "is_filled_with_evens_64" isFilledWithEvens64 :: Ptr Int64 -> CInt -> IO CInt
+foreign import ccall memcpy :: Ptr a -> Ptr b -> Int -> IO ()
+
diff --git a/tests/io/Vector.hs b/tests/io/Vector.hs
new file mode 100644
--- /dev/null
+++ b/tests/io/Vector.hs
@@ -0,0 +1,75 @@
+{-# LANGUAGE TemplateHaskell, TypeFamilies #-}
+
+module Vector where
+
+import Control.Applicative
+
+import Data.Array.Accelerate hiding (fromList)
+import Data.Array.Accelerate.Array.Sugar (EltRepr)
+import Data.Array.Accelerate.IO
+
+import Data.Vector.Storable
+
+import Foreign
+
+import Test.QuickCheck
+import Test.QuickCheck.All
+import Test.QuickCheck.Monadic
+
+roundtrip :: ( Arbitrary a
+             , Eq a
+             , Elt a
+             , Storable a
+             , BlockPtrs (EltRepr a) ~ ((), Ptr a) )
+          => [a] -> Property
+roundtrip xs = monadicIO $ do
+  let xsv = fromList xs
+  xsarr <- run $ fromVectorIO xsv
+  xsv'  <- run $ toVectorIO xsarr
+  assert (xsv == xsv')
+
+unsaferoundtrip :: ( Arbitrary a
+                   , Eq a
+                   , Elt a
+                   , Storable a
+                   , BlockPtrs (EltRepr a) ~ ((), Ptr a) )
+                => [a] -> Bool
+unsaferoundtrip xs = xsv == (toVector (fromVector xsv))
+  where xsv = fromList xs
+
+prop_Int8_roundtrip :: [Int8] -> Property
+prop_Int8_roundtrip = roundtrip
+prop_Int8_unsaferoundtrip :: [Int8] -> Bool
+prop_Int8_unsaferoundtrip = unsaferoundtrip
+
+prop_Int16_roundtrip :: [Int16] -> Property
+prop_Int16_roundtrip = roundtrip
+prop_Int16_unsaferoundtrip :: [Int16] -> Bool
+prop_Int16_unsaferoundtrip = unsaferoundtrip
+
+prop_Int32_roundtrip :: [Int32] -> Property
+prop_Int32_roundtrip = roundtrip
+prop_Int32_unsaferoundtrip :: [Int32] -> Bool
+prop_Int32_unsaferoundtrip = unsaferoundtrip
+
+prop_Int64_roundtrip :: [Int64] -> Property
+prop_Int64_roundtrip = roundtrip
+prop_Int64_unsaferoundtrip :: [Int64] -> Bool
+prop_Int64_unsaferoundtrip = unsaferoundtrip
+
+prop_Int_roundtrip :: [Int] -> Property
+prop_Int_roundtrip = roundtrip
+prop_Int_unsaferoundtrip :: [Int] -> Bool
+prop_Int_unsaferoundtrip = unsaferoundtrip
+
+prop_Float_roundtrip :: [Float] -> Property
+prop_Float_roundtrip = roundtrip
+prop_Float_unsaferoundtrip :: [Float] -> Bool
+prop_Float_unsaferoundtrip = unsaferoundtrip
+
+prop_Double_roundtrip :: [Double] -> Property
+prop_Double_roundtrip = roundtrip
+prop_Double_unsaferoundtrip :: [Double] -> Bool
+prop_Double_unsaferoundtrip = unsaferoundtrip
+
+test = $quickCheckAll
diff --git a/tests/io/fill_with_values.cpp b/tests/io/fill_with_values.cpp
new file mode 100644
--- /dev/null
+++ b/tests/io/fill_with_values.cpp
@@ -0,0 +1,81 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdint.h>
+
+/* Returns one if it's filled with even values starting at 0 */
+template <typename T>
+int is_filled_with_evens(T *p, int size) {
+  T   prev   = 0;
+  int result = 1; // default to true
+  int i;
+
+  if (p[0] != 0) {
+    result = 0;
+  }
+
+  for (i=1; result && i < size; i++) {
+      if (p[i] != prev + 2) {
+          result = 0;
+      }
+      else {
+          prev = p[i];
+      }
+  }
+
+  return result;
+}
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+int32_t *one_to_ten() {
+  int32_t *p = (int32_t*) malloc(sizeof(int32_t) * 10);
+  int i;
+  for (i=0; i<10; i++) {
+    p[i] = i+1;
+  }
+  return p;
+}
+
+double *ten_to_one() {
+  double *p = (double*) malloc(sizeof(double) * 10);
+  int i;
+  for (i=0; i< 10; i++) {
+    p[i] = 10.0 - (double) i;
+  }
+  return p;
+}
+
+int32_t *n_int_32s (int n) {
+  return (int32_t*) malloc(sizeof(int32_t) * n);
+}
+
+int16_t *n_int_16s(int n) {
+  return (int16_t*) malloc(sizeof(int16_t) * n);
+}
+
+int64_t *n_int_64s(int n) {
+  return (int64_t*) malloc(sizeof(int64_t) * n);
+}
+
+int is_filled_with_evens_16(int16_t *p, int size)
+{
+    return is_filled_with_evens(p, size);
+}
+
+int is_filled_with_evens_32(int32_t *p, int size)
+{
+    return is_filled_with_evens(p, size);
+}
+
+int is_filled_with_evens_64(int64_t *p, int size)
+{
+    return is_filled_with_evens(p, size);
+}
+
+#ifdef __cplusplus
+}
+#endif
+
diff --git a/tests/primitives/Backpermute.hs b/tests/primitives/Backpermute.hs
new file mode 100644
--- /dev/null
+++ b/tests/primitives/Backpermute.hs
@@ -0,0 +1,69 @@
+{-# LANGUAGE TypeOperators #-}
+
+module Backpermute where
+
+import Random
+
+import Control.Monad
+import Control.Exception
+import System.Random.MWC
+import Data.Array.Unboxed
+import Data.Array.Accelerate as Acc
+import Prelude               as P
+
+
+-- Tests
+-- -----
+
+reverseAcc :: Vector Float -> Acc (Vector Float)
+reverseAcc xs =
+  let xs' = use xs
+      len = unindex1 (shape xs')
+  in
+  backpermute (shape xs') (\ix -> index1 $ len - (unindex1 ix) - 1) xs'
+
+reverseRef :: UArray Int Float -> UArray Int Float
+reverseRef xs = listArray (bounds xs) (reverse (elems xs))
+
+
+transposeAcc :: Acc.Array DIM2 Float -> Acc (Acc.Array DIM2 Float)
+transposeAcc mat =
+  let mat' = use mat
+      swap = lift1 $ \(Z:.x:.y) -> Z:.y:.x :: Z:.Exp Int:.Exp Int
+  in
+  backpermute (swap $ shape mat') swap mat'
+
+transposeRef :: UArray (Int,Int) Float -> UArray (Int,Int) Float
+transposeRef mat =
+  let swap (x,y) = (y,x)
+      (u,v)      = bounds mat
+  in
+  array (swap u, swap v) [(swap ix, e) | (ix, e) <- assocs mat]
+
+
+-- Main
+-- ----
+
+run :: String -> Int -> IO (() -> UArray Int Float, () -> Acc (Vector Float))
+run alg n = withSystemRandom $ \gen -> do
+  vec  <- randomUArrayR (-1,1) gen n
+  vec' <- convertUArray vec
+  --
+  let go f g = return (\() -> f vec, \() -> g vec')
+  case alg of
+    "reverse" -> go reverseRef reverseAcc
+    _         -> error $ "unknown variant: " ++ alg
+
+run2d :: String -> Int -> IO (() -> UArray (Int,Int) Float, () -> Acc (Acc.Array DIM2 Float))
+run2d alg n = withSystemRandom $ \gen -> do
+  let n'    = P.round $ sqrt (P.fromIntegral n :: Double)
+      (u,v) = (n'*2, n'`div`2)
+  mat  <- listArray ((0,0), (u-1,v-1)) `fmap` replicateM (u*v) (uniformR (-1,1) gen)
+  mat' <- let m = fromIArray mat :: Acc.Array DIM2 Float
+          in  evaluate (m `Acc.indexArray` (Z:.0:.0)) >> return m
+  --
+  let go f g = return (\() -> f mat, \() -> g mat')
+  case alg of
+    "transpose" -> go transposeRef transposeAcc
+    _           -> error $ "unknown variant: " ++ alg
+
diff --git a/tests/primitives/Fold.hs b/tests/primitives/Fold.hs
new file mode 100644
--- /dev/null
+++ b/tests/primitives/Fold.hs
@@ -0,0 +1,84 @@
+{-# LANGUAGE FlexibleContexts, TypeOperators #-}
+
+module Fold where
+
+import Random
+
+import Control.Monad
+import Control.Exception
+import System.Random.MWC
+import Data.Array.Unboxed    hiding (Array)
+import Data.Array.Accelerate as Acc
+import Prelude		     as P
+
+
+-- one-dimension ah-ha-ha
+-- ----------------------
+
+toUA :: (IArray UArray a, IArray UArray b) => ([a] -> b) -> UArray Int a -> UArray () b
+toUA f = listArray ((),()) . return . f . elems
+
+
+sumAcc, prodAcc, maxAcc, minAcc :: Shape ix => Array (ix:.Int) Float -> Acc (Array ix Float)
+sumAcc  = Acc.fold (+) 0 . Acc.use
+prodAcc = Acc.fold (*) 1 . Acc.use
+maxAcc  = Acc.fold1 Acc.max . Acc.use
+minAcc  = Acc.fold1 Acc.min . Acc.use
+
+sumRef, prodRef, maxRef, minRef :: UArray Int Float -> UArray () Float
+sumRef  = toUA sum
+prodRef = toUA product
+maxRef  = toUA maximum
+minRef  = toUA minimum
+
+
+-- two-dimensions ah-ha-ha
+-- -----------------------
+
+foldU2D :: IArray UArray a => (a -> a -> a) -> a -> UArray (Int,Int) a -> UArray Int a
+foldU2D f z arr =
+  let (_,(m,_)) = bounds arr
+  in  accumArray f z (0,m) [ (i,e) | ((i,_),e) <- assocs arr ]
+
+sum2DRef, prod2DRef :: UArray (Int,Int) Float -> UArray Int Float
+sum2DRef  = foldU2D (+) 0
+prod2DRef = foldU2D (*) 1
+
+
+-- Main
+-- ----
+run :: String -> Int -> IO (() -> UArray () Float, () -> Acc (Scalar Float))
+run alg n = withSystemRandom $ \gen -> do
+  vec  <- randomUArrayR (-1,1) gen n
+  vec' <- convertUArray vec
+  --
+  let go  f g = return (run_ref f vec, run_acc g vec')
+  case alg of
+    "sum"        -> go sumRef sumAcc
+    "product"    -> go prodRef prodAcc
+    "maximum"    -> go maxRef maxAcc
+    "minimum"    -> go minRef minAcc
+    x            -> error $ "unknown variant: " ++ x
+  where
+    {-# NOINLINE run_ref #-}
+    run_ref f xs () = f xs
+    run_acc f xs () = f xs
+
+run2d :: String -> Int -> IO (() -> UArray Int Float, () -> Acc (Vector Float))
+run2d alg n = withSystemRandom $ \gen -> do
+  let u = P.floor . sqrt $ (P.fromIntegral n :: Double)
+      v = 2*u+1 :: Int
+  mat  <- listArray ((0,0), (u-1,v-1)) `fmap` replicateM (u*v) (uniformR (-1,1) gen)
+  mat' <- let m = fromIArray mat :: Array DIM2 Float
+          in  evaluate (m `Acc.indexArray` (Z:.0:.0)) >> return m
+  --
+  let go f g = return (run_ref f mat, run_acc g mat')
+  case alg of
+    "sum-2d"     -> go sum2DRef sumAcc
+    "product-2d" -> go prod2DRef prodAcc
+    x            -> error $ "unknown variant: " ++ x
+  where
+    {-# NOINLINE run_ref #-}
+    run_ref f xs () = f xs
+    run_acc f xs () = f xs
+
diff --git a/tests/primitives/Map.hs b/tests/primitives/Map.hs
new file mode 100644
--- /dev/null
+++ b/tests/primitives/Map.hs
@@ -0,0 +1,52 @@
+{-# LANGUAGE FlexibleContexts #-}
+
+module Map where
+
+import Random
+
+import System.Random.MWC
+import Data.Array.Unboxed
+import Data.Array.Accelerate as Acc
+
+
+-- Tests
+-- -----
+sqAcc, absAcc :: Vector Float -> Acc (Vector Float)
+absAcc = Acc.map abs . Acc.use
+sqAcc  = Acc.map (\x -> x * x) . Acc.use
+
+plusAcc :: Exp Float -> Vector Float -> Acc (Vector Float)
+plusAcc alpha = Acc.map (+ alpha) . Acc.use
+
+
+toUA :: (IArray UArray a, IArray UArray b) => ([a] -> [b]) -> UArray Int a -> UArray Int b
+toUA f xs = listArray (bounds xs) $ f (elems xs)
+
+sqRef, absRef :: UArray Int Float -> UArray Int Float
+absRef = toUA $ Prelude.map abs
+sqRef  = toUA $ Prelude.map (\x -> x*x)
+
+plusRef :: Float -> UArray Int Float -> UArray Int Float
+plusRef alpha = toUA $ Prelude.map (+alpha)
+
+
+-- Main
+-- ----
+run :: String -> Int -> IO (() -> UArray Int Float, () -> Acc (Vector Float))
+run alg n = withSystemRandom $ \gen -> do
+  vec   <- randomUArrayR (-1,1) gen n
+  vec'  <- convertUArray vec
+  alpha <- uniform gen
+  --
+  let go f g = return (run_ref f vec, run_acc g vec')
+  case alg of
+    "abs"    -> go absRef absAcc
+    "plus"   -> go (plusRef alpha) (plusAcc $ constant alpha)
+    "square" -> go sqRef sqAcc
+    x        -> error $ "unknown variant: " ++ x
+
+  where
+    {-# NOINLINE run_ref #-}
+    run_ref f xs () = f xs
+    run_acc f xs () = f xs
+
diff --git a/tests/primitives/Permute.hs b/tests/primitives/Permute.hs
new file mode 100644
--- /dev/null
+++ b/tests/primitives/Permute.hs
@@ -0,0 +1,39 @@
+module Permute where
+
+import Random
+
+import System.Random.MWC
+import Data.Array.Unboxed
+import Data.Array.Accelerate as Acc
+import Prelude               as P
+
+
+-- Tests
+-- -----
+
+histogramAcc :: (Int,Int) -> Vector Float -> Acc (Vector Int)
+histogramAcc (m,n) vec =
+  let vec'  = use vec
+      zeros = generate (constant (Z:. n-m)) (const 0)
+      ones  = generate (shape vec') (const 1)
+  in
+  permute (+) zeros (\ix -> index1 $ Acc.floor (vec' Acc.! ix)) ones
+
+histogramRef :: (Int,Int) -> UArray Int Float -> UArray Int Int
+histogramRef (m,n) vec =
+  accumArray (+) 0 (0,n-m-1) [(P.floor e, 1) | e <- elems vec]
+
+
+-- Main
+-- ----
+run :: String -> Int -> IO (() -> UArray Int Int, () -> Acc (Vector Int))
+run alg n = withSystemRandom $ \gen -> do
+  vec  <- randomUArrayR (0,100::Float) gen n
+  vec' <- convertUArray vec
+  --
+  let go f g = return (\() -> f vec, \() -> g vec')
+  case alg of
+    "histogram" -> go (histogramRef (0,100)) (histogramAcc (0,100))
+    _           -> error $ "unknown variant: " ++ alg
+
+
diff --git a/tests/primitives/ScanSeg.hs b/tests/primitives/ScanSeg.hs
new file mode 100644
--- /dev/null
+++ b/tests/primitives/ScanSeg.hs
@@ -0,0 +1,58 @@
+
+module ScanSeg where
+
+import Random
+
+import System.IO
+import System.Random.MWC
+import Data.Array.Unboxed
+import Data.Array.Accelerate as Acc
+import Prelude               as P
+
+
+-- Segmented prefix-sum
+-- --------------------
+prefixSumSegAcc :: Vector Float -> Segments -> Acc (Vector Float)
+prefixSumSegAcc xs seg
+  = let
+      xs'  = use xs
+      seg' = use seg
+    in
+    prescanlSeg (+) 0 xs' seg'
+
+
+prefixSumSegRef :: UArray Int Float -> UArray Int Int -> UArray Int Float
+prefixSumSegRef xs seg
+  = listArray (bounds xs)
+  $ list_prescanlSeg (+) 0 (elems xs) (elems seg)
+
+list_prescanlSeg :: (a -> a -> a) -> a -> [a] -> [Int] -> [a]
+list_prescanlSeg f x xs seg = concatMap (init . P.scanl f x) (split seg xs)
+  where
+    split [] _      = []
+    split _  []     = []
+    split (i:is) vs =
+      let (h,t) = splitAt i vs
+      in  h : split is t
+
+
+-- Main
+-- ----
+run :: String -> Int -> IO (() -> UArray Int Float, () -> Acc (Vector Float))
+run alg m = withSystemRandom $ \gen -> do
+  let n = P.round . sqrt $ (P.fromIntegral m :: Double)
+  seg  <- randomUArrayR (0,n) gen n
+  seg' <- convertUArray seg
+  let ne = sum (elems seg)
+  vec  <- randomUArrayR (-1,1) gen ne
+  vec' <- convertUArray vec
+  --
+  let go f g = return (run_ref f vec seg, run_acc g vec' seg')
+  case alg of
+    "sum" -> go prefixSumSegRef prefixSumSegAcc
+    x     -> error $ "unknown variant: " ++ x
+  where
+    {-# NOINLINE run_ref #-}
+    run_ref f xs seg () = f xs seg
+    run_acc f xs seg () = f xs seg
+
diff --git a/tests/primitives/Stencil.hs b/tests/primitives/Stencil.hs
new file mode 100644
--- /dev/null
+++ b/tests/primitives/Stencil.hs
@@ -0,0 +1,229 @@
+{-# LANGUAGE FlexibleContexts #-}
+
+module Stencil where
+
+import Random
+
+import Control.Monad
+import Control.Exception
+import System.Random.MWC
+
+import Data.Array.Unboxed hiding (Array)
+import qualified Data.Array.IArray as IArray
+
+import Data.Array.Accelerate hiding (min, max, round, fromIntegral)
+import qualified Data.Array.Accelerate as Acc
+
+
+
+-- Stencil Tests
+-- -------------
+
+-- 1D --------------------------------------------------------------------------
+
+stencil1D :: Floating a
+          => (a, a, a) -> a
+stencil1D (x, y, z) = (x + z - 2 * y) / 2
+
+test_stencil1D :: Int -> IO (() -> UArray Int Float, () -> Acc (Vector Float))
+test_stencil1D n = withSystemRandom $ \gen -> do
+  vec  <- randomUArrayR (-1,1) gen n
+  vec' <- convertUArray vec
+  return (\() -> run_ref vec, \() -> run_acc vec')
+  where
+    run_acc   = stencil stencil1D Clamp . use
+    run_ref v =
+      let (minx,maxx) = bounds v
+          clamp x     = minx `max` x `min` maxx
+
+          f ix = let x = v IArray.! clamp (ix-1)
+                     y = v IArray.! ix
+                     z = v IArray.! clamp (ix+1)
+                 in
+                 (x + z - 2 * y) / 2
+      in
+      array (bounds v) [(ix, f ix) | ix <- indices v]
+
+
+-- 2D --------------------------------------------------------------------------
+
+stencil2D :: Floating (Exp a)
+          => Stencil3x3 a -> Exp a
+stencil2D ( (t1, t2, t3)
+          , (l , m,  r )
+          , (b1, b2, b3)
+          )
+          = (t1/2 + t2 + t3/2 + l + r + b1/2 + b2 + b3/2 - 4 * m) / 4
+
+test_stencil2D :: Int -> IO (() -> IArray.Array (Int,Int) Float, () -> Acc (Array DIM2 Float))
+test_stencil2D n2 = withSystemRandom $ \gen -> do
+  let n = round . (/3) . sqrt $ (fromIntegral n2 :: Double)
+      m = n * 4
+  mat  <- listArray ((0,0),(n-1,m-1)) `fmap` replicateM (n*m) (uniformR (-1,1) gen) :: IO (IArray.Array (Int,Int) Float)
+  mat' <- let v = fromIArray mat                                                    :: Array DIM2 Float
+          in  evaluate (v `indexArray` (Z:.0:.0)) >> return v
+  --
+  return (\() -> run_ref mat, \() -> run_acc mat')
+  where
+    run_acc     = stencil stencil2D (Constant 0) . use
+    run_ref arr =
+      let get ix
+            | inRange (bounds arr) ix = arr IArray.! ix
+            | otherwise               = 0
+
+          f (x,y) = let t1 = get (x-1,y-1)
+                        t2 = get (x,  y-1)
+                        t3 = get (x+1,y-1)
+                        l  = get (x-1,y)
+                        m  = get (x,  y)
+                        r  = get (x+1,y)
+                        b1 = get (x-1,y+1)
+                        b2 = get (x,  y+1)
+                        b3 = get (x+1,y+1)
+                    in
+                    (t1/2 + t2 + t3/2 + l + r + b1/2 + b2 + b3/2 - 4 * m) / 4
+      in
+      array (bounds arr) [(ix, f ix) | ix <- indices arr]
+
+
+
+stencil2D5 :: Floating (Exp a)
+           => Stencil3x3 a -> Exp a
+stencil2D5 ( (_, t, _)
+           , (l, m, r)
+           , (_, b, _)
+           )
+           = (t + l + r + b - 4 * m) / 4
+
+test_stencil2D5 :: Int -> IO (() -> IArray.Array (Int,Int) Float, () -> Acc (Array DIM2 Float))
+test_stencil2D5 n2 = withSystemRandom $ \gen -> do
+  let n = round . sqrt $ (fromIntegral n2 :: Double)
+  mat  <- listArray ((0,0),(n-1,n-1)) `fmap` replicateM (n*n) (uniformR (-1,1) gen) :: IO (IArray.Array (Int,Int) Float)
+  mat' <- let m = fromIArray mat                                                    :: Array DIM2 Float
+          in  evaluate (m `indexArray` (Z:.0:.0)) >> return m
+  --
+  return (\() -> run_ref mat, \() -> run_acc mat')
+  where
+    run_acc     = stencil stencil2D5 Clamp . use
+    run_ref arr =
+      let ((minx,miny),(maxx,maxy)) = bounds arr
+          clamp (x,y) = (minx `max` x `min` maxx
+                        ,miny `max` y `min` maxy)
+          f (x,y)     = let t = arr IArray.! clamp (x,y-1)
+                            b = arr IArray.! clamp (x,y+1)
+                            l = arr IArray.! clamp (x-1,y)
+                            r = arr IArray.! clamp (x+1,y)
+                            m = arr IArray.! (x,y)
+                        in
+                        (t + l + r + b - 4 * m) / 4
+      in
+      array (bounds arr) [(ix, f ix) | ix <- indices arr]
+
+
+
+stencil2Dpair :: Stencil3x3 (Int,Float) -> Exp Float
+stencil2Dpair ( (_, _, _)
+              , (x, _, _)
+              , (y, _, z)
+              )
+  = let (x1,x2) = unlift x :: (Exp Int, Exp Float)
+        (y1,y2) = unlift y
+        (z1,z2) = unlift z
+    in
+    (x2 * Acc.fromIntegral x1 + y2 * Acc.fromIntegral z1 - z2 * Acc.fromIntegral y1)
+
+test_stencil2Dpair :: Int -> IO (() -> IArray.Array (Int,Int) Float, () -> Acc (Array DIM2 Float))
+test_stencil2Dpair n2 = withSystemRandom $ \gen -> do
+  let n = round (fromIntegral n2 ** 0.5 :: Double)
+      m = 2 * n
+  mat  <- listArray ((0,0),(n-1,m-1)) `fmap` replicateM (n*m) (uniformR ((-100,0), (100,0)) gen) :: IO (IArray.Array (Int,Int) (Int,Float))
+  mat' <- let a = fromIArray mat
+          in  evaluate (a `indexArray` (Z:.0:.0)) >> return a
+  --
+  return (\() -> run_ref mat, \() -> run_acc mat')
+  where
+    run_acc     = stencil stencil2Dpair Wrap . use
+    run_ref arr =
+      let ((minx,miny),(maxx,maxy)) = bounds arr
+          wrap (m,n) i
+            | i < m     = n + (i-m)
+            | i > n     = i - n + m
+            | otherwise = i
+
+          get (x,y) = arr IArray.! ( wrap (minx,maxx) x, wrap (miny,maxy) y)
+          f   (x,y) = let (a1,a2) = get (x-1,y)
+                          (b1,b2) = get (x-1,y+1)
+                          (c1,c2) = get (x+1,y+1)
+                      in
+                      (a2 * fromIntegral a1 + b2 * fromIntegral c1 - c2 * fromIntegral b1)
+      in
+      array (bounds arr) [(ix, f ix) | ix <- indices arr]
+
+
+-- 3D --------------------------------------------------------------------------
+
+stencil3D :: Num (Exp a)
+          => Stencil3x3x3 a -> Exp a
+stencil3D (front, back, _) =      -- 'b4' is the focal point
+  let ((f1, f2, _),
+       (f3, f4, _),
+       _          ) = front
+      ((b1, b2, _),
+       (b3, b4, _),
+       _          ) = back
+  in
+  f1 + f2 + f3 + f4 + b1 + b2 + b3 + b4
+
+test_stencil3D :: Int -> IO (() -> UArray (Int,Int,Int) Float, () -> Acc (Array DIM3 Float))
+test_stencil3D n3 = withSystemRandom $ \gen -> do
+  let u = round (fromIntegral n3 ** (1/3) :: Double)
+      v = u `div` 2
+      w = u * 3
+  arr  <- listArray ((0,0,0), (u-1, v-1, w-1)) `fmap` replicateM (u*v*w) (uniformR (-1,1) gen) :: IO (UArray (Int,Int,Int) Float)
+  arr' <- let a = fromIArray arr
+          in  evaluate (a `indexArray` (Z:.0:.0:.0)) >> return a
+  --
+  return (\() -> run_ref arr, \() -> run_acc arr')
+  where
+    run_acc     = stencil stencil3D Mirror . use
+    run_ref arr =
+      let ((minx,miny,minz),(maxx,maxy,maxz)) = bounds arr
+          mirror (m,n) i
+            | i < m     = -i + m
+            | i > n     = n - (i-n+2)
+            | otherwise = i
+
+          get (x,y,z) = arr IArray.! ( mirror (minx,maxx) x, mirror (miny,maxy) y, mirror (minz,maxz) z)
+          f   (x,y,z) = let f1 = get (x-1,y-1,z-1)
+                            f2 = get (x,  y-1,z-1)
+                            f3 = get (x-1,y,  z-1)
+                            f4 = get (x,  y,  z-1)
+                            b1 = get (x-1,y-1,z  )
+                            b2 = get (x,  y-1,z  )
+                            b3 = get (x-1,y,  z  )
+                            b4 = get (x,  y,  z  )
+                        in
+                        f1 + f2 + f3 + f4 + b1 + b2 + b3 + b4
+      in
+      array (bounds arr) [(ix, f ix) | ix <- indices arr]
+
+
+-- Main
+-- ----
+
+run :: String -> Int -> IO (() -> UArray Int Float, () -> Acc (Vector Float))
+run "1D" = test_stencil1D
+run x    = error $ "unknown variant: " ++ x
+
+
+run2D :: String -> Int -> IO (() -> IArray.Array (Int,Int) Float, () -> Acc (Array DIM2 Float))
+run2D "2D"        = test_stencil2D
+run2D "3x3-cross" = test_stencil2D5
+run2D "3x3-pair"  = test_stencil2Dpair
+run2D x    = error $ "unknown variant: " ++ x
+
+
+run3D :: String -> Int -> IO (() -> UArray (Int,Int,Int) Float, () -> Acc (Array DIM3 Float))
+run3D "3D" = test_stencil3D
+run3D x    = error $ "unknown variant: " ++ x
+
diff --git a/tests/primitives/Stencil2.hs b/tests/primitives/Stencil2.hs
new file mode 100644
--- /dev/null
+++ b/tests/primitives/Stencil2.hs
@@ -0,0 +1,70 @@
+{-# LANGUAGE FlexibleContexts #-}
+
+module Stencil2 where
+
+import Control.Monad
+import Control.Exception
+import System.Random.MWC
+import Data.Array.Unboxed           hiding (Array)
+import Data.Array.Accelerate        hiding (round, min, max, fromIntegral)
+import qualified Data.Array.IArray  as IArray
+
+
+
+stencil2D2 :: Floating (Exp a) => Stencil3x3 a -> Stencil3x3 a -> Exp a
+stencil2D2 ((_,t,_), (_,x,_), (_,b,_))
+           ((_,_,_), (l,y,r), (_,_,_)) = t + b + l + r - ((x+y) / 2)
+
+
+test_stencil2_2D :: Int -> IO (() -> UArray (Int,Int) Float, () -> Acc (Array DIM2 Float))
+test_stencil2_2D n2 = withSystemRandom $ \gen -> do
+  let n = round $ sqrt (fromIntegral n2 :: Double)
+      m = n * 2
+      u = m `div` 3
+      v = n + m
+  m1  <- listArray ((0,0),(n-1,m-1)) `fmap` replicateM (n*m) (uniformR (-1,1) gen) :: IO (UArray (Int,Int) Float)
+  m2  <- listArray ((0,0),(u-1,v-1)) `fmap` replicateM (u*v) (uniformR (-1,1) gen) :: IO (UArray (Int,Int) Float)
+  m1' <- let m1' = fromIArray m1 in evaluate (m1' `indexArray` (Z:.0:.0)) >> return m1'
+  m2' <- let m2' = fromIArray m2 in evaluate (m2' `indexArray` (Z:.0:.0)) >> return m2'
+  --
+  return (\() -> run_ref m1 m2, \() -> run_acc m1' m2')
+  where
+    run_acc xs ys = stencil2 stencil2D2 Mirror (use xs) Wrap (use ys)
+    run_ref xs ys =
+      let (_,(n,m)) = bounds xs
+          (_,(u,v)) = bounds ys
+          sh        = ((0,0), (n `min` u, m `min` v))
+
+          -- boundary conditions are placed on the *source* arrays
+          --
+          get1 (x,y) = xs IArray.! (mirror n x, mirror m y)
+          get2 (x,y) = ys IArray.! (wrap   u x, wrap   v y)
+
+          mirror sz i
+            | i < 0     = -i
+            | i > sz    = sz - (i-sz)
+            | otherwise = i
+
+          wrap sz i
+            | i < 0     = sz + i + 1
+            | i > sz    = i - sz - 1
+            | otherwise = i
+
+          f (ix,iy) = let t = get1 (ix,  iy-1)
+                          b = get1 (ix,  iy+1)
+                          x = get1 (ix,  iy)
+                          l = get2 (ix-1,iy)
+                          r = get2 (ix+1,iy)
+                          y = get2 (ix,  iy)
+                    in
+                    t + b + l + r - ((x+y) / 2)
+      in
+      array sh [(ix, f ix) | ix <- range sh]
+
+-- Main
+-- ----
+
+run2D :: String -> Int -> IO (() -> UArray (Int,Int) Float, () -> Acc (Array DIM2 Float))
+run2D "2D" = test_stencil2_2D
+run2D x    = error $ "unknown variant: " ++ x
+
diff --git a/tests/primitives/Zip.hs b/tests/primitives/Zip.hs
new file mode 100644
--- /dev/null
+++ b/tests/primitives/Zip.hs
@@ -0,0 +1,35 @@
+
+module Zip where
+
+import Random
+
+import System.Random.MWC
+import Data.Array.Unboxed       as IArray
+import Data.Array.Accelerate    as Acc hiding (min)
+
+
+-- Tests
+-- -----
+
+zipAcc :: Vector Float -> Vector Int -> Acc (Vector (Float,Int))
+zipAcc xs ys = Acc.zip (use xs) (use ys)
+
+
+zipRef :: UArray Int Float -> UArray Int Int -> IArray.Array Int (Float,Int)
+zipRef xs ys =
+  let mn      = bounds xs
+      uv      = bounds ys
+      newSize = (0, (rangeSize mn `min` rangeSize uv) - 1)
+  in
+  listArray newSize $ Prelude.zip (elems xs) (elems ys)
+
+-- Main
+-- ----
+run :: Int -> IO (() -> IArray.Array Int (Float,Int), () -> Acc (Vector (Float,Int)))
+run n = withSystemRandom $ \gen -> do
+  xs  <- randomUArrayR (-1,1) gen n
+  ys  <- randomUArrayR (-1,1) gen n
+  xs' <- convertUArray xs
+  ys' <- convertUArray ys
+  return $ (\() -> zipRef xs ys, \() -> zipAcc xs' ys')
+
diff --git a/tests/primitives/ZipWith.hs b/tests/primitives/ZipWith.hs
new file mode 100644
--- /dev/null
+++ b/tests/primitives/ZipWith.hs
@@ -0,0 +1,42 @@
+{-# LANGUAGE ParallelListComp #-}
+
+module ZipWith where
+
+import Random
+
+import System.Random.MWC
+import Data.Array.Unboxed
+import Data.Array.Accelerate as Acc hiding (min)
+
+
+-- Tests
+-- -----
+
+plusAcc :: Vector Float -> Vector Float -> Acc (Vector Float)
+plusAcc xs ys = Acc.zipWith (+) (use xs) (use ys)
+
+plusRef :: UArray Int Float -> UArray Int Float -> UArray Int Float
+plusRef = zipWithRef (+)
+
+zipWithRef :: (IArray array a, IArray array b, IArray array c)
+           => (a -> b -> c) -> array Int a -> array Int b -> array Int c
+zipWithRef f xs ys =
+  let mn      = bounds xs
+      uv      = bounds ys
+      newSize = (0, (rangeSize mn `min` rangeSize uv) - 1)
+  in
+  listArray newSize [f x y | x <- elems xs | y <- elems ys]
+
+-- Main
+-- ----
+run :: String -> Int -> IO (() -> UArray Int Float, () -> Acc (Vector Float))
+run alg n = withSystemRandom $ \gen -> do
+  xs  <- randomUArrayR (-1,1) gen n
+  ys  <- randomUArrayR (-1,1) gen n
+  xs' <- convertUArray xs
+  ys' <- convertUArray ys
+  let go f g = return (\() -> f xs ys, \() -> g xs' ys')
+  case alg of
+    "plus" -> go plusRef plusAcc
+    _      -> error $ "unknown variant: " ++ alg
+
diff --git a/tests/simple/BlackScholes.hs b/tests/simple/BlackScholes.hs
new file mode 100644
--- /dev/null
+++ b/tests/simple/BlackScholes.hs
@@ -0,0 +1,88 @@
+
+module BlackScholes where
+
+import Random
+
+import Prelude               as P
+import System.Random.MWC
+import Data.Array.IArray     as IArray
+import Data.Array.Accelerate as Acc
+
+
+riskfree, volatility :: Float
+riskfree   = 0.02
+volatility = 0.30
+
+-- Black-Scholes option pricing
+-------------------------------
+
+horner :: Num a => [a] -> a -> a
+horner coeff x = foldr1 madd coeff
+  where
+    madd a b = b*x + a
+
+cnd' :: Floating a => a -> a
+cnd' d =
+  let poly     = horner coeff
+      coeff    = [0.0,0.31938153,-0.356563782,1.781477937,-1.821255978,1.330274429]
+      rsqrt2pi = 0.39894228040143267793994605993438
+      k        = 1.0 / (1.0 + 0.2316419 * abs d)
+  in
+  rsqrt2pi * exp (-0.5*d*d) * poly k
+
+
+blackscholesAcc :: Vector (Float, Float, Float) -> Acc (Vector (Float, Float))
+blackscholesAcc xs = Acc.map go (Acc.use xs)
+  where
+  go x =
+    let (price, strike, years) = Acc.unlift x
+        r     = Acc.constant riskfree
+        v     = Acc.constant volatility
+        sqrtT = sqrt years
+        d1    = (log (price / strike) + (r + 0.5 * v * v) * years) / (v * sqrtT)
+        d2    = d1 - v * sqrtT
+        cnd d = d >* 0 ? (1.0 - cnd' d, cnd' d)
+        cndD1 = cnd d1
+        cndD2 = cnd d2
+        expRT = exp (-r * years)
+    in
+    Acc.lift ( price * cndD1 - strike * expRT * cndD2
+             , strike * expRT * (1.0 - cndD2) - price * (1.0 - cndD1))
+
+
+blackscholesRef :: IArray.Array Int (Float,Float,Float) -> IArray.Array Int (Float,Float)
+blackscholesRef xs = listArray (bounds xs) [ go x | x <- elems xs ]
+  where
+  go (price, strike, years) =
+    let r     = riskfree
+        v     = volatility
+        sqrtT = sqrt years
+        d1    = (log (price / strike) + (r + 0.5 * v * v) * years) / (v * sqrtT)
+        d2    = d1 - v * sqrtT
+        cnd d = if d > 0 then 1.0 - cnd' d else cnd' d
+        cndD1 = cnd d1
+        cndD2 = cnd d2
+        expRT = exp (-r * years)
+    in
+    ( price * cndD1 - strike * expRT * cndD2
+    , strike * expRT * (1.0 - cndD2) - price * (1.0 - cndD1))
+
+
+-- Main
+-- ----
+
+run :: Int -> IO (() -> IArray.Array Int (Float,Float), () -> Acc (Vector (Float,Float)))
+run n = withSystemRandom $ \gen -> do
+  v_sp <- randomUArrayR (5,30)    gen n
+  v_os <- randomUArrayR (1,100)   gen n
+  v_oy <- randomUArrayR (0.25,10) gen n
+
+  let v_psy = listArray (0,n-1) $ P.zip3 (elems v_sp) (elems v_os) (elems v_oy)
+      a_psy = Acc.fromIArray v_psy
+  --
+  return (run_ref v_psy, run_acc a_psy)
+  where
+    {-# NOINLINE run_ref #-}
+    run_ref psy () = blackscholesRef psy
+    run_acc psy () = blackscholesAcc psy
+
diff --git a/tests/simple/DotP.hs b/tests/simple/DotP.hs
new file mode 100644
--- /dev/null
+++ b/tests/simple/DotP.hs
@@ -0,0 +1,44 @@
+{-# LANGUAGE ParallelListComp #-}
+
+module DotP where
+
+import Random
+
+import System.Random.MWC
+import Data.Array.Unboxed
+import Data.Array.Accelerate as Acc
+
+
+-- Dot product
+-- -----------
+dotpAcc :: Vector Float -> Vector Float -> Acc (Scalar Float)
+dotpAcc xs ys
+  = let
+      xs' = use xs
+      ys' = use ys
+    in
+    Acc.fold (+) 0 (Acc.zipWith (*) xs' ys')
+
+dotpRef :: UArray Int Float
+        -> UArray Int Float
+        -> UArray ()  Float
+dotpRef xs ys
+  = listArray ((), ()) [sum [x * y | x <- elems xs | y <- elems ys]]
+
+
+-- Main
+-- ----
+
+run :: Int -> IO (() -> UArray () Float, () -> Acc (Scalar Float))
+run n = withSystemRandom $ \gen -> do
+  v1  <- randomUArrayR (-1,1) gen n
+  v2  <- randomUArrayR (-1,1) gen n
+  v1' <- convertUArray v1
+  v2' <- convertUArray v2
+  --
+  return (run_ref v1 v2, run_acc v1' v2')
+  where
+    {-# NOINLINE run_ref #-}
+    run_ref xs ys () = dotpRef xs ys
+    run_acc xs ys () = dotpAcc xs ys
+
diff --git a/tests/simple/Filter.hs b/tests/simple/Filter.hs
new file mode 100644
--- /dev/null
+++ b/tests/simple/Filter.hs
@@ -0,0 +1,55 @@
+{-# LANGUAGE FlexibleContexts #-}
+
+module Filter where
+
+import Random
+
+import System.Random.MWC
+import Data.Array.Unboxed       (IArray, UArray, elems, listArray)
+import Data.Array.Accelerate    as Acc
+
+
+-- Filter
+-- ------
+filterAcc :: Elt a
+          => (Exp a -> Exp Bool)
+          -> Vector a
+          -> Acc (Vector a)
+filterAcc p vec
+  = let arr              = Acc.use vec
+        flags            = Acc.map (boolToInt . p) arr
+        (targetIdx, len) = Acc.scanl' (+) 0 flags
+        arr'             = Acc.backpermute (index1 $ the len) id arr
+    in
+    Acc.permute const arr' (\ix -> flags!ix ==* 0 ? (ignore, index1 $ targetIdx!ix)) arr
+    -- FIXME: This is abusing 'permute' in that the first two arguments are
+    --        only justified because we know the permutation function will
+    --        write to each location in the target exactly once.
+    --        Instead, we should have a primitive that directly encodes the
+    --        compaction pattern of the permutation function.
+
+
+filterRef :: IArray UArray e
+          => (e -> Bool)
+          -> UArray Int e
+          -> UArray Int e
+filterRef p xs
+  = let xs' = Prelude.filter p (elems xs)
+    in
+    listArray (0, Prelude.length xs' - 1) xs'
+
+
+-- Main
+-- ----
+
+run :: Int -> IO (() -> UArray Int Float, () -> Acc (Vector Float))
+run n = withSystemRandom $ \gen -> do
+  vec  <- randomUArrayR (-1,1) gen n
+  vec' <- convertUArray vec
+  --
+  return (run_ref vec, run_acc vec')
+  where
+    {-# NOINLINE run_ref #-}
+    run_ref xs () = filterRef (> 0) xs
+    run_acc xs () = filterAcc (>*0) xs
+
diff --git a/tests/simple/Radix.hs b/tests/simple/Radix.hs
new file mode 100644
--- /dev/null
+++ b/tests/simple/Radix.hs
@@ -0,0 +1,86 @@
+--
+-- Radix sort for a subclass of element types
+--
+
+module Radix where
+
+import Random
+
+import qualified Prelude
+import Prelude               hiding (zip, map, scanl, scanr, zipWith, fst)
+import Data.Bits             hiding (shiftL, shiftR, bit, testBit)
+import Data.Array.Accelerate as Acc
+
+import Data.List             (sort)
+import Data.Array.Unboxed    (IArray, UArray, listArray, bounds, elems)
+import System.Random.MWC
+import Unsafe.Coerce
+
+
+-- Radix sort
+-- ----------
+
+class Elt e => Radix e where
+    passes :: Exp e   -> Int                -- Haskell-side control needs to know this
+    radix  :: Exp Int -> Exp e -> Exp Int
+
+instance Radix Int where                    -- may be 32- or 64-bit
+    passes    = bitSize . (undefined :: Exp t -> t)
+    radix i e = i ==* (passes' e - 1) ? (radix' (e `xor` minBound), radix' e)
+      where
+        radix' x = (x `shiftR` i) .&. 1
+        passes'  = constant . passes
+
+-- For IEEE-754 floating-point representation. Unsafe, but widely supported.
+--
+instance Radix Float where
+    passes _   = 32
+    radix i e  = let x = (unsafeCoerce e :: Exp Int32)
+                 in  i ==* 31 ? (radix' (x `xor` minBound), radix' (floatFlip x))
+      where
+        floatFlip x = x `testBit` 31 ? (complement x, x)  -- twos-complement negative numbers
+        radix'    x = x `testBit` i  ? (1,0)
+
+
+--
+-- A simple (parallel) radix sort implementation [1].
+--
+-- [1] G. E. Blelloch. "Prefix sums and their applications." Technical Report
+--     CMU-CS-90-190. Carnegie Mellon University. 1990.
+--
+sortAcc :: Radix a => Vector a -> Acc (Vector a)
+sortAcc = sortAccBy id
+
+sortAccBy :: (Elt a, Radix r) => (Exp a -> Exp r) -> Vector a -> Acc (Vector a)
+sortAccBy rdx arr = foldr1 (>->) (Prelude.map radixPass [0..p-1]) (use arr)
+  where
+    n = constant $ (arraySize $ arrayShape arr) - 1
+    p = passes . rdx . (undefined :: Vector e -> Exp e) $ arr
+
+    deal f x      = let (a,b) = unlift x in (f ==* 0) ? (a,b)
+    radixPass k v = let flags = map (radix (constant k) . rdx) v
+                        idown = prescanl (+) 0 . map (xor 1) $ flags
+                        iup   = map (n-) . prescanr (+) 0    $ flags
+                        index = zipWith deal flags (zip idown iup)
+                    in
+                    permute const v (\ix -> index1 (index!ix)) v
+
+
+sortRef :: UArray Int Int -> UArray Int Int
+sortRef xs = listArray (bounds xs) $ sort (elems xs)
+
+
+-- Main
+-- ----
+
+run :: Int -> IO (() -> UArray Int Int, () -> Acc (Vector Int))
+run n = withSystemRandom $ \gen -> do
+  vec  <- randomUArrayR (minBound,maxBound) gen n
+  vec' <- convertUArray vec
+  --
+  return (run_ref vec, run_acc vec')
+  where
+    {-# NOINLINE run_ref #-}
+    run_ref xs () = sortRef xs
+    run_acc xs () = sortAcc xs
+
diff --git a/tests/simple/SASUM.hs b/tests/simple/SASUM.hs
new file mode 100644
--- /dev/null
+++ b/tests/simple/SASUM.hs
@@ -0,0 +1,35 @@
+
+module SASUM where
+
+import Random
+
+import System.Random.MWC
+import Data.Array.Unboxed
+import Data.Array.Accelerate as Acc
+
+
+-- Sum of absolute values
+-- ----------------------
+sasumAcc :: Vector Float -> Acc (Scalar Float)
+sasumAcc xs
+  = Acc.fold (+) 0 . Acc.map abs $ Acc.use xs
+
+sasumRef :: UArray Int Float -> UArray () Float
+sasumRef xs
+  = listArray ((), ()) [Prelude.sum . Prelude.map abs $ elems xs]
+
+
+-- Main
+-- ----
+
+run :: Int -> IO (() -> UArray () Float, () -> Acc (Scalar Float))
+run n = withSystemRandom $ \gen -> do
+  vec  <- randomUArrayR (-1,1) gen n
+  vec' <- convertUArray vec
+  --
+  return (run_ref vec, run_acc vec')
+  where
+    {-# NOINLINE run_ref #-}
+    run_ref xs () = sasumRef xs
+    run_acc xs () = sasumAcc xs
+
diff --git a/tests/simple/SAXPY.hs b/tests/simple/SAXPY.hs
new file mode 100644
--- /dev/null
+++ b/tests/simple/SAXPY.hs
@@ -0,0 +1,42 @@
+{-# LANGUAGE ParallelListComp #-}
+
+module SAXPY where
+
+import Random
+
+import System.Random.MWC
+import Data.Array.Unboxed
+import Data.Array.Accelerate as Acc
+
+-- SAXPY
+-- -----
+saxpyAcc :: Float -> Vector Float -> Vector Float -> Acc (Vector Float)
+saxpyAcc alpha xs ys
+  = let
+      xs' = use xs
+      ys' = use ys
+    in
+    Acc.zipWith (\x y -> constant alpha * x + y) xs' ys'
+
+saxpyRef :: Float -> UArray Int Float -> UArray Int Float -> UArray Int Float
+saxpyRef alpha xs ys
+  = listArray (bounds xs) [alpha * x + y | x <- elems xs | y <- elems ys]
+
+
+-- Main
+-- ----
+
+run :: Int -> IO (() -> UArray Int Float, () -> Acc (Vector Float))
+run nelements = withSystemRandom $ \gen -> do
+  v1    <- randomUArrayR (-1,1) gen nelements
+  v2    <- randomUArrayR (-1,1) gen nelements
+  v1'   <- convertUArray v1
+  v2'   <- convertUArray v2
+  alpha <- uniform gen
+  --
+  return (run_ref alpha v1 v2, run_acc alpha v1' v2')
+  where
+    {-# NOINLINE run_ref #-}
+    run_ref alpha xs ys () = saxpyRef alpha xs ys
+    run_acc alpha xs ys () = saxpyAcc alpha xs ys
+
diff --git a/tests/simple/SMVM.hs b/tests/simple/SMVM.hs
new file mode 100644
--- /dev/null
+++ b/tests/simple/SMVM.hs
@@ -0,0 +1,78 @@
+{-# LANGUAGE FlexibleContexts #-}
+
+module SMVM where
+
+import Random
+import SMVM.Matrix
+
+import System.Random.MWC
+import Data.Array.Unboxed
+import Data.Array.Accelerate           (Vector, Segments, Acc)
+import qualified Data.Array.Accelerate as Acc
+import qualified Data.Vector.Unboxed   as V
+
+
+-- Sparse-matrix vector multiplication
+-- -----------------------------------
+
+type SparseVector a = (Vector Int, Vector a)
+type SparseMatrix a = (Segments, SparseVector a)
+
+smvmAcc :: SparseMatrix Float -> Vector Float -> Acc (Vector Float)
+smvmAcc (segd', (inds', vals')) vec'
+  = let
+      segd     = Acc.use segd'
+      inds     = Acc.use inds'
+      vals     = Acc.use vals'
+      vec      = Acc.use vec'
+      ---
+      vecVals  = Acc.backpermute (Acc.shape inds) (\i -> Acc.index1 $ inds Acc.! i) vec
+      products = Acc.zipWith (*) vecVals vals
+    in
+    Acc.foldSeg (+) 0 products segd
+
+
+-- The reference version will be slow, with many conversions between
+-- array/vector/list representations. This will likely skew heap usage
+-- calculations, but oh well...
+--
+type USparseMatrix a = (UArray Int Int, (UArray Int Int, UArray Int a))
+
+smvmRef :: USparseMatrix Float -> UArray Int Float -> UArray Int Float
+smvmRef (segd, (inds, values)) vec
+  = listArray (0, rangeSize (bounds segd) - 1)
+    [sum [ values!i * vec!(inds!i) | i <- range seg] | seg <- segd' ]
+  where
+    segbegin = scanl  (+) 0 $ elems segd
+    segend   = scanl1 (+)   $ elems segd
+    segd'    = zipWith (\x y -> (x,y-1)) segbegin segend
+
+
+-- Main
+-- ----
+
+run :: Maybe FilePath -> IO (() -> UArray Int Float, () -> Acc (Vector Float))
+run f = withSystemRandom $ \gen ->  do
+  -- sparse-matrix
+  (segd', smat') <- maybe (randomCSRMatrix gen 512 512) (readCSRMatrix gen) f
+  let (ind',val') = V.unzip smat'
+
+  segd <- convertVector segd'
+  ind  <- convertVector ind'
+  val  <- convertVector val'
+  let smat = (segd, (ind,val))
+
+  -- vector
+  vec' <- uniformVector gen (V.length segd')
+  vec  <- convertVector vec'
+
+  -- multiply!
+  return (run_ref (v2a segd', (v2a ind',v2a val')) (v2a vec'), run_acc smat vec)
+  where
+    {-# NOINLINE run_ref #-}
+    run_ref smat vec () = smvmRef smat vec
+    run_acc smat vec () = smvmAcc smat vec
+    --
+    v2a :: (V.Unbox a, IArray UArray a) => V.Vector a -> UArray Int a
+    v2a vec = listArray (0, V.length vec - 1) $ V.toList vec
+
diff --git a/tests/simple/SMVM/Matrix.hs b/tests/simple/SMVM/Matrix.hs
new file mode 100644
--- /dev/null
+++ b/tests/simple/SMVM/Matrix.hs
@@ -0,0 +1,75 @@
+{-# LANGUAGE BangPatterns, TupleSections #-}
+
+module SMVM.Matrix (readCSRMatrix, randomCSRMatrix) where
+
+import Random
+import SMVM.MatrixMarket
+import System.Random.MWC
+import System.IO.Unsafe
+
+import Data.Vector.Unboxed (Vector)
+import qualified Data.Vector.Unboxed            as V
+import qualified Data.Vector.Unboxed.Mutable    as M
+import qualified Data.Vector.Algorithms.Intro   as V
+
+type CSRMatrix a = (Vector Int, Vector (Int,a))
+
+
+-- Read a sparse matrix from a MatrixMarket file. Pattern matrices are filled
+-- with random numbers in the range (-1,1).
+--
+{-# INLINE readCSRMatrix #-}
+readCSRMatrix :: GenIO -> FilePath -> IO (CSRMatrix Float)
+readCSRMatrix gen file = do
+  mtx <- readMatrix file
+  case mtx of
+    (RealMatrix    dim l vals) -> csr dim l vals
+    (PatternMatrix dim l ix)   -> csr dim l =<< mapM' (\(a,b) -> (a,b,) `fmap` uniformR (-1,1) gen) ix
+    (IntMatrix _ _ _)          -> error "IntMatrix type not supported"
+    (ComplexMatrix _ _ _)      -> error "ComplexMatrix type not supported"
+
+
+-- A randomly generated matrix of given size
+--
+{-# INLINE randomCSRMatrix #-}
+randomCSRMatrix :: GenIO -> Int -> Int -> IO (CSRMatrix Float)
+randomCSRMatrix gen rows cols = do
+  segd <- randomVectorR ( 0,cols`div`3) gen rows
+  let nnz = V.sum segd
+  inds <- randomVectorR ( 0,cols-1) gen nnz
+  vals <- randomVectorR (-1,1) gen nnz
+  return (segd, V.zip inds vals)
+
+
+-- Read elements into unboxed arrays, convert to zero-indexed compress sparse
+-- row format.
+--
+{-# INLINE csr #-}
+csr :: (Int,Int) -> Int -> [(Int,Int,Float)] -> IO (Vector Int, Vector (Int,Float))
+csr (m,_) l elems = do
+  mu <- M.new l
+  let goe  _ []     = return ()
+      goe !n (x:xs) = let (i,j,v) = x in M.unsafeWrite mu n (i-1,j-1,v) >> goe (n+1) xs
+  goe 0 elems
+
+  let cmp (x1,y1,_) (x2,y2,_) | x1 == x2  = compare y1 y2
+                              | otherwise = compare x1 x2
+  V.sortBy cmp mu
+
+  (i,j,v) <- V.unzip3 `fmap` V.unsafeFreeze mu
+  mseg    <- M.new m
+  let gos !n rows | n < m     = let (s,ss) = V.span (==n) rows in M.unsafeWrite mseg n (V.length s) >> gos (n+1) ss
+                  | otherwise = V.unsafeFreeze mseg
+  seg <- gos 0 i
+  return (seg , V.zip j v)
+
+
+-- Lazier versions of things in Control.Monad
+--
+sequence' :: [IO a] -> IO [a]
+sequence' ms = foldr k (return []) ms
+    where k m m' = do { x <- m; xs <- unsafeInterleaveIO m'; return (x:xs) }
+
+mapM' :: (a -> IO b) -> [a] -> IO [b]
+mapM' f as = sequence' (map f as)
+
diff --git a/tests/simple/SMVM/MatrixMarket.hs b/tests/simple/SMVM/MatrixMarket.hs
new file mode 100644
--- /dev/null
+++ b/tests/simple/SMVM/MatrixMarket.hs
@@ -0,0 +1,133 @@
+{-# LANGUAGE OverloadedStrings, GADTs, StandaloneDeriving #-}
+
+module SMVM.MatrixMarket (Matrix(..), readMatrix) where
+
+import Control.Applicative                      hiding (many)
+
+import Data.Complex
+import Data.Attoparsec.Char8                    hiding (parse, Result(..))
+import Data.Attoparsec.Lazy                     (parse, Result(..))
+import Data.ByteString.Lex.Double
+import qualified Data.ByteString.Lazy           as L
+
+
+-- | Specifies the element type.  Pattern matrices do not have any elements,
+-- only indices, and only make sense for coordinate matrices and vectors.
+--
+data Field  = Real | Complex | Integer | Pattern
+    deriving (Eq, Show)
+
+-- | Specifies either sparse or dense storage.  In sparse (\"coordinate\")
+-- storage, elements are given in (i,j,x) triplets for matrices (or (i,x) for
+-- vectors).  Indices are 1-based, so that A(1,1) is the first element of a
+-- matrix, and x(1) is the first element of a vector.
+--
+-- In dense (\"array\") storage, elements are given in column-major order.
+--
+-- In both cases, each element is given on a separate line.
+--
+data Format = Coordinate | Array 
+    deriving (Eq, Show)
+
+-- | Specifies any special structure in the matrix.  For symmetric and hermition
+-- matrices, only the lower-triangular part of the matrix is given. For skew
+-- matrices, only the entries below the diagonal are stored.
+--
+data Structure = General | Symmetric | Hermitian | Skew 
+    deriving (Eq, Show)
+
+
+-- We really want a type parameter to Matrix, but I think that requires some
+-- kind of dynamic typing so that we can determine (a ~ Integral) or (a ~
+-- RealFloat), and so forth, depending on the file being read. This will do for
+-- our purposes...
+--
+-- Format is: (rows,columns) nnz [(row,column,value)]
+--
+data Matrix where
+  PatternMatrix :: (Int,Int) -> Int -> [(Int,Int)]               -> Matrix
+  IntMatrix     :: (Int,Int) -> Int -> [(Int,Int,Int)]           -> Matrix
+  RealMatrix    :: (Int,Int) -> Int -> [(Int,Int,Float)]         -> Matrix
+  ComplexMatrix :: (Int,Int) -> Int -> [(Int,Int,Complex Float)] -> Matrix
+
+deriving instance Show Matrix
+
+
+--------------------------------------------------------------------------------
+-- Combinators
+--------------------------------------------------------------------------------
+
+comment :: Parser ()
+comment = char '%' *> skipWhile (not . eol) *> endOfLine
+  where eol w = w `elem` "\n\r"
+
+floating :: Fractional a => Parser a
+floating = do
+  mv <- readDouble <$> (skipSpace *> takeTill isSpace)  -- readDouble does the fancy stuff
+  case mv of
+       Just (v,_) -> return . realToFrac $ v
+       Nothing    -> fail "floating-point number"
+
+integral :: Integral a => Parser a
+integral = skipSpace *> decimal
+
+format :: Parser Format
+format =  string "coordinate" *> pure Coordinate
+      <|> string "array"      *> pure Array
+      <?> "matrix format"
+
+field :: Parser Field
+field =  string "real"    *> pure Real
+     <|> string "complex" *> pure Complex
+     <|> string "integer" *> pure Integer
+     <|> string "pattern" *> pure Pattern
+     <?> "matrix field"
+
+structure :: Parser Structure
+structure =  string "general"        *> pure General
+         <|> string "symmetric"      *> pure Symmetric
+         <|> string "hermitian"      *> pure Hermitian
+         <|> string "skew-symmetric" *> pure Skew
+         <?> "matrix structure"
+
+header :: Parser (Format,Field,Structure)
+header =  string "%%MatrixMarket matrix"
+       >> (,,) <$> (skipSpace *> format)
+               <*> (skipSpace *> field)
+               <*> (skipSpace *> structure)
+               <*  endOfLine
+               <?> "MatrixMarket header"
+
+extent :: Parser (Int,Int,Int)
+extent = do
+  [m,n,l] <- skipWhile isSpace *> count 3 integral <* endOfLine
+  return (m,n,l)
+
+line :: Parser a -> Parser (Int,Int,a)
+line f = (,,) <$> integral
+              <*> integral
+              <*> f
+              <*  endOfLine
+
+--------------------------------------------------------------------------------
+-- Matrix Market
+--------------------------------------------------------------------------------
+
+matrix :: Parser Matrix
+matrix = do
+  (_,t,_) <- header
+  (m,n,l) <- skipMany comment *> extent
+  case t of
+    Real    -> RealMatrix    (m,n) l `fmap` many (line floating)
+    Complex -> ComplexMatrix (m,n) l `fmap` many (line ((:+) <$> floating <*> floating))
+    Integer -> IntMatrix     (m,n) l `fmap` many (line integral)
+    Pattern -> PatternMatrix (m,n) l `fmap` many ((,) <$> integral <*> integral)
+
+
+readMatrix :: FilePath -> IO Matrix
+readMatrix file = do
+  chunks <- L.readFile file
+  case parse matrix chunks of
+       Fail _ _ msg -> error $ file ++ ": " ++ msg
+       Done _ mtx   -> return mtx
+
diff --git a/tests/simple/SharingRecovery.hs b/tests/simple/SharingRecovery.hs
new file mode 100644
--- /dev/null
+++ b/tests/simple/SharingRecovery.hs
@@ -0,0 +1,177 @@
+{-# LANGUAGE TypeOperators, ScopedTypeVariables #-}
+
+
+--
+-- Some tests to make sure that sharing recovery is working.
+--
+module SharingRecovery where
+
+import Prelude hiding (zip3)
+
+import Data.Array.Accelerate as Acc
+
+
+mkArray :: Int -> Acc (Array DIM1 Int)
+mkArray n = use $ fromList (Z:.1) [n]
+
+muchSharing :: Int -> Acc (Array DIM1 Int)
+muchSharing 0 = (mkArray 0)
+muchSharing n = Acc.map (\_ -> newArr ! (lift (Z:.(0::Int))) +
+                               newArr ! (lift (Z:.(1::Int)))) (mkArray n)
+  where
+    newArr = muchSharing (n-1)
+
+idx :: Int -> Exp DIM1
+idx i = lift (Z:.i)
+
+bfsFail :: Acc (Array DIM1 Int)
+bfsFail = Acc.map (\x -> (map2 ! (idx 1)) +  (map1 ! (idx 2)) + x) arr
+  where
+    map1 :: Acc (Array DIM1 Int)
+    map1 =  Acc.map (\y -> (map2 ! (idx 3)) + y) arr
+    map2 :: Acc (Array DIM1 Int)
+    map2 =  Acc.map (\z -> z + 1) arr
+    arr  :: Acc (Array DIM1 Int)
+    arr =  mkArray 666
+
+twoLetsSameLevel :: Acc (Array DIM1 Int)
+twoLetsSameLevel =
+  let arr1 = mkArray 1
+  in let arr2 = mkArray 2
+     in  Acc.map (\_ -> arr1!(idx 1) + arr1!(idx 2) + arr2!(idx 3) + arr2!(idx 4)) (mkArray 3)
+
+twoLetsSameLevel2 :: Acc (Array DIM1 Int)
+twoLetsSameLevel2 =
+ let arr2 = mkArray 2
+ in let arr1 = mkArray 1
+    in  Acc.map (\_ -> arr1!(idx 1) + arr1!(idx 2) + arr2!(idx 3) + arr2!(idx 4)) (mkArray 3)
+
+--
+-- These two programs test that lets can be introduced not just at the top of a AST
+-- but in intermediate nodes.
+--
+noLetAtTop :: Acc (Array DIM1 Int)
+noLetAtTop = Acc.map (\x -> x + 1) bfsFail
+
+noLetAtTop2 :: Acc (Array DIM1 Int)
+noLetAtTop2 = Acc.map (\x -> x + 2) $ Acc.map (\x -> x + 1) bfsFail
+
+--
+--
+--
+simple :: Acc (Array DIM1 (Int,Int))
+simple = Acc.map (\_ -> a ! (idx 1))  d
+  where
+    c = use $ Acc.fromList (Z :. 3) [1..]
+    d = Acc.map (+1) c
+    a = Acc.zip d c
+
+--------------------------------------------------------------------------------
+--
+-- sortKey is a real program that Ben Lever wrote. It has some pretty interesting
+-- sharing going on.
+--
+sortKey :: (Elt e)
+        => (Exp e -> Exp Int)         -- ^mapping function to produce key array from input array
+        -> Acc (Vector e)
+        -> Acc (Vector e)
+sortKey keyFun arr =  foldl sortOneBit arr (Prelude.map lift ([0..31] :: [Int]))
+  where
+    sortOneBit inArr bitNum = outArr
+      where
+        keys    = Acc.map keyFun inArr
+
+        bits    = Acc.map (\a -> (Acc.testBit a bitNum) ? (1, 0)) keys
+        bitsInv = Acc.map (\b -> (b ==* 0) ? (1, 0)) bits
+
+        (falses, numZeroes) = Acc.scanl' (+) 0 bitsInv
+        trues               = Acc.map (\x -> (Acc.the numZeroes) + (Acc.fst x) - (Acc.snd x)) $
+                               Acc.zip ixs falses
+
+        dstIxs = Acc.map (\x -> let (b, t, f) = unlift x  in (b ==* (constant (0::Int))) ? (f, t)) $
+                   zip3 bits trues falses
+        outArr = scatter dstIxs inArr inArr -- just use input as default array
+                                            --(we're writing over everything anyway)
+    --
+    ixs   = enumeratedArray (shape arr)
+
+-- | Copy elements from source array to destination array according to a map. For
+--   example:
+--
+--    default = [0, 0, 0, 0, 0, 0, 0, 0, 0]
+--    map     = [1, 3, 7, 2, 5, 8]
+--    input   = [1, 9, 6, 4, 4, 2, 5]
+--
+--    output  = [0, 1, 4, 9, 0, 4, 0, 6, 2]
+--
+--   Note if the same index appears in the map more than once, the result is
+--   undefined. The map vector cannot be larger than the input vector.
+--
+scatter :: (Elt e)
+        => Acc (Vector Int)      -- ^map
+        -> Acc (Vector e)        -- ^default
+        -> Acc (Vector e)        -- ^input
+        -> Acc (Vector e)        -- ^output
+scatter mapV defaultV inputV = Acc.permute (const) defaultV pF inputV
+  where
+    pF ix = lift (Z :. (mapV ! ix))
+
+
+-- | Create an array where each element is the value of its corresponding row-major
+--   index.
+--
+--enumeratedArray :: (Shape sh) => Exp sh -> Acc (Array sh Int)
+--enumeratedArray sh = Acc.reshape sh
+--                   $ Acc.generate (index1 $ shapeSize sh) unindex1
+
+enumeratedArray :: Exp DIM1 -> Acc (Array DIM1 Int)
+enumeratedArray sh = Acc.generate sh unindex1
+
+unzip3 :: forall sh. forall e1. forall e2. forall e3. (Shape sh, Elt e1, Elt e2, Elt e3)
+       => Acc (Array sh (e1, e2, e3))
+       -> (Acc (Array sh e1), Acc (Array sh e2), Acc (Array sh e3))
+unzip3 abcs = (as, bs, cs)
+  where
+    (bs, cs)  = Acc.unzip bcs
+    (as, bcs) = Acc.unzip
+              $ Acc.map (\abc -> let (a, b, c) = unlift abc :: (Exp e1, Exp e2, Exp e3)
+                                 in lift (a, lift (b, c))) abcs
+
+testSort :: Acc (Vector Int)
+testSort = sortKey id $ use $ fromList (Z:.10) [9,8,7,6,5,4,3,2,1,0]
+
+----------------------------------------------------------------------
+
+--
+-- map1 has children map3 and map2.
+-- map2 has child map3.
+-- Back when we still used a list for the NodeCounts data structure this mean that
+-- you would be merging [1,3,2] with [2,3] which violated precondition of (+++).
+-- This tests that the new algorithm works just fine on this.
+--
+orderFail :: Acc (Array DIM1 Int)
+orderFail = Acc.map (\_ -> map1 ! (idx 1) + map2 ! (idx 1)) arr
+  where
+    map1 = Acc.map (\_ -> map3 ! (idx 1) + map2 ! (idx 2)) arr
+    map2 = Acc.map (\_ -> map3 ! (idx 3)) arr
+    map3 = Acc.map (+1) arr
+    arr = mkArray 42
+
+----------------------------------------------------------------------
+
+-- Tests array-valued lambdas in conjunction with sharing recovery.
+--
+pipe :: Acc (Vector Int)
+pipe = (acc1 >-> acc2) xs
+  where
+    z :: Acc (Scalar Int)
+    z = unit 0
+
+    xs :: Acc (Vector Int)
+    xs = use $ fromList (Z:.10) [0..]
+
+    acc1 :: Acc (Vector Int) -> Acc (Vector Int)
+    acc1 = Acc.map (\_ -> the z)
+
+    acc2 :: Acc (Vector Int) -> Acc (Vector Int)
+    acc2 arr = let arr2 = use $ fromList (Z:.10) [10..] in Acc.map (\_ -> arr2!constant (Z:.(0::Int))) (Acc.zip arr arr2)
diff --git a/tests/simple/SliceExamples.hs b/tests/simple/SliceExamples.hs
new file mode 100644
--- /dev/null
+++ b/tests/simple/SliceExamples.hs
@@ -0,0 +1,114 @@
+{-# LANGUAGE TypeOperators, ScopedTypeVariables #-}
+module SliceExamples where
+
+import Data.Array.Accelerate as Acc
+import qualified Data.Array.Unboxed as UA
+
+--    y
+--    ^
+--    |  3   4 
+--    |  1   2
+--     -------> x
+--
+arr :: Acc (Array DIM2 Int)
+arr = use $ fromList (Z:.2:.2) [1,2,3,4]
+
+slice1 :: Exp (Z:.Int:.All:.All)
+slice1 = lift $ Z:.(2::Int):.All:.All
+
+slice2 :: Exp (Z:.All:.Int:.All)
+slice2 = lift $ Z:.All:.(2::Int):.All
+
+slice3 :: Exp (Z:.All:.All:.Int)
+slice3 = lift $ Z:.All:.All:.(2::Int)
+
+-- Replicate into z-axis
+-- should produce [1,2,3,4,1,2,3,4]
+test1 :: Acc (Array DIM3 Int)
+test1 = Acc.replicate slice1 arr
+
+-- Replicate into y-axis
+-- should produce [1,2,1,2,3,4,3,4]
+test2 :: Acc (Array DIM3 Int)
+test2 =  Acc.replicate slice2 arr
+
+-- Replicate into x-axis
+-- should produce [1,1,2,2,3,3,4,4]
+test3 :: Acc (Array DIM3 Int)
+test3 =  Acc.replicate slice3 arr
+
+--
+-- repN. Replicates an array into the rightmost dimension of
+-- the result array.
+--
+repN :: forall sh e. (Shape sh, Elt e)
+     => Int 
+     -> Acc (Array sh e)
+     -> Acc (Array (sh:.Int) e)
+repN n a = Acc.replicate (lift (Any:.n :: Any sh:.Int)) a
+
+repExample :: Acc (Array DIM2 Int) -> Acc (Array DIM3 Int)
+repExample = repN 2
+
+repExample' :: Acc (Array DIM2 Int) -> Acc (Array DIM3 Int)
+repExample' = Acc.replicate (lift (Z:.All:.All:.(2::Int)))
+
+slice1' :: Any (Z:.Int:.Int) :. Int
+slice1' = Any:.2
+
+slice2' :: Z:.All:.All:.Int
+slice2' = Z:.All:.All:.2
+
+run1 :: IO (() -> UA.UArray (Int,Int,Int) Int, () -> Acc (Array DIM3 Int))
+run1 = return (\() -> UA.array ((0,0,0),(1,1,1)) [ ((0,0,0), 1)
+                                                 , ((0,0,1), 2)
+                                                 , ((0,1,0), 3)
+                                                 , ((0,1,1), 4)
+                                                 , ((1,0,0), 1)
+                                                 , ((1,0,1), 2)
+                                                 , ((1,1,0), 3)
+                                                 , ((1,1,1), 4) ]
+              ,\() -> test1)
+
+run2 :: IO (() -> UA.UArray (Int,Int,Int) Int, () -> Acc (Array DIM3 Int))
+run2 = return (\() -> UA.array ((0,0,0),(1,1,1)) [ ((0,0,0), 1)
+                                                 , ((0,0,1), 2)
+                                                 , ((0,1,0), 1)
+                                                 , ((0,1,1), 2)
+                                                 , ((1,0,0), 3)
+                                                 , ((1,0,1), 4)
+                                                 , ((1,1,0), 3)
+                                                 , ((1,1,1), 4) ]
+              ,\() -> test2)
+
+run3 :: IO (() -> UA.UArray (Int,Int,Int) Int, () -> Acc (Array DIM3 Int))
+run3 = return (\() -> UA.array ((0,0,0),(1,1,1)) [ ((0,0,0), 1)
+                                                 , ((0,0,1), 1)
+                                                 , ((0,1,0), 2)
+                                                 , ((0,1,1), 2)
+                                                 , ((1,0,0), 3)
+                                                 , ((1,0,1), 3)
+                                                 , ((1,1,0), 4)
+                                                 , ((1,1,1), 4) ]
+              ,\() -> test3)
+run4 :: IO (() -> UA.UArray (Int,Int,Int) Int, () -> Acc (Array DIM3 Int))
+run4 = return (\() -> UA.array ((0,0,0),(1,1,1)) [ ((0,0,0), 1)
+                                                 , ((0,0,1), 1)
+                                                 , ((0,1,0), 2)
+                                                 , ((0,1,1), 2)
+                                                 , ((1,0,0), 3)
+                                                 , ((1,0,1), 3)
+                                                 , ((1,1,0), 4)
+                                                 , ((1,1,1), 4) ]
+              ,\() -> repExample arr)
+
+run5 :: IO (() -> UA.UArray (Int,Int,Int) Int, () -> Acc (Array DIM3 Int))
+run5 = return (\() -> UA.array ((0,0,0),(1,1,1)) [ ((0,0,0), 1)
+                                                 , ((0,0,1), 1)
+                                                 , ((0,1,0), 2)
+                                                 , ((0,1,1), 2)
+                                                 , ((1,0,0), 3)
+                                                 , ((1,0,1), 3)
+                                                 , ((1,1,0), 4)
+                                                 , ((1,1,1), 4) ]
+              ,\() -> repExample' arr)