diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,30 +1,23 @@ -Copyright (c)2011, Trevor L. McDonell - -All rights reserved. +Copyright (c) [2007..2012] The Accelerate Team. 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. + * 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 names of the contributors nor of their affiliations 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. +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ''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 COPYRIGHT HOLDERS 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/accelerate-examples.cabal b/accelerate-examples.cabal --- a/accelerate-examples.cabal +++ b/accelerate-examples.cabal @@ -1,69 +1,118 @@ -Name: accelerate-examples -Version: 0.2.0.1 -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 -Category: Compilers/Interpreters, -Build-type: Simple -Cabal-version: >=1.6 -Homepage: http://www.cse.unsw.edu.au/~chak/project/accelerate/ - -Category: Compilers/Interpreters, Concurrency, Data, Parallelism -Stability: Experimental +Name: accelerate-examples +Version: 0.12.0.0 +License: BSD3 +License-file: LICENSE +Author: The Accelerate Team +Maintainer: Trevor L. McDonell +Category: Compilers/Interpreters, Concurrency, Data, Parallelism +Stability: Experimental +Build-type: Simple +Cabal-version: >=1.6 +Homepage: http://www.cse.unsw.edu.au/~chak/project/accelerate/ +Synopsis: Examples using the Accelerate library +Description: This package defines a number of example programs for verification and + performance testing of the Accelerate language and associated backend + implementations. By default the package attempts to build with all available + backends, but this might not be desirable or even possible depending on your + available hardware. To disable a specific component, install with the extra + cabal flag(s): + . + . > cabal install accelerate-examples -f-opt + . + Where the available options are: + . + * cuda: A parallel CUDA backend for NVIDIA GPUs + . + * opencl: A parallel OpenCL backend + . + * io: Extra tests for reading and writing arrays in various formats + . Flag cuda Description: Enable the CUDA parallel backend for NVIDIA GPUs Default: True +Flag opencl + Description: Enable the OpenCL parallel backend + Default: True + Flag io Description: Provide access to the block copy I/O functionality - Default: False + Default: True +-- Run some basic randomised tests of the primitive library operations, and unit +-- tests of known (hopefully prior) failures. +-- +Executable accelerate-quickcheck + Main-is: Main.hs + hs-source-dirs: examples/quickcheck/ + + other-modules: Config + Arbitrary.Array + Arbitrary.Shape + Test.Base + Test.IndexSpace + Test.Mapping + Test.PrefixSum + Test.Reduction + + ghc-options: -Wall -O2 + -threaded + -fpedantic-bottoms + -fno-warn-orphans + -fno-full-laziness + -fno-excess-precision + if impl(ghc >= 7.0) + ghc-options: -rtsopts + + if flag(cuda) + CPP-options: -DACCELERATE_CUDA_BACKEND + build-depends: accelerate-cuda >= 0.12, + cuda >= 0.4 + + if flag(opencl) + CPP-options: -DACCELERATE_OPENCL_BACKEND + build-depends: accelerate-opencl == 0.1.* + + build-depends: accelerate >= 0.12, + QuickCheck >= 2.0, + test-framework >= 0.5, + test-framework-quickcheck2 >= 0.2, + random + + +-- The main examples program including verification and timing tests for several +-- simple accelerate programs +-- 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 + hs-source-dirs: examples/tests + examples/tests/io + examples/tests/image-processing + examples/tests/primitives + examples/tests/simple - c-sources: tests/io/fill_with_values.cpp + other-modules: Benchmark, Config, Random, Test, Util, Validate, + + -- tests: + -- image-processing + Canny, IntegralImage, PGM, + + -- io + BlockCopy, VectorCopy, + + -- primitives + Backpermute, Fold, FoldSeg, Gather, Map, Permute, ScanSeg, + Scatter, Stencil, Stencil2, Vector, Zip, ZipWith, + + -- simple + BlackScholes, DotP, Filter, Radix, SASUM, SAXPY, + SharingRecovery, SliceExamples, SMVM, SMVM.Matrix, SMVM.MatrixMarket + + + c-sources: examples/tests/io/fill_with_values.cpp extra-libraries: stdc++ ghc-options: -Wall -O2 @@ -72,25 +121,90 @@ if flag(cuda) CPP-options: -DACCELERATE_CUDA_BACKEND + build-depends: accelerate-cuda >= 0.12 + if flag(opencl) + CPP-options: -DACCELERATE_OPENCL_BACKEND + build-depends: accelerate-opencl == 0.1.* + if flag(io) CPP-options: -DACCELERATE_IO + build-depends: accelerate-io >= 0.12 - build-depends: accelerate == 0.10.*, - array >= 0.3 && < 0.5, - attoparsec == 0.8.*, + build-depends: accelerate >= 0.12, + array >= 0.3, + attoparsec >= 0.10, 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.* + bytestring >= 0.9, + bytestring-lexing >= 0.2, + cmdargs >= 0.6, + criterion >= 0.5, + deepseq >= 1.1, + directory >= 1.0, + filepath >= 1.0, + hashtables >= 1.0.1, + mtl >= 1.1, + mwc-random >= 0.8, + pgm >= 0.1, + pretty >= 1.0, + vector >= 0.7, + vector-algorithms >= 0.4 + +-- A quasicrystal demo as the sum of waves in a plane +-- +Executable accelerate-crystal + hs-source-dirs: examples/crystal + Main-is: Main.hs + other-modules: Config + ghc-options: -O2 -Wall + if impl(ghc >= 7.0) + ghc-options: -rtsopts + + if flag(cuda) + CPP-options: -DACCELERATE_CUDA_BACKEND + build-depends: accelerate-cuda >= 0.12 + + if flag(opencl) + CPP-options: -DACCELERATE_OPENCL_BACKEND + build-depends: accelerate-opencl == 0.1.* + + build-depends: accelerate >= 0.12, + base == 4.*, + criterion >= 0.5, + fclabels >= 1.0, + gloss >= 1.5 + + +-- A stable fluid simulation +-- +Executable accelerate-fluid + Main-is: Main.hs + hs-source-dirs: examples/fluid/src + other-modules: Config + Event + Fluid + Type + World + + ghc-options: -O2 -Wall + + if flag(cuda) + cpp-options: -DACCELERATE_CUDA_BACKEND + build-depends: accelerate-cuda >= 0.12 + + Build-depends: accelerate >= 0.12, + accelerate-io >= 0.12, + base == 4.*, + bmp >= 1.2, + criterion >= 0.5, + fclabels >= 1.0, + gloss >= 1.7 + + +source-repository head + type: git + location: https://github.com/AccelerateHS/accelerate-examples + +-- vim: nospell +-- diff --git a/examples/crystal/Config.hs b/examples/crystal/Config.hs new file mode 100644 --- /dev/null +++ b/examples/crystal/Config.hs @@ -0,0 +1,89 @@ +{-# LANGUAGE CPP, TemplateHaskell, PatternGuards #-} + +module Config ( + + Options, optBackend, optSize, optZoom, optScale, optDegree, optBench, + processArgs, run + +) where + +import Data.Label +import System.Exit +import System.Console.GetOpt +import Data.Array.Accelerate ( Arrays, Acc ) +import qualified Data.Array.Accelerate.Interpreter as Interp +#ifdef ACCELERATE_CUDA_BACKEND +import qualified Data.Array.Accelerate.CUDA as CUDA +#endif + +data Backend = Interpreter +#ifdef ACCELERATE_CUDA_BACKEND + | CUDA +#endif + deriving (Bounded, Show) + +data Options = Options + { + _optBackend :: Backend + , _optSize :: Int + , _optZoom :: Int + , _optScale :: Float + , _optDegree :: Int + , _optBench :: Bool + , _optHelp :: Bool + } + deriving Show + +$(mkLabels [''Options]) + +defaultOptions :: Options +defaultOptions = Options + { _optBackend = maxBound + , _optSize = 200 + , _optZoom = 3 + , _optScale = 30 + , _optDegree = 5 + , _optBench = False + , _optHelp = False + } + + +run :: (Arrays a, Arrays b) => Options -> (Acc a -> Acc b) -> a -> b +run opts f = case _optBackend opts of + Interpreter -> head . Interp.stream f . return +#ifdef ACCELERATE_CUDA_BACKEND + CUDA -> CUDA.run1 f +#endif + + +options :: [OptDescr (Options -> Options)] +options = + [ Option [] ["interpreter"] (NoArg (set optBackend Interpreter)) "reference implementation (sequential)" +#ifdef ACCELERATE_CUDA_BACKEND + , Option [] ["cuda"] (NoArg (set optBackend CUDA)) "implementation for NVIDIA GPUs (parallel)" +#endif + , Option [] ["size"] (ReqArg (set optSize . read) "INT") "visualisation size (200)" + , Option [] ["zoom"] (ReqArg (set optZoom . read) "INT") "pixel replication factor (3)" + , Option [] ["scale"] (ReqArg (set optScale . read) "FLOAT") "feature size of visualisation (30)" + , Option [] ["degree"] (ReqArg (set optDegree . read) "INT") "number of waves to sum for each point (5)" + , Option [] ["benchmark"] (NoArg (set optBench True)) "benchmark instead of displaying animation (False)" + , Option "h?" ["help"] (NoArg (set optHelp True)) "show help message" + ] + + +processArgs :: [String] -> IO (Options, [String]) +processArgs argv = + case getOpt' Permute options argv of + (o,_,n,[]) -> case foldl (flip id) defaultOptions o of + opts | False <- get optHelp opts -> return (opts, n) + opts | True <- get optBench opts -> return (opts, "--help":n) + _ -> putStrLn (helpMsg []) >> exitSuccess + (_,_,_,err) -> error (helpMsg err) + where + helpMsg err = concat err ++ usageInfo header options + header = unlines + [ "accelerate-crystal (c) [2011..2012] The Accelerate Team" + , "" + , "Usage: accelerate-crystal [OPTIONS]" + ] + diff --git a/examples/crystal/Main.hs b/examples/crystal/Main.hs new file mode 100644 --- /dev/null +++ b/examples/crystal/Main.hs @@ -0,0 +1,196 @@ +{-# LANGUAGE CPP #-} +-- Quasicrystals demo. +-- +-- Based on code from: +-- http://hackage.haskell.org/package/repa-examples +-- http://mainisusuallyafunction.blogspot.com/2011/10/quasicrystals-as-sums-of-waves-in-plane.html +-- + +module Main where + +import Config + +import Data.Word +import Data.Label +import Foreign.Ptr +import Control.Monad +import Control.Exception +import Criterion.Main ( defaultMain, bench, whnf ) +import Foreign.ForeignPtr +import System.Environment +import System.IO.Unsafe +import Data.Array.Accelerate ( Array, Scalar, Exp, Acc, DIM2, Z(..), (:.)(..) ) +import qualified Data.Array.Accelerate as A +import qualified Graphics.Gloss as G + +import Data.Array.Accelerate.Array.Data ( ptrsOfArrayData ) +import Data.Array.Accelerate.Array.Sugar ( Array(..) ) + + +-- Types ---------------------------------------------------------------------- +-- | Real value +type R = Float + +-- | Point on the 2D plane. +type R2 = (R, R) + +-- | Angle in radians. +type Angle = R + +-- | Angle offset used for animation. +type Phi = Float + +-- | Number of waves to sum for each pixel. +type Degree = Int + +-- | Feature size of visualisation. +type Scale = Float + +-- | Size of image to render. +type Size = Int + +-- | How many times to duplicate each pixel / image zoom. +type Zoom = Int + +-- | Type of the generated image data +type RGBA = Word32 +type Bitmap = Array DIM2 RGBA + +-- | Action to render a frame +type Render = Scalar Phi -> Bitmap + + +-- Point ---------------------------------------------------------------------- +-- | Compute a single point of the visualisation. +quasicrystal :: Size -> Scale -> Degree -> Acc (Scalar Phi) -> Exp DIM2 -> Exp R +quasicrystal size scale degree phi p + = waves degree phi $ point size scale p + + +-- | Sum up all the waves at a particular point. +waves :: Degree -> Acc (Scalar Phi) -> Exp R2 -> Exp R +waves degree phi x = wrap $ waver degree 0 + where + waver :: Int -> Exp Float -> Exp Float + waver n acc + | n == 0 = acc + | otherwise = waver (n - 1) (acc + wave (A.constant (fromIntegral n) * phi') x) + + phi' + = A.the phi + + wrap n + = let n_ = A.truncate n :: Exp Int + n' = n - A.fromIntegral n_ + in + (n_ `rem` 2 A./=* 0) A.? (1-n', n') + + +-- | Generate the value for a single wave. +wave :: Exp Angle -> Exp R2 -> Exp R +wave th pt = (cos (cth*x + sth*y) + 1) / 2 + where + (x,y) = A.unlift pt + cth = cos th + sth = sin th + +-- | Convert an image point to a point on our wave plane. +point :: Size -> Scale -> Exp DIM2 -> Exp R2 +point size scale ix = A.lift (adj x, adj y) + where + (Z:.x:.y) = A.unlift ix + denom = A.constant (fromIntegral size - 1) + adj n = A.constant scale * ((2 * A.fromIntegral n / denom) - 1) + + +-- Computation ---------------------------------------------------------------- +-- | Compute a single frame +makeImage :: Size -> Scale -> Degree -> Acc (Scalar Phi) -> Acc Bitmap +makeImage size scale degree phi = arrPixels + where + -- Compute values for the wave density at each point in the range [0..1] + arrVals :: Acc (Array DIM2 Float) + arrVals = A.generate + (A.constant $ Z :. size :. size) + (quasicrystal size scale degree phi) + + -- Convert the [0..1] values of wave density to an RGBA flat image + arrPixels :: Acc Bitmap + arrPixels = A.map rampColour arrVals + + +-- | Colour ramp from red to white, convert into RGBA +rampColour :: Exp Float -> Exp RGBA +rampColour v = ra + g + b + where + u = 0 `A.max` v `A.min` 1 + ra = 0xFF0000FF + g = A.truncate ((0.4 + (u * 0.6)) * 0xFF) * 0x10000 + b = A.truncate (u * 0xFF) * 0x100 + + +-- Rendering ------------------------------------------------------------------ +-- | Compute a single frame of the animation as a Gloss picture. +--frame :: Size -> Scale -> Zoom -> Degree -> Float -> G.Picture + +frame :: Render -> Size -> Zoom -> Float -> G.Picture +frame render size zoom time = G.scale zoom' zoom' pic + where + -- Scale the time to be the phi value of the animation. The action seems to + -- slow down at increasing phi values, so we increase phi faster as time + -- moves on. + x = 1 + (time ** 1.5) * 0.005 + + -- lift to a singleton array, else we would generate new code with the + -- constant embedded at every frame + phi = A.fromList Z [pi / x] + + -- Compute the image + arrPixels = render phi + + -- Wrap the array data in a Foreign pointer and turn into a Gloss picture + {-# NOINLINE rawData #-} + rawData = let (Array _ adata) = arrPixels + ((),ptr) = ptrsOfArrayData adata + in + unsafePerformIO $ newForeignPtr_ (castPtr ptr) + + pic = G.bitmapOfForeignPtr + size size -- raw image size + rawData -- the image data + False -- don't cache this in texture memory + + -- Zoom the image so we get a bigger window. + zoom' = fromIntegral zoom + + +-- Main ----------------------------------------------------------------------- +main :: IO () +main + = do (config, nops) <- processArgs =<< getArgs + let size = get optSize config + zoom = get optZoom config + scale = get optScale config + degree = get optDegree config + render = run config $ makeImage size scale degree + + void . evaluate $ render (A.fromList Z [0]) + + if get optBench config + then withArgs nops $ defaultMain + [ bench "crystal" $ whnf render (A.fromList Z [1.0]) ] + +#if MIN_VERSION_gloss(1,6,0) + else G.animate + (G.InWindow "Quasicrystals" (size * zoom, size * zoom) (10, 10)) + (G.black) + (frame render size zoom) +#else + else G.animateInWindow + "Quasicrystals" + (size * zoom, size * zoom) + (10, 10) + (G.black) + (frame render size zoom) +#endif + diff --git a/examples/fluid/src/Config.hs b/examples/fluid/src/Config.hs new file mode 100644 --- /dev/null +++ b/examples/fluid/src/Config.hs @@ -0,0 +1,164 @@ +{-# LANGUAGE CPP, TemplateHaskell, PatternGuards #-} +-- +-- Configuration parameters +-- + +module Config ( + + Options, + viscosity, diffusion, densityBMP, velocityBMP, simulationWidth, + simulationHeight, displayScale, displayFramerate, optBench, + + processArgs, run, run1 + +) where + +import Data.Label +import System.Exit +import System.Console.GetOpt +import Data.Array.Accelerate ( Acc, Arrays ) +import qualified Data.Array.Accelerate.Interpreter as I +#ifdef ACCELERATE_CUDA_BACKEND +import qualified Data.Array.Accelerate.CUDA as CUDA +#endif + +data Backend + = Interpreter +#ifdef ACCELERATE_CUDA_BACKEND + | CUDA +#endif + deriving (Show, Bounded) + + +data Options = Options + { + -- simulation + _viscosity :: Float + , _diffusion :: Float + + -- visualisation + , _densityBMP :: Maybe FilePath + , _velocityBMP :: Maybe FilePath + , _simulationWidth :: Int + , _simulationHeight :: Int + , _displayScale :: Int + , _displayFramerate :: Int + + -- misc + , _optBackend :: Backend + , _optBench :: Bool + , _optHelp :: Bool + } + deriving Show + +$(mkLabels [''Options]) + +defaultOptions :: Options +defaultOptions = Options + { _viscosity = 0 + , _diffusion = 0 + + , _densityBMP = Nothing + , _velocityBMP = Nothing + , _simulationWidth = 250 + , _simulationHeight = 250 + , _displayScale = 2 + , _displayFramerate = 20 + + , _optBackend = maxBound + , _optBench = False + , _optHelp = False + } + +-- Execute an Accelerate expression using the selected backend +-- +run :: Arrays a => Options -> Acc a -> a +run opts = case _optBackend opts of + Interpreter -> I.run +#ifdef ACCELERATE_CUDA_BACKEND + CUDA -> CUDA.run +#endif + +run1 :: (Arrays a, Arrays b) => Options -> (Acc a -> Acc b) -> a -> b +run1 opts f = case _optBackend opts of + Interpreter -> head . I.stream f . return +#ifdef ACCELERATE_CUDA_BACKEND + CUDA -> CUDA.run1 f +#endif + + +-- Process command line arguments +-- +options :: [OptDescr (Options -> Options)] +options = + [ Option [] ["interpreter"] (NoArg (set optBackend Interpreter)) + "reference implementation (sequential)" + +#ifdef ACCELERATE_CUDA_BACKEND + , Option [] ["cuda"] (NoArg (set optBackend CUDA)) + "implementation for NVIDIA GPUs (parallel)" + +#endif + , Option [] ["viscosity"] (ReqArg (set viscosity . read) "FLOAT") + $ "viscosity for velocity dampening " ++ def viscosity + + , Option [] ["diffusion"] (ReqArg (set diffusion . read) "FLOAT") + $ "diffusion rate for mass dispersion " ++ def diffusion + + , Option [] ["width"] (ReqArg (set simulationWidth . read) "INT") + $ "grid width for simulation " ++ def simulationWidth + + , Option [] ["height"] (ReqArg (set simulationHeight . read) "INT") + $ "grid height for simulation " ++ def simulationHeight + + , Option [] ["scale"] (ReqArg (set displayScale . read) "INT") + $ "feature size of visualisation " ++ def displayScale + + , Option [] ["framerate"] (ReqArg (set displayFramerate . read) "INT") + $ "frame rate for visualisation " ++ def displayFramerate + + , Option [] ["bmp-density"] (ReqArg (set densityBMP . Just) "FILE.bmp") + "file for initial fluid density" + + , Option [] ["bmp-velocity"] (ReqArg (set velocityBMP . Just) "FILE.bmp") + "file for initial fluid velocity" + + -- + , Option [] ["benchmark"] (NoArg (set optBench True)) + "benchmark instead of displaying animation (False)" + + , Option "h?" ["help"] (NoArg (set optHelp True)) + "show help message" + ] + where + parens s = "(" ++ s ++ ")" + def f = parens (show (get f defaultOptions)) + + +processArgs :: [String] -> IO (Options, [String]) +processArgs argv = + case getOpt' Permute options argv of + (o,_,n,[]) -> case foldl (flip id) defaultOptions o of + opts | False <- get optHelp opts -> return (opts, n) + opts | True <- get optBench opts -> return (opts, "--help":n) + _ -> putStrLn (helpMsg []) >> exitSuccess + (_,_,_,err) -> error (helpMsg err) + where + helpMsg err = concat err ++ usageInfo header options ++ footer + header = unlines + [ "accelerate-fluid (c) 2011 Trevor L. McDonell" + , "" + , "Usage: accelerate-fluid [OPTIONS]" + ] + + footer = unlines + [ "" + , "Runtime usage:" + , "" + , " click add density sources to the image" + , " shift-click add velocity sources" + , " r reset the image" + , " d toggle display of density field" + , " v toggle display of velocity field lines" + ] + diff --git a/examples/fluid/src/Event.hs b/examples/fluid/src/Event.hs new file mode 100644 --- /dev/null +++ b/examples/fluid/src/Event.hs @@ -0,0 +1,81 @@ +{-# LANGUAGE PatternGuards #-} +-- +-- Event handling +-- + +module Event where + +import Config +import World +import Data.Label +import Graphics.Gloss.Interface.Pure.Game +import Data.Array.Accelerate ( Z(..), (:.)(..) ) + + +-- Event locations are returned as window coordinates, where the origin is in +-- the centre of the window and increases to the right and up. If the simulation +-- size is (100,100) with scale factor of 4, then the event coordinates are +-- returned in the range [-200,200]. +-- +react :: Options -> Event -> World -> IO World +react opt event world = + case event of + EventKey (Char c) s m _ -> keyboard c s m + EventKey (MouseButton LeftButton) s m uv -> return $ mouse m s (coord uv) + EventMotion uv -> return $ motion (coord uv) + _ -> return $ world + where + -- Inject a new density source when the left button is clicked. + -- + -- If the shift key is held, remember the location and add a new velocity + -- source between the old and new points as the mouse moves. + -- + mouse key button xy + | Up <- shift key, Down <- button + = world { currentSource = Density xy + , densitySource = addDensity xy } + + | Down <- shift key, Down <- button + = world { currentSource = Velocity xy } + + | Down <- shift key, Up <- button + = case currentSource world of + Velocity x0y0 -> world { currentSource = None + , velocitySource = addVelocity x0y0 xy } + _ -> world { currentSource = None} + + | otherwise + = world { currentSource = None } + + -- Handle key presses + -- + keyboard 'r' Down _ = initialise opt + keyboard 'd' Down _ = return $ world { displayDensity = not (displayDensity world) } + keyboard 'v' Down _ = return $ world { displayVelocity = not (displayVelocity world) } + keyboard _ _ _ = return world + + -- As the mouse moves, keep inserting density sources, or adding source + -- velocities + -- + motion xy = + case currentSource world of + Density _ -> world { currentSource = Density xy + , densitySource = addDensity xy } + Velocity x0y0 -> world { currentSource = Velocity xy + , velocitySource = addVelocity x0y0 xy } + _ -> world + -- + addDensity (x,y) = (Z:.y:.x, 1) : densitySource world + + addVelocity (x0,y0) (x1,y1) = let u = fromIntegral (x1-x0) + v = fromIntegral (y1-y0) + in (Z:.y0:.x0, (u * width, v * height)) : velocitySource world + -- + zoom = fromIntegral $ get displayScale opt + width = fromIntegral $ get simulationWidth opt + height = fromIntegral $ get simulationHeight opt + scaleX = width / (width * zoom + 1) + scaleY = height / (height * zoom + 1) + coord (u,v) = ( truncate $ u * scaleX + width /2 + , truncate $ v * scaleY + height/2) + diff --git a/examples/fluid/src/Fluid.hs b/examples/fluid/src/Fluid.hs new file mode 100644 --- /dev/null +++ b/examples/fluid/src/Fluid.hs @@ -0,0 +1,169 @@ +-- +-- Fluid simulation +-- + +module Fluid ( + + Simulation, fluid + +) where + +import Type + +import Data.Array.Accelerate ( Z(..), (:.)(..), Exp, Acc, Scalar, Vector, (?|), (==*) ) +import qualified Data.Array.Accelerate as A + + +type Simulation + = Acc ( Scalar Timestep -- time to evolve the simulation + , DensitySource -- locations to add density sources + , VelocitySource -- locations to add velocity sources + , DensityField -- the current density field + , VelocityField ) -- the current velocity field + -> Acc ( DensityField, VelocityField ) + + +-- A fluid simulation +-- +fluid :: Viscosity -> Diffusion -> Simulation +fluid dp dn inputs = + let (dt, ds, vs, df, vf) = A.unlift inputs + df' = density dn dt ds vf df + vf' = velocity dp dt vs vf + in + A.lift (df', vf') + +-- The velocity over a timestep evolves due to three causes: +-- 1. the addition of forces +-- 2. viscous diffusion +-- 3. self-advection +-- +velocity + :: Viscosity + -> Acc (Scalar Timestep) + -> Acc VelocitySource + -> Acc VelocityField + -> Acc VelocityField +velocity dp dt vs vf0 + = project + . advect dt vf0 + . project + . diffuse dp dt + $ inject vs vf0 + + +-- Ensure the velocity field conserves mass +-- +project :: Acc VelocityField -> Acc VelocityField +project vf = A.stencil2 poisson A.Mirror vf A.Mirror p + where + steps = 20 + grad = A.stencil divF A.Mirror vf + p1 = A.stencil2 pF (A.Constant 0) grad A.Mirror + p = foldl1 (.) (replicate steps p1) grad + + poisson :: A.Stencil3x3 Velocity -> A.Stencil3x3 Float -> Exp Velocity + poisson (_,(_,uv,_),_) ((_,t,_), (l,_,r), (_,b,_)) = uv .-. 0.5 .*. A.lift (r-l, t-b) + + divF :: A.Stencil3x3 Velocity -> Exp Float + divF ((_,t,_), (l,_,r), (_,b,_)) = -0.5 * (A.fst r - A.fst l + A.snd t - A.snd b) + + pF :: A.Stencil3x3 Float -> A.Stencil3x3 Float -> Exp Float + pF (_,(_,x,_),_) ((_,t,_), (l,_,r), (_,b,_)) = 0.25 * (x + l + t + r + b) + + + +-- The density over a timestep evolves due to three causes: +-- 1. the addition of source particles +-- 2. self-diffusion +-- 3. motion through the velocity field +-- +density + :: Diffusion + -> Acc (Scalar Timestep) + -> Acc DensitySource + -> Acc VelocityField + -> Acc DensityField + -> Acc DensityField +density dn dt ds vf + = advect dt vf + . diffuse dn dt + . inject ds + + +-- Inject sources into the field +-- +-- TLM: sources should be a vector of (index, value) pairs, but no fusion means +-- that extracting the components for permute (via unzip) is extra work. +-- +inject + :: FieldElt e + => Acc (Vector Index, Vector e) + -> Acc (Field e) + -> Acc (Field e) +inject source field = + let (is, ps) = A.unlift source + in A.size ps ==* 0 + ?| ( field, A.permute (.+.) field (is A.!) ps ) + + +diffuse + :: FieldElt e + => Diffusion + -> Acc (Scalar Timestep) + -> Acc (Field e) + -> Acc (Field e) +diffuse dn dt df0 = + a ==* 0 + ?| ( df0 , foldl1 (.) (replicate steps diffuse1) df0 ) + where + steps = 20 + a = A.the dt * A.constant dn * (A.fromIntegral (A.size df0)) + c = 1 + 4*a + + diffuse1 df = A.stencil2 relax (A.Constant zero) df0 A.Mirror df + + relax :: FieldElt e => A.Stencil3x3 e -> A.Stencil3x3 e -> Exp e + relax (_,(_,x0,_),_) ((_,t,_), (l,_,r), (_,b,_)) = (x0 .+. a .*. (l.+.t.+.r.+.b)) ./. c + + +advect + :: FieldElt e + => Acc (Scalar Timestep) + -> Acc VelocityField + -> Acc (Field e) + -> Acc (Field e) +advect dt' vf df = A.generate sh backtrace + where + dt = A.the dt' + sh = A.shape vf + Z :. h :. w = A.unlift sh + + backtrace ix = s0.*.(t0.*.d00 .+. t1.*.d10) .+. s1.*.(t0.*.d01 .+. t1.*.d11) + where + Z:.j:.i = A.unlift ix + (u, v) = A.unlift (vf A.! ix) + + -- backtrack densities based on velocity field + clamp z = A.max 0.5 . A.min (A.fromIntegral z - 1.5) + x = w `clamp` (A.fromIntegral i - dt * u) + y = h `clamp` (A.fromIntegral j - dt * v) + + -- discrete locations surrounding point + i0 = A.truncate x + j0 = A.truncate y + i1 = i0 + 1 + j1 = j0 + 1 + + -- weighting based on location between the discrete points + s1 = x - A.fromIntegral i0 + t1 = y - A.fromIntegral j0 + s0 = 1 - s1 + t0 = 1 - t1 + + -- read the density values surrounding the calculated advection point + d00 = df A.! A.lift (Z :. j0 :. i0) + d10 = df A.! A.lift (Z :. j1 :. i0) + d01 = df A.! A.lift (Z :. j0 :. i1) + d11 = df A.! A.lift (Z :. j1 :. i1) + diff --git a/examples/fluid/src/Main.hs b/examples/fluid/src/Main.hs new file mode 100644 --- /dev/null +++ b/examples/fluid/src/Main.hs @@ -0,0 +1,72 @@ +-- +-- A stable fluid simulation +-- +-- Jos Stam, "Real-time Fluid Dynamics for Games" +-- + +module Main where + +import Config +import World +import Fluid +import Event +import Data.Label +import Criterion.Main +import System.Environment +import Graphics.Gloss.Interface.IO.Game + +import Prelude as P +import Data.Array.Accelerate as A + + +main :: IO () +main = do + (opt,noms) <- processArgs =<< getArgs + let -- configuration parameters + -- + width = get simulationWidth opt * get displayScale opt + height = get simulationHeight opt * get displayScale opt + fps = get displayFramerate opt + dp = get viscosity opt + dn = get diffusion opt + + -- Prepare user-input density and velocity sources + -- + sources s = let (ix, ss) = P.unzip s + sh = Z :. length ix + in ( A.fromList sh ix, A.fromList sh ss ) + + -- Prepare to execute the next step of the simulation. + -- + -- Critically, we use the run1 execution form to ensure we bypass all + -- front-end conversion phases. + -- + simulate timestep world = + let step = run1 opt (fluid dp dn) + dt = A.fromList Z [timestep] + ds = sources (densitySource world) + vs = sources (velocitySource world) + (df', vf') = step ( dt, ds, vs, densityField world, velocityField world ) + in + return $ world { densityField = df', velocityField = vf' + , densitySource = [], velocitySource = [] } + + -- warming up... + initialWorld <- initialise opt + _ <- simulate 0.1 initialWorld + + if get optBench opt + -- benchmark + then withArgs noms $ defaultMain + [ bench "fluid" $ whnfIO (simulate 1.0 initialWorld) ] + + -- simulate + else playIO + (InWindow "accelerate-fluid" (width, height) (10, 20)) + black -- background colour + fps -- display framerate + initialWorld -- initial state of the simulation + (render opt) -- render world state into a picture + (react opt) -- handle user events + simulate -- one step of the simulation + diff --git a/examples/fluid/src/Type.hs b/examples/fluid/src/Type.hs new file mode 100644 --- /dev/null +++ b/examples/fluid/src/Type.hs @@ -0,0 +1,67 @@ +{-# LANGUAGE FlexibleInstances #-} +-- +-- Types used throughout the simulation +-- + +module Type ( + + Timestep, Viscosity, Diffusion, Index, Density, Velocity, + Field, FieldElt(..), DensityField, VelocityField, DensitySource, VelocitySource, + + RGBA, Image + +) where + +import Data.Word +import Data.Array.Accelerate + +type Timestep = Float +type Viscosity = Float +type Diffusion = Float +type Index = DIM2 +type Density = Float +type Velocity = (Float, Float) + +type Field a = Array DIM2 a +type DensityField = Field Density +type VelocityField = Field Velocity + +type Source a = (Vector Index, Vector a) +type DensitySource = Source Density +type VelocitySource = Source Velocity + +type RGBA = Word32 +type Image a = Array DIM2 a + + +infixl 6 .+. +infixl 6 .-. +infixl 7 .*. +infixl 7 ./. + +class Elt e => FieldElt e where + zero :: e + (.+.) :: Exp e -> Exp e -> Exp e + (.-.) :: Exp e -> Exp e -> Exp e + (.*.) :: Exp Float -> Exp e -> Exp e + (./.) :: Exp e -> Exp Float -> Exp e + +instance FieldElt Density where + zero = 0 + (.+.) = (+) + (.-.) = (-) + (.*.) = (*) + (./.) = (/) + +instance FieldElt Velocity where + zero = (0, 0) + (.+.) = app2 (+) + (.-.) = app2 (-) + c .*. xy = let (x,y) = unlift xy in lift (c*x, c*y) + xy ./. c = let (x,y) = unlift xy in lift (x/c, y/c) + +app2 :: Elt e => (Exp e -> Exp e -> Exp e) -> Exp (e,e) -> Exp (e,e) -> Exp (e,e) +app2 f xu yv = let (x,u) = unlift xu + (y,v) = unlift yv + in lift (f x y, f u v) + diff --git a/examples/fluid/src/World.hs b/examples/fluid/src/World.hs new file mode 100644 --- /dev/null +++ b/examples/fluid/src/World.hs @@ -0,0 +1,198 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE PatternGuards #-} +-- +-- Displaying the world state +-- + +module World (World(..), Source(..), initialise, render) where + +import Type +import Config + +import Codec.BMP +import Data.Bits +import Data.Int +import Data.Word +import Data.Label +import Control.Monad +import Graphics.Gloss.Interface.IO.Game +import Foreign.Ptr +import Foreign.ForeignPtr +import Foreign.Storable +import Foreign.Marshal.Alloc +import Data.Array.Accelerate ( Z(..), (:.)(..), Exp ) +import Data.Array.Accelerate.Array.Sugar ( Array(..) ) + +import qualified Data.Array.Accelerate.Array.Data as A +import qualified Data.Array.Accelerate.IO as A +import qualified Data.Array.Accelerate as A + + +data World = World + { + -- current state of the simulation + densityField :: !DensityField + , velocityField :: !VelocityField + + -- user input + , densitySource :: [(Index, Density)] + , velocitySource :: [(Index, Velocity)] + , currentSource :: Source + , displayDensity :: Bool + , displayVelocity :: Bool + } + deriving Show + +data Source = Density (Int, Int) + | Velocity (Int, Int) + | None + deriving (Eq, Show) + + +-- Initialisation -------------------------------------------------------------- +-- -------------- -- + +initialise :: Options -> IO World +initialise opt = + let width = get simulationWidth opt + height = get simulationHeight opt + in do + density <- initialDensity opt width height + velocity <- initialVelocity opt width height + + return $ World + { densityField = density + , velocityField = velocity + , densitySource = [] + , velocitySource = [] + , currentSource = None + , displayDensity = True + , displayVelocity = False + } + + +initialDensity :: Options -> Int -> Int -> IO DensityField +initialDensity opt width height + -- Load the file, and use the luminance value as scalar density + -- + | Just file <- get densityBMP opt + = do + arr <- readImageFromBMP file + let (Z:.h:.w) = A.arrayShape arr + + when (w /= width || h /= height) + $ error "fluid: density-bmp does not match width x height" + + return . run opt $ A.map densityOfRGBA (A.use arr) + + -- No density file given, just set the field to zero + -- + | otherwise + = return $ A.fromList (Z :. height :. width) (repeat 0) + + +initialVelocity :: Options -> Int -> Int -> IO VelocityField +initialVelocity opt width height + -- Load the file, and use the red and green channels for x- and y- velocity + -- values respectively + -- + | Just file <- get velocityBMP opt + = do + arr <- readImageFromBMP file + let (Z:.h:.w) = A.arrayShape arr + + when (w /= width || h /= height) + $ error "fluid: velocity-bmp does not match width x height" + + return . run opt $ A.map velocityOfRGBA (A.use arr) + + -- No density file given, just set to zero + -- + | otherwise + = return $ A.fromList (Z :. height :. width) (repeat (0,0)) + + +readImageFromBMP :: FilePath -> IO (Image RGBA) +readImageFromBMP file = do + bmp <- either (error . show) id `fmap` readBMP file + let (w,h) = bmpDimensions bmp + -- + A.fromByteString (Z :. h :. w) ((), unpackBMPToRGBA32 bmp) + + +densityOfRGBA :: Exp RGBA -> Exp Density +densityOfRGBA rgba = + let b = (0.11 / 255) * A.fromIntegral ((rgba `div` 0x100) .&. 0xFF) + g = (0.59 / 255) * A.fromIntegral ((rgba `div` 0x10000) .&. 0xFF) + r = (0.3 / 255) * A.fromIntegral ((rgba `div` 0x1000000) .&. 0xFF) + in + r + g + b + +velocityOfRGBA :: Exp RGBA -> Exp Velocity +velocityOfRGBA rgba = + let g = A.fromIntegral (-128 + A.fromIntegral ((rgba `div` 0x10000) .&. 0xFF) :: Exp Int32) + r = A.fromIntegral (-128 + A.fromIntegral ((rgba `div` 0x1000000) .&. 0xFF) :: Exp Int32) + in + A.lift (r * 0.001, g * 0.001) + + +-- Rendering ------------------------------------------------------------------- +-- --------- -- + +render :: Options -> World -> IO Picture +render opt world = do + den <- if displayDensity world + then renderDensity $ densityField world + else return blank + + vel <- if displayVelocity world + then renderVelocity $ velocityField world + else return blank + -- + return $ Scale zoom zoom $ Pictures [ den, vel ] + where + zoom = fromIntegral $ get displayScale opt + + +renderDensity :: DensityField -> IO Picture +renderDensity df@(Array _ ad) = do + dst <- mallocBytes (n*4) + fill 0 src dst + fptr <- newForeignPtr finalizerFree dst + return $ bitmapOfForeignPtr w h fptr False + where + ((),src) = A.ptrsOfArrayData ad + Z:.h:.w = A.arrayShape df + n = h*w + colour !f = let c = 0 `max` f `min` 1 + in floatToWord8 (255*c) + + fill !i !s !d | i >= n = return () + | otherwise = do c <- colour `fmap` peek s + poke d 0xFF -- A + pokeByteOff d 1 c -- B + pokeByteOff d 2 c -- G + pokeByteOff d 3 c -- R + fill (i+1) (plusPtr s 4) (plusPtr d 4) + + +renderVelocity :: VelocityField -> IO Picture +renderVelocity vf + = return + $ Translate (fromIntegral $ -w `div` 2) (fromIntegral $ -h `div` 2) + $ Pictures [ field (x,y) | y <- [2,7..h], x <- [2,7..w] ] + where + Z:.h:.w = A.arrayShape vf + field (x0,y0) = + let x = fromIntegral x0 + y = fromIntegral y0 + (u,v) = A.indexArray vf (Z:.y0:.x0) + in Color red $ Line [ (x,y), (x+u, y+v) ] + + +-- Float to Word8 conversion because the one in the GHC libraries doesn't have +-- enough specialisations and goes via Integer. +{-# INLINE floatToWord8 #-} +floatToWord8 :: Float -> Word8 +floatToWord8 f = fromIntegral (truncate f :: Int) + diff --git a/examples/quickcheck/Arbitrary/Array.hs b/examples/quickcheck/Arbitrary/Array.hs new file mode 100644 --- /dev/null +++ b/examples/quickcheck/Arbitrary/Array.hs @@ -0,0 +1,78 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE TypeOperators #-} + +module Arbitrary.Array where + +import Arbitrary.Shape + +import Data.List +import Test.QuickCheck +import Data.Array.Accelerate.Array.Sugar +import System.Random ( Random ) + + +instance (Elt e, Arbitrary e) => Arbitrary (Array DIM0 e) where + arbitrary = do + e <- arbitrary + return $! fromList Z [e] + -- + shrink arr = + [ fromList Z [x] | x <- shrink (arr ! Z) ] + + +instance (Elt e, Arbitrary e) => Arbitrary (Array DIM1 e) where + arbitrary = do + sh <- sized arbitrarySmallShape + adata <- vectorOf (size sh) arbitrary + return $! fromList sh adata + -- + shrink arr = + let (Z :. n) = shape arr + indices = [ map (Z:.) (nub sz) | sz <- shrink [0 .. n-1] ] + in + [ fromList (Z :. length sl) (map (arr!) sl) | sl <- indices ] + + +instance (Elt e, Arbitrary e) => Arbitrary (Array DIM2 e) where + arbitrary = do + sh <- sized arbitrarySmallShape + adata <- vectorOf (size sh) arbitrary + return $! fromList sh adata + -- + shrink arr = + let (Z :. width :. height) = shape arr + in + [ fromList (Z :. length slx :. length sly) [ arr ! (Z:.x:.y) | x <- slx, y <- sly ] + | slx <- map nub $ shrink [0 .. width - 1] + , sly <- map nub $ shrink [0 .. height - 1] + ] + + +arbitraryArrayOfShape + :: (Shape sh, Elt e, Arbitrary e) + => sh + -> Gen (Array sh e) +arbitraryArrayOfShape sh = do + adata <- vectorOf (size sh) arbitrary + return $! fromList sh adata + +arbitrarySegmentedArray + :: (Integral i, Shape sh, Elt e, Arbitrary sh, Arbitrary e) + => Segments i -> Gen (Array (sh :. Int) e) +arbitrarySegmentedArray segs = do + let sz = fromIntegral . sum $ toList segs + sh <- sized $ \n -> arbitrarySmallShape (n `div` 2) + adata <- vectorOf (size sh * sz) arbitrary + return $! fromList (sh :. sz) adata + +arbitrarySegments :: (Elt i, Integral i, Arbitrary i, Random i) => Gen (Segments i) +arbitrarySegments = do + seg <- listOf (sized $ \n -> choose (0, fromIntegral n)) + return $! fromList (Z :. length seg) seg + +arbitrarySegments1 :: (Elt i, Integral i, Arbitrary i, Random i) => Gen (Segments i) +arbitrarySegments1 = do + seg <- listOf (sized $ \n -> choose (1, fromIntegral n)) + return $! fromList (Z :. length seg) seg + diff --git a/examples/quickcheck/Arbitrary/Shape.hs b/examples/quickcheck/Arbitrary/Shape.hs new file mode 100644 --- /dev/null +++ b/examples/quickcheck/Arbitrary/Shape.hs @@ -0,0 +1,51 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE TypeOperators #-} + +module Arbitrary.Shape where + +import Test.QuickCheck +import Data.Array.Accelerate.Array.Sugar + + +instance Arbitrary Z where + arbitrary = return Z + shrink = return + +-- Keep shapes relatively small, because we are generating arrays via lists, but +-- scale the 'sized' parameter so we still get non-trivial test vectors. +-- +instance (Shape sh, Arbitrary sh) => Arbitrary (sh :. Int) where + arbitrary = do + sh <- arbitrary + sz <- sized $ \n -> + let nMax = 2048 `div` max 1 (size sh) + nMaxed = (n*4) `mod` nMax + in + choose (0, nMaxed) + -- + return (sh :. sz) + + shrink (sh :. sz) = [ sh' :. sz' | sh' <- shrink sh, sz' <- shrink sz ] + + +-- Generate an arbitrary shape, where each dimension is less than some specific +-- value. This will generate zero-sized shapes, but will take care not to +-- generate too many of them. +-- +arbitrarySmallShape :: forall sh. (Shape sh, Arbitrary sh) => Int -> Gen sh +arbitrarySmallShape maxDim = resize maxDim arbitrary + +-- Generate a shape as above, but also restrict the dimensions to be either all +-- zero, or all non-zero. +-- +arbitrarySmallNonEmptyShape :: forall sh. (Shape sh, Arbitrary sh) => Int -> Gen sh +arbitrarySmallNonEmptyShape maxDim = do + sh <- resize maxDim arbitrary :: Gen sh + let ix = map (`mod` max 1 maxDim) (shapeToList sh) + sh' + | all (== 0) ix = ix + | otherwise = map (max 1) ix + -- + return (listToShape sh' :: sh) + diff --git a/examples/quickcheck/Config.hs b/examples/quickcheck/Config.hs new file mode 100644 --- /dev/null +++ b/examples/quickcheck/Config.hs @@ -0,0 +1,150 @@ +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE CPP #-} + +module Config ( + + -- options & test configuration + Options, + processArgs, + optBackend, double, float, int64, int32, int16, int8, + + -- running tests + run, run1 + +) where + +import Data.Char +import Data.List +import Data.Label +import Control.Monad +import Test.Framework +import System.Console.GetOpt ( OptDescr(..), ArgDescr(..) ) + +import Data.Array.Accelerate ( Arrays, Acc ) +import qualified Data.Array.Accelerate.Interpreter as Interp + +#ifdef ACCELERATE_CUDA_BACKEND +import qualified Foreign.CUDA.Driver as CUDA +import qualified Data.Array.Accelerate.CUDA as CUDA +#endif + +data Backend = Interpreter +#ifdef ACCELERATE_CUDA_BACKEND + | CUDA +#endif + deriving (Eq, Enum, Bounded, Show) + + +data Options = Options + { + _optBackend :: !Backend, + _double :: !Bool, + _float :: !Bool, + _int64 :: !Bool, + _int32 :: !Bool, + _int16 :: !Bool, + _int8 :: !Bool + } + deriving Show + +$( mkLabels [''Options] ) + +defaultOptions :: Options +defaultOptions = Options + { _optBackend = maxBound + , _double = True + , _float = True + , _int64 = True + , _int32 = True + , _int16 = True + , _int8 = True + } + + +backends :: [OptDescr (Options -> Options)] +backends = + [ Option [] ["interpreter"] (NoArg (set optBackend Interpreter)) "reference implementation (sequential)" +#ifdef ACCELERATE_CUDA_BACKEND + , Option [] ["cuda"] (NoArg (set optBackend CUDA)) "implementation for NVIDIA GPUs (parallel)" +#endif + ] + + +run :: Arrays a => Options -> Acc a -> a +run opts = case _optBackend opts of + Interpreter -> Interp.run +#ifdef ACCELERATE_CUDA_BACKEND + CUDA -> CUDA.run +#endif + +run1 :: (Arrays a, Arrays b) => Options -> (Acc a -> Acc b) -> a -> b +run1 opts f = case _optBackend opts of + Interpreter -> head . Interp.stream f . return +#ifdef ACCELERATE_CUDA_BACKEND + CUDA -> CUDA.run1 f +#endif + + +-- Display a very basic usage message info: specifically, list the available +-- backends and highlight the active one. +-- +usageInfo :: Options -> String +usageInfo opts = unlines (header : table) + where + active this = if this == map toLower (show $ get optBackend opts) then "*" else "" + (ss,bs,ds) = unzip3 $ map (\(b,d) -> (active b, b, d)) $ concatMap extract backends + table = zipWith3 paste (sameLen ss) (sameLen bs) ds + paste x y z = " " ++ x ++ " " ++ y ++ " " ++ z + sameLen xs = flushLeft ((maximum . map length) xs) xs + flushLeft n xs = [ take n (x ++ repeat ' ') | x <- xs ] + -- + extract (Option _ los _ descr) = + let losFmt = intercalate ", " los + in case lines descr of + [] -> [(losFmt, "")] + (x:xs) -> (losFmt, x) : [ ("",x') | x' <- xs ] + -- + header = intercalate "\n" $ + [ "accelerate-quickcheck (c) 2012 The Accelerate Team" + , "" + , "Usage: accelerate-quickcheck [BACKEND] [OPTIONS]" + , "" + , "Available backends:" + ] + +-- Once the backend has been selected, we might need to enable or disable +-- certain classes of tests. +-- +-- CUDA: double precision is only supported on certain kinds of hardware. We +-- can't get the exact device the backend will choose to run on, but we do +-- know that it will pick the most capable device available. +-- +configureBackend :: Options -> IO Options +configureBackend opts = case _optBackend opts of + Interpreter -> return opts +#ifdef ACCELERATE_CUDA_BACKEND + CUDA -> do + CUDA.initialise [] + n <- CUDA.count + devs <- mapM CUDA.device [0 .. n-1] + props <- mapM CUDA.props devs + return $! opts { _double = any (\dev -> CUDA.computeCapability dev >= 1.3) props } +#endif + + +processArgs :: [String] -> IO (Options, RunnerOptions) +processArgs argv = do + args <- interpretArgs argv + (options, runner) <- case args of + Left msg -> error msg + Right (opts, rest) -> (,opts) `liftM` configureBackend (foldl parse1 defaultOptions rest) + -- + putStrLn $ usageInfo options + return $ (options, runner) + where + parse1 opts x = + case filter (\(Option _ [f] _ _) -> map toLower x `isPrefixOf` f) backends of + [Option _ _ (NoArg go) _] -> go opts + _ -> opts + diff --git a/examples/quickcheck/Main.hs b/examples/quickcheck/Main.hs new file mode 100644 --- /dev/null +++ b/examples/quickcheck/Main.hs @@ -0,0 +1,34 @@ + +module Main where + +import Test.Framework +import System.Environment + +import Config +import Test.Mapping +import Test.Reduction +import Test.PrefixSum +import Test.IndexSpace + + +main :: IO () +main = do + -- process command line args, and print a brief usage message + -- + (options, runner) <- processArgs =<< getArgs + + -- the default execution order uses some knowledge of what functionality is + -- required for each operation in the CUDA backend. + -- + defaultMainWithOpts + [ test_map options + , test_zipWith options + , test_foldAll options + , test_fold options + , test_prefixsum options -- requires fold + , test_foldSeg options -- requires scan + , test_permute options + , test_backpermute options + ] + runner + diff --git a/examples/quickcheck/Test/Base.hs b/examples/quickcheck/Test/Base.hs new file mode 100644 --- /dev/null +++ b/examples/quickcheck/Test/Base.hs @@ -0,0 +1,121 @@ +{-# LANGUAGE DefaultSignatures #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE TypeOperators #-} + +module Test.Base where + +import Prelude as P +import Test.QuickCheck +import Data.Array.Accelerate +import Data.Array.Accelerate.Array.Sugar as Sugar + + +-- A class of things that support almost-equality, so that we can disregard +-- small amounts of floating-point round-off error. +-- +class Similar a where + {-# INLINE (~=) #-} + (~=) :: a -> a -> Bool + default (~=) :: Eq a => a -> a -> Bool + (~=) = (==) + +infix 4 ~= + +instance Similar a => Similar [a] where + [] ~= [] = True + (x:xs) ~= (y:ys) = x ~= y && xs ~= ys + _ ~= _ = False + +instance (Similar a, Similar b) => Similar (a, b) where + (x1, y1) ~= (x2, y2) = x1 ~= x2 && y1 ~= y2 + + +instance Similar Int +instance Similar Int8 +instance Similar Int16 +instance Similar Int32 +instance Similar Int64 +instance Similar Word +instance Similar Word8 +instance Similar Word16 +instance Similar Word32 +instance Similar Word64 + +instance Similar Float where (~=) = absRelTol +instance Similar Double where (~=) = absRelTol + +{-# INLINE relTol #-} +relTol :: (Fractional a, Ord a) => a -> a -> a -> Bool +relTol epsilon x y = abs ((x-y) / (x+y+epsilon)) < epsilon + +{-# INLINE absRelTol #-} +absRelTol :: (Fractional a, Ord a) => a -> a -> Bool +absRelTol 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 + + +-- some "missing" instances +-- +instance Eq Z where + Z == Z = True + +instance (Eq sh, Shape sh) => Eq (sh :. Int) where + (sh1 :. n1) == (sh2 :. n2) = n1 == n2 && sh1 == sh2 + (sh1 :. n1) /= (sh2 :. n2) = n1 /= n2 || sh1 /= sh2 + +instance (Eq e, Eq sh, Shape sh) => Eq (Array sh e) where + a1 == a2 = arrayShape a1 == arrayShape a2 + && toList a1 == toList a2 + + a1 /= a2 = arrayShape a1 /= arrayShape a2 + || toList a1 /= toList a2 + +instance (Similar e, Eq sh, Shape sh) => Similar (Array sh e) where + a1 ~= a2 = arrayShape a1 == arrayShape a2 + && toList a1 ~= toList a2 + + +-- Print expected/received message on inequality +-- +infix 4 .==. +(.==.) :: (Similar a, Show a) => a -> a -> Property +(.==.) ans ref = printTestCase message (ref ~= ans) + where + message = unlines ["*** Expected:", show ref + ,"*** Received:", show ans ] + + +-- Miscellaneous +-- + +indexHead :: Shape sh => (sh:.Int) -> Int +indexHead (_ :. sz) = sz + +indexTail :: Shape sh => (sh:.Int) -> sh +indexTail (sh :. _) = sh + +isEmptyArray :: Shape sh => Array sh e -> Bool +isEmptyArray arr = arraySize (arrayShape arr) == 0 + +mkDim :: Shape sh => Int -> sh +mkDim n = listToShape (P.replicate n 0) + +dim0 :: DIM0 +dim0 = mkDim 0 + +dim1 :: DIM1 +dim1 = mkDim 1 + +dim2 :: DIM2 +dim2 = mkDim 2 + +splitEvery :: Int -> [a] -> [[a]] +splitEvery _ [] = cycle [[]] +splitEvery n xs = let (h,t) = splitAt n xs + in h : splitEvery n t + diff --git a/examples/quickcheck/Test/IndexSpace.hs b/examples/quickcheck/Test/IndexSpace.hs new file mode 100644 --- /dev/null +++ b/examples/quickcheck/Test/IndexSpace.hs @@ -0,0 +1,121 @@ +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ScopedTypeVariables #-} + +module Test.IndexSpace where + +import Prelude as P +import Data.Label +import Data.Maybe +import Data.Typeable +import Test.QuickCheck +import Test.Framework +import Test.Framework.Providers.QuickCheck2 + +import Config +import Test.Base +import Arbitrary.Array ( ) +import Data.Array.Unboxed as IArray hiding ( Array ) +import Data.Array.Accelerate as Acc +import Data.Array.Accelerate.Array.Sugar ( newArray ) + + +-- +-- Forward permutation --------------------------------------------------------- +-- + +test_permute :: Options -> Test +test_permute opt = testGroup "permute" $ catMaybes + [ testIntegralElt int8 (undefined :: Int8) + , testIntegralElt int16 (undefined :: Int16) + , testIntegralElt int32 (undefined :: Int32) + , testIntegralElt int64 (undefined :: Int64) + , testIntegralElt int8 (undefined :: Word8) + , testIntegralElt int16 (undefined :: Word16) + , testIntegralElt int32 (undefined :: Word32) + , testIntegralElt int64 (undefined :: Word64) + , testFloatingElt float (undefined :: Float) + , testFloatingElt double (undefined :: Double) + ] + where + testIntegralElt :: forall e. (Elt e, Integral e, IsIntegral e, Arbitrary e) => (Options :-> Bool) -> e -> Maybe Test + testIntegralElt ok _ + | P.not (get ok opt) = Nothing + | otherwise = Just $ testGroup (show (typeOf (undefined :: e))) + [ testProperty "histogram" (test_histogram Acc.fromIntegral P.fromIntegral :: Vector e -> Property) + ] + + testFloatingElt :: forall e. (Elt e, RealFrac e, IsFloating e, Arbitrary e) => (Options :-> Bool) -> e -> Maybe Test + testFloatingElt ok _ + | P.not (get ok opt) = Nothing + | otherwise = Just $ testGroup (show (typeOf (undefined :: e))) + [ testProperty "histogram" (test_histogram Acc.floor P.floor :: Vector e -> Property) + ] + + test_histogram f g xs = + sized $ \n -> run opt (histogramAcc n f xs) .==. histogramRef n g xs + + histogramAcc n f xs = + let n' = unit (constant n) + xs' = use xs + zeros = generate (constant (Z :. n)) (const 0) + ones = generate (shape xs') (const 1) + in + permute (+) zeros (\ix -> index1 $ f (xs' Acc.! ix) `mod` the n') ones + + histogramRef n f xs = + let arr :: UArray Int Int32 + arr = accumArray (+) 0 (0, n-1) [ (f e `mod` n, 1) | e <- toList xs ] + in + fromIArray arr + + +-- +-- Backward permutation -------------------------------------------------------- +-- + +test_backpermute :: Options -> Test +test_backpermute opt = testGroup "backpermute" $ catMaybes + [ testElt int8 (undefined :: Int8) + , testElt int16 (undefined :: Int16) + , testElt int32 (undefined :: Int32) + , testElt int64 (undefined :: Int64) + , testElt int8 (undefined :: Word8) + , testElt int16 (undefined :: Word16) + , testElt int32 (undefined :: Word32) + , testElt int64 (undefined :: Word64) + , testElt float (undefined :: Float) + , testElt double (undefined :: Double) + ] + where + testElt :: forall e. (Elt e, Similar e, Arbitrary e) => (Options :-> Bool) -> e -> Maybe Test + testElt ok _ + | P.not (get ok opt) = Nothing + | otherwise = Just $ testGroup (show (typeOf (undefined::e))) + [ testProperty "reverse" (test_reverse :: Array DIM1 e -> Property) + , testProperty "transpose" (test_transpose :: Array DIM2 e -> Property) + ] + + test_reverse xs = run opt (reverseAcc xs) .==. reverseRef xs + test_transpose xs = run opt (transposeAcc xs) .==. transposeRef xs + + -- Reverse a vector + -- + reverseRef xs = fromList (arrayShape xs) (reverse $ toList xs) + reverseAcc xs = + let xs' = use xs + n = unindex1 $ shape xs' + in + backpermute (shape xs') (\ix -> index1 $ n - unindex1 ix - 1) xs' + + -- Transpose a 2D matrix + -- + transposeAcc xs = + let xs' = use xs + swap = lift1 $ \(Z:.x:.y) -> Z:.y:.x :: Z :. Exp Int :. Exp Int + in + backpermute (swap (shape xs')) swap xs' + + transposeRef xs = + let swap (Z:.x:.y) = Z :. y :. x + in newArray (swap (arrayShape xs)) (\ix -> indexArray xs (swap ix)) + diff --git a/examples/quickcheck/Test/Mapping.hs b/examples/quickcheck/Test/Mapping.hs new file mode 100644 --- /dev/null +++ b/examples/quickcheck/Test/Mapping.hs @@ -0,0 +1,121 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE RankNTypes #-} + +module Test.Mapping where + +import Prelude as P +import Data.Label +import Data.Maybe +import Data.Typeable +import Test.QuickCheck +import Test.Framework +import Test.Framework.Providers.QuickCheck2 + +import Config +import Test.Base +import Arbitrary.Array () +import Data.Array.Accelerate as Acc +import Data.Array.Accelerate.Array.Sugar as Sugar +import qualified Data.Array.Accelerate.Array.Representation as R + +-- +-- Map ------------------------------------------------------------------------- +-- + +test_map :: Options -> Test +test_map opt = testGroup "map" $ catMaybes + [ testElt int8 (undefined :: Int8) + , testElt int16 (undefined :: Int16) + , testElt int32 (undefined :: Int32) + , testElt int64 (undefined :: Int64) + , testElt int8 (undefined :: Word8) + , testElt int16 (undefined :: Word16) + , testElt int32 (undefined :: Word32) + , testElt int64 (undefined :: Word64) + , testElt float (undefined :: Float) + , testElt double (undefined :: Double) + ] + where + test_square xs = run opt (Acc.map (\x -> x*x) (use xs)) .==. mapRef (\x -> x*x) xs + test_abs xs = run opt (Acc.map abs (use xs)) .==. mapRef abs xs + test_plus z xs = + let z' = unit (constant z) + in run opt (Acc.map (+ the z') (use xs)) .==. mapRef (+ z) xs + + testElt :: forall a. (Elt a, IsNum a, Similar a, Arbitrary a) + => (Options :-> Bool) + -> a + -> Maybe Test + testElt ok _ + | P.not (get ok opt) = Nothing + | otherwise = Just $ testGroup (show (typeOf (undefined :: a))) + [ testDim dim0 + , testDim dim1 + , testDim dim2 + ] + where + testDim :: forall sh. (Shape sh, Eq sh, Arbitrary sh, Arbitrary (Array sh a)) => sh -> Test + testDim sh = testGroup ("DIM" ++ show (dim sh)) + [ testProperty "abs" (test_abs :: Array sh a -> Property) + , testProperty "plus" (test_plus :: a -> Array sh a -> Property) + , testProperty "square" (test_square :: Array sh a -> Property) + ] + + +test_zipWith :: Options -> Test +test_zipWith opt = testGroup "zipWith" $ catMaybes + [ testElt int8 (undefined :: Int8) + , testElt int16 (undefined :: Int16) + , testElt int32 (undefined :: Int32) + , testElt int64 (undefined :: Int64) + , testElt int8 (undefined :: Word8) + , testElt int16 (undefined :: Word16) + , testElt int32 (undefined :: Word32) + , testElt int64 (undefined :: Word64) + , testElt float (undefined :: Float) + , testElt double (undefined :: Double) + ] + where + test_zip xs ys = run opt (Acc.zip (use xs) (use ys)) .==. zipWithRef (,) xs ys + test_plus xs ys = run opt (Acc.zipWith (+) (use xs) (use ys)) .==. zipWithRef (+) xs ys + test_min xs ys = run opt (Acc.zipWith Acc.min (use xs) (use ys)) .==. zipWithRef P.min xs ys + + testElt :: forall a. (Elt a, IsNum a, Ord a, Similar a, Arbitrary a) + => (Options :-> Bool) + -> a + -> Maybe Test + testElt ok _ + | P.not (get ok opt) = Nothing + | otherwise = Just $ testGroup (show (typeOf (undefined :: a))) + [ testDim dim0 + , testDim dim1 + , testDim dim2 + ] + where + testDim :: forall sh. (Shape sh, Eq sh, Arbitrary sh, Arbitrary (Array sh a)) => sh -> Test + testDim sh = testGroup ("DIM" ++ show (dim sh)) + [ testProperty "zip" (test_zip :: Array sh a -> Array sh a -> Property) + , testProperty "plus" (test_plus :: Array sh a -> Array sh a -> Property) + , testProperty "min" (test_min :: Array sh a -> Array sh a -> Property) + ] + + +-- Reference Implementation +-- ------------------------ + +mapRef :: (Shape sh, Elt b) => (a -> b) -> Array sh a -> Array sh b +mapRef f xs + = fromList (arrayShape xs) + $ P.map f + $ toList xs + +zipWithRef :: (Shape sh, Elt c) => (a -> b -> c) -> Array sh a -> Array sh b -> Array sh c +zipWithRef f xs ys = + let shx = fromElt (arrayShape xs) + shy = fromElt (arrayShape ys) + sh = toElt (R.intersect shx shy) + in + newArray sh (\ix -> f (xs Sugar.! ix) (ys Sugar.! ix)) + diff --git a/examples/quickcheck/Test/PrefixSum.hs b/examples/quickcheck/Test/PrefixSum.hs new file mode 100644 --- /dev/null +++ b/examples/quickcheck/Test/PrefixSum.hs @@ -0,0 +1,105 @@ +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ScopedTypeVariables #-} + +module Test.PrefixSum where + +import Prelude as P hiding ( sum ) +import Test.QuickCheck +import Data.Label +import Data.Maybe +import Data.Typeable +import Test.Framework +import Test.Framework.Providers.QuickCheck2 + +import Config +import Test.Base +import Arbitrary.Array () +import Data.Array.Accelerate as Acc + + +-- +-- prefix sum ------------------------------------------------------------------ +-- + +test_prefixsum :: Options -> Test +test_prefixsum opt = testGroup "prefix sum" $ catMaybes + [ testElt int8 (undefined :: Int8) + , testElt int16 (undefined :: Int16) + , testElt int32 (undefined :: Int32) + , testElt int64 (undefined :: Int64) + , testElt int8 (undefined :: Word8) + , testElt int16 (undefined :: Word16) + , testElt int32 (undefined :: Word32) + , testElt int64 (undefined :: Word64) + , testElt float (undefined :: Float) + , testElt double (undefined :: Double) + ] + where + testElt :: forall e. (Elt e, IsNum e, Ord e, Similar e, Arbitrary e) => (Options :-> Bool) -> e -> Maybe Test + testElt ok _ + | P.not (get ok opt) = Nothing + | otherwise = Just $ testGroup (show (typeOf (undefined :: e))) + [ testProperty "scanl" (test_scanl :: Vector e -> Property) + , testProperty "scanl'" (test_scanl' :: Vector e -> Property) + , testProperty "scanl1" (test_scanl1 :: Vector e -> Property) + , testProperty "scanr" (test_scanr :: Vector e -> Property) + , testProperty "scanr'" (test_scanr' :: Vector e -> Property) + , testProperty "scanr1" (test_scanr1 :: Vector e -> Property) + ] + + -- left scan + -- + test_scanl xs = run opt (Acc.scanl (+) 0 (use xs)) .==. scanlRef (+) 0 xs + test_scanl1 xs = run opt (Acc.scanl1 Acc.min (use xs)) .==. scanl1Ref P.min xs + test_scanl' xs = + let (vec, sum) = Acc.scanl' (+) 0 (use xs) + in (run opt vec, run opt sum) .==. scanl'Ref (+) 0 xs + + -- right scan + -- + test_scanr xs = run opt (Acc.scanr (+) 0 (use xs)) .==. scanrRef (+) 0 xs + test_scanr1 xs = run opt (Acc.scanr1 Acc.max (use xs)) .==. scanr1Ref P.max xs + test_scanr' xs = + let (vec, sum) = Acc.scanr' (+) 0 (use xs) + in (run opt vec, run opt sum) .==. scanr'Ref (+) 0 xs + + + +-- Reference implementation +-- ------------------------ + +scanlRef :: Elt e => (e -> e -> e) -> e -> Vector e -> Vector e +scanlRef f z vec = + let (Z :. n) = arrayShape vec + in Acc.fromList (Z :. n+1) . P.scanl f z . Acc.toList $ vec + +scanl'Ref :: Elt e => (e -> e -> e) -> e -> Vector e -> (Vector e, Scalar e) +scanl'Ref f z vec = + let (Z :. n) = arrayShape vec + result = P.scanl f z (Acc.toList vec) + in (Acc.fromList (Z :. n) result, Acc.fromList Z (P.drop n result)) + +scanl1Ref :: Elt e => (e -> e -> e) -> Vector e -> Vector e +scanl1Ref f vec + = Acc.fromList (arrayShape vec) + . P.scanl1 f + . Acc.toList $ vec + +scanrRef :: Elt e => (e -> e -> e) -> e -> Vector e -> Vector e +scanrRef f z vec = + let (Z :. n) = arrayShape vec + in Acc.fromList (Z :. n+1) . P.scanr f z . Acc.toList $ vec + +scanr'Ref :: Elt e => (e -> e -> e) -> e -> Vector e -> (Vector e, Scalar e) +scanr'Ref f z vec = + let (Z :. n) = arrayShape vec + result = P.scanr f z (Acc.toList vec) + in (Acc.fromList (Z :. n) (P.tail result), Acc.fromList Z result) + +scanr1Ref :: Elt e => (e -> e -> e) -> Vector e -> Vector e +scanr1Ref f vec + = Acc.fromList (arrayShape vec) + . P.scanr1 f + . Acc.toList $ vec + diff --git a/examples/quickcheck/Test/Reduction.hs b/examples/quickcheck/Test/Reduction.hs new file mode 100644 --- /dev/null +++ b/examples/quickcheck/Test/Reduction.hs @@ -0,0 +1,225 @@ +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ScopedTypeVariables #-} + +module Test.Reduction where + +import Prelude as P +import Data.List +import Data.Label +import Data.Maybe +import Data.Typeable +import Test.QuickCheck +import Test.Framework +import Test.Framework.Providers.QuickCheck2 + +import Config +import Test.Base +import Arbitrary.Array +import Data.Array.Accelerate as Acc +import Data.Array.Accelerate.Array.Sugar as Sugar + + + +-- +-- Reduction ------------------------------------------------------------------- +-- + +-- foldAll +-- ------- + +test_foldAll :: Options -> Test +test_foldAll opt = testGroup "foldAll" $ catMaybes + [ testElt int8 (undefined :: Int8) + , testElt int16 (undefined :: Int16) + , testElt int32 (undefined :: Int32) + , testElt int64 (undefined :: Int64) + , testElt int8 (undefined :: Word8) + , testElt int16 (undefined :: Word16) + , testElt int32 (undefined :: Word32) + , testElt int64 (undefined :: Word64) + , testElt float (undefined :: Float) + , testElt double (undefined :: Double) + ] + where + testElt :: forall e. (Elt e, IsNum e, Ord e, Similar e, Arbitrary e) => (Options :-> Bool) -> e -> Maybe Test + testElt ok _ + | P.not (get ok opt) = Nothing + | otherwise = Just $ testGroup (show (typeOf (undefined :: e))) + [ testDim dim0 + , testDim dim1 + , testDim dim2 + ] + where + testDim :: forall sh. (Shape sh, Arbitrary sh, Arbitrary (Array sh e)) => sh -> Test + testDim sh = testGroup ("DIM" ++ show (dim sh)) + [ + testProperty "sum" (test_sum :: Array sh e -> Property) + , testProperty "non-neutral sum" (test_sum' :: Array sh e -> e -> Property) + , testProperty "minimum" (test_min :: Array sh e -> Property) + , testProperty "maximum" (test_max :: Array sh e -> Property) + ] + -- + -- The tests + -- + test_min xs + = arraySize (arrayShape xs) > 0 + ==> run opt (Acc.fold1All Acc.min (use xs)) .==. fold1AllRef P.min xs + + test_max xs + = arraySize (arrayShape xs) > 0 + ==> run opt (Acc.fold1All Acc.max (use xs)) .==. fold1AllRef P.max xs + + test_sum xs = run opt (Acc.foldAll (+) 0 (use xs)) .==. foldAllRef (+) 0 xs + test_sum' xs z = + let z' = unit (constant z) + in run opt (Acc.foldAll (+) (the z') (use xs)) .==. foldAllRef (+) z xs + + + +-- multidimensional fold +-- --------------------- + +test_fold :: Options -> Test +test_fold opt = testGroup "fold" $ catMaybes + [ testElt int8 (undefined :: Int8) + , testElt int16 (undefined :: Int16) + , testElt int32 (undefined :: Int32) + , testElt int64 (undefined :: Int64) + , testElt int8 (undefined :: Word8) + , testElt int16 (undefined :: Word16) + , testElt int32 (undefined :: Word32) + , testElt int64 (undefined :: Word64) + , testElt float (undefined :: Float) + , testElt double (undefined :: Double) + ] + where + testElt :: forall e. (Elt e, IsNum e, Ord e, Similar e, Arbitrary e) => (Options :-> Bool) -> e -> Maybe Test + testElt ok _ + | P.not (get ok opt) = Nothing + | otherwise = Just $ testGroup (show (typeOf (undefined :: e))) + [ testDim dim1 + , testDim dim2 + ] + where + testDim :: forall sh. (Shape sh, Eq sh, Arbitrary sh, Arbitrary (Array (sh:.Int) e)) => (sh:.Int) -> Test + testDim sh = testGroup ("DIM" ++ show (dim sh)) + [ + testProperty "sum" (test_sum :: Array (sh :. Int) e -> Property) + , testProperty "non-neutral sum" (test_sum' :: Array (sh :. Int) e -> e -> Property) + , testProperty "minimum" (test_min :: Array (sh :. Int) e -> Property) + , testProperty "maximum" (test_max :: Array (sh :. Int) e -> Property) + ] + -- + -- The tests + -- + test_min xs + = indexHead (arrayShape xs) > 0 + ==> run opt (Acc.fold1 Acc.min (use xs)) .==. fold1Ref P.min xs + + test_max xs + = indexHead (arrayShape xs) > 0 + ==> run opt (Acc.fold1 Acc.max (use xs)) .==. fold1Ref P.max xs + + test_sum xs = run opt (Acc.fold (+) 0 (use xs)) .==. foldRef (+) 0 xs + test_sum' xs z = + let z' = unit (constant z) + in run opt (Acc.fold (+) (the z') (use xs)) .==. foldRef (+) z xs + + +-- segmented fold +-- -------------- + +test_foldSeg :: Options -> Test +test_foldSeg opt = testGroup "foldSeg" $ catMaybes + [ testElt int8 (undefined :: Int8) + , testElt int16 (undefined :: Int16) + , testElt int32 (undefined :: Int32) + , testElt int64 (undefined :: Int64) + , testElt int8 (undefined :: Word8) + , testElt int16 (undefined :: Word16) + , testElt int32 (undefined :: Word32) + , testElt int64 (undefined :: Word64) + , testElt float (undefined :: Float) + , testElt double (undefined :: Double) + ] + where + testElt :: forall e. (Elt e, IsNum e, Ord e, Similar e, Arbitrary e) => (Options :-> Bool) -> e -> Maybe Test + testElt ok _ + | P.not (get ok opt) = Nothing + | otherwise = Just $ testGroup (show (typeOf (undefined :: e))) + [ testDim dim1 + , testDim dim2 + ] + where + testDim :: forall sh. (Shape sh, Eq sh, Arbitrary sh, Arbitrary (Array (sh:.Int) e)) => (sh:.Int) -> Test + testDim sh = testGroup ("DIM" ++ show (dim sh)) + [ + testProperty "sum" + $ forAll arbitrarySegments $ \(seg :: Segments Int32) -> + forAll (arbitrarySegmentedArray seg) $ \(xs :: Array (sh:.Int) e) -> + run opt (Acc.foldSeg (+) 0 (use xs) (use seg)) .==. foldSegRef (+) 0 xs seg + + , testProperty "non-neutral sum" + $ forAll arbitrarySegments $ \(seg :: Segments Int32) -> + forAll (arbitrarySegmentedArray seg) $ \(xs :: Array (sh:.Int) e) -> + forAll arbitrary $ \z -> + let z' = unit (constant z) + in run opt (Acc.foldSeg (+) (the z') (use xs) (use seg)) .==. foldSegRef (+) z xs seg + + , testProperty "minimum" + $ forAll arbitrarySegments1 $ \(seg :: Segments Int32) -> + forAll (arbitrarySegmentedArray seg) $ \(xs :: Array (sh:.Int) e) -> + run opt (Acc.fold1Seg Acc.min (use xs) (use seg)) .==. fold1SegRef P.min xs seg + ] + + +-- Reference implementation +-- ------------------------ + +foldAllRef :: Elt e => (e -> e -> e) -> e -> Array sh e -> Array Z e +foldAllRef f z + = Acc.fromList Z + . return + . foldl f z + . Acc.toList + +fold1AllRef :: Elt e => (e -> e -> e) -> Array sh e -> Array Z e +fold1AllRef f + = Acc.fromList Z + . return + . foldl1 f + . Acc.toList + +foldRef :: (Shape sh, Elt e) => (e -> e -> e) -> e -> Array (sh :. Int) e -> Array sh e +foldRef f z arr = + let (sh :. n) = arrayShape arr + in fromList sh [ foldl f z sub | sub <- splitEvery n (toList arr) ] + +fold1Ref :: (Shape sh, Elt e) => (e -> e -> e) -> Array (sh :. Int) e -> Array sh e +fold1Ref f arr = + let (sh :. n) = arrayShape arr + in fromList sh [ foldl1 f sub | sub <- splitEvery n (toList arr) ] + +foldSegRef :: (Shape sh, Elt e, Elt i, Integral i) => (e -> e -> e) -> e -> Array (sh :. Int) e -> Segments i -> Array (sh :. Int) e +foldSegRef f z arr seg = fromList (sh :. sz) $ concat [ foldseg sub | sub <- splitEvery n (toList arr) ] + where + (sh :. n) = arrayShape arr + (Z :. sz) = arrayShape seg + seg' = toList seg + foldseg xs = P.map (foldl' f z) (split seg' xs) + +fold1SegRef :: (Shape sh, Elt e, Elt i, Integral i) => (e -> e -> e) -> Array (sh :. Int) e -> Segments i -> Array (sh :. Int) e +fold1SegRef f arr seg = fromList (sh :. sz) $ concat [ foldseg sub | sub <- splitEvery n (toList arr) ] + where + (sh :. n) = arrayShape arr + (Z :. sz) = arrayShape seg + seg' = toList seg + foldseg xs = P.map (foldl1' f) (split seg' xs) + +split :: Integral i => [i] -> [a] -> [[a]] +split [] _ = [] +split (i:is) vs = + let (h,t) = splitAt (P.fromIntegral i) vs + in h : split is t + diff --git a/examples/tests/Benchmark.hs b/examples/tests/Benchmark.hs new file mode 100644 --- /dev/null +++ b/examples/tests/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/examples/tests/Config.hs b/examples/tests/Config.hs new file mode 100644 --- /dev/null +++ b/examples/tests/Config.hs @@ -0,0 +1,191 @@ +{-# LANGUAGE CPP, DeriveDataTypeable #-} + +module Config where + +import Data.Version +import Text.PrettyPrint +import System.Console.CmdArgs +import Paths_accelerate_examples + +import Data.Array.Accelerate ( Arrays, Acc ) +import qualified Data.Array.Accelerate.Interpreter as Interpreter + +-- ----------------------------------------------------------------------------- +-- [NOTE: Adding backend support] +-- +-- To add support for additional Accelerate backends: +-- +-- 1. Edit the cabal file so that the backend can be optionally included via +-- CPP macros, as has been done for the CUDA backend. Add the appropriate +-- import statement below. +-- +-- 2. Add a new constructor to the 'Backend' data type, and add this case to +-- the 'backend' dispatch function. This function is called to evaluate an +-- Accelerate expression. +-- +-- 3. Add the new constructor to the 'cfgBackend' list in 'defaultConfig'. +-- This will make the backend appear in the command-line options. +-- +#ifdef ACCELERATE_CUDA_BACKEND +import qualified Data.Array.Accelerate.CUDA as CUDA +#endif + +#ifdef ACCELERATE_OPENCL_BACKEND +import qualified Data.Array.Accelerate.OpenCL as OpenCL +#endif + +-- 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 +#ifdef ACCELERATE_OPENCL_BACKEND + | OpenCL +#endif + deriving (Show, Data, Typeable) + + +-- How to evaluate Accelerate programs with the chosen backend? +-- +backend :: Arrays a => Config -> Acc a -> a +backend cfg = + case cfgBackend cfg of + Interpreter -> Interpreter.run +#ifdef ACCELERATE_CUDA_BACKEND + CUDA -> CUDA.run +#endif +#ifdef ACCELERATE_OPENCL_BACKEND + OpenCL -> OpenCL.run +#endif + +-- +-- ----------------------------------------------------------------------------- +-- + +-- 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 = + let numElements = 100000 + in Config + { + cfgBackend = enum + [ Interpreter + &= help "Reference implementation (sequential)" +#ifdef ACCELERATE_CUDA_BACKEND + , CUDA + &= explicit + &= name "cuda" + &= help "Implementation for NVIDIA GPUs (parallel)" +#endif +#ifdef ACCELERATE_OPENCL_BACKEND + , OpenCL + &= explicit + &= name "opencl" + &= help "Implementation for OpenCL (parallel)" +#endif + + ] + + , cfgVerify = def + &= explicit + &= name "k" + &= name "verify" + &= help "Only verify examples, do not run timing tests" + + , cfgElements = numElements + &= explicit + &= name "n" + &= name "size" + &= help ("Canonical test data size (" ++ shows numElements ")") + + , 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:" + , " 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/examples/tests/Main.hs b/examples/tests/Main.hs new file mode 100644 --- /dev/null +++ b/examples/tests/Main.hs @@ -0,0 +1,78 @@ +{-# 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, whenNormal, 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 + whenNormal . putStrLn $ + unlines ["running with " ++ shows (cfgBackend config) " backend" + ,"to display available options, rerun with '--help'" + ] + valid <- runVerify config tests + -- + unless (null valid || cfgVerify config) $ runTiming config valid diff --git a/examples/tests/Random.hs b/examples/tests/Random.hs new file mode 100644 --- /dev/null +++ b/examples/tests/Random.hs @@ -0,0 +1,73 @@ +{-# 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.Unsafe as U +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 = U.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 = U.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/examples/tests/Test.hs b/examples/tests/Test.hs new file mode 100644 --- /dev/null +++ b/examples/tests/Test.hs @@ -0,0 +1,237 @@ +{-# 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 FoldSeg +import qualified ScanSeg +import qualified Stencil +import qualified Stencil2 +import qualified Permute +import qualified Backpermute +import qualified Vector +import qualified Gather +import qualified Scatter + +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 +import qualified VectorCopy () -- FIXME: Need to add tests from here +#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.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 + + +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 "fold-2d-maximum" "maximum along innermost matrix dimension" $ Fold.run2d "maximum-2d" n + , mkTest "fold-2d-minimum" "minimum along innermost matrix dimension" $ Fold.run2d "minimum-2d" n + , mkTest "foldseg-sum" "segmented vector reduction" $ FoldSeg.run "sum" n + , mkTest "scanseg-sum" "segmented vector prescanl (+) 0" $ 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 + , mkTest "init" "vector init" $ Vector.run "init" n + , mkTest "tail" "vector tail" $ Vector.run "tail" n + , mkTest "take" "vector take" $ Vector.run "take" n + , mkTest "drop" "vector drop" $ Vector.run "drop" n + , mkTest "slit" "vector slit" $ Vector.run "slit" n + , mkTest "gather" "backpermute via index mapping vector" $ Gather.run "gather" n + , mkTest "gather-if" "cond. backpermute via index mapping vector" $ Gather.run "gather-if" n + , mkTest "scatter" "permute via index mapping vector" $ Scatter.run "scatter" n + , mkTest "scatter-if" "cond. permute via index mapping vector" $ Scatter.run "scatter-if" 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) + -- + , mkIO "bound-variables" "stencil2" $ return (show Stencil2.varUse) + ] + 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 + + +-- 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 + putStr (title test ++ ": ") >> hFlush stdout + result <- case test of + Test _ _ ref acc cvt -> + return $ case validate (ref ()) (cvt $ run acc) of + [] -> Ok + errs -> Failed $ + if quiet then ("(" ++ show (length errs) ++ " differences)") + else unlines . ("":) $ map (\(i,v) -> ">>> " ++ shows i " : " ++ show v) errs + TestNoRef _ _ acc -> return $ run acc `seq` Ok + TestIO _ _ act -> act >>= \v -> v `seq` return Ok + -- + 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/examples/tests/Util.hs b/examples/tests/Util.hs new file mode 100644 --- /dev/null +++ b/examples/tests/Util.hs @@ -0,0 +1,20 @@ + +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' + +replicateM' :: Int -> IO a -> IO [a] +replicateM' n x = sequence' (replicate n x) + diff --git a/examples/tests/Validate.hs b/examples/tests/Validate.hs new file mode 100644 --- /dev/null +++ b/examples/tests/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/examples/tests/image-processing/Canny.hs b/examples/tests/image-processing/Canny.hs new file mode 100644 --- /dev/null +++ b/examples/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 ! clamp (x+offsetx, y+offsety) + rev = gradM ! 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/examples/tests/image-processing/IntegralImage.hs b/examples/tests/image-processing/IntegralImage.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/tests/image-processing/PGM.hs b/examples/tests/image-processing/PGM.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/tests/io/BlockCopy.hs b/examples/tests/io/BlockCopy.hs new file mode 100644 --- /dev/null +++ b/examples/tests/io/BlockCopy.hs @@ -0,0 +1,91 @@ +{-# LANGUAGE ForeignFunctionInterface, ScopedTypeVariables, TypeOperators #-} + +module BlockCopy where + +-- standard libraries +import Prelude as P +import Foreign.Ptr +import Foreign.C +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/examples/tests/io/VectorCopy.hs b/examples/tests/io/VectorCopy.hs new file mode 100644 --- /dev/null +++ b/examples/tests/io/VectorCopy.hs @@ -0,0 +1,74 @@ +{-# LANGUAGE TemplateHaskell, TypeFamilies #-} + +module VectorCopy where + +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 :: IO Bool +test = $quickCheckAll diff --git a/examples/tests/io/fill_with_values.cpp b/examples/tests/io/fill_with_values.cpp new file mode 100644 --- /dev/null +++ b/examples/tests/io/fill_with_values.cpp @@ -0,0 +1,81 @@ +#include +#include +#include + +/* Returns one if it's filled with even values starting at 0 */ +template +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/examples/tests/primitives/Backpermute.hs b/examples/tests/primitives/Backpermute.hs new file mode 100644 --- /dev/null +++ b/examples/tests/primitives/Backpermute.hs @@ -0,0 +1,70 @@ +{-# LANGUAGE TypeOperators #-} + +module Backpermute where + +import Util +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/examples/tests/primitives/Fold.hs b/examples/tests/primitives/Fold.hs new file mode 100644 --- /dev/null +++ b/examples/tests/primitives/Fold.hs @@ -0,0 +1,90 @@ +{-# LANGUAGE FlexibleContexts, TypeOperators #-} + +module Fold where + +import Util +import Random + +import Control.Monad +import Control.Exception +import System.Random.MWC +import Data.List +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 (foldl' (+) 0) +prodRef = toUA (foldl' (*) 1) +maxRef = toUA (foldl1' P.max) +minRef = toUA (foldl1' P.min) + + +-- 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, min2DRef, max2DRef :: UArray (Int,Int) Float -> UArray Int Float +sum2DRef = foldU2D (+) 0 +prod2DRef = foldU2D (*) 1 +min2DRef = foldU2D P.min ( 1/0) +max2DRef = foldU2D P.max (-1/0) + + +-- 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 + "maximum-2d" -> go max2DRef maxAcc + "minimum-2d" -> go min2DRef minAcc + x -> error $ "unknown variant: " ++ x + where + {-# NOINLINE run_ref #-} + run_ref f xs () = f xs + run_acc f xs () = f xs + diff --git a/examples/tests/primitives/FoldSeg.hs b/examples/tests/primitives/FoldSeg.hs new file mode 100644 --- /dev/null +++ b/examples/tests/primitives/FoldSeg.hs @@ -0,0 +1,62 @@ + +module FoldSeg where + +import Random + +import Data.List +import System.IO +import System.Random.MWC +import Data.Array.Unboxed +import Data.Array.Accelerate as A +import Prelude as P + + +-- segmented reduction +-- ------------------- + +sumSegAcc :: Vector Float -> Segments Int32 -> Acc (Vector Float) +sumSegAcc xs seg + = let xs' = use xs + seg' = use seg + in + A.foldSeg (+) 0 xs' seg' + +sumSegRef :: UArray Int Float -> UArray Int Int32 -> UArray Int Float +sumSegRef xs seg + = listArray (bounds seg) + $ list_foldSeg (+) 0 (elems xs) (elems seg) + +list_foldSeg :: (a -> a -> a) -> a -> [a] -> [Int32] -> [a] +list_foldSeg f s xs seg = P.map (foldl' f s) (split seg xs) + where + split [] _ = [] + split _ [] = [] + split (i:is) vs = + let (h,t) = splitAt (P.fromIntegral i) vs + in h : split is t + + +-- main +-- ---- + +run :: String -> Int -> IO (() -> UArray Int Float, () -> Acc (Vector Float)) +run alg m = withSystemRandom $ \gen -> do + -- generate segments + -- + let n = P.round $ sqrt (P.fromIntegral m :: Double) + seg <- randomUArrayR (0, 2*n) gen (P.fromIntegral n) + seg' <- convertUArray seg + + -- generate elements + -- + let x = P.fromIntegral $ sum (elems seg) + vec <- randomUArrayR (-1,1) gen x + vec' <- convertUArray vec + + -- super-happy-fun-times + -- + let go f g = return (\() -> f vec seg, \() -> g vec' seg') + case alg of + "sum" -> go sumSegRef sumSegAcc + unknown -> error $ "unknown variant: " ++ unknown + diff --git a/examples/tests/primitives/Gather.hs b/examples/tests/primitives/Gather.hs new file mode 100644 --- /dev/null +++ b/examples/tests/primitives/Gather.hs @@ -0,0 +1,73 @@ +{-# LANGUAGE FlexibleContexts #-} + +module Gather where + +import Random + +import System.Random.MWC +import Data.Array.Unboxed +import Data.Array.Accelerate as Acc hiding ((!)) +import Prelude as P + + +-- Tests +-- ----- + +gatherAcc :: Vector Int -> Vector Float -> Acc (Vector Float) +gatherAcc mapV inputV = Acc.gather (use mapV) (use inputV) + +gatherIfAcc :: Vector Int -> Vector Int -> Vector Float -> Vector Float -> Acc (Vector Float) +gatherIfAcc mapV maskV defaultV inputV + = Acc.gatherIf (use mapV) (use maskV) evenAcc (use defaultV) (use inputV) + +evenAcc :: Exp Int -> Exp Bool +evenAcc v = (v `mod` 2) ==* 0 + + +gatherRef :: UArray Int Int -> UArray Int Float -> UArray Int Float +gatherRef mapV inputV = amap (\ix -> inputV ! ix) mapV + +gatherIfRef :: UArray Int Int -> UArray Int Int -> UArray Int Float -> UArray Int Float -> UArray Int Float +gatherIfRef mapV maskV defaultV inputV + = listArray (bounds mapV) + $ P.map (\(mIx, mV, dV) -> if evenRef mV then (inputV ! mIx) else dV) + $ P.zip3 mapL maskL defaultL + where + mapL = elems mapV + maskL = elems maskV + defaultL = elems defaultV + +evenRef :: Int -> Bool +evenRef = even + + +-- 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 + + mapV <- randomUArrayR (0, n - 1) gen n + mapV' <- convertUArray mapV + + maskV <- randomUArrayR (0, n) gen n + maskV' <- convertUArray maskV + + defaultV <- randomUArrayR (-1, 1) gen n + defaultV' <- convertUArray defaultV + + -- + let go f g = return (run_ref f vec, run_acc g vec') + + case alg of + "gather" -> go (gatherRef mapV) (gatherAcc mapV') + "gather-if" -> go (gatherIfRef mapV maskV defaultV) (gatherIfAcc mapV' maskV' defaultV') + x -> error $ "unknown variant: " ++ x + + where + {-# NOINLINE run_ref #-} + run_ref f xs () = f xs + run_acc f xs () = f xs + + diff --git a/examples/tests/primitives/Map.hs b/examples/tests/primitives/Map.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/tests/primitives/Permute.hs b/examples/tests/primitives/Permute.hs new file mode 100644 --- /dev/null +++ b/examples/tests/primitives/Permute.hs @@ -0,0 +1,40 @@ +module Permute where + +import Random + +import Data.Int +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 Int32) +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 Int32 +histogramRef (m,n) vec = + accumArray (+) 0 (0,n-m-1) [(P.floor e, 1) | e <- elems vec] + + +-- Main +-- ---- +run :: String -> Int -> IO (() -> UArray Int Int32, () -> Acc (Vector Int32)) +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/examples/tests/primitives/ScanSeg.hs b/examples/tests/primitives/ScanSeg.hs new file mode 100644 --- /dev/null +++ b/examples/tests/primitives/ScanSeg.hs @@ -0,0 +1,63 @@ + +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 Int32 -> Acc (Vector Float) +prefixSumSegAcc xs seg + = let + xs' = use xs + seg' = use seg + in + prescanlSeg (+) 0 xs' seg' + + +prefixSumSegRef :: UArray Int Float -> UArray Int Int32 -> UArray Int Float +prefixSumSegRef xs seg + = listArray (bounds xs) + $ list_prescanlSeg (+) 0 (elems xs) (elems seg) + +list_prescanlSeg :: (a -> a -> a) -> a -> [a] -> [Int32] -> [a] +list_prescanlSeg f x xs seg = concatMap (P.init . P.scanl f x) (split seg xs) + where + split [] _ = [] + split _ [] = [] + split (i:is) vs = + let (h,t) = splitAt (P.fromIntegral i) vs + in h : split is t + + +-- Main +-- ---- +run :: String -> Int -> IO (() -> UArray Int Float, () -> Acc (Vector Float)) +run alg m = withSystemRandom $ \gen -> do + -- generate segments + -- + let n = P.round . sqrt $ (P.fromIntegral m :: Double) + seg <- randomUArrayR (0,n) gen (P.fromIntegral n) + seg' <- convertUArray seg + + -- generate elements + -- + let ne = P.fromIntegral $ 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/examples/tests/primitives/Scatter.hs b/examples/tests/primitives/Scatter.hs new file mode 100644 --- /dev/null +++ b/examples/tests/primitives/Scatter.hs @@ -0,0 +1,107 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE FlexibleContexts #-} + +module Scatter where + +import Random + +import Data.List +import Data.Maybe +import Data.Array.ST +import Control.Monad +import Control.Applicative +import System.Random.MWC +import Data.Array.Unboxed + +import Prelude as P +import Data.Array.Accelerate as Acc hiding ((!)) +import qualified Data.Array.MArray as M +import qualified Data.HashTable.IO as Hash + + +-- Tests +-- ----- + +scatterAcc :: Vector Int -> Vector Float -> Vector Float -> Acc (Vector Float) +scatterAcc mapV defaultV inputV = Acc.scatter (use mapV) (use defaultV) (use inputV) + +scatterIfAcc :: Vector Int -> Vector Int -> Vector Float -> Vector Float -> Acc (Vector Float) +scatterIfAcc mapV maskV defaultV inputV + = Acc.scatterIf (use mapV) (use maskV) evenAcc (use defaultV) (use inputV) + +evenAcc :: Exp Int -> Exp Bool +evenAcc v = (v `mod` 2) ==* 0 + + +scatterRef :: UArray Int Int -> UArray Int Float -> UArray Int Float -> UArray Int Float +scatterRef mapV defaultV inputV = runSTUArray $ do + mu <- M.thaw defaultV + forM_ (P.zip [0..] $ elems mapV) $ \(inIx, outIx) -> do + writeArray mu outIx (inputV ! inIx) + return mu + +scatterIfRef :: UArray Int Int -> UArray Int Int -> UArray Int Float -> UArray Int Float -> UArray Int Float +scatterIfRef mapV maskV defaultV inputV = runSTUArray $ do + mu <- M.thaw defaultV + forM_ (P.zip [0..] $ elems mapV) $ \(inIx, outIx) -> do + when (evenRef (maskV ! inIx)) $ do + writeArray mu outIx (inputV ! inIx) + return mu + +evenRef :: Int -> Bool +evenRef = even + + +-- Random +-- ------ + +uniqueRandomUArrayR :: GenIO -> (Int,Int) -> Int -> IO (UArray Int Int) +uniqueRandomUArrayR gen lim n = do + set <- Hash.new :: IO (Hash.BasicHashTable Int ()) + + let go !i !m | i >= n = return m + | otherwise = do + v <- uniformR lim gen + exists <- isJust <$> Hash.lookup set v + if exists + then go (i+1) m + else Hash.insert set v () >> go (i+1) (m+1) + + n' <- go 0 0 + listArray (0, n'-1) . P.map P.fst <$> Hash.toList set + + +-- Main +-- ---- +run :: String -> Int -> IO (() -> UArray Int Float, () -> Acc (Vector Float)) +run alg n = withSystemRandom $ \gen -> do + let m = 2 * n + + mapV <- uniqueRandomUArrayR gen (0, m-1) n + mapV' <- convertUArray mapV + let n' = rangeSize (bounds mapV) + + vec <- randomUArrayR (-1, 1) gen n' + vec' <- convertUArray vec + + maskV <- randomUArrayR (0, n') gen n' + maskV' <- convertUArray maskV + + defaultV <- randomUArrayR (-1, 1) gen m + defaultV' <- convertUArray defaultV + + -- + let go f g = return (run_ref f vec, run_acc g vec') + + case alg of + "scatter" -> go (scatterRef mapV defaultV) (scatterAcc mapV' defaultV') + "scatter-if" -> go (scatterIfRef mapV maskV defaultV) (scatterIfAcc mapV' maskV' defaultV') + x -> error $ "unknown variant: " ++ x + + where + {-# NOINLINE run_ref #-} + run_ref f xs () = f xs + run_acc f xs () = f xs + + + diff --git a/examples/tests/primitives/Stencil.hs b/examples/tests/primitives/Stencil.hs new file mode 100644 --- /dev/null +++ b/examples/tests/primitives/Stencil.hs @@ -0,0 +1,229 @@ +{-# LANGUAGE FlexibleContexts #-} + +module Stencil where + +import Util +import Random + +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/examples/tests/primitives/Stencil2.hs b/examples/tests/primitives/Stencil2.hs new file mode 100644 --- /dev/null +++ b/examples/tests/primitives/Stencil2.hs @@ -0,0 +1,109 @@ +{-# LANGUAGE FlexibleContexts #-} + +module Stencil2 where + +import Util + +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.Accelerate as A +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) + + +stencil2D2Ref + :: (Floating a, IArray UArray a) + => UArray (Int,Int) a + -> UArray (Int,Int) a + -> UArray (Int,Int) a +stencil2D2Ref xs ys = array sh [(ix, f ix) | ix <- range sh] + where + (_,(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) + get0 (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-1, iy ) + b = get1 (ix+1, iy ) + x = get1 (ix, iy ) + l = get0 (ix, iy-1) + r = get0 (ix, iy+1) + y = get0 (ix, iy ) + in + 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 = stencil2D2Ref xs ys + + +varUse :: (Acc (Array DIM2 Int), Acc (Array DIM2 Float), Acc (Array DIM2 Float)) +varUse = (first, both, second) + where + is :: Array DIM2 Int + is = fromList (Z:.10:.10) [0..] + + fs :: Array DIM2 Float + fs = fromList (Z:.10:.10) [0..] + + -- Ignoring the first parameter + first = stencil2 centre Clamp (use fs) Clamp (use is) + where + centre :: Stencil3x3 Float -> Stencil3x3 Int -> Exp Int + centre _ (_,(_,y,_),_) = y + + -- Using both + both = stencil2 centre Clamp (use fs) Clamp (use is) + where + centre :: Stencil3x3 Float -> Stencil3x3 Int -> Exp Float + centre (_,(_,x,_),_) (_,(_,y,_),_) = x + A.fromIntegral y + + -- Not using the second parameter + second = stencil2 centre Clamp (use fs) Clamp (use is) + where + centre :: Stencil3x3 Float -> Stencil3x3 Int -> Exp Float + centre (_,(_,x,_),_) _ = x + + +-- 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/examples/tests/primitives/Vector.hs b/examples/tests/primitives/Vector.hs new file mode 100644 --- /dev/null +++ b/examples/tests/primitives/Vector.hs @@ -0,0 +1,67 @@ +{-# LANGUAGE FlexibleContexts #-} + +module Vector where + +import Random + +import System.Random.MWC +import Data.Array.Unboxed +import Data.Array.Accelerate as Acc + + +-- Tests +-- ----- + +initAcc, tailAcc :: Vector Float -> Acc (Vector Float) +initAcc = Acc.init . Acc.use +tailAcc = Acc.tail . Acc.use + +takeAcc, dropAcc :: Int -> Vector Float -> Acc (Vector Float) +takeAcc n = (Acc.take $ constant n) . Acc.use +dropAcc n = (Acc.drop $ constant n) . Acc.use + +slitAcc :: Int -> Int -> Vector Float -> Acc (Vector Float) +slitAcc i n = (Acc.slit (constant i) (constant n)) . 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) + +initRef, tailRef :: UArray Int Float -> UArray Int Float +initRef = toUA $ Prelude.init +tailRef = toUA $ Prelude.tail + +takeRef, dropRef :: Int -> UArray Int Float -> UArray Int Float +takeRef n = toUA $ Prelude.take n +dropRef n = toUA $ Prelude.drop n + +slitRef :: Int -> Int -> UArray Int Float -> UArray Int Float +slitRef i n = toUA $ (Prelude.take n . Prelude.drop i) + + +-- 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 + ri0 <- uniform gen -- ri = random int + ri1 <- uniform gen + + -- + let go f g = return (run_ref f vec, run_acc g vec') + m = (abs ri0) `mod` n + len = (abs ri1) `mod` (n - m) + + case alg of + "init" -> go initRef initAcc + "tail" -> go tailRef tailAcc + "take" -> go (takeRef m) (takeAcc m) + "drop" -> go (dropRef m) (dropAcc m) + "slit" -> go (slitRef m len) (slitAcc m len) + x -> error $ "unknown variant: " ++ x + + where + {-# NOINLINE run_ref #-} + run_ref f xs () = f xs + run_acc f xs () = f xs + diff --git a/examples/tests/primitives/Zip.hs b/examples/tests/primitives/Zip.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/tests/primitives/ZipWith.hs b/examples/tests/primitives/ZipWith.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/tests/simple/BlackScholes.hs b/examples/tests/simple/BlackScholes.hs new file mode 100644 --- /dev/null +++ b/examples/tests/simple/BlackScholes.hs @@ -0,0 +1,88 @@ + +module BlackScholes where + +import Random + +import System.Random.MWC +import Data.Array.IArray as IArray +import Data.Array.Accelerate as Acc +import Prelude as P + + +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/examples/tests/simple/DotP.hs b/examples/tests/simple/DotP.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/tests/simple/Filter.hs b/examples/tests/simple/Filter.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/tests/simple/Radix.hs b/examples/tests/simple/Radix.hs new file mode 100644 --- /dev/null +++ b/examples/tests/simple/Radix.hs @@ -0,0 +1,90 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ScopedTypeVariables #-} +-- +-- Radix sort for a subclass of element types +-- + +module Radix where + +import Random + +import Prelude as P +import Data.Array.Accelerate as A + +import Data.Bits +import Data.List ( sort ) +import Data.Array.Unboxed ( IArray, UArray, listArray, bounds, elems ) +import System.Random.MWC + + +-- Radix sort +-- ---------- + +class Elt e => Radix e where + passes :: e {- dummy -} -> Int + radix :: Exp Int -> Exp e -> Exp Int + +instance Radix Int32 where + passes = bitSize + radix i e = i ==* (passes' - 1) ? (radix' (e `xor` minBound), radix' e) + where + radix' x = A.fromIntegral $ (x `A.shiftR` i) .&. 1 + passes' = constant (passes (undefined :: Int32)) + +-- For IEEE-754 floating-point representation. Unsafe, but widely supported. +-- TLM: unsafeCoerce does not work in the CUDA backend. +-- +-- 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 => Acc (Vector a) -> Acc (Vector a) +sortAcc = sortAccBy id + +sortAccBy + :: forall a r. (Elt a, Radix r) + => (Exp a -> Exp r) + -> Acc (Vector a) + -> Acc (Vector a) +sortAccBy rdx arr = foldr1 (>->) (P.map radixPass [0..p-1]) arr + where + p = passes (undefined :: r) + -- + deal f x = let (a,b) = unlift x in (f ==* 0) ? (a,b) + radixPass k v = let k' = unit (constant k) + flags = A.map (radix (the k') . rdx) v + idown = prescanl (+) 0 . A.map (xor 1) $ flags + iup = A.map (size v - 1 -) . prescanr (+) 0 $ flags + index = A.zipWith deal flags (A.zip idown iup) + in + permute const v (\ix -> index1 (index!ix)) v + + +sortRef :: (Ord a, IArray UArray a) => UArray Int a -> UArray Int a +sortRef xs = listArray (bounds xs) $ sort (elems xs) + + +-- Main +-- ---- + +run :: Int -> IO (() -> UArray Int Int32, () -> Acc (Vector Int32)) +run n = withSystemRandom $ \gen -> do + vec <- randomUArrayR (minBound,maxBound) gen n + vec' <- use `fmap` 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/examples/tests/simple/SASUM.hs b/examples/tests/simple/SASUM.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/tests/simple/SAXPY.hs b/examples/tests/simple/SAXPY.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/tests/simple/SMVM.hs b/examples/tests/simple/SMVM.hs new file mode 100644 --- /dev/null +++ b/examples/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 Int, 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/examples/tests/simple/SMVM/Matrix.hs b/examples/tests/simple/SMVM/Matrix.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/tests/simple/SMVM/MatrixMarket.hs b/examples/tests/simple/SMVM/MatrixMarket.hs new file mode 100644 --- /dev/null +++ b/examples/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 +import Data.ByteString.Lex.Double +import qualified Data.Attoparsec.Lazy as L +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` many1 (line floating) + Complex -> ComplexMatrix (m,n) l `fmap` many1 (line ((:+) <$> floating <*> floating)) + Integer -> IntMatrix (m,n) l `fmap` many1 (line integral) + Pattern -> PatternMatrix (m,n) l `fmap` many1 ((,) <$> integral <*> integral) + + +readMatrix :: FilePath -> IO Matrix +readMatrix file = do + chunks <- L.readFile file + case L.parse matrix chunks of + L.Fail _ _ msg -> error $ file ++ ": " ++ msg + L.Done _ mtx -> return mtx + diff --git a/examples/tests/simple/SharingRecovery.hs b/examples/tests/simple/SharingRecovery.hs new file mode 100644 --- /dev/null +++ b/examples/tests/simple/SharingRecovery.hs @@ -0,0 +1,145 @@ +{-# 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) + +-- | 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 + +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/examples/tests/simple/SliceExamples.hs b/examples/tests/simple/SliceExamples.hs new file mode 100644 --- /dev/null +++ b/examples/tests/simple/SliceExamples.hs @@ -0,0 +1,115 @@ +{-# 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) diff --git a/src/Benchmark.hs b/src/Benchmark.hs deleted file mode 100644 --- a/src/Benchmark.hs +++ /dev/null @@ -1,58 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/src/Config.hs +++ /dev/null @@ -1,135 +0,0 @@ -{-# 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:" - , " 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 deleted file mode 100644 --- a/src/Main.hs +++ /dev/null @@ -1,74 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/src/Random.hs +++ /dev/null @@ -1,72 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/src/Test.hs +++ /dev/null @@ -1,234 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/src/Util.hs +++ /dev/null @@ -1,17 +0,0 @@ - -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 deleted file mode 100644 --- a/src/Validate.hs +++ /dev/null @@ -1,102 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/image-processing/Canny.hs +++ /dev/null @@ -1,117 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/image-processing/IntegralImage.hs +++ /dev/null @@ -1,49 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/image-processing/PGM.hs +++ /dev/null @@ -1,32 +0,0 @@ --- --- 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 deleted file mode 100644 --- a/tests/io/BlockCopy.hs +++ /dev/null @@ -1,90 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/io/Vector.hs +++ /dev/null @@ -1,75 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/io/fill_with_values.cpp +++ /dev/null @@ -1,81 +0,0 @@ -#include -#include -#include - -/* Returns one if it's filled with even values starting at 0 */ -template -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 deleted file mode 100644 --- a/tests/primitives/Backpermute.hs +++ /dev/null @@ -1,69 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/primitives/Fold.hs +++ /dev/null @@ -1,84 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/primitives/Map.hs +++ /dev/null @@ -1,52 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/primitives/Permute.hs +++ /dev/null @@ -1,39 +0,0 @@ -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 deleted file mode 100644 --- a/tests/primitives/ScanSeg.hs +++ /dev/null @@ -1,58 +0,0 @@ - -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 deleted file mode 100644 --- a/tests/primitives/Stencil.hs +++ /dev/null @@ -1,229 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/primitives/Stencil2.hs +++ /dev/null @@ -1,70 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/primitives/Zip.hs +++ /dev/null @@ -1,35 +0,0 @@ - -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 deleted file mode 100644 --- a/tests/primitives/ZipWith.hs +++ /dev/null @@ -1,42 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/simple/BlackScholes.hs +++ /dev/null @@ -1,88 +0,0 @@ - -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 deleted file mode 100644 --- a/tests/simple/DotP.hs +++ /dev/null @@ -1,44 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/simple/Filter.hs +++ /dev/null @@ -1,55 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/simple/Radix.hs +++ /dev/null @@ -1,86 +0,0 @@ --- --- 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 deleted file mode 100644 --- a/tests/simple/SASUM.hs +++ /dev/null @@ -1,35 +0,0 @@ - -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 deleted file mode 100644 --- a/tests/simple/SAXPY.hs +++ /dev/null @@ -1,42 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/simple/SMVM.hs +++ /dev/null @@ -1,78 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/simple/SMVM/Matrix.hs +++ /dev/null @@ -1,75 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/simple/SMVM/MatrixMarket.hs +++ /dev/null @@ -1,133 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/simple/SharingRecovery.hs +++ /dev/null @@ -1,177 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/tests/simple/SliceExamples.hs +++ /dev/null @@ -1,114 +0,0 @@ -{-# 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)