diff --git a/Data/Array/Repa.hs b/Data/Array/Repa.hs --- a/Data/Array/Repa.hs +++ b/Data/Array/Repa.hs @@ -1,209 +1,147 @@ -{-# LANGUAGE PatternGuards, PackageImports, ScopedTypeVariables, RankNTypes #-} -{-# LANGUAGE TypeOperators, FlexibleContexts, NoMonomorphismRestriction, FlexibleInstances, UndecidableInstances #-} -{-# OPTIONS -fno-warn-orphans #-} --- | See the repa-examples package for examples. +-- | Repa arrays are wrappers around a linear structure that holds the element +-- data. The representation tag determines what structure holds the data. -- --- More information at . +-- Delayed Representations (functions that compute elements) -- --- There is a draft tutorial at +-- * `D` -- Functions from indices to elements. -- --- @Release Notes: --- For 2.2.0.1: --- * Added unsafeFromForeignPtr, which helps use foreign source --- arrays without intermediate copying. --- * Added forceWith and forceWith2, which can be used to force --- arrays into foreign result buffers without intermediate copying. +-- * `C` -- Cursor functions. -- --- For 2.1.0.1: --- * The fold and foldAll functions now run in parallel and require the --- starting element to be neutral with respect to the reduction operator. --- -- thanks to Trevor McDonell --- * Added (\/\/) update function. -- thanks to Trevor McDonell --- * Dropped unneeded Elt constraints from traverse functions. --- @ +-- Manifest Representations (real data) +-- +-- * `U` -- Adaptive unboxed vectors. +-- +-- * `V` -- Boxed vectors. +-- +-- * `B` -- Strict ByteStrings. +-- +-- * `F` -- Foreign memory buffers. +-- +-- Meta Representations +-- +-- * `P` -- Arrays that are partitioned into several representations. +-- +-- * `X` -- Arrays whose elements are all undefined. +-- +-- Array fusion is achieved via the delayed (`D`) and cursored (`C`) +-- representations. At compile time, the GHC simplifier combines the functions +-- contained within `D` and `C` arrays without needing to create manifest +-- intermediate arrays. +-- +-- Converting between the parallel manifest representations (eg `U` and `B`) +-- is either constant time or parallel copy, depending on the compatability +-- of the physical representation. +-- +-- /Writing fast code:/ +-- +-- 1. Repa does not support nested parallellism. +-- This means that you cannot `map` a parallel worker function across +-- an array and then call `computeP` to evaluate it, or pass a parallel +-- worker to parallel reductions such as `foldP`. If you do then you will +-- get a run-time warning and the code will run very slowly. +-- +-- 2. Arrays of type @(Array D sh a)@ or @(Array C sh a)@ are /not real arrays/. +-- They are represented as functions that compute each element on demand. +-- You need to use a function like `computeS`, `computeP`, `computeUnboxedP` +-- and so on to actually evaluate the elements. +-- +-- 3. You should add @INLINE@ pragmas to all leaf-functions in your code, +-- expecially ones that compute numberic results. This ensures they are +-- specialised at the appropriate element types. +-- +-- 4. Scheduling a parallel computation takes about 200us on an OSX machine. +-- You should sequential computation for small arrays in inner loops, +-- or a the bottom of a divide-and-conquer algorithm. +-- module Data.Array.Repa - ( module Data.Array.Repa.Shape - , module Data.Array.Repa.Index - , module Data.Array.Repa.Slice - - -- from Data.Array.Repa.Internals.Elt ----------------------- - , Elt(..) - - -- from Data.Array.Repa.Internals.Base ---------------------- - , Array(..) - , Region(..) - , Range(..) - , Rect(..) - , Generator(..) - , deepSeqArray, deepSeqArrays - , singleton, toScalar - , extent, delay - - -- - , withManifest, withManifest' - - -- * Indexing - , (!), index - , (!?), safeIndex - , unsafeIndex + ( -- * Abstract array representation + Array(..) + , module Data.Array.Repa.Shape + , module Data.Array.Repa.Index + , Repr(..), (!), toList + , deepSeqArrays - -- * Construction - , fromFunction - , fromVector - , fromList - , unsafeFromForeignPtr + -- * Converting between array representations + , computeP, computeS + , copyP, copyS + , now - -- from Data.Array.Repa.Interlals.Forcing ------------------- - -- * Forcing - , force, forceWith - , force2, forceWith2 - , toVector - , toList + -- * Concrete array representations + -- ** Delayed representation + , D, fromFunction, toFunction + , delay + -- ** Unboxed vector representation + , U + , computeUnboxedP, computeUnboxedS + , fromListUnboxed + , fromUnboxed + , toUnboxed + -- from Data.Array.Repa.Operators.IndexSpace ---------------- - -- * Index space transformations + -- * Operators + -- ** Index space transformations , reshape , append, (++) , transpose , extend - , slice - , backpermute + , backpermute, unsafeBackpermute , backpermuteDft + , module Data.Array.Repa.Slice + , slice + -- from Data.Array.Repa.Operators.Mapping ------------------- - -- * Structure preserving operations + -- ** Structure preserving operations , map , zipWith , (+^), (-^), (*^), (/^) - - -- from Data.Array.Repa.Operations.Modify ------------------- - -- * Bulk updates - , (//) - - -- from Data.Array.Repa.Operators.Reduction ----------------- - -- * Reductions - , fold, foldAll - , sum, sumAll - - -- from Data.Array.Repa.Operators.Traverse ------------------ - -- * Generic Traversal - , traverse - , traverse2 - , traverse3 - , traverse4 - , unsafeTraverse - , unsafeTraverse2 - , unsafeTraverse3 - , unsafeTraverse4 + , Combine(..) - -- from Data.Array.Repa.Operators.Interleave ---------------- - -- * Interleaving + -- from Data.Array.Repa.Operators.Traversal ------------------ + -- ** Generic traversal + , traverse, unsafeTraverse + , traverse2, unsafeTraverse2 + , traverse3, unsafeTraverse3 + , traverse4, unsafeTraverse4 + + -- from Data.Array.Repa.Operators.Interleave ----------------- + -- ** Interleaving , interleave2 , interleave3 , interleave4 - - -- from Data.Array.Repa.Operators.Select -------------------- - -- * Selection + + -- from Data.Array.Repa.Operators.Reduction ------------------ + -- ** Reduction + , foldP, foldS + , foldAllP, foldAllS + , sumP, sumS + , sumAllP, sumAllS + + -- from Data.Array.Repa.Operators.Selection ------------------ , select) - where +import Data.Array.Repa.Base +import Data.Array.Repa.Shape import Data.Array.Repa.Index import Data.Array.Repa.Slice -import Data.Array.Repa.Shape -import Data.Array.Repa.Internals.Elt -import Data.Array.Repa.Internals.Base -import Data.Array.Repa.Internals.Forcing -import Data.Array.Repa.Operators.Traverse +import Data.Array.Repa.Eval +import Data.Array.Repa.Repr.Delayed +import Data.Array.Repa.Repr.Vector +import Data.Array.Repa.Repr.Unboxed +import Data.Array.Repa.Repr.ByteString +import Data.Array.Repa.Repr.ForeignPtr +import Data.Array.Repa.Repr.Cursored +import Data.Array.Repa.Repr.Partitioned +import Data.Array.Repa.Repr.Undefined () +import Data.Array.Repa.Operators.Mapping +import Data.Array.Repa.Operators.Traversal import Data.Array.Repa.Operators.IndexSpace import Data.Array.Repa.Operators.Interleave -import Data.Array.Repa.Operators.Mapping -import Data.Array.Repa.Operators.Modify import Data.Array.Repa.Operators.Reduction -import Data.Array.Repa.Operators.Select -import qualified Data.Array.Repa.Shape as S - -import Prelude hiding (sum, map, zipWith, (++)) -import qualified Prelude as P - -stage = "Data.Array.Repa" - - --- Instances -------------------------------------------------------------------------------------- --- Show -instance (Shape sh, Elt a, Show a) => Show (Array sh a) where - show arr = - let shape = showShape (extent arr) - elems = show (toList arr) - in - "Array (" P.++ shape P.++ ") " P.++ elems - - --- Eq -instance (Shape sh, Elt a, Eq a) => Eq (Array sh a) where - - {-# INLINE (==) #-} - (==) arr1 arr2 - = foldAll (&&) True - $ reshape (Z :. (S.size $ extent arr1)) - $ zipWith (==) arr1 arr2 - - {-# INLINE (/=) #-} - (/=) a1 a2 - = not $ (==) a1 a2 - --- Num --- All operators apply elementwise. -instance (Shape sh, Elt a, Num a) => Num (Array sh a) where - {-# INLINE (+) #-} - (+) = zipWith (+) - - {-# INLINE (-) #-} - (-) = zipWith (-) - - {-# INLINE (*) #-} - (*) = zipWith (*) - - {-# INLINE negate #-} - negate = map negate - - {-# INLINE abs #-} - abs = map abs - - {-# INLINE signum #-} - signum = map signum - - {-# INLINE fromInteger #-} - fromInteger n = fromFunction failShape (\_ -> fromInteger n) - where failShape = error $ stage P.++ ".fromInteger: Constructed array has no shape." - - --- | Force an array before passing it to a function. -withManifest - :: (Shape sh, Elt a) - => (Array sh a -> b) -> Array sh a -> b - -{-# INLINE withManifest #-} -withManifest f arr - = case arr of - Array sh [Region RangeAll (GenManifest vec)] - -> vec `seq` f (Array sh [Region RangeAll (GenManifest vec)]) - - _ -> f (force arr) - - --- | Force an array before passing it to a function. -withManifest' - :: (Shape sh, Elt a) - => Array sh a -> (Array sh a -> b) -> b - -{-# INLINE withManifest' #-} -withManifest' arr f - = case arr of - Array sh [Region RangeAll (GenManifest vec)] - -> vec `seq` f (Array sh [Region RangeAll (GenManifest vec)]) - - _ -> f (force arr) +import Data.Array.Repa.Operators.Selection +import Prelude () diff --git a/Data/Array/Repa/Arbitrary.hs b/Data/Array/Repa/Arbitrary.hs deleted file mode 100644 --- a/Data/Array/Repa/Arbitrary.hs +++ /dev/null @@ -1,99 +0,0 @@ -{-# LANGUAGE TypeOperators, FlexibleInstances #-} -{-# OPTIONS -fno-warn-orphans #-} - --- Utils to help with testing. Not exported. -module Data.Array.Repa.Arbitrary - ( arbitraryShape - , arbitrarySmallShape - , arbitraryListOfLength - , arbitrarySmallArray) -where -import Data.Array.Repa.Internals.Elt -import Data.Array.Repa.Internals.Base -import Data.Array.Repa.Index -import Data.Array.Repa.Shape as S -import Control.Monad -import Test.QuickCheck - - --- Arbitrary -------------------------------------------------------------------------------------- -instance Arbitrary Z where - arbitrary = return Z - --- | Generate an arbitrary index, which may have 0's for some components. -instance (Shape sh, Arbitrary sh) => Arbitrary (sh :. Int) where - arbitrary - = do sh1 <- arbitrary - let sh1Unit = if size sh1 == 0 then unitDim else sh1 - - -- Make sure not to create an index so big that we get - -- integer overflow when converting it to the linear form. - n <- liftM abs $ arbitrary - let nMax = maxBound `div` (size sh1Unit) - let nMaxed = n `mod` nMax - - return $ sh1 :. nMaxed - --- | Generate an aribrary shape that does not have 0's for any component. -arbitraryShape - :: (Shape sh, Arbitrary sh) - => Gen (sh :. Int) - -arbitraryShape - = do sh1 <- arbitrary - let sh1Unit = if size sh1 == 0 then unitDim else sh1 - - -- Make sure not to create an index so big that we get - -- integer overflow when converting it to the linear form. - n <- liftM abs $ arbitrary - let nMax = maxBound `div` size sh1Unit - let nMaxed = n `mod` nMax - let nClamped = if nMaxed == 0 then 1 else nMaxed - - return $ sh1Unit :. nClamped - - --- | Generate an arbitrary shape where each dimension is more than zero, --- but less than a specific value. -arbitrarySmallShape - :: (Shape sh, Arbitrary sh) - => Int - -> Gen (sh :. Int) - -arbitrarySmallShape maxDim - = do sh <- arbitraryShape - let dims = listOfShape sh - - let clamp x - = case x `mod` maxDim of - 0 -> 1 - n -> n - - return $ if True - then shapeOfList $ map clamp dims - else sh - - -arbitraryListOfLength - :: Arbitrary a - => Int -> Gen [a] - -arbitraryListOfLength n - | n == 0 = return [] - | otherwise - = do i <- arbitrary - rest <- arbitraryListOfLength (n - 1) - return $ i : rest - --- | Create an arbitrary small array, restricting the size of each of the --- dimensions to some value. -arbitrarySmallArray - :: (Shape sh, Elt a, Arbitrary sh, Arbitrary a) - => Int - -> Gen (Array (sh :. Int) a) - -arbitrarySmallArray maxDim - = do sh <- arbitrarySmallShape maxDim - xx <- arbitraryListOfLength (S.size sh) - return $ fromList sh xx - diff --git a/Data/Array/Repa/Base.hs b/Data/Array/Repa/Base.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Base.hs @@ -0,0 +1,85 @@ + +module Data.Array.Repa.Base + ( Array + , Repr (..), (!), toList + , deepSeqArrays) +where +import Data.Array.Repa.Shape + +-- | Arrays with a representation tag, shape, and element type. +-- Use one of the type tags like `D`, `U` and so on for @r@, +-- one of `DIM1`, `DIM2` ... for @sh@. +data family Array r sh e + + +-- | Class of array representations that we can read elements from. +-- +class Repr r e where + -- | O(1). Take the extent of an array. + extent :: Shape sh => Array r sh e -> sh + + -- | O(1). Shape polymorphic indexing. + index, unsafeIndex + :: Shape sh => Array r sh e -> sh -> e + + {-# INLINE index #-} + index arr ix = arr `linearIndex` toIndex (extent arr) ix + + {-# INLINE unsafeIndex #-} + unsafeIndex arr ix = arr `unsafeLinearIndex` toIndex (extent arr) ix + + -- | O(1). Linear indexing into underlying, row-major, array representation. + linearIndex, unsafeLinearIndex + :: Shape sh => Array r sh e -> Int -> e + + {-# INLINE unsafeLinearIndex #-} + unsafeLinearIndex = linearIndex + + -- | Ensure an array's data structure is fully evaluated. + deepSeqArray :: Shape sh => Array r sh e -> b -> b + + +-- | O(1). Alias for `index` +(!) :: (Repr r e, Shape sh) => Array r sh e -> sh -> e +(!) = index + + +-- | O(n). Convert an array to a list. +toList :: (Shape sh, Repr r e) + => Array r sh e -> [e] +{-# INLINE toList #-} +toList arr + = go 0 + where len = size (extent arr) + go ix + | ix == len = [] + | otherwise = unsafeLinearIndex arr ix : go (ix + 1) + + +-- | Apply `deepSeqArray` to up to four arrays. +-- +-- The implementation of this function has been hand-unwound to work for up to +-- four arrays. Putting more in the list yields `error`. +-- +deepSeqArrays + :: (Shape sh, Repr r e) + => [Array r sh e] -> b -> b +{-# INLINE deepSeqArrays #-} +deepSeqArrays arrs x + = case arrs of + [] -> x + + [a1] + -> a1 `deepSeqArray` x + + [a1, a2] + -> a1 `deepSeqArray` a2 `deepSeqArray` x + + [a1, a2, a3] + -> a1 `deepSeqArray` a2 `deepSeqArray` a3 `deepSeqArray` x + + [a1, a2, a3, a4] + -> a1 `deepSeqArray` a2 `deepSeqArray` a3 `deepSeqArray` a4 `deepSeqArray` x + + _ -> error "deepSeqArrays: only works for up to four arrays" + diff --git a/Data/Array/Repa/Eval.hs b/Data/Array/Repa/Eval.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Eval.hs @@ -0,0 +1,130 @@ +{-# LANGUAGE UndecidableInstances #-} + +-- | Low level interface to parallel array filling operators. +module Data.Array.Repa.Eval + ( -- * Element types + Elt (..) + + -- * Parallel array filling + , Fillable (..) + , Fill (..) + , FillRange (..) + , fromList + + -- * Converting between representations + , computeP, computeS + , copyP, copyS + , now + + -- * Chunked filling + , fillChunkedS + , fillChunkedP + , fillChunkedIOP + + -- * Blockwise filling + , fillBlock2P + , fillBlock2S + + -- * Cursored blockwise filling + , fillCursoredBlock2S + , fillCursoredBlock2P + + -- * Chunked selection + , selectChunkedS + , selectChunkedP) +where +import Data.Array.Repa.Eval.Elt +import Data.Array.Repa.Eval.Fill +import Data.Array.Repa.Eval.Chunked +import Data.Array.Repa.Eval.Cursored +import Data.Array.Repa.Eval.Selection +import Data.Array.Repa.Repr.Delayed +import Data.Array.Repa.Base +import Data.Array.Repa.Shape +import System.IO.Unsafe + + +-- | Parallel computation of array elements. +-- +-- * The `Fill` class is defined so that the source array must have a +-- delayed representation (`D` or `C`) +-- +-- * If you want to copy data between manifest representations then use +-- `copyP` instead. +-- +-- * If you want to convert a manifest array back to a delayed representation +-- then use `delay` instead. +-- +computeP :: Fill r1 r2 sh e + => Array r1 sh e -> Array r2 sh e +{-# INLINE [4] computeP #-} +computeP arr1 + = arr1 `deepSeqArray` + unsafePerformIO + $ do marr2 <- newMArr (size $ extent arr1) + fillP arr1 marr2 + unsafeFreezeMArr (extent arr1) marr2 + + +-- | Sequential computation of array elements. +computeS + :: Fill r1 r2 sh e + => Array r1 sh e -> Array r2 sh e +{-# INLINE [4] computeS #-} +computeS arr1 + = arr1 `deepSeqArray` + unsafePerformIO + $ do marr2 <- newMArr (size $ extent arr1) + fillS arr1 marr2 + unsafeFreezeMArr (extent arr1) marr2 + + + + +-- | Parallel copying of arrays. +-- +-- * This is a wrapper that delays an array before calling `computeP`. +-- +-- * You can use it to copy manifest arrays between representations. +-- +-- * You can also use it to compute elements, but doing this may not be as +-- efficient. This is because delaying it the second time can hide +-- information about the structure of the original computation. +-- +copyP :: (Repr r1 e, Fill D r2 sh e) + => Array r1 sh e -> Array r2 sh e +{-# INLINE [4] copyP #-} +copyP arr1 = computeP $ delay arr1 + + +-- | Sequential copying of arrays. +copyS :: (Repr r1 e, Fill D r2 sh e) + => Array r1 sh e -> Array r2 sh e +{-# INLINE [4] copyS #-} +copyS arr1 = computeS $ delay arr1 + + + + +-- | Apply `deepSeqArray` to an array so the result is actually constructed +-- at this point in a monadic computation. +-- +-- * Haskell's laziness means that applications of `computeP` and `copyP` are +-- automatically suspended. +-- +-- * Laziness can be problematic for data parallel programs, because we want +-- each array to be constructed in parallel before moving onto the next one. +-- +-- For example: +-- +-- @ do arr2 <- now $ computeP $ map f arr1 +-- arr3 <- now $ computeP $ zipWith arr2 arr1 +-- return arr3 +-- @ +-- +now :: (Shape sh, Repr r e, Monad m) + => Array r sh e -> m (Array r sh e) +{-# INLINE [4] now #-} +now arr + = do arr `deepSeqArray` return () + return arr diff --git a/Data/Array/Repa/Eval/Chunked.hs b/Data/Array/Repa/Eval/Chunked.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Eval/Chunked.hs @@ -0,0 +1,138 @@ +{-# LANGUAGE MagicHash #-} +-- | Evaluate an array by breaking it up into linear chunks and filling +-- each chunk in parallel. +module Data.Array.Repa.Eval.Chunked + ( fillChunkedP + , fillChunkedS + , fillChunkedS' + , fillChunkedIOP) +where +import Data.Array.Repa.Eval.Gang +import GHC.Exts +import Prelude as P + +-- | Fill something sequentially. +-- +-- * The array is filled linearly from start to finish. +-- +fillChunkedS + :: Int -- ^ Number of elements. + -> (Int -> a -> IO ()) -- ^ Update function to write into result buffer. + -> (Int -> a) -- ^ Fn to get the value at a given index. + -> IO () + +{-# INLINE [0] fillChunkedS #-} +fillChunkedS !(I# len) !write !getElem + = fill 0# + where fill !ix + | ix >=# len = return () + | otherwise + = do write (I# ix) (getElem (I# ix)) + fill (ix +# 1#) + +fillChunkedS' + :: Int + -> (Int -> IO ()) + -> IO () + +fillChunkedS' !(I# len) eat + = fill 0# + where fill !ix + | ix >=# len = return () + | otherwise + = do eat (I# ix) + fill (ix +# 1#) + + + + +-- | Fill something in parallel. +-- +-- * The array is split into linear chunks and each thread fills one chunk. +-- +fillChunkedP + :: Int -- ^ Number of elements. + -> (Int -> a -> IO ()) -- ^ Update function to write into result buffer. + -> (Int -> a) -- ^ Fn to get the value at a given index. + -> IO () + +{-# INLINE [0] fillChunkedP #-} +fillChunkedP !(I# len) !write !getElem + = gangIO theGang + $ \(I# thread) -> + let !start = splitIx thread + !end = splitIx (thread +# 1#) + in fill start end + + where + -- Decide now to split the work across the threads. + -- If the length of the vector doesn't divide evenly among the threads, + -- then the first few get an extra element. + !(I# threads) = gangSize theGang + !chunkLen = len `quotInt#` threads + !chunkLeftover = len `remInt#` threads + + {-# INLINE splitIx #-} + splitIx thread + | thread <# chunkLeftover = thread *# (chunkLen +# 1#) + | otherwise = thread *# chunkLen +# chunkLeftover + + -- Evaluate the elements of a single chunk. + {-# INLINE fill #-} + fill !ix !end + | ix >=# end = return () + | otherwise + = do write (I# ix) (getElem (I# ix)) + fill (ix +# 1#) end + + +-- | Fill something in parallel, using a separate IO action for each thread. +fillChunkedIOP + :: Int -- ^ Number of elements. + -> (Int -> a -> IO ()) -- ^ Update fn to write into result buffer. + -> (Int -> IO (Int -> IO a)) -- ^ Create a fn to get the value at a given index. + -- The first `Int` is the thread number, so you can do some + -- per-thread initialisation. + -> IO () + +{-# INLINE [0] fillChunkedIOP #-} +fillChunkedIOP !(I# len) !write !mkGetElem + = gangIO theGang + $ \(I# thread) -> + let !start = splitIx thread + !end = splitIx (thread +# 1#) + in fillChunk thread start end + + where + -- Decide now to split the work across the threads. + -- If the length of the vector doesn't divide evenly among the threads, + -- then the first few get an extra element. + !(I# threads) = gangSize theGang + !chunkLen = len `quotInt#` threads + !chunkLeftover = len `remInt#` threads + + {-# INLINE splitIx #-} + splitIx thread + | thread <# chunkLeftover = thread *# (chunkLen +# 1#) + | otherwise = thread *# chunkLen +# chunkLeftover + + + -- Given the threadId, starting and ending indices. + -- Make a function to get each element for this chunk + -- and call it for every index. + {-# INLINE fillChunk #-} + fillChunk !thread !ixStart !ixEnd + = do getElem <- mkGetElem (I# thread) + fill getElem ixStart ixEnd + + -- Call the provided getElem function for every element + -- in a chunk, and feed the result to the write function. + {-# INLINE fill #-} + fill !getElem !ix0 !end + = go ix0 + where go !ix + | ix >=# end = return () + | otherwise + = do x <- getElem (I# ix) + write (I# ix) x + go (ix +# 1#) diff --git a/Data/Array/Repa/Eval/Cursored.hs b/Data/Array/Repa/Eval/Cursored.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Eval/Cursored.hs @@ -0,0 +1,217 @@ +{-# LANGUAGE MagicHash #-} +-- | Evaluate an array by dividing it into rectangular blocks and filling +-- each block in parallel. +module Data.Array.Repa.Eval.Cursored + ( fillBlock2P + , fillBlock2S + , fillCursoredBlock2P + , fillCursoredBlock2S ) +where +import Data.Array.Repa.Index +import Data.Array.Repa.Shape +import Data.Array.Repa.Eval.Elt +import Data.Array.Repa.Eval.Gang +import GHC.Base (remInt, quotInt) +import Prelude as P +import GHC.Exts + +-- Non-cursored interface ----------------------------------------------------- +-- | Fill a block in a rank-2 array in parallel. +-- +-- * Blockwise filling can be more cache-efficient than linear filling for +-- rank-2 arrays. +-- +-- * Coordinates given are of the filled edges of the block. +-- +-- * We divide the block into columns, and give one column to each thread. +-- +-- * Each column is filled in row major order from top to bottom. +-- +fillBlock2P + :: Elt a + => (Int -> a -> IO ()) -- ^ Update function to write into result buffer. + -> (DIM2 -> a) -- ^ Function to evaluate the element at an index. + -> Int -- ^ Width of the whole array. + -> Int -- ^ x0 lower left corner of block to fill + -> Int -- ^ y0 + -> Int -- ^ x1 upper right corner of block to fill + -> Int -- ^ y1 + -> IO () + +{-# INLINE [0] fillBlock2P #-} +fillBlock2P !write !getElem !imageWidth !x0 !y0 !x1 !y1 + = fillCursoredBlock2P + write id addDim getElem + imageWidth x0 y0 x1 y1 + + +-- | Fill a block in a rank-2 array sequentially. +-- +-- * Blockwise filling can be more cache-efficient than linear filling for +-- rank-2 arrays. +-- +-- * Coordinates given are of the filled edges of the block. +-- +-- * The block is filled in row major order from top to bottom. +-- +fillBlock2S + :: Elt a + => (Int -> a -> IO ()) -- ^ Update function to write into result buffer. + -> (DIM2 -> a) -- ^ Function to evaluate the element at an index. + -> Int# -- ^ Width of the whole array. + -> Int# -- ^ x0 lower left corner of block to fill + -> Int# -- ^ y0 + -> Int# -- ^ x1 upper right corner of block to fill + -> Int# -- ^ y1 + -> IO () + +{-# INLINE [0] fillBlock2S #-} +fillBlock2S !write !getElem imageWidth x0 y0 x1 y1 + = fillCursoredBlock2S + write id addDim getElem + imageWidth x0 y0 x1 y1 + + +-- Block filling ---------------------------------------------------------------------------------- +-- | Fill a block in a rank-2 array in parallel. +-- +-- * Blockwise filling can be more cache-efficient than linear filling for rank-2 arrays. +-- +-- * Using cursor functions can help to expose inter-element indexing computations to +-- the GHC and LLVM optimisers. +-- +-- * Coordinates given are of the filled edges of the block. +-- +-- * We divide the block into columns, and give one column to each thread. +-- +-- * Each column is filled in row major order from top to bottom. +-- +fillCursoredBlock2P + :: Elt a + => (Int -> a -> IO ()) -- ^ Update function to write into result buffer. + -> (DIM2 -> cursor) -- ^ Make a cursor to a particular element. + -> (DIM2 -> cursor -> cursor) -- ^ Shift the cursor by an offset. + -> (cursor -> a) -- ^ Function to evaluate the element at an index. + -> Int -- ^ Width of the whole array. + -> Int -- ^ x0 lower left corner of block to fill + -> Int -- ^ y0 + -> Int -- ^ x1 upper right corner of block to fill + -> Int -- ^ y1 + -> IO () + +{-# INLINE [0] fillCursoredBlock2P #-} +fillCursoredBlock2P + !write + !makeCursorFCB !shiftCursorFCB !getElemFCB + !(I# imageWidth) !x0 !y0 !x1 !y1 + = gangIO theGang fillBlock + where !threads = gangSize theGang + !blockWidth = x1 - x0 + 1 + + -- All columns have at least this many pixels. + !colChunkLen = blockWidth `quotInt` threads + + -- Extra pixels that we have to divide between some of the threads. + !colChunkSlack = blockWidth `remInt` threads + + -- Get the starting pixel of a column in the image. + {-# INLINE colIx #-} + colIx !ix + | ix < colChunkSlack = x0 + ix * (colChunkLen + 1) + | otherwise = x0 + ix * colChunkLen + colChunkSlack + + -- Give one column to each thread + {-# INLINE fillBlock #-} + fillBlock :: Int -> IO () + fillBlock !ix + = let !(I# x0') = colIx ix + !(I# x1') = colIx (ix + 1) - 1 + !(I# y0') = y0 + !(I# y1') = y1 + in fillCursoredBlock2S + write + makeCursorFCB shiftCursorFCB getElemFCB + imageWidth x0' y0' x1' y1' + + +-- | Fill a block in a rank-2 array, sequentially. +-- +-- * Blockwise filling can be more cache-efficient than linear filling for rank-2 arrays. +-- +-- * Using cursor functions can help to expose inter-element indexing computations to +-- the GHC and LLVM optimisers. +-- +-- * Coordinates given are of the filled edges of the block. +-- +-- * The block is filled in row major order from top to bottom. +-- +fillCursoredBlock2S + :: Elt a + => (Int -> a -> IO ()) -- ^ Update function to write into result buffer. + -> (DIM2 -> cursor) -- ^ Make a cursor to a particular element. + -> (DIM2 -> cursor -> cursor) -- ^ Shift the cursor by an offset. + -> (cursor -> a) -- ^ Function to evaluate an element at the given index. + -> Int# -- ^ Width of the whole array. + -> Int# -- ^ x0 lower left corner of block to fill. + -> Int# -- ^ y0 + -> Int# -- ^ x1 upper right corner of block to fill. + -> Int# -- ^ y1 + -> IO () + +{-# INLINE [0] fillCursoredBlock2S #-} +fillCursoredBlock2S + !write + !makeCursor !shiftCursor !getElem + !imageWidth !x0 !y0 !x1 !y1 + + = fillBlock y0 + + where {-# INLINE fillBlock #-} + fillBlock !y + | y ># y1 = return () + | otherwise + = do fillLine4 x0 + fillBlock (y +# 1#) + + where {-# INLINE fillLine4 #-} + fillLine4 !x + | x +# 4# ># x1 = fillLine1 x + | otherwise + = do -- Compute each source cursor based on the previous one so that + -- the variable live ranges in the generated code are shorter. + let srcCur0 = makeCursor (Z :. (I# y) :. (I# x)) + let srcCur1 = shiftCursor (Z :. 0 :. 1) srcCur0 + let srcCur2 = shiftCursor (Z :. 0 :. 1) srcCur1 + let srcCur3 = shiftCursor (Z :. 0 :. 1) srcCur2 + + -- Get the result value for each cursor. + let val0 = getElem srcCur0 + let val1 = getElem srcCur1 + let val2 = getElem srcCur2 + let val3 = getElem srcCur3 + + -- Ensure that we've computed each of the result values before we + -- write into the array. If the backend code generator can't tell + -- our destination array doesn't alias with the source then writing + -- to it can prevent the sharing of intermediate computations. + touch val0 + touch val1 + touch val2 + touch val3 + + -- Compute cursor into destination array. + let !dstCur0 = x +# (y *# imageWidth) + write (I# dstCur0) val0 + write (I# (dstCur0 +# 1#)) val1 + write (I# (dstCur0 +# 2#)) val2 + write (I# (dstCur0 +# 3#)) val3 + fillLine4 (x +# 4#) + + {-# INLINE fillLine1 #-} + fillLine1 !x + | x ># x1 = return () + | otherwise + = do write (I# (x +# (y *# imageWidth))) + (getElem $ makeCursor (Z :. (I# y) :. (I# x))) + fillLine1 (x +# 1#) + diff --git a/Data/Array/Repa/Eval/Elt.hs b/Data/Array/Repa/Eval/Elt.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Eval/Elt.hs @@ -0,0 +1,282 @@ +-- | Values that can be stored in Repa Arrays. +{-# LANGUAGE MagicHash, UnboxedTuples, TypeSynonymInstances, FlexibleInstances #-} +module Data.Array.Repa.Eval.Elt + (Elt (..)) +where +import GHC.Prim +import GHC.Exts +import GHC.Types +import GHC.Word +import GHC.Int + + +-- Note that the touch# function is special because we can pass it boxed or unboxed +-- values. The argument type has kind ?, not just * or #. + +-- | Element types that can be used with the blockwise filling functions. +-- +-- This class is mainly used to define the `touch` method. This is used internally +-- in the imeplementation of Repa to prevent let-binding from being floated +-- inappropriately by the GHC simplifier. Doing a `seq` sometimes isn't enough, +-- because the GHC simplifier can erase these, and still move around the bindings. +-- +class Elt a where + + -- | Place a demand on a value at a particular point in an IO computation. + touch :: a -> IO () + + -- | Generic zero value, helpful for debugging. + zero :: a + + -- | Generic one value, helpful for debugging. + one :: a + + +-- Bool ----------------------------------------------------------------------- +instance Elt Bool where + {-# INLINE touch #-} + touch b + = IO (\state -> case touch# b state of + state' -> (# state', () #)) + + {-# INLINE zero #-} + zero = False + + {-# INLINE one #-} + one = True + + +-- Floating ------------------------------------------------------------------- +instance Elt Float where + {-# INLINE touch #-} + touch (F# f) + = IO (\state -> case touch# f state of + state' -> (# state', () #)) + + {-# INLINE zero #-} + zero = 0 + + {-# INLINE one #-} + one = 1 + + +instance Elt Double where + {-# INLINE touch #-} + touch (D# d) + = IO (\state -> case touch# d state of + state' -> (# state', () #)) + + {-# INLINE zero #-} + zero = 0 + + {-# INLINE one #-} + one = 1 + + +-- Int ------------------------------------------------------------------------ +instance Elt Int where + {-# INLINE touch #-} + touch (I# i) + = IO (\state -> case touch# i state of + state' -> (# state', () #)) + + {-# INLINE zero #-} + zero = 0 + + {-# INLINE one #-} + one = 1 + +instance Elt Int8 where + {-# INLINE touch #-} + touch (I8# w) + = IO (\state -> case touch# w state of + state' -> (# state', () #)) + + {-# INLINE zero #-} + zero = 0 + + {-# INLINE one #-} + one = 1 + + +instance Elt Int16 where + {-# INLINE touch #-} + touch (I16# w) + = IO (\state -> case touch# w state of + state' -> (# state', () #)) + + {-# INLINE zero #-} + zero = 0 + + {-# INLINE one #-} + one = 1 + + +instance Elt Int32 where + {-# INLINE touch #-} + touch (I32# w) + = IO (\state -> case touch# w state of + state' -> (# state', () #)) + + {-# INLINE zero #-} + zero = 0 + + {-# INLINE one #-} + one = 1 + + +instance Elt Int64 where + {-# INLINE touch #-} + touch (I64# w) + = IO (\state -> case touch# w state of + state' -> (# state', () #)) + + {-# INLINE zero #-} + zero = 0 + + {-# INLINE one #-} + one = 1 + + +-- Word ----------------------------------------------------------------------- +instance Elt Word where + {-# INLINE touch #-} + touch (W# i) + = IO (\state -> case touch# i state of + state' -> (# state', () #)) + + {-# INLINE zero #-} + zero = 0 + + {-# INLINE one #-} + one = 1 + + +instance Elt Word8 where + {-# INLINE touch #-} + touch (W8# w) + = IO (\state -> case touch# w state of + state' -> (# state', () #)) + + {-# INLINE zero #-} + zero = 0 + + {-# INLINE one #-} + one = 1 + + +instance Elt Word16 where + {-# INLINE touch #-} + touch (W16# w) + = IO (\state -> case touch# w state of + state' -> (# state', () #)) + + {-# INLINE zero #-} + zero = 0 + + {-# INLINE one #-} + one = 1 + + +instance Elt Word32 where + {-# INLINE touch #-} + touch (W32# w) + = IO (\state -> case touch# w state of + state' -> (# state', () #)) + + {-# INLINE zero #-} + zero = 0 + + {-# INLINE one #-} + one = 1 + + +instance Elt Word64 where + {-# INLINE touch #-} + touch (W64# w) + = IO (\state -> case touch# w state of + state' -> (# state', () #)) + + {-# INLINE zero #-} + zero = 0 + + {-# INLINE one #-} + one = 1 + + +-- Tuple ---------------------------------------------------------------------- +instance (Elt a, Elt b) => Elt (a, b) where + {-# INLINE touch #-} + touch (a, b) + = do touch a + touch b + + {-# INLINE zero #-} + zero = (zero, zero) + + {-# INLINE one #-} + one = (one, one) + + +instance (Elt a, Elt b, Elt c) => Elt (a, b, c) where + {-# INLINE touch #-} + touch (a, b, c) + = do touch a + touch b + touch c + + {-# INLINE zero #-} + zero = (zero, zero, zero) + + {-# INLINE one #-} + one = (one, one, one) + + +instance (Elt a, Elt b, Elt c, Elt d) => Elt (a, b, c, d) where + {-# INLINE touch #-} + touch (a, b, c, d) + = do touch a + touch b + touch c + touch d + + {-# INLINE zero #-} + zero = (zero, zero, zero, zero) + + {-# INLINE one #-} + one = (one, one, one, one) + + +instance (Elt a, Elt b, Elt c, Elt d, Elt e) => Elt (a, b, c, d, e) where + {-# INLINE touch #-} + touch (a, b, c, d, e) + = do touch a + touch b + touch c + touch d + touch e + + {-# INLINE zero #-} + zero = (zero, zero, zero, zero, zero) + + {-# INLINE one #-} + one = (one, one, one, one, one) + + +instance (Elt a, Elt b, Elt c, Elt d, Elt e, Elt f) => Elt (a, b, c, d, e, f) where + {-# INLINE touch #-} + touch (a, b, c, d, e, f) + = do touch a + touch b + touch c + touch d + touch e + touch f + + {-# INLINE zero #-} + zero = (zero, zero, zero, zero, zero, zero) + + {-# INLINE one #-} + one = (one, one, one, one, one, one) + + diff --git a/Data/Array/Repa/Eval/Fill.hs b/Data/Array/Repa/Eval/Fill.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Eval/Fill.hs @@ -0,0 +1,71 @@ + +module Data.Array.Repa.Eval.Fill + ( Fillable (..), fromList + , Fill (..) + , FillRange (..)) +where +import Data.Array.Repa.Base +import Data.Array.Repa.Shape +import Control.Monad +import System.IO.Unsafe + +-- Fillable ------------------------------------------------------------------- +-- | Class of manifest array representations that can be filled in parallel +-- and then frozen into immutable Repa arrays. +class Fillable r e where + + -- | Mutable version of the representation. + data MArr r e + + -- | Allocate a new mutable array of the given size. + newMArr :: Int -> IO (MArr r e) + + -- | Write an element into the mutable array. + unsafeWriteMArr :: MArr r e -> Int -> e -> IO () + + -- | Freeze the mutable array into an immutable Repa array. + unsafeFreezeMArr :: sh -> MArr r e -> IO (Array r sh e) + + +-- | O(n). Construct a manifest array from a list. +fromList + :: (Shape sh, Fillable r e) + => sh -> [e] -> Array r sh e +fromList sh xx + = unsafePerformIO + $ do let len = length xx + if len /= size sh + then error "Data.Array.Repa.Eval.Fill.fromList: provide array shape does not match list length" + else do + marr <- newMArr len + zipWithM_ (unsafeWriteMArr marr) [0..] xx + unsafeFreezeMArr sh marr + + +-- Fill ----------------------------------------------------------------------- +-- | Compute all elements defined by an array and write them to a fillable +-- representation. +-- +-- Note that instances require that the source array to have a delayed +-- representation such as `D` or `C`. If you want to use a pre-existing +-- manifest array as the source then `delay` it first. +class (Shape sh, Repr r1 e, Fillable r2 e) => Fill r1 r2 sh e where + -- | Fill an entire array sequentially. + fillS :: Array r1 sh e -> MArr r2 e -> IO () + + -- | Fill an entire array in parallel. + fillP :: Array r1 sh e -> MArr r2 e -> IO () + + +-- FillRange ------------------------------------------------------------------ +-- | Compute a range of elements defined by an array and write them to a fillable +-- representation. +class (Shape sh, Repr r1 e, Fillable r2 e) => FillRange r1 r2 sh e where + -- | Fill a range of an array sequentially. + fillRangeS :: Array r1 sh e -> MArr r2 e -> sh -> sh -> IO () + + -- | Fill a range of an array in parallel. + fillRangeP :: Array r1 sh e -> MArr r2 e -> sh -> sh -> IO () + + + diff --git a/Data/Array/Repa/Eval/Gang.hs b/Data/Array/Repa/Eval/Gang.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Eval/Gang.hs @@ -0,0 +1,218 @@ +{-# LANGUAGE CPP #-} + +-- | Gang Primitives. +module Data.Array.Repa.Eval.Gang + ( theGang + , Gang, forkGang, gangSize, gangIO, gangST) +where +import GHC.IO +import GHC.ST +import GHC.Conc (forkOn) +import Control.Concurrent.MVar +import Control.Exception (assert) +import Control.Monad +import GHC.Conc (numCapabilities) +import System.IO + + +-- TheGang -------------------------------------------------------------------- +-- | This globally shared gang is auto-initialised at startup and shared by all +-- Repa computations. +-- +-- In a data parallel setting, it does not help to have multiple gangs +-- running at the same time. This is because a single data parallel +-- computation should already be able to keep all threads busy. If we had +-- multiple gangs running at the same time, then the system as a whole would +-- run slower as the gangs would contend for cache and thrash the scheduler. +-- +-- If, due to laziness or otherwise, you try to start multiple parallel +-- Repa computations at the same time, then you will get the following +-- warning on stderr at runtime: +-- +-- @Data.Array.Repa: Performing nested parallel computation sequentially. +-- You've probably called the 'compute' or 'copy' function while another +-- instance was already running. This can happen if the second version +-- was suspended due to lazy evaluation. Use 'deepSeqArray' to ensure that +-- each array is fully evaluated before you 'compute' the next one. +-- @ +-- +theGang :: Gang +{-# NOINLINE theGang #-} +theGang + = unsafePerformIO + $ do let caps = numCapabilities + forkGang caps + + +-- Requests ------------------------------------------------------------------- +-- | The 'Req' type encapsulates work requests for individual members of a gang. +data Req + -- | Instruct the worker to run the given action. + = ReqDo (Int -> IO ()) + + -- | Tell the worker that we're shutting the gang down. + -- The worker should signal that it's receieved the request by + -- writing to its result var before returning to the caller (forkGang). + | ReqShutdown + + +-- Gang ----------------------------------------------------------------------- +-- | A 'Gang' is a group of threads that execute arbitrary work requests. +data Gang + = Gang + { -- | Number of threads in the gang. + _gangThreads :: !Int + + -- | Workers listen for requests on these vars. + , _gangRequestVars :: [MVar Req] + + -- | Workers put their results in these vars. + , _gangResultVars :: [MVar ()] + + -- | Indicates that the gang is busy. + , _gangBusy :: MVar Bool + } + +instance Show Gang where + showsPrec p (Gang n _ _ _) + = showString "<<" + . showsPrec p n + . showString " threads>>" + + +-- | O(1). Yield the number of threads in the 'Gang'. +gangSize :: Gang -> Int +gangSize (Gang n _ _ _) + = n + + +-- | Fork a 'Gang' with the given number of threads (at least 1). +forkGang :: Int -> IO Gang +forkGang n + = assert (n > 0) + $ do + -- Create the vars we'll use to issue work requests. + mvsRequest <- sequence $ replicate n $ newEmptyMVar + + -- Create the vars we'll use to signal that threads are done. + mvsDone <- sequence $ replicate n $ newEmptyMVar + + -- Add finalisers so we can shut the workers down cleanly if they + -- become unreachable. + zipWithM_ (\varReq varDone + -> addMVarFinalizer varReq (finaliseWorker varReq varDone)) + mvsRequest + mvsDone + + -- Create all the worker threads + zipWithM_ forkOn [0..] + $ zipWith3 gangWorker + [0 .. n-1] mvsRequest mvsDone + + -- The gang is currently idle. + busy <- newMVar False + + return $ Gang n mvsRequest mvsDone busy + + + +-- | The worker thread of a 'Gang'. +-- The threads blocks on the MVar waiting for a work request. +gangWorker :: Int -> MVar Req -> MVar () -> IO () +gangWorker threadId varRequest varDone + = do + -- Wait for a request + req <- takeMVar varRequest + + case req of + ReqDo action + -> do -- Run the action we were given. + action threadId + + -- Signal that the action is complete. + putMVar varDone () + + -- Wait for more requests. + gangWorker threadId varRequest varDone + + ReqShutdown + -> putMVar varDone () + + +-- | Finaliser for worker threads. +-- We want to shutdown the corresponding thread when it's MVar becomes +-- unreachable. +-- Without this Repa programs can complain about "Blocked indefinitely +-- on an MVar" because worker threads are still blocked on the request +-- MVars when the program ends. Whether the finalizer is called or not +-- is very racey. It happens about 1 in 10 runs when for the +-- repa-edgedetect benchmark, and less often with the others. +-- +-- We're relying on the comment in System.Mem.Weak that says +-- "If there are no other threads to run, the runtime system will +-- check for runnablefinalizers before declaring the system to be +-- deadlocked." +-- +-- If we were creating and destroying the gang cleanly we wouldn't need +-- this, but theGang is created with a top-level unsafePerformIO. +-- Hacks beget hacks beget hacks... +-- +finaliseWorker :: MVar Req -> MVar () -> IO () +finaliseWorker varReq varDone + = do putMVar varReq ReqShutdown + takeMVar varDone + return () + + +-- | Issue work requests for the 'Gang' and wait until they complete. +-- +-- If the gang is already busy then print a warning to `stderr` and just +-- run the actions sequentially in the requesting thread. +gangIO :: Gang + -> (Int -> IO ()) + -> IO () + +{-# NOINLINE gangIO #-} +gangIO gang@(Gang _ _ _ busy) action + = do b <- swapMVar busy True + if b + then do + seqIO gang action + + else do + parIO gang action + _ <- swapMVar busy False + return () + + +-- | Run an action on the gang sequentially. +seqIO :: Gang -> (Int -> IO ()) -> IO () +seqIO (Gang n _ _ _) action + = do hPutStr stderr + $ unlines + [ "Data.Array.Repa: Performing nested parallel computation sequentially." + , " You've probably called the 'compute' or 'copy' function while another" + , " instance was already running. This can happen if the second version" + , " was suspended due to lazy evaluation. Use 'deepSeqArray' to ensure" + , " that each array is fully evaluated before you 'compute' the next one." + , "" ] + + mapM_ action [0 .. n-1] + +-- | Run an action on the gang in parallel. +parIO :: Gang -> (Int -> IO ()) -> IO () +parIO (Gang _ mvsRequest mvsResult _) action + = do + -- Send requests to all the threads. + mapM_ (\v -> putMVar v (ReqDo action)) mvsRequest + + -- Wait for all the requests to complete. + mapM_ takeMVar mvsResult + + +-- | Same as 'gangIO' but in the 'ST' monad. +gangST :: Gang -> (Int -> ST s ()) -> ST s () +gangST g p = unsafeIOToST . gangIO g $ unsafeSTToIO . p + + + diff --git a/Data/Array/Repa/Eval/Reduction.hs b/Data/Array/Repa/Eval/Reduction.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Eval/Reduction.hs @@ -0,0 +1,256 @@ +{-# LANGUAGE BangPatterns, MagicHash #-} +module Data.Array.Repa.Eval.Reduction + ( foldS, foldP + , foldAllS, foldAllP) +where +import Data.Array.Repa.Eval.Elt +import Data.Array.Repa.Eval.Gang +import qualified Data.Vector.Unboxed as V +import qualified Data.Vector.Unboxed.Mutable as M +import GHC.Base ( quotInt, divInt ) +import GHC.Exts + +-- | Sequential reduction of a multidimensional array along the innermost dimension. +foldS :: (Elt a, V.Unbox a) + => M.IOVector a -- ^ vector to write elements into + -> (Int# -> a) -- ^ function to get an element from the given index + -> (a -> a -> a) -- ^ binary associative combination function + -> a -- ^ starting value (typically an identity) + -> Int# -- ^ inner dimension (length to fold over) + -> IO () +{-# INLINE [1] foldS #-} +foldS vec !get !c !r !n + = iter 0# 0# + where + !(I# end) = M.length vec + + {-# INLINE iter #-} + iter !sh !sz + | sh >=# end = return () + | otherwise + = do let !next = sz +# n + M.unsafeWrite vec (I# sh) (reduceAny get c r sz next) + iter (sh +# 1#) next + + +-- | Parallel reduction of a multidimensional array along the innermost dimension. +-- Each output value is computed by a single thread, with the output values +-- distributed evenly amongst the available threads. +foldP :: (Elt a, V.Unbox a) + => M.IOVector a -- ^ vector to write elements into + -> (Int -> a) -- ^ function to get an element from the given index + -> (a -> a -> a) -- ^ binary associative combination operator + -> a -- ^ starting value. Must be neutral with respect + -- ^ to the operator. eg @0 + a = a@. + -> Int -- ^ inner dimension (length to fold over) + -> IO () +{-# INLINE [1] foldP #-} +foldP vec !f !c !r !(I# n) + = gangIO theGang + $ \(I# tid) -> fill (split tid) (split (tid +# 1#)) + where + !(I# threads) = gangSize theGang + !(I# len) = M.length vec + !step = (len +# threads -# 1#) `quotInt#` threads + + {-# INLINE split #-} + split !ix + = let !ix' = ix *# step + in if len <# ix' + then len + else ix' + + {-# INLINE fill #-} + fill !start !end + = iter start (start *# n) + where + {-# INLINE iter #-} + iter !sh !sz + | sh >=# end = return () + | otherwise + = do let !next = sz +# n + M.unsafeWrite vec (I# sh) (reduce f c r (I# sz) (I# next)) + iter (sh +# 1#) next + + +-- | Sequential reduction of all the elements in an array. +foldAllS :: (Elt a, V.Unbox a) + => (Int# -> a) -- ^ function to get an element from the given index + -> (a -> a -> a) -- ^ binary associative combining function + -> a -- ^ starting value + -> Int# -- ^ number of elements + -> a + +{-# INLINE [1] foldAllS #-} +foldAllS !f !c !r !len + = reduceAny (\i -> f i) c r 0# len + + + +-- | Parallel tree reduction of an array to a single value. Each thread takes an +-- equally sized chunk of the data and computes a partial sum. The main thread +-- then reduces the array of partial sums to the final result. +-- +-- We don't require that the initial value be a neutral element, so each thread +-- computes a fold1 on its chunk of the data, and the seed element is only +-- applied in the final reduction step. +-- +foldAllP :: (Elt a, V.Unbox a) + => (Int -> a) -- ^ function to get an element from the given index + -> (a -> a -> a) -- ^ binary associative combining function + -> a -- ^ starting value + -> Int -- ^ number of elements + -> IO a +{-# INLINE [1] foldAllP #-} + +foldAllP !f !c !r !len + | len == 0 = return r + | otherwise = do + mvec <- M.unsafeNew chunks + gangIO theGang $ \tid -> fill mvec tid (split tid) (split (tid+1)) + vec <- V.unsafeFreeze mvec + return $! V.foldl' c r vec + where + !threads = gangSize theGang + !step = (len + threads - 1) `quotInt` threads + chunks = ((len + step - 1) `divInt` step) `min` threads + + {-# INLINE split #-} + split !ix = len `min` (ix * step) + + {-# INLINE fill #-} + fill !mvec !tid !start !end + | start >= end = return () + | otherwise = M.unsafeWrite mvec tid (reduce f c (f start) (start+1) end) + + + +-- Reduce --------------------------------------------------------------------- +-- | This is the primitive reduction function. +-- We use manual specialisations and rewrite rules to avoid the result +-- being boxed up in the final iteration. +{-# INLINE [0] reduce #-} +reduce :: (Int -> a) -- ^ Get data from the array. + -> (a -> a -> a) -- ^ Function to combine elements. + -> a -- ^ Starting value. + -> Int -- ^ Starting index in array. + -> Int -- ^ Ending index in array. + -> a -- ^ Result. +reduce f c r (I# start) (I# end) + = reduceAny (\i -> f (I# i)) c r start end + + +-- | Sequentially reduce values between the given indices +{-# INLINE [0] reduceAny #-} +reduceAny :: (Int# -> a) -> (a -> a -> a) -> a -> Int# -> Int# -> a +reduceAny !f !c !r !start !end + = iter start r + where + {-# INLINE iter #-} + iter !i !z + | i >=# end = z + | otherwise = iter (i +# 1#) (f i `c` z) + + +{-# INLINE [0] reduceInt #-} +reduceInt + :: (Int# -> Int#) + -> (Int# -> Int# -> Int#) + -> Int# + -> Int# -> Int# + -> Int# + +reduceInt !f !c !r !start !end + = iter start r + where + {-# INLINE iter #-} + iter !i !z + | i >=# end = z + | otherwise = iter (i +# 1#) (f i `c` z) + + +{-# INLINE [0] reduceFloat #-} +reduceFloat + :: (Int# -> Float#) + -> (Float# -> Float# -> Float#) + -> Float# + -> Int# -> Int# + -> Float# + +reduceFloat !f !c !r !start !end + = iter start r + where + {-# INLINE iter #-} + iter !i !z + | i >=# end = z + | otherwise = iter (i +# 1#) (f i `c` z) + + +{-# INLINE [0] reduceDouble #-} +reduceDouble + :: (Int# -> Double#) + -> (Double# -> Double# -> Double#) + -> Double# + -> Int# -> Int# + -> Double# + +reduceDouble !f !c !r !start !end + = iter start r + where + {-# INLINE iter #-} + iter !i !z + | i >=# end = z + | otherwise = iter (i +# 1#) (f i `c` z) + + +{-# INLINE unboxInt #-} +unboxInt :: Int -> Int# +unboxInt (I# i) = i + + +{-# INLINE unboxFloat #-} +unboxFloat :: Float -> Float# +unboxFloat (F# f) = f + + +{-# INLINE unboxDouble #-} +unboxDouble :: Double -> Double# +unboxDouble (D# d) = d + + +{-# RULES "reduceInt" + forall (get :: Int# -> Int) f r start end + . reduceAny get f r start end + = I# (reduceInt + (\i -> unboxInt (get i)) + (\d1 d2 -> unboxInt (f (I# d1) (I# d2))) + (unboxInt r) + start + end) + #-} + + +{-# RULES "reduceFloat" + forall (get :: Int# -> Float) f r start end + . reduceAny get f r start end + = F# (reduceFloat + (\i -> unboxFloat (get i)) + (\d1 d2 -> unboxFloat (f (F# d1) (F# d2))) + (unboxFloat r) + start + end) + #-} + + +{-# RULES "reduceDouble" + forall (get :: Int# -> Double) f r start end + . reduceAny get f r start end + = D# (reduceDouble + (\i -> unboxDouble (get i)) + (\d1 d2 -> unboxDouble (f (D# d1) (D# d2))) + (unboxDouble r) + start + end) + #-} + + diff --git a/Data/Array/Repa/Eval/Selection.hs b/Data/Array/Repa/Eval/Selection.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Eval/Selection.hs @@ -0,0 +1,126 @@ +{-# LANGUAGE BangPatterns, ExplicitForAll, ScopedTypeVariables, PatternGuards #-} +module Data.Array.Repa.Eval.Selection + (selectChunkedS, selectChunkedP) +where +import Data.Array.Repa.Eval.Gang +import Data.Array.Repa.Shape +import Data.Vector.Unboxed as V +import Data.Vector.Unboxed.Mutable as VM +import GHC.Base (remInt, quotInt) +import Prelude as P +import Control.Monad as P +import Data.IORef + + +-- | Select indices matching a predicate. +-- +-- * This primitive can be useful for writing filtering functions. +-- +selectChunkedS + :: Shape sh + => (sh -> a -> IO ()) -- ^ Update function to write into result. + -> (sh -> Bool) -- ^ See if this predicate matches. + -> (sh -> a) -- ^ .. and apply fn to the matching index + -> sh -- ^ Extent of indices to apply to predicate. + -> IO Int -- ^ Number of elements written to destination array. + +{-# INLINE selectChunkedS #-} +selectChunkedS !fnWrite !fnMatch !fnProduce !shSize + = fill 0 0 + where lenSrc = size shSize + + fill !nSrc !nDst + | nSrc >= lenSrc = return nDst + + | ixSrc <- fromIndex shSize nSrc + , fnMatch ixSrc + = do fnWrite ixSrc (fnProduce ixSrc) + fill (nSrc + 1) (nDst + 1) + + | otherwise + = fill (nSrc + 1) nDst + + +-- | Select indices matching a predicate, in parallel. +-- +-- * This primitive can be useful for writing filtering functions. +-- +-- * The array is split into linear chunks, with one chunk being given to +-- each thread. +-- +-- * The number of elements in the result array depends on how many threads +-- you're running the program with. +-- +selectChunkedP + :: forall a + . Unbox a + => (Int -> Bool) -- ^ See if this predicate matches. + -> (Int -> a) -- .. and apply fn to the matching index + -> Int -- Extent of indices to apply to predicate. + -> IO [IOVector a] -- Chunks containing array elements. + +{-# INLINE selectChunkedP #-} +selectChunkedP !fnMatch !fnProduce !len + = do + -- Make IORefs that the threads will write their result chunks to. + -- We start with a chunk size proportial to the number of threads we have, + -- but the threads themselves can grow the chunks if they run out of space. + refs <- P.replicateM threads + $ do vec <- VM.new $ len `div` threads + newIORef vec + + -- Fire off a thread to fill each chunk. + gangIO theGang + $ \thread -> makeChunk (refs !! thread) + (splitIx thread) + (splitIx (thread + 1) - 1) + + -- Read the result chunks back from the IORefs. + -- If a thread had to grow a chunk, then these might not be the same ones + -- we created back in the first step. + P.mapM readIORef refs + + where -- See how many threads we have available. + !threads = gangSize theGang + !chunkLen = len `quotInt` threads + !chunkLeftover = len `remInt` threads + + + -- Decide where to split the source array. + {-# INLINE splitIx #-} + splitIx thread + | thread < chunkLeftover = thread * (chunkLen + 1) + | otherwise = thread * chunkLen + chunkLeftover + + + -- Fill the given chunk with elements selected from this range of indices. + makeChunk :: IORef (IOVector a) -> Int -> Int -> IO () + makeChunk !ref !ixSrc !ixSrcEnd + = do vecDst <- VM.new (len `div` threads) + vecDst' <- fillChunk ixSrc ixSrcEnd vecDst 0 (VM.length vecDst - 1) + writeIORef ref vecDst' + + + -- The main filling loop. + fillChunk :: Int -> Int -> IOVector a -> Int -> Int -> IO (IOVector a) + fillChunk !ixSrc !ixSrcEnd !vecDst !ixDst !ixDstEnd + -- If we've finished selecting elements, then slice the vector down + -- so it doesn't have any empty space at the end. + | ixSrc >= ixSrcEnd + = return $ VM.slice 0 ixDst vecDst + + -- If we've run out of space in the chunk then grow it some more. + | ixDst >= ixDstEnd + = do let ixDstEnd' = VM.length vecDst * 2 - 1 + vecDst' <- VM.grow vecDst (ixDstEnd + 1) + fillChunk (ixSrc + 1) ixSrcEnd vecDst' (ixDst + 1) ixDstEnd' + + -- We've got a maching element, so add it to the chunk. + | fnMatch ixSrc + = do VM.unsafeWrite vecDst ixDst (fnProduce ixSrc) + fillChunk (ixSrc + 1) ixSrcEnd vecDst (ixDst + 1) ixDstEnd + + -- The element doesnt match, so keep going. + | otherwise + = fillChunk (ixSrc + 1) ixSrcEnd vecDst ixDst ixDstEnd + diff --git a/Data/Array/Repa/Index.hs b/Data/Array/Repa/Index.hs --- a/Data/Array/Repa/Index.hs +++ b/Data/Array/Repa/Index.hs @@ -27,7 +27,7 @@ -- | Our index type, used for both shapes and indices. infixl 3 :. data tail :. head - = tail :. head + = !tail :. !head deriving (Show, Eq, Ord) -- Common dimensions @@ -41,39 +41,42 @@ -- Shape ------------------------------------------------------------------------------------------ instance Shape Z where - {-# INLINE rank #-} + {-# INLINE [1] rank #-} rank _ = 0 - {-# INLINE zeroDim #-} - zeroDim = Z + {-# INLINE [1] zeroDim #-} + zeroDim = Z - {-# INLINE unitDim #-} + {-# INLINE [1] unitDim #-} unitDim = Z - {-# INLINE intersectDim #-} + {-# INLINE [1] intersectDim #-} intersectDim _ _ = Z - {-# INLINE addDim #-} + {-# INLINE [1] addDim #-} addDim _ _ = Z - {-# INLINE size #-} + {-# INLINE [1] size #-} size _ = 1 - {-# INLINE sizeIsValid #-} + {-# INLINE [1] sizeIsValid #-} sizeIsValid _ = True - {-# INLINE toIndex #-} + {-# INLINE [1] toIndex #-} toIndex _ _ = 0 - {-# INLINE fromIndex #-} + {-# INLINE [1] fromIndex #-} fromIndex _ _ = Z - {-# INLINE inShapeRange #-} + {-# INLINE [1] inShapeRange #-} inShapeRange Z Z Z = True + {-# NOINLINE listOfShape #-} listOfShape _ = [] + + {-# NOINLINE shapeOfList #-} shapeOfList [] = Z shapeOfList _ = error $ stage ++ ".fromList: non-empty list when converting to Z." @@ -82,29 +85,29 @@ instance Shape sh => Shape (sh :. Int) where - {-# INLINE rank #-} + {-# INLINE [1] rank #-} rank (sh :. _) = rank sh + 1 - {-# INLINE zeroDim #-} + {-# INLINE [1] zeroDim #-} zeroDim = zeroDim :. 0 - {-# INLINE unitDim #-} + {-# INLINE [1] unitDim #-} unitDim = unitDim :. 1 - {-# INLINE intersectDim #-} + {-# INLINE [1] intersectDim #-} intersectDim (sh1 :. n1) (sh2 :. n2) = (intersectDim sh1 sh2 :. (min n1 n2)) - {-# INLINE addDim #-} + {-# INLINE [1] addDim #-} addDim (sh1 :. n1) (sh2 :. n2) = addDim sh1 sh2 :. (n1 + n2) - {-# INLINE size #-} + {-# INLINE [1] size #-} size (sh1 :. n) = size sh1 * n - {-# INLINE sizeIsValid #-} + {-# INLINE [1] sizeIsValid #-} sizeIsValid (sh1 :. n) | size sh1 > 0 = n <= maxBound `div` size sh1 @@ -112,29 +115,30 @@ | otherwise = False - {-# INLINE toIndex #-} + {-# INLINE [1] toIndex #-} toIndex (sh1 :. sh2) (sh1' :. sh2') = toIndex sh1 sh1' * sh2 + sh2' - {-# INLINE fromIndex #-} - fromIndex (ds :. d) n - = fromIndex ds (n `quotInt` d) :. r - where - -- If we assume that the index is in range, there is no point - -- in computing the remainder for the highest dimension since - -- n < d must hold. This saves one remInt per element access which - -- is quite a big deal. - r | rank ds == 0 = n - | otherwise = n `remInt` d + {-# INLINE [1] fromIndex #-} + fromIndex (ds :. d) n + = fromIndex ds (n `quotInt` d) :. r + where + -- If we assume that the index is in range, there is no point + -- in computing the remainder for the highest dimension since + -- n < d must hold. This saves one remInt per element access which + -- is quite a big deal. + r | rank ds == 0 = n + | otherwise = n `remInt` d - {-# INLINE inShapeRange #-} + {-# INLINE [1] inShapeRange #-} inShapeRange (zs :. z) (sh1 :. n1) (sh2 :. n2) = (n2 >= z) && (n2 < n1) && (inShapeRange zs sh1 sh2) - + {-# NOINLINE listOfShape #-} listOfShape (sh :. n) = n : listOfShape sh + {-# NOINLINE shapeOfList #-} shapeOfList xx = case xx of [] -> error $ stage ++ ".toList: empty list when converting to (_ :. Int)" diff --git a/Data/Array/Repa/Internals/Base.hs b/Data/Array/Repa/Internals/Base.hs deleted file mode 100644 --- a/Data/Array/Repa/Internals/Base.hs +++ /dev/null @@ -1,410 +0,0 @@ -{-# LANGUAGE ExplicitForAll, TypeOperators, FlexibleInstances, UndecidableInstances, BangPatterns, - ExistentialQuantification #-} -module Data.Array.Repa.Internals.Base - ( Array (..) - , Region(..) - , Range (..) - , Rect (..) - , Generator(..) - , deepSeqArray, deepSeqArrays - - , singleton, toScalar - , extent, delay - - -- * Predicates - , inRange - - -- * Indexing - , (!), index - , (!?), safeIndex - , unsafeIndex - - -- * Construction - , fromFunction - , fromVector - , fromList - , unsafeFromForeignPtr) -where -import Data.Array.Repa.Index -import Data.Array.Repa.Internals.Elt -import Data.Array.Repa.Shape as S -import qualified Data.Vector.Unboxed as V -import Data.Vector.Unboxed (Vector) -import Foreign.ForeignPtr -import Foreign.Storable -import System.IO.Unsafe - -stage = "Data.Array.Repa.Array" - --- Array ---------------------------------------------------------------------- --- | Repa arrays. -data Array sh a - = Array - { -- | The entire extent of the array. - arrayExtent :: sh - - -- | Arrays can be partitioned into several regions. - , arrayRegions :: [Region sh a] } - - --- | Defines the values in a region of the array. -data Region sh a - = Region - { -- | The range of elements this region applies to. - regionRange :: Range sh - - -- | How to compute the array elements in this region. - , regionGenerator :: Generator sh a } - - --- | Represents a range of elements in the array. -data Range sh - -- | Covers the entire array. - = RangeAll - - -- | The union of a possibly disjoint set of rectangles. - | RangeRects - { rangeMatch :: sh -> Bool - , rangeRects :: [Rect sh] } - - --- | A rectangle\/cube of arbitrary dimension. --- The indices are of the minimum and maximim elements to fill. -data Rect sh - = Rect sh sh - --- | Generates array elements for a particular region in the array. -data Generator sh a - -- | Elements are already computed and sitting in this vector. - = GenManifest (Vector a) - -- NOTE: Don't make the vector field strict. If you do then deepSeqing arrays - -- outside of loops won't cause the unboxings to be floated out. - - -- | Elements can be computed using these cursor functions. - | forall cursor - . GenCursor - { -- | Make a cursor to a particular element. - genMakeCursor :: sh -> cursor - - -- | Shift the cursor by an offset, to get to another element. - , genShiftCursor :: sh -> cursor -> cursor - - -- | Load\/compute the element at the given cursor. - , genLoadElem :: cursor -> a } - - --- DeepSeqs ------------------------------------------------------------------- --- | Ensure the structure for an array is fully evaluated. --- As we are in a lazy language, applying the @force@ function to a delayed array doesn't --- actually compute it at that point. Rather, Haskell builds a suspension representing the --- appliction of the @force@ function to that array. Use @deepSeqArray@ to ensure the array --- is actually computed at a particular point in the program. -infixr 0 `deepSeqArray` -deepSeqArray :: Shape sh => Array sh a -> b -> b -{-# INLINE deepSeqArray #-} -deepSeqArray (Array ex rgns) x - = ex `S.deepSeq` rgns `deepSeqRegions` x - --- | Like `deepSeqArray` but seqs all the arrays in a list. --- This is specialised up to lists of 4 arrays. Using more in the list will break fusion. -infixr 0 `deepSeqArrays` -deepSeqArrays :: Shape sh => [Array sh a] -> b -> b -{-# INLINE deepSeqArrays #-} -deepSeqArrays as y - = case as of - [] -> y - [a] -> a `deepSeqArray` y - [a1, a2] -> a1 `deepSeqArray` a2 `deepSeqArray` y - [a1, a2, a3] -> a1 `deepSeqArray` a2 `deepSeqArray` a3 `deepSeqArray` y - [a1, a2, a3, a4]-> a1 `deepSeqArray` a2 `deepSeqArray` a3 `deepSeqArray` a4 `deepSeqArray` y - _ -> deepSeqArrays' as y - -deepSeqArrays' as' y - = case as' of - [] -> y - x : xs -> x `deepSeqArray` xs `deepSeqArrays` y - --- | Ensure the structure for a region is fully evaluated. -infixr 0 `deepSeqRegion` -deepSeqRegion :: Shape sh => Region sh a -> b -> b -{-# INLINE deepSeqRegion #-} -deepSeqRegion (Region range gen) x - = range `deepSeqRange` gen `deepSeqGen` x - - --- | Ensure the structure for some regions are fully evaluated. -infixr 0 `deepSeqRegions` -deepSeqRegions :: Shape sh => [Region sh a] -> b -> b -{-# INLINE deepSeqRegions #-} -deepSeqRegions rs y - = case rs of - [] -> y - [r] -> r `deepSeqRegion` y - [r1, r2] -> r1 `deepSeqRegion` r2 `deepSeqRegion` y - rs' -> deepSeqRegions' rs' y - -deepSeqRegions' rs' y - = case rs' of - [] -> y - x : xs -> x `deepSeqRegion` xs `deepSeqRegions'` y - - --- | Ensure a range is fully evaluated. -infixr 0 `deepSeqRange` -deepSeqRange :: Shape sh => Range sh -> b -> b -{-# INLINE deepSeqRange #-} -deepSeqRange range x - = case range of - RangeAll -> x - RangeRects f rects -> f `seq` rects `seq` x - - --- | Ensure a Generator's structure is fully evaluated. -infixr 0 `deepSeqGen` -deepSeqGen :: Shape sh => Generator sh a -> b -> b -{-# INLINE deepSeqGen #-} -deepSeqGen gen x - = case gen of - GenManifest vec -> vec `seq` x - GenCursor{} -> x - - --- Predicates ------------------------------------------------------------------------------------- -inRange :: Shape sh => Range sh -> sh -> Bool -{-# INLINE inRange #-} -inRange RangeAll _ = True -inRange (RangeRects fn _) ix = fn ix - - --- Singletons ------------------------------------------------------------------------------------- --- | Wrap a scalar into a singleton array. -singleton :: Elt a => a -> Array Z a -{-# INLINE singleton #-} -singleton = fromFunction Z . const - --- | Take the scalar value from a singleton array. -toScalar :: Elt a => Array Z a -> a -{-# INLINE toScalar #-} -toScalar arr = arr ! Z - - --- Projections ------------------------------------------------------------------------------------ --- | Take the extent of an array. -extent :: Array sh a -> sh -{-# INLINE extent #-} -extent arr = arrayExtent arr - - --- | Unpack an array into delayed form. -delay :: (Shape sh, Elt a) - => Array sh a - -> (sh, sh -> a) - -{-# INLINE delay #-} -delay arr@(Array sh _) - = (sh, (arr !)) - - --- Indexing --------------------------------------------------------------------------------------- --- | Get an indexed element from an array. --- This uses the same level of bounds checking as your Data.Vector installation. -(!), index - :: forall sh a - . (Shape sh, Elt a) - => Array sh a - -> sh - -> a - -{-# INLINE (!) #-} -(!) arr ix = index arr ix - -{-# INLINE index #-} -index arr ix - = case arr of - Array _ [] - -> zero - - Array sh [Region _ gen1] - -> indexGen sh gen1 ix - - Array sh [Region r1 gen1, Region _ gen2] - | inRange r1 ix -> indexGen sh gen1 ix - | otherwise -> indexGen sh gen2 ix - - _ -> index' arr ix - - - where {-# INLINE indexGen #-} - indexGen sh gen ix' - = case gen of - GenManifest vec - -> vec V.! (S.toIndex sh ix') - - GenCursor makeCursor _ loadElem - -> loadElem $ makeCursor ix' - - index' (Array sh (Region range gen : rs)) ix' - | inRange range ix = indexGen sh gen ix' - | otherwise = index' (Array sh rs) ix' - - index' (Array _ []) _ - = zero - - - --- | Get an indexed element from an array. --- If the element is out of range then `Nothing`. -(!?), safeIndex - :: forall sh a - . (Shape sh, Elt a) - => Array sh a - -> sh - -> Maybe a - -{-# INLINE (!?) #-} -(!?) arr ix = safeIndex arr ix - - -{-# INLINE safeIndex #-} -safeIndex arr ix - = case arr of - Array _ [] - -> Nothing - - Array sh [Region _ gen1] - -> indexGen sh gen1 ix - - Array sh [Region r1 gen1, Region r2 gen2] - | inRange r1 ix -> indexGen sh gen1 ix - | inRange r2 ix -> indexGen sh gen2 ix - | otherwise -> Nothing - - _ -> index' arr ix - - - where {-# INLINE indexGen #-} - indexGen sh gen ix' - = case gen of - GenManifest vec - -> vec V.!? (S.toIndex sh ix') - - GenCursor makeCursor _ loadElem - -> Just (loadElem $ makeCursor ix') - - index' (Array sh (Region range gen : rs)) ix' - | inRange range ix = indexGen sh gen ix' - | otherwise = index' (Array sh rs) ix' - - index' (Array _ []) _ - = Nothing - - --- | Get an indexed element from an array, without bounds checking. --- This assumes that the regions in the array give full coverage. --- An array with no regions gets zero for every element. -unsafeIndex - :: forall sh a - . (Shape sh, Elt a) - => Array sh a - -> sh - -> a - -{-# INLINE unsafeIndex #-} -unsafeIndex arr ix - = case arr of - Array _ [] - -> zero - - Array sh [Region _ gen1] - -> unsafeIndexGen sh gen1 ix - - Array sh [Region r1 gen1, Region _ gen2] - | inRange r1 ix -> unsafeIndexGen sh gen1 ix - | otherwise -> unsafeIndexGen sh gen2 ix - - _ -> unsafeIndex' arr ix - - where {-# INLINE unsafeIndexGen #-} - unsafeIndexGen sh gen ix' - = case gen of - GenManifest vec - -> vec `V.unsafeIndex` (S.toIndex sh ix') - - GenCursor makeCursor _ loadElem - -> loadElem $ makeCursor ix' - - unsafeIndex' (Array sh (Region range gen : rs)) ix' - | inRange range ix = unsafeIndexGen sh gen ix' - | otherwise = unsafeIndex' (Array sh rs) ix' - - unsafeIndex' (Array _ []) _ - = zero - - --- Conversions ------------------------------------------------------------------------------------ --- | Create a `Delayed` array from a function. -fromFunction - :: Shape sh - => sh - -> (sh -> a) - -> Array sh a - -{-# INLINE fromFunction #-} -fromFunction sh fnElems - = sh `S.deepSeq` - Array sh [Region - RangeAll - (GenCursor id addDim fnElems)] - --- | Create a `Manifest` array from an unboxed `Vector`. --- The elements are in row-major order. -fromVector - :: Shape sh - => sh - -> Vector a - -> Array sh a - -{-# INLINE fromVector #-} -fromVector sh vec - = sh `S.deepSeq` vec `seq` - Array sh [Region RangeAll (GenManifest vec)] - - --- | Convert a list to an array. --- The length of the list must be exactly the `size` of the extent given, else `error`. -fromList - :: (Shape sh, Elt a) - => sh - -> [a] - -> Array sh a - -{-# INLINE fromList #-} -fromList sh xx - | V.length vec /= S.size sh - = error $ unlines - [ stage ++ ".fromList: size of array shape does not match size of list" - , " size of shape = " ++ (show $ S.size sh) ++ "\n" - , " size of list = " ++ (show $ V.length vec) ++ "\n" ] - - | otherwise - = Array sh [Region RangeAll (GenManifest vec)] - - where vec = V.fromList xx - - --- | Convert a `Ptr` to an `Array`. --- The data is used directly, and not copied. --- You promise not to modify the pointed-to data any further. --- -unsafeFromForeignPtr - :: (Shape sh, Elt a, Storable a) - => sh - -> ForeignPtr a - -> Array sh a - -unsafeFromForeignPtr sh fptr - = fromFunction sh - (\ix -> unsafePerformIO - $ withForeignPtr fptr - (\ptr -> peekElemOff ptr $ toIndex sh ix)) - diff --git a/Data/Array/Repa/Internals/Elt.hs b/Data/Array/Repa/Internals/Elt.hs deleted file mode 100644 --- a/Data/Array/Repa/Internals/Elt.hs +++ /dev/null @@ -1,283 +0,0 @@ --- | Values that can be stored in Repa Arrays. -{-# LANGUAGE MagicHash, UnboxedTuples, TypeSynonymInstances, FlexibleInstances #-} -module Data.Array.Repa.Internals.Elt - (Elt (..)) -where -import GHC.Prim -import GHC.Exts -import GHC.Types -import GHC.Word -import GHC.Int -import Data.Vector.Unboxed - - --- Note that the touch# function is special because we can pass it boxed or unboxed --- values. The argument type has kind ?, not just * or #. - --- | Element types that can be stored in Repa arrays. --- Repa uses `Data.Vector.Unboxed` to store the actual data. The implementation --- of this library is based on type families and picks an efficient, specialised --- representation for every element type. In particular, unboxed vectors of pairs --- are represented as pairs of unboxed vectors. -class (Show a, Unbox a) => Elt a where - - -- | We use this to prevent bindings from being floated inappropriatey. - -- Doing a `seq` sometimes isn't enough, because the GHC simplifier can - -- erase these, and/or still move around the bindings. - touch :: a -> IO () - - -- | Generic zero value, helpful for debugging. - zero :: a - - -- | Generic one value, helpful for debugging. - one :: a - - --- Bool ----------------------------------------------------------------------- -instance Elt Bool where - {-# INLINE touch #-} - touch b - = IO (\state -> case touch# b state of - state' -> (# state', () #)) - - {-# INLINE zero #-} - zero = False - - {-# INLINE one #-} - one = True - - --- Floating ------------------------------------------------------------------- -instance Elt Float where - {-# INLINE touch #-} - touch (F# f) - = IO (\state -> case touch# f state of - state' -> (# state', () #)) - - {-# INLINE zero #-} - zero = 0 - - {-# INLINE one #-} - one = 1 - - -instance Elt Double where - {-# INLINE touch #-} - touch (D# d) - = IO (\state -> case touch# d state of - state' -> (# state', () #)) - - {-# INLINE zero #-} - zero = 0 - - {-# INLINE one #-} - one = 1 - - --- Int ------------------------------------------------------------------------ -instance Elt Int where - {-# INLINE touch #-} - touch (I# i) - = IO (\state -> case touch# i state of - state' -> (# state', () #)) - - {-# INLINE zero #-} - zero = 0 - - {-# INLINE one #-} - one = 1 - -instance Elt Int8 where - {-# INLINE touch #-} - touch (I8# w) - = IO (\state -> case touch# w state of - state' -> (# state', () #)) - - {-# INLINE zero #-} - zero = 0 - - {-# INLINE one #-} - one = 1 - - -instance Elt Int16 where - {-# INLINE touch #-} - touch (I16# w) - = IO (\state -> case touch# w state of - state' -> (# state', () #)) - - {-# INLINE zero #-} - zero = 0 - - {-# INLINE one #-} - one = 1 - - -instance Elt Int32 where - {-# INLINE touch #-} - touch (I32# w) - = IO (\state -> case touch# w state of - state' -> (# state', () #)) - - {-# INLINE zero #-} - zero = 0 - - {-# INLINE one #-} - one = 1 - - -instance Elt Int64 where - {-# INLINE touch #-} - touch (I64# w) - = IO (\state -> case touch# w state of - state' -> (# state', () #)) - - {-# INLINE zero #-} - zero = 0 - - {-# INLINE one #-} - one = 1 - - --- Word ----------------------------------------------------------------------- -instance Elt Word where - {-# INLINE touch #-} - touch (W# i) - = IO (\state -> case touch# i state of - state' -> (# state', () #)) - - {-# INLINE zero #-} - zero = 0 - - {-# INLINE one #-} - one = 1 - - -instance Elt Word8 where - {-# INLINE touch #-} - touch (W8# w) - = IO (\state -> case touch# w state of - state' -> (# state', () #)) - - {-# INLINE zero #-} - zero = 0 - - {-# INLINE one #-} - one = 1 - - -instance Elt Word16 where - {-# INLINE touch #-} - touch (W16# w) - = IO (\state -> case touch# w state of - state' -> (# state', () #)) - - {-# INLINE zero #-} - zero = 0 - - {-# INLINE one #-} - one = 1 - - -instance Elt Word32 where - {-# INLINE touch #-} - touch (W32# w) - = IO (\state -> case touch# w state of - state' -> (# state', () #)) - - {-# INLINE zero #-} - zero = 0 - - {-# INLINE one #-} - one = 1 - - -instance Elt Word64 where - {-# INLINE touch #-} - touch (W64# w) - = IO (\state -> case touch# w state of - state' -> (# state', () #)) - - {-# INLINE zero #-} - zero = 0 - - {-# INLINE one #-} - one = 1 - - --- Tuple ---------------------------------------------------------------------- -instance (Elt a, Elt b) => Elt (a, b) where - {-# INLINE touch #-} - touch (a, b) - = do touch a - touch b - - {-# INLINE zero #-} - zero = (zero, zero) - - {-# INLINE one #-} - one = (one, one) - - -instance (Elt a, Elt b, Elt c) => Elt (a, b, c) where - {-# INLINE touch #-} - touch (a, b, c) - = do touch a - touch b - touch c - - {-# INLINE zero #-} - zero = (zero, zero, zero) - - {-# INLINE one #-} - one = (one, one, one) - - -instance (Elt a, Elt b, Elt c, Elt d) => Elt (a, b, c, d) where - {-# INLINE touch #-} - touch (a, b, c, d) - = do touch a - touch b - touch c - touch d - - {-# INLINE zero #-} - zero = (zero, zero, zero, zero) - - {-# INLINE one #-} - one = (one, one, one, one) - - -instance (Elt a, Elt b, Elt c, Elt d, Elt e) => Elt (a, b, c, d, e) where - {-# INLINE touch #-} - touch (a, b, c, d, e) - = do touch a - touch b - touch c - touch d - touch e - - {-# INLINE zero #-} - zero = (zero, zero, zero, zero, zero) - - {-# INLINE one #-} - one = (one, one, one, one, one) - - -instance (Elt a, Elt b, Elt c, Elt d, Elt e, Elt f) => Elt (a, b, c, d, e, f) where - {-# INLINE touch #-} - touch (a, b, c, d, e, f) - = do touch a - touch b - touch c - touch d - touch e - touch f - - {-# INLINE zero #-} - zero = (zero, zero, zero, zero, zero, zero) - - {-# INLINE one #-} - one = (one, one, one, one, one, one) - - diff --git a/Data/Array/Repa/Internals/EvalBlockwise.hs b/Data/Array/Repa/Internals/EvalBlockwise.hs deleted file mode 100644 --- a/Data/Array/Repa/Internals/EvalBlockwise.hs +++ /dev/null @@ -1,153 +0,0 @@ -{-# LANGUAGE BangPatterns #-} - --- | Old non-cursored, blockwise filling functions. --- NOTE: this isn't currently used. -module Data.Array.Repa.Internals.EvalBlockwise - ( fillVectorBlockwiseP - , fillVectorBlock - , fillVectorBlockP) -where -import Data.Array.Repa.Internals.Elt -import Data.Array.Repa.Internals.Gang -import Data.Vector.Unboxed.Mutable as VM -import GHC.Base (remInt, quotInt) -import Prelude as P - - --- Blockwise filling ------------------------------------------------------------------------------ -fillVectorBlockwiseP - :: Elt a - => IOVector a -- ^ vector to write elements into - -> (Int -> a) -- ^ fn to evaluate an element at the given index - -> Int -- ^ width of image. - -> IO () - -{-# INLINE [0] fillVectorBlockwiseP #-} -fillVectorBlockwiseP !vec !getElemFVBP !imageWidth - = gangIO theGang fillBlock - - where !threads = gangSize theGang - !vecLen = VM.length vec - !imageHeight = vecLen `div` imageWidth - !colChunkLen = imageWidth `quotInt` threads - !colChunkSlack = imageWidth `remInt` threads - - - {-# INLINE colIx #-} - colIx !ix - | ix < colChunkSlack = ix * (colChunkLen + 1) - | otherwise = ix * colChunkLen + colChunkSlack - - - -- just give one column to each thread - {-# INLINE fillBlock #-} - fillBlock :: Int -> IO () - fillBlock !ix - = let !x0 = colIx ix - !x1 = colIx (ix + 1) - !y0 = 0 - !y1 = imageHeight - in fillVectorBlock vec getElemFVBP imageWidth x0 y0 x1 y1 - - --- Block filling ---------------------------------------------------------------------------------- --- | Fill a block in a 2D image, in parallel. --- Coordinates given are of the filled edges of the block. --- We divide the block into columns, and give one column to each thread. -fillVectorBlockP - :: Elt a - => IOVector a -- ^ vector to write elements into - -> (Int -> a) -- ^ fn to evaluate an element at the given index. - -> Int -- ^ width of whole image - -> Int -- ^ x0 lower left corner of block to fill - -> Int -- ^ y0 (low x and y value) - -> Int -- ^ x1 upper right corner of block - -> Int -- ^ y1 (high x and y value, last index to fill) - -> IO () - -{-# INLINE [0] fillVectorBlockP #-} -fillVectorBlockP !vec !getElem !imageWidth !x0 !y0 !x1 !y1 - = gangIO theGang fillBlock - where !threads = gangSize theGang - !blockWidth = x1 - x0 + 1 - - -- All columns have at least this many pixels. - !colChunkLen = blockWidth `quotInt` threads - - -- Extra pixels that we have to divide between some of the threads. - !colChunkSlack = blockWidth `remInt` threads - - -- Get the starting pixel of a column in the image. - {-# INLINE colIx #-} - colIx !ix - | ix < colChunkSlack = x0 + ix * (colChunkLen + 1) - | otherwise = x0 + ix * colChunkLen + colChunkSlack - - -- Give one column to each thread - {-# INLINE fillBlock #-} - fillBlock :: Int -> IO () - fillBlock !ix - = let !x0' = colIx ix - !x1' = colIx (ix + 1) - 1 - !y0' = y0 - !y1' = y1 - in fillVectorBlock vec getElem imageWidth x0' y0' x1' y1' - - --- | Fill a block in a 2D image. --- Coordinates given are of the filled edges of the block. -fillVectorBlock - :: Elt a - => IOVector a -- ^ vector to write elements into. - -> (Int -> a) -- ^ fn to evaluate an element at the given index. - -> Int -- ^ width of whole image - -> Int -- ^ x0 lower left corner of block to fill - -> Int -- ^ y0 (low x and y value) - -> Int -- ^ x1 upper right corner of block - -> Int -- ^ y1 (high x and y value, last index to fill) - -> IO () - -{-# INLINE [0] fillVectorBlock #-} -fillVectorBlock !vec !getElemFVB !imageWidth !x0 !y0 !x1 !y1 - = do -- putStrLn $ "fillVectorBlock: " P.++ show (x0, y0, x1, y1) - fillBlock ixStart (ixStart + (x1 - x0)) - where - -- offset from end of one line to the start of the next. - !ixStart = x0 + y0 * imageWidth - !ixFinal = x1 + y1 * imageWidth - - {-# INLINE fillBlock #-} - fillBlock !ixLineStart !ixLineEnd - | ixLineStart > ixFinal = return () - | otherwise - = do fillLine4 ixLineStart - fillBlock (ixLineStart + imageWidth) (ixLineEnd + imageWidth) - - where {-# INLINE fillLine4 #-} - fillLine4 !ix - | ix + 4 > ixLineEnd = fillLine1 ix - | otherwise - = do - let d0 = getElemFVB (ix + 0) - let d1 = getElemFVB (ix + 1) - let d2 = getElemFVB (ix + 2) - let d3 = getElemFVB (ix + 3) - - touch d0 - touch d1 - touch d2 - touch d3 - - VM.unsafeWrite vec (ix + 0) d0 - VM.unsafeWrite vec (ix + 1) d1 - VM.unsafeWrite vec (ix + 2) d2 - VM.unsafeWrite vec (ix + 3) d3 - fillLine4 (ix + 4) - - {-# INLINE fillLine1 #-} - fillLine1 !ix - | ix > ixLineEnd = return () - | otherwise - = do VM.unsafeWrite vec ix (getElemFVB ix) - fillLine1 (ix + 1) - diff --git a/Data/Array/Repa/Internals/EvalChunked.hs b/Data/Array/Repa/Internals/EvalChunked.hs deleted file mode 100644 --- a/Data/Array/Repa/Internals/EvalChunked.hs +++ /dev/null @@ -1,62 +0,0 @@ --- | Evaluate a vector by breaking it up into linear chunks and filling each chunk --- in parallel. -{-# LANGUAGE BangPatterns #-} -module Data.Array.Repa.Internals.EvalChunked - ( fillChunkedS - , fillChunkedP) -where -import Data.Array.Repa.Internals.Gang -import GHC.Base (remInt, quotInt) -import Prelude as P - - --- | Fill something sequentially. -fillChunkedS - :: Int -- ^ Number of elements - -> (Int -> a -> IO ()) -- ^ Update function to write into result buffer - -> (Int -> a) -- ^ Fn to get the value at a given index. - -> IO () - -{-# INLINE [0] fillChunkedS #-} -fillChunkedS !len !write !getElem - = fill 0 - where fill !ix - | ix >= len = return () - | otherwise - = do write ix (getElem ix) - fill (ix + 1) - - --- | Fill something in parallel. -fillChunkedP - :: Int -- ^ Number of elements - -> (Int -> a -> IO ()) -- ^ Update function to write into result buffer - -> (Int -> a) -- ^ Fn to get the value at a given index. - -> IO () - -{-# INLINE [0] fillChunkedP #-} -fillChunkedP !len !write !getElem - = gangIO theGang - $ \thread -> fill (splitIx thread) (splitIx (thread + 1)) - - where - -- Decide now to split the work across the threads. - -- If the length of the vector doesn't divide evenly among the threads, - -- then the first few get an extra element. - !threads = gangSize theGang - !chunkLen = len `quotInt` threads - !chunkLeftover = len `remInt` threads - - {-# INLINE splitIx #-} - splitIx thread - | thread < chunkLeftover = thread * (chunkLen + 1) - | otherwise = thread * chunkLen + chunkLeftover - - -- Evaluate the elements of a single chunk. - {-# INLINE fill #-} - fill !ix !end - | ix >= end = return () - | otherwise - = do write ix (getElem ix) - fill (ix + 1) end - diff --git a/Data/Array/Repa/Internals/EvalCursored.hs b/Data/Array/Repa/Internals/EvalCursored.hs deleted file mode 100644 --- a/Data/Array/Repa/Internals/EvalCursored.hs +++ /dev/null @@ -1,136 +0,0 @@ - -{-# LANGUAGE BangPatterns, UnboxedTuples #-} -module Data.Array.Repa.Internals.EvalCursored - ( fillCursoredBlock2P - , fillCursoredBlock2 ) -where -import Data.Array.Repa.Index -import Data.Array.Repa.Internals.Elt -import Data.Array.Repa.Internals.Gang -import GHC.Base (remInt, quotInt) -import Prelude as P - - --- Block filling ---------------------------------------------------------------------------------- --- | Fill a block in a 2D image, in parallel. --- Coordinates given are of the filled edges of the block. --- We divide the block into columns, and give one column to each thread. -fillCursoredBlock2P - :: Elt a - => (Int -> a -> IO ()) -- ^ Update function to write into result buffer - -> (DIM2 -> cursor) -- ^ make a cursor to a particular element - -> (DIM2 -> cursor -> cursor) -- ^ shift the cursor by an offset - -> (cursor -> a) -- ^ fn to evaluate an element at the given index. - -> Int -- ^ width of whole image - -> Int -- ^ x0 lower left corner of block to fill - -> Int -- ^ y0 (low x and y value) - -> Int -- ^ x1 upper right corner of block to fill - -> Int -- ^ y1 (high x and y value, index of last elem to fill) - -> IO () - -{-# INLINE [0] fillCursoredBlock2P #-} -fillCursoredBlock2P - !write - !makeCursorFCB !shiftCursorFCB !getElemFCB - !imageWidth !x0 !y0 !x1 !y1 - = gangIO theGang fillBlock - where !threads = gangSize theGang - !blockWidth = x1 - x0 + 1 - - -- All columns have at least this many pixels. - !colChunkLen = blockWidth `quotInt` threads - - -- Extra pixels that we have to divide between some of the threads. - !colChunkSlack = blockWidth `remInt` threads - - -- Get the starting pixel of a column in the image. - {-# INLINE colIx #-} - colIx !ix - | ix < colChunkSlack = x0 + ix * (colChunkLen + 1) - | otherwise = x0 + ix * colChunkLen + colChunkSlack - - -- Give one column to each thread - {-# INLINE fillBlock #-} - fillBlock :: Int -> IO () - fillBlock !ix - = let !x0' = colIx ix - !x1' = colIx (ix + 1) - 1 - !y0' = y0 - !y1' = y1 - in fillCursoredBlock2 - write - makeCursorFCB shiftCursorFCB getElemFCB - imageWidth x0' y0' x1' y1' - - --- | Fill a block in a 2D image. --- Coordinates given are of the filled edges of the block. -fillCursoredBlock2 - :: Elt a - => (Int -> a -> IO ()) -- ^ Update function to write into result buffer - -> (DIM2 -> cursor) -- ^ make a cursor to a particular element - -> (DIM2 -> cursor -> cursor) -- ^ shift the cursor by an offset - -> (cursor -> a) -- ^ fn to evaluate an element at the given index. - -> Int -- ^ width of whole image - -> Int -- ^ x0 lower left corner of block to fill - -> Int -- ^ y0 (low x and y value) - -> Int -- ^ x1 upper right corner of block to fill - -> Int -- ^ y1 (high x and y value, index of last elem to fill) - -> IO () - -{-# INLINE [0] fillCursoredBlock2 #-} -fillCursoredBlock2 - !write - !makeCursor !shiftCursor !getElem - !imageWidth !x0 !y0 !x1 !y1 - - = fillBlock y0 - - where {-# INLINE fillBlock #-} - fillBlock !y - | y > y1 = return () - | otherwise - = do fillLine4 x0 - fillBlock (y + 1) - - where {-# INLINE fillLine4 #-} - fillLine4 !x - | x + 4 > x1 = fillLine1 x - | otherwise - = do -- Compute each source cursor based on the previous one so that - -- the variable live ranges in the generated code are shorter. - let srcCur0 = makeCursor (Z :. y :. x) - let srcCur1 = shiftCursor (Z :. 0 :. 1) srcCur0 - let srcCur2 = shiftCursor (Z :. 0 :. 1) srcCur1 - let srcCur3 = shiftCursor (Z :. 0 :. 1) srcCur2 - - -- Get the result value for each cursor. - let val0 = getElem srcCur0 - let val1 = getElem srcCur1 - let val2 = getElem srcCur2 - let val3 = getElem srcCur3 - - -- Ensure that we've computed each of the result values before we - -- write into the array. If the backend code generator can't tell - -- our destination array doesn't alias with the source then writing - -- to it can prevent the sharing of intermediate computations. - touch val0 - touch val1 - touch val2 - touch val3 - - -- Compute cursor into destination array. - let !dstCur0 = x + y * imageWidth - write (dstCur0) val0 - write (dstCur0 + 1) val1 - write (dstCur0 + 2) val2 - write (dstCur0 + 3) val3 - fillLine4 (x + 4) - - {-# INLINE fillLine1 #-} - fillLine1 !x - | x > x1 = return () - | otherwise - = do write (x + y * imageWidth) (getElem $ makeCursor (Z :. y :. x)) - fillLine1 (x + 1) - diff --git a/Data/Array/Repa/Internals/EvalReduction.hs b/Data/Array/Repa/Internals/EvalReduction.hs deleted file mode 100644 --- a/Data/Array/Repa/Internals/EvalReduction.hs +++ /dev/null @@ -1,121 +0,0 @@ -{-# LANGUAGE BangPatterns #-} -module Data.Array.Repa.Internals.EvalReduction - ( foldS, foldP - , foldAllS, foldAllP) -where -import Data.Array.Repa.Internals.Elt -import Data.Array.Repa.Internals.Gang -import qualified Data.Vector.Unboxed as V -import qualified Data.Vector.Unboxed.Mutable as M -import GHC.Base ( quotInt, divInt ) - - --- | Sequential reduction of a multidimensional array along the innermost dimension. -foldS :: Elt a - => M.IOVector a -- ^ vector to write elements into - -> (Int -> a) -- ^ function to get an element from the given index - -> (a -> a -> a) -- ^ binary associative combination function - -> a -- ^ starting value (typically an identity) - -> Int -- ^ inner dimension (length to fold over) - -> IO () -{-# INLINE foldS #-} -foldS vec !f !c !r !n = iter 0 0 - where - !end = M.length vec - - {-# INLINE iter #-} - iter !sh !sz | sh >= end = return () - | otherwise = - let !next = sz + n - in M.unsafeWrite vec sh (reduce f c r sz next) >> iter (sh+1) next - - --- | Parallel reduction of a multidimensional array along the innermost dimension. --- Each output value is computed by a single thread, with the output values --- distributed evenly amongst the available threads. -foldP :: Elt a - => M.IOVector a -- ^ vector to write elements into - -> (Int -> a) -- ^ function to get an element from the given index - -> (a -> a -> a) -- ^ binary associative combination operator - -> a -- ^ starting value. Must be neutral with respect - -- ^ to the operator. eg @0 + a = a@. - -> Int -- ^ inner dimension (length to fold over) - -> IO () -{-# INLINE foldP #-} -foldP vec !f !c !r !n - = gangIO theGang - $ \tid -> fill (split tid) (split (tid+1)) - where - !threads = gangSize theGang - !len = M.length vec - !step = (len + threads - 1) `quotInt` threads - - {-# INLINE split #-} - split !ix = len `min` (ix * step) - - {-# INLINE fill #-} - fill !start !end = iter start (start * n) - where - {-# INLINE iter #-} - iter !sh !sz | sh >= end = return () - | otherwise = - let !next = sz + n - in M.unsafeWrite vec sh (reduce f c r sz next) >> iter (sh+1) next - - --- | Sequential reduction of all the elements in an array. -foldAllS :: Elt a - => (Int -> a) -- ^ function to get an element from the given index - -> (a -> a -> a) -- ^ binary associative combining function - -> a -- ^ starting value - -> Int -- ^ number of elements - -> IO a -{-# INLINE foldAllS #-} -foldAllS !f !c !r !len = return $! reduce f c r 0 len - - --- | Parallel tree reduction of an array to a single value. Each thread takes an --- equally sized chunk of the data and computes a partial sum. The main thread --- then reduces the array of partial sums to the final result. --- --- We don't require that the initial value be a neutral element, so each thread --- computes a fold1 on its chunk of the data, and the seed element is only --- applied in the final reduction step. --- -foldAllP :: Elt a - => (Int -> a) -- ^ function to get an element from the given index - -> (a -> a -> a) -- ^ binary associative combining function - -> a -- ^ starting value - -> Int -- ^ number of elements - -> IO a -{-# INLINE foldAllP #-} -foldAllP !f !c !r !len - | len == 0 = return r - | otherwise = do - mvec <- M.unsafeNew chunks - gangIO theGang $ \tid -> fill mvec tid (split tid) (split (tid+1)) - vec <- V.unsafeFreeze mvec - return $! V.foldl' c r vec - where - !threads = gangSize theGang - !step = (len + threads - 1) `quotInt` threads - chunks = ((len + step - 1) `divInt` step) `min` threads - - {-# INLINE split #-} - split !ix = len `min` (ix * step) - - {-# INLINE fill #-} - fill !mvec !tid !start !end - | start >= end = return () - | otherwise = M.unsafeWrite mvec tid (reduce f c (f start) (start+1) end) - - --- | Sequentially reduce values between the given indices -{-# INLINE reduce #-} -reduce :: (Int -> a) -> (a -> a -> a) -> a -> Int -> Int -> a -reduce !f !c !r !start !end = iter start r - where - {-# INLINE iter #-} - iter !i !z | i >= end = z - | otherwise = iter (i+1) (f i `c` z) - diff --git a/Data/Array/Repa/Internals/Forcing.hs b/Data/Array/Repa/Internals/Forcing.hs deleted file mode 100644 --- a/Data/Array/Repa/Internals/Forcing.hs +++ /dev/null @@ -1,215 +0,0 @@ -{-# LANGUAGE BangPatterns #-} -module Data.Array.Repa.Internals.Forcing - ( toVector - , toList - , force, forceWith - , force2, forceWith2) -where -import Data.Array.Repa.Internals.EvalChunked -import Data.Array.Repa.Internals.EvalCursored -import Data.Array.Repa.Internals.Elt -import Data.Array.Repa.Internals.Base -import Data.Array.Repa.Index -import Data.Array.Repa.Shape as S -import qualified Data.Vector.Unboxed as V -import qualified Data.Vector.Unboxed.Mutable as VM -import Data.Vector.Unboxed (Vector) -import System.IO.Unsafe - -stage = "Data.Array.Repa.Internals.Forcing" - - --- Conversions that also force the array ---------------------------------------------------------- --- | Convert an array to an unboxed `Data.Vector`, forcing it if required. --- The elements come out in row-major order. -toVector - :: (Shape sh, Elt a) - => Array sh a - -> Vector a -{-# INLINE toVector #-} -toVector arr - = case force arr of - Array _ [Region _ (GenManifest vec)] -> vec - _ -> error $ stage ++ ".toVector: force failed" - - --- | Convert an array to a list, forcing it if required. -toList :: (Shape sh, Elt a) - => Array sh a - -> [a] - -{-# INLINE toList #-} -toList arr - = V.toList $ toVector arr - - --- Forcing ---------------------------------------------------------------------------------------- --- | Force an array, so that it becomes `Manifest`. --- The array is split into linear chunks and each chunk evaluated in parallel. -force :: (Shape sh, Elt a) - => Array sh a -> Array sh a - -{-# INLINE [2] force #-} -force arr - = unsafePerformIO - $ do (sh, vec) <- forceIO arr - return $ sh `seq` vec `seq` - Array sh [Region RangeAll (GenManifest vec)] - - where forceIO arr' - = case arr' of - -- Don't force an already forced array. - Array sh [Region RangeAll (GenManifest vec)] - -> return (sh, vec) - - Array sh _ - -> do mvec <- VM.unsafeNew (S.size sh) - forceWith (VM.unsafeWrite mvec) arr' - vec <- V.unsafeFreeze mvec - return (sh, vec) - - --- | Force an array, passing elements to the provided update function. --- Provide something like @(Foreign.Ptr.pokeElemOff ptr)@ to write elements into a buffer. --- The array is split into linear chunks and each chunk is evaluated in parallel. -forceWith - :: (Shape sh, Elt a) - => (Int -> a -> IO ()) - -> Array sh a - -> IO () - -{-# INLINE [2] forceWith #-} -forceWith !update arr@(Array sh _) - = fillChunkedP - (S.size sh) - update - (\ix -> arr `unsafeIndex` fromIndex sh ix) - - --- | Force an array, so that it becomes `Manifest`. --- This forcing function is specialised for DIM2 arrays, and does blockwise filling. -force2 :: Elt a => Array DIM2 a -> Array DIM2 a -{-# INLINE [2] force2 #-} -force2 arr - = unsafePerformIO - $ do (sh, vec) <- forceIO2 arr - return $ sh `seq` vec `seq` - Array sh [Region RangeAll (GenManifest vec)] - - where forceIO2 arr' - = arr' `deepSeqArray` - case arr' of - -- Don't force an already forced array. - Array sh [Region RangeAll (GenManifest vec)] - -> return (sh, vec) - - -- Create a vector to hold the new array and load in the regions. - Array sh _ - -> do mvec <- VM.new (S.size sh) - forceWith2 (VM.unsafeWrite mvec) arr' - vec <- V.unsafeFreeze mvec - return (sh, vec) - - --- | Force an array, passing elements to the provided update function. --- Provide something like @(Foreign.Ptr.pokeElemOff ptr)@ to write elements into a buffer. --- This forcing function is specialised for DIM2 arrays, and does blockwise filling. -forceWith2 - :: Elt a - => (Int -> a -> IO ()) - -> Array DIM2 a - -> IO () - -{-# INLINE [2] forceWith2 #-} -forceWith2 !write arr - = arr `deepSeqArray` - case arr of - -- If the array is already manifest then copy it into the buffer. - -- We don't need a particular traversal order just for a copy. - Array _ [Region RangeAll (GenManifest _)] - -> forceWith write arr - - -- NOTE We must specialise this for common numbers of regions so that - -- we get fusion for them. If we just have the last case (arbitrary - -- region list) then the worker won't fuse with the filling / - -- evaluation code. - Array sh [r1] - -> do fillRegion2P write sh r1 - - Array sh [r1, r2] - -> do fillRegion2P write sh r1 - fillRegion2P write sh r2 - - Array sh regions - -> do mapM_ (fillRegion2P write sh) regions - - --- FillRegion2P ----------------------------------------------------------------------------------- --- | Fill an array region into a vector. --- This is specialised for DIM2 regions. --- The region is evaluated in parallel in a blockwise manner, where each block is --- evaluated independently and in a separate thread. For delayed or cursored regions --- access their source elements from the local neighbourhood, this specialised version --- should given better cache performance than plain `fillRegionP`. --- -fillRegion2P - :: Elt a - => (Int -> a -> IO ()) -- ^ Update function to write into result buffer - -> DIM2 -- ^ Extent of entire array. - -> Region DIM2 a -- ^ Region to fill. - -> IO () - -{-# INLINE [1] fillRegion2P #-} -fillRegion2P write sh@(_ :. height :. width) (Region range gen) - = write `seq` height `seq` width `seq` - case range of - RangeAll - -> fillRect2 write sh gen - (Rect (Z :. 0 :. 0) - (Z :. height - 1 :. width - 1)) - - RangeRects _ [r1] - -> do fillRect2 write sh gen r1 - - RangeRects _ [r1, r2] - -> do fillRect2 write sh gen r1 - fillRect2 write sh gen r2 - - RangeRects _ [r1, r2, r3] - -> do fillRect2 write sh gen r1 - fillRect2 write sh gen r2 - fillRect2 write sh gen r3 - - RangeRects _ [r1, r2, r3, r4] - -> do fillRect2 write sh gen r1 - fillRect2 write sh gen r2 - fillRect2 write sh gen r3 - fillRect2 write sh gen r4 - - RangeRects _ rects - -> mapM_ (fillRect2 write sh gen) rects - - --- | Fill a rectangle in a vector. -fillRect2 - :: Elt a - => (Int -> a -> IO ()) -- ^ Update function to write into result buffer - -> DIM2 -- ^ Extent of entire array. - -> Generator DIM2 a -- ^ Generator for array elements. - -> Rect DIM2 -- ^ Rectangle to fill. - -> IO () - -{-# INLINE fillRect2 #-} -fillRect2 write sh@(_ :. _ :. width) gen (Rect (Z :. y0 :. x0) (Z :. y1 :. x1)) - = write `seq` width `seq` y0 `seq` x0 `seq` y1 `seq` x1 `seq` - case gen of - GenManifest vec - -> fillCursoredBlock2P write - id addDim (\ix -> vec `V.unsafeIndex` toIndex sh ix) - width x0 y0 x1 y1 - - -- Cursor based arrays. - GenCursor makeCursor shiftCursor loadElem - -> fillCursoredBlock2P write - makeCursor shiftCursor loadElem - width x0 y0 x1 y1 diff --git a/Data/Array/Repa/Internals/Gang.hs b/Data/Array/Repa/Internals/Gang.hs deleted file mode 100644 --- a/Data/Array/Repa/Internals/Gang.hs +++ /dev/null @@ -1,249 +0,0 @@ -{-# LANGUAGE CPP #-} - --- | Gang Primitives. --- Based on DPH code by Roman Leshchinskiy --- --- Gang primitives. --- -#define TRACE_GANG 0 - -module Data.Array.Repa.Internals.Gang - ( Gang, seqGang, forkGang, gangSize, gangIO, gangST, traceGang, traceGangST - , theGang) -where -import GHC.IO -import GHC.ST -import GHC.Conc (forkOn) - -import Control.Concurrent.MVar -import Control.Exception (assert) - -import Control.Monad (zipWithM, zipWithM_) -import GHC.Conc (numCapabilities) -import System.IO - -#if TRACE_GANG -import GHC.Exts (traceEvent) -import System.Time ( ClockTime(..), getClockTime ) -#endif - --- TheGang ---------------------------------------------------------------------------------------- --- | The gang is shared by all computations. -theGang :: Gang -{-# NOINLINE theGang #-} -theGang = unsafePerformIO $ forkGang numCapabilities - - --- Requests --------------------------------------------------------------------------------------- --- | The 'Req' type encapsulates work requests for individual members of a gang. -data Req - -- | Instruct the worker to run the given action then signal it's done - -- by writing to the MVar. - = ReqDo (Int -> IO ()) (MVar ()) - - -- | Tell the worker that we're shutting the gang down. The worker should - -- signal that it's received the request down by writing to the MVar - -- before returning to its caller (forkGang) - | ReqShutdown (MVar ()) - - --- | Create a new request for the given action. -newReq :: (Int -> IO ()) -> IO Req -newReq p - = do mv <- newEmptyMVar - return $ ReqDo p mv - - --- | Block until a thread request has been executed. --- NOTE: only one thread can wait for the request. -waitReq :: Req -> IO () -waitReq req - = case req of - ReqDo _ varDone -> takeMVar varDone - ReqShutdown varDone -> takeMVar varDone - - --- Gang ------------------------------------------------------------------------------------------ --- | A 'Gang' is a group of threads which execute arbitrary work requests. --- To get the gang to do work, write Req-uest values to its MVars -data Gang - = Gang !Int -- Number of 'Gang' threads - [MVar Req] -- One 'MVar' per thread - (MVar Bool) -- Indicates whether the 'Gang' is busy - - -instance Show Gang where - showsPrec p (Gang n _ _) - = showString "<<" - . showsPrec p n - . showString " threads>>" - - --- | A sequential gang has no threads. -seqGang :: Gang -> Gang -seqGang (Gang n _ mv) = Gang n [] mv - - --- | The worker thread of a 'Gang'. --- The threads blocks on the MVar waiting for a work request. -gangWorker :: Int -> MVar Req -> IO () -gangWorker threadId varReq - = do traceGang $ "Worker " ++ show threadId ++ " waiting for request." - req <- takeMVar varReq - - case req of - ReqDo action varDone - -> do traceGang $ "Worker " ++ show threadId ++ " begin" - start <- getGangTime - action threadId - end <- getGangTime - traceGang $ "Worker " ++ show threadId ++ " end (" ++ diffTime start end ++ ")" - - putMVar varDone () - gangWorker threadId varReq - - ReqShutdown varDone - -> do traceGang $ "Worker " ++ show threadId ++ " shutting down." - putMVar varDone () - - --- | Finaliser for worker threads. --- We want to shutdown the corresponding thread when it's MVar becomes unreachable. --- Without this Repa programs can complain about "Blocked indefinitely on an MVar" --- because worker threads are still blocked on the request MVars when the program ends. --- Whether the finalizer is called or not is very racey. It happens about 1 in 10 runs --- when for the repa-edgedetect benchmark, and less often with the others. --- --- We're relying on the comment in System.Mem.Weak that says --- "If there are no other threads to run, the runtime system will check for runnable --- finalizers before declaring the system to be deadlocked." --- --- If we were creating and destroying the gang cleanly we wouldn't need this, but theGang --- is created with a top-level unsafePerformIO. Hacks beget hacks beget hacks... --- -finaliseWorker :: MVar Req -> IO () -finaliseWorker varReq - = do varDone <- newEmptyMVar - putMVar varReq (ReqShutdown varDone) - takeMVar varDone - return () - - --- | Fork a 'Gang' with the given number of threads (at least 1). -forkGang :: Int -> IO Gang -forkGang n - = assert (n > 0) - $ do - -- Create the vars we'll use to issue work requests. - mvs <- sequence . replicate n $ newEmptyMVar - - -- Add finalisers so we can shut the workers down cleanly if they become unreachable. - mapM_ (\var -> addMVarFinalizer var (finaliseWorker var)) mvs - - -- Create all the worker threads - zipWithM_ forkOn [0..] - $ zipWith gangWorker [0 .. n-1] mvs - - -- The gang is currently idle. - busy <- newMVar False - - return $ Gang n mvs busy - - --- | The number of threads in the 'Gang'. -gangSize :: Gang -> Int -gangSize (Gang n _ _) = n - - --- | Issue work requests for the 'Gang' and wait until they have been executed. --- If the gang is already busy then just run the action in the --- requesting thread. --- --- TODO: We might want to print a configurable warning that this is happening. --- -gangIO :: Gang - -> (Int -> IO ()) - -> IO () - -{-# NOINLINE gangIO #-} -gangIO (Gang n mvs busy) p - = do traceGang "gangIO: issuing work requests (SEQ_IF_GANG_BUSY)" - b <- swapMVar busy True - - traceGang $ "gangIO: gang is currently " ++ (if b then "busy" else "idle") - if b - then do - hPutStr stderr - $ unlines [ "Data.Array.Repa: Performing nested parallel computation sequentially." - , " You've probably called the 'force' function while another instance was" - , " already running. This can happen if the second version was suspended due" - , " to lazy evaluation. Use 'deepSeqArray' to ensure that each array is fully" - , " evaluated before you 'force' the next one." - , "" ] - - mapM_ p [0 .. n-1] - - else do - parIO n mvs p - _ <- swapMVar busy False - return () - - --- | Issue some requests to the worker threads and wait for them to complete. -parIO :: Int -- ^ Number of threads in the gang. - -> [MVar Req] -- ^ Request vars for worker threads. - -> (Int -> IO ()) -- ^ Action to run in all the workers, it's given the ix of - -- the particular worker thread it's running on. - -> IO () - -parIO n mvs p - = do traceGang "parIO: begin" - - start <- getGangTime - reqs <- sequence . replicate n $ newReq p - - traceGang "parIO: issuing requests" - _ <- zipWithM putMVar mvs reqs - - traceGang "parIO: waiting for requests to complete" - mapM_ waitReq reqs - end <- getGangTime - - traceGang $ "parIO: end " ++ diffTime start end - - --- | Same as 'gangIO' but in the 'ST' monad. -gangST :: Gang -> (Int -> ST s ()) -> ST s () -gangST g p = unsafeIOToST . gangIO g $ unsafeSTToIO . p - - --- Tracing ---------------------------------------------------------------------------------------- -#if TRACE_GANG -getGangTime :: IO Integer -getGangTime - = do TOD sec pico <- getClockTime - return (pico + sec * 1000000000000) - -diffTime :: Integer -> Integer -> String -diffTime x y = show (y-x) - -traceGang :: String -> IO () -traceGang s - = do t <- getGangTime - traceEvent $ show t ++ " @ " ++ s - -#else -getGangTime :: IO () -getGangTime = return () - -diffTime :: () -> () -> String -diffTime _ _ = "" - -traceGang :: String -> IO () -traceGang _ = return () - -#endif - -traceGangST :: String -> ST s () -traceGangST s = unsafeIOToST (traceGang s) - diff --git a/Data/Array/Repa/Internals/Select.hs b/Data/Array/Repa/Internals/Select.hs deleted file mode 100644 --- a/Data/Array/Repa/Internals/Select.hs +++ /dev/null @@ -1,118 +0,0 @@ -{-# LANGUAGE BangPatterns, ExplicitForAll, ScopedTypeVariables, PatternGuards #-} -module Data.Array.Repa.Internals.Select - (selectChunkedS, selectChunkedP) -where -import Data.Array.Repa.Internals.Gang -import Data.Array.Repa.Shape -import Data.Vector.Unboxed as V -import Data.Vector.Unboxed.Mutable as VM -import GHC.Base (remInt, quotInt) -import Prelude as P -import Control.Monad as P -import Data.IORef - --- | Select indices matching a predicate -selectChunkedS - :: (Shape sh, Unbox a) - => (sh -> Bool) -- ^ See if this predicate matches. - -> (sh -> a) -- ^ .. and apply fn to the matching index - -> IOVector a -- ^ .. then write the result into the vector. - -> sh -- ^ Extent of indices to apply to predicate. - -> IO Int -- ^ Number of elements written to destination array. - -{-# INLINE selectChunkedS #-} -selectChunkedS match produce !vDst !shSize - = fill 0 0 - where lenSrc = size shSize - lenDst = VM.length vDst - - fill !nSrc !nDst - | nSrc >= lenSrc = return nDst - | nDst >= lenDst = return nDst - - | ixSrc <- fromIndex shSize nSrc - , match ixSrc - = do VM.unsafeWrite vDst nDst (produce ixSrc) - fill (nSrc + 1) (nDst + 1) - - | otherwise - = fill (nSrc + 1) nDst - - --- | Select indices matching a predicate, in parallel. --- The array is chunked up, with one chunk being given to each thread. --- The number of elements in the result array depends on how many threads --- you're running the program with. -selectChunkedP - :: forall a - . Unbox a - => (Int -> Bool) -- ^ See if this predicate matches. - -> (Int -> a) -- .. and apply fn to the matching index - -> Int -- Extent of indices to apply to predicate. - -> IO [IOVector a] -- Chunks containing array elements. - -{-# INLINE selectChunkedP #-} -selectChunkedP !match !produce !len - = do - -- Make IORefs that the threads will write their result chunks to. - -- We start with a chunk size proportial to the number of threads we have, - -- but the threads themselves can grow the chunks if they run out of space. - refs <- P.replicateM threads - $ do vec <- VM.new $ len `div` threads - newIORef vec - - -- Fire off a thread to fill each chunk. - gangIO theGang - $ \thread -> makeChunk (refs !! thread) - (splitIx thread) - (splitIx (thread + 1) - 1) - - -- Read the result chunks back from the IORefs. - -- If a thread had to grow a chunk, then these might not be the same ones - -- we created back in the first step. - P.mapM readIORef refs - - where -- See how many threads we have available. - !threads = gangSize theGang - !chunkLen = len `quotInt` threads - !chunkLeftover = len `remInt` threads - - - -- Decide where to split the source array. - {-# INLINE splitIx #-} - splitIx thread - | thread < chunkLeftover = thread * (chunkLen + 1) - | otherwise = thread * chunkLen + chunkLeftover - - - -- Fill the given chunk with elements selected from this range of indices. - makeChunk :: IORef (IOVector a) -> Int -> Int -> IO () - makeChunk !ref !ixSrc !ixSrcEnd - = do vecDst <- VM.new (len `div` threads) - vecDst' <- fillChunk ixSrc ixSrcEnd vecDst 0 (VM.length vecDst - 1) - writeIORef ref vecDst' - - - -- The main filling loop. - fillChunk :: Int -> Int -> IOVector a -> Int -> Int -> IO (IOVector a) - fillChunk !ixSrc !ixSrcEnd !vecDst !ixDst !ixDstEnd - -- If we've finished selecting elements, then slice the vector down - -- so it doesn't have any empty space at the end. - | ixSrc >= ixSrcEnd - = return $ VM.slice 0 ixDst vecDst - - -- If we've run out of space in the chunk then grow it some more. - | ixDst >= ixDstEnd - = do let ixDstEnd' = VM.length vecDst * 2 - 1 - vecDst' <- VM.grow vecDst (ixDstEnd + 1) - fillChunk (ixSrc + 1) ixSrcEnd vecDst' (ixDst + 1) ixDstEnd' - - -- We've got a maching element, so add it to the chunk. - | match ixSrc - = do VM.unsafeWrite vecDst ixDst (produce ixSrc) - fillChunk (ixSrc + 1) ixSrcEnd vecDst (ixDst + 1) ixDstEnd - - -- The element doesnt match, so keep going. - | otherwise - = fillChunk (ixSrc + 1) ixSrcEnd vecDst ixDst ixDstEnd - diff --git a/Data/Array/Repa/Operators/IndexSpace.hs b/Data/Array/Repa/Operators/IndexSpace.hs --- a/Data/Array/Repa/Operators/IndexSpace.hs +++ b/Data/Array/Repa/Operators/IndexSpace.hs @@ -1,4 +1,3 @@ -{-# OPTIONS_HADDOCK hide #-} {-# LANGUAGE TypeOperators, ExplicitForAll, FlexibleContexts #-} module Data.Array.Repa.Operators.IndexSpace @@ -7,61 +6,49 @@ , transpose , extend , slice - , backpermute - , backpermuteDft) + , backpermute, unsafeBackpermute + , backpermuteDft, unsafeBackpermuteDft) where import Data.Array.Repa.Index import Data.Array.Repa.Slice -import Data.Array.Repa.Internals.Elt -import Data.Array.Repa.Internals.Base -import Data.Array.Repa.Operators.Traverse +import Data.Array.Repa.Base +import Data.Array.Repa.Repr.Delayed +import Data.Array.Repa.Operators.Traversal import Data.Array.Repa.Shape as S import Prelude hiding ((++)) import qualified Prelude as P stage = "Data.Array.Repa.Operators.IndexSpace" --- Index space transformations -------------------------------------------------------------------- +-- Index space transformations ------------------------------------------------ -- | Impose a new shape on the elements of an array. -- The new extent must be the same size as the original, else `error`. --- --- TODO: This only works for arrays with a single region. --- -reshape :: (Shape sh, Shape sh', Elt a) - => sh' - -> Array sh a - -> Array sh' a - -{-# INLINE reshape #-} -reshape sh' arr - | not $ S.size sh' == S.size (extent arr) - = error $ stage P.++ ".reshape: reshaped array will not match size of the original" - -reshape sh' (Array sh [Region RangeAll gen]) - = Array sh' [Region RangeAll gen'] - where gen' = case gen of - GenManifest vec - -> GenManifest vec - - GenCursor makeCursor _ loadElem - -> GenCursor - id - addDim - (loadElem . makeCursor . fromIndex sh . toIndex sh') +reshape :: (Shape sh2, Shape sh1 + , Repr r1 e) + => sh2 + -> Array r1 sh1 e + -> Array D sh2 e -reshape _ _ - = error $ stage P.++ ".reshape: can't reshape a partitioned array" +{-# INLINE [3] reshape #-} +reshape sh2 arr + | not $ S.size sh2 == S.size (extent arr) + = error + $ stage P.++ ".reshape: reshaped array will not match size of the original" +reshape sh2 arr + = fromFunction sh2 + $ unsafeIndex arr . fromIndex (extent arr) . toIndex sh2 + -- | Append two arrays. --- append, (++) - :: (Shape sh, Elt a) - => Array (sh :. Int) a - -> Array (sh :. Int) a - -> Array (sh :. Int) a + :: ( Shape sh + , Repr r1 e, Repr r2 e) + => Array r1 (sh :. Int) e + -> Array r2 (sh :. Int) e + -> Array D (sh :. Int) e -{-# INLINE append #-} +{-# INLINE [3] append #-} append arr1 arr2 = unsafeTraverse2 arr1 arr2 fnExtent fnElem where @@ -81,11 +68,12 @@ -- | Transpose the lowest two dimensions of an array. -- Transposing an array twice yields the original. transpose - :: (Shape sh, Elt a) - => Array (sh :. Int :. Int) a - -> Array (sh :. Int :. Int) a + :: ( Shape sh + , Repr r e) + => Array r (sh :. Int :. Int) e + -> Array D (sh :. Int :. Int) e -{-# INLINE transpose #-} +{-# INLINE [3] transpose #-} transpose arr = unsafeTraverse arr (\(sh :. m :. n) -> (sh :. n :.m)) @@ -93,19 +81,18 @@ -- | Extend an array, according to a given slice specification. --- (used to be called replicate). extend :: ( Slice sl , Shape (FullShape sl) , Shape (SliceShape sl) - , Elt e) + , Repr r e) => sl - -> Array (SliceShape sl) e - -> Array (FullShape sl) e + -> Array r (SliceShape sl) e + -> Array D (FullShape sl) e -{-# INLINE extend #-} +{-# INLINE [3] extend #-} extend sl arr - = backpermute + = unsafeBackpermute (fullOfSlice sl (extent arr)) (sliceOfFull sl) arr @@ -114,14 +101,14 @@ slice :: ( Slice sl , Shape (FullShape sl) , Shape (SliceShape sl) - , Elt e) - => Array (FullShape sl) e + , Repr r e) + => Array r (FullShape sl) e -> sl - -> Array (SliceShape sl) e + -> Array D (SliceShape sl) e -{-# INLINE slice #-} +{-# INLINE [3] slice #-} slice arr sl - = backpermute + = unsafeBackpermute (sliceOfFull sl (extent arr)) (fullOfSlice sl) arr @@ -129,37 +116,51 @@ -- | Backwards permutation of an array's elements. -- The result array has the same extent as the original. -backpermute - :: forall sh sh' a - . (Shape sh, Shape sh', Elt a) - => sh' -- ^ Extent of result array. - -> (sh' -> sh) -- ^ Function mapping each index in the result array - -- to an index of the source array. - -> Array sh a -- ^ Source array. - -> Array sh' a +backpermute, unsafeBackpermute + :: forall r sh1 sh2 e + . ( Shape sh1, Shape sh2 + , Repr r e) + => sh2 -- ^ Extent of result array. + -> (sh2 -> sh1) -- ^ Function mapping each index in the result array + -- to an index of the source array. + -> Array r sh1 e -- ^ Source array. + -> Array D sh2 e -{-# INLINE backpermute #-} +{-# INLINE [3] backpermute #-} backpermute newExtent perm arr = traverse arr (const newExtent) (. perm) +{-# INLINE [3] unsafeBackpermute #-} +unsafeBackpermute newExtent perm arr + = unsafeTraverse arr (const newExtent) (. perm) + -- | Default backwards permutation of an array's elements. -- If the function returns `Nothing` then the value at that index is taken -- from the default array (@arrDft@) -backpermuteDft - :: forall sh sh' a - . (Shape sh, Shape sh', Elt a) - => Array sh' a -- ^ Default values (@arrDft@) - -> (sh' -> Maybe sh) -- ^ Function mapping each index in the result array - -- to an index in the source array. - -> Array sh a -- ^ Source array. - -> Array sh' a +backpermuteDft, unsafeBackpermuteDft + :: forall r0 r1 sh1 sh2 e + . ( Shape sh1, Shape sh2 + , Repr r0 e, Repr r1 e) + => Array r0 sh2 e -- ^ Default values (@arrDft@) + -> (sh2 -> Maybe sh1) -- ^ Function mapping each index in the result array + -- to an index in the source array. + -> Array r1 sh1 e -- ^ Source array. + -> Array D sh2 e -{-# INLINE backpermuteDft #-} +{-# INLINE [3] backpermuteDft #-} backpermuteDft arrDft fnIndex arrSrc = fromFunction (extent arrDft) fnElem where fnElem ix = case fnIndex ix of - Just ix' -> arrSrc ! ix' - Nothing -> arrDft ! ix + Just ix' -> arrSrc `index` ix' + Nothing -> arrDft `index` ix + +{-# INLINE [3] unsafeBackpermuteDft #-} +unsafeBackpermuteDft arrDft fnIndex arrSrc + = fromFunction (extent arrDft) fnElem + where fnElem ix + = case fnIndex ix of + Just ix' -> arrSrc `unsafeIndex` ix' + Nothing -> arrDft `unsafeIndex` ix diff --git a/Data/Array/Repa/Operators/Interleave.hs b/Data/Array/Repa/Operators/Interleave.hs --- a/Data/Array/Repa/Operators/Interleave.hs +++ b/Data/Array/Repa/Operators/Interleave.hs @@ -1,5 +1,4 @@ -{-# OPTIONS_HADDOCK hide #-} -{-# LANGUAGE TypeOperators, PatternGuards #-} +{-# LANGUAGE TypeOperators, ExplicitForAll, FlexibleContexts #-} module Data.Array.Repa.Operators.Interleave ( interleave2 @@ -7,11 +6,13 @@ , interleave4) where import Data.Array.Repa.Index -import Data.Array.Repa.Internals.Elt -import Data.Array.Repa.Internals.Base -import Data.Array.Repa.Operators.Traverse -import Data.Array.Repa.Shape as S +import Data.Array.Repa.Base +import Data.Array.Repa.Repr.Delayed +import Data.Array.Repa.Operators.Traversal +import Data.Array.Repa.Shape as S +import Prelude hiding ((++)) +-- Interleave ----------------------------------------------------------------- -- | Interleave the elements of two arrays. -- All the input arrays must have the same extent, else `error`. -- The lowest dimension of the result array is twice the size of the inputs. @@ -22,12 +23,13 @@ -- @ -- interleave2 - :: (Shape sh, Elt a) - => Array (sh :. Int) a - -> Array (sh :. Int) a - -> Array (sh :. Int) a + :: (Shape sh + , Repr r1 a, Repr r2 a) + => Array r1 (sh :. Int) a + -> Array r2 (sh :. Int) a + -> Array D (sh :. Int) a -{-# INLINE interleave2 #-} +{-# INLINE [3] interleave2 #-} interleave2 arr1 arr2 = arr1 `deepSeqArray` arr2 `deepSeqArray` unsafeTraverse2 arr1 arr2 shapeFn elemFn @@ -49,13 +51,14 @@ -- | Interleave the elements of three arrays. interleave3 - :: (Shape sh, Elt a) - => Array (sh :. Int) a - -> Array (sh :. Int) a - -> Array (sh :. Int) a - -> Array (sh :. Int) a + :: ( Shape sh + , Repr r1 a, Repr r2 a, Repr r3 a) + => Array r1 (sh :. Int) a + -> Array r2 (sh :. Int) a + -> Array r3 (sh :. Int) a + -> Array D (sh :. Int) a -{-# INLINE interleave3 #-} +{-# INLINE [3] interleave3 #-} interleave3 arr1 arr2 arr3 = arr1 `deepSeqArray` arr2 `deepSeqArray` arr3 `deepSeqArray` unsafeTraverse3 arr1 arr2 arr3 shapeFn elemFn @@ -79,14 +82,15 @@ -- | Interleave the elements of four arrays. interleave4 - :: (Shape sh, Elt a) - => Array (sh :. Int) a - -> Array (sh :. Int) a - -> Array (sh :. Int) a - -> Array (sh :. Int) a - -> Array (sh :. Int) a + :: ( Shape sh + , Repr r1 a, Repr r2 a, Repr r3 a, Repr r4 a) + => Array r1 (sh :. Int) a + -> Array r2 (sh :. Int) a + -> Array r3 (sh :. Int) a + -> Array r4 (sh :. Int) a + -> Array D (sh :. Int) a -{-# INLINE interleave4 #-} +{-# INLINE [3] interleave4 #-} interleave4 arr1 arr2 arr3 arr4 = arr1 `deepSeqArray` arr2 `deepSeqArray` arr3 `deepSeqArray` arr4 `deepSeqArray` unsafeTraverse4 arr1 arr2 arr3 arr4 shapeFn elemFn @@ -108,4 +112,3 @@ 2 -> get3 (sh :. ix `div` 4) 3 -> get4 (sh :. ix `div` 4) _ -> error "Data.Array.Repa.interleave4: this never happens :-P" - diff --git a/Data/Array/Repa/Operators/Mapping.hs b/Data/Array/Repa/Operators/Mapping.hs --- a/Data/Array/Repa/Operators/Mapping.hs +++ b/Data/Array/Repa/Operators/Mapping.hs @@ -1,98 +1,57 @@ -{-# OPTIONS_HADDOCK hide #-} -{-# LANGUAGE NoMonomorphismRestriction, PatternGuards #-} +{-# LANGUAGE FunctionalDependencies, UndecidableInstances #-} module Data.Array.Repa.Operators.Mapping - ( map - , zipWith - , (+^) - , (-^) - , (*^) - , (/^)) + ( -- * Generic maps + map + , zipWith + , (+^), (-^), (*^), (/^) + + -- * Combining maps + , Combine(..)) where -import Data.Array.Repa.Internals.Elt -import Data.Array.Repa.Internals.Base -import Data.Array.Repa.Shape as S -import qualified Data.Vector.Unboxed as V -import qualified Prelude as P -import Prelude (($), (.), (+), (*), (+), (/), (-)) +import Data.Array.Repa.Shape +import Data.Array.Repa.Base +import Data.Array.Repa.Repr.ByteString +import Data.Array.Repa.Repr.Cursored +import Data.Array.Repa.Repr.Delayed +import Data.Array.Repa.Repr.ForeignPtr +import Data.Array.Repa.Repr.Partitioned +import Data.Array.Repa.Repr.Unboxed +import Data.Array.Repa.Repr.Undefined +import Prelude hiding (map, zipWith) +import Foreign.Storable +import Data.Word --- | Apply a worker function to each element of an array, yielding a new array with the same extent. --- --- This is specialised for arrays of up to four regions, using more breaks fusion. +-- | Apply a worker function to each element of an array, +-- yielding a new array with the same extent. -- -map :: (Shape sh, Elt a, Elt b) - => (a -> b) - -> Array sh a - -> Array sh b - -{-# INLINE map #-} -map f (Array sh regions) - = Array sh (mapRegions regions) - - where {-# INLINE mapRegions #-} - mapRegions rs - = case rs of - [] -> [] - [r] -> [mapRegion r] - [r1, r2] -> [mapRegion r1, mapRegion r2] - [r1, r2, r3] -> [mapRegion r1, mapRegion r2, mapRegion r3] - [r1, r2, r3, r4] -> [mapRegion r1, mapRegion r2, mapRegion r3, mapRegion r4] - _ -> mapRegions' rs - - mapRegions' rs - = case rs of - [] -> [] - (r : rs') -> mapRegion r : mapRegions' rs' - - {-# INLINE mapRegion #-} - mapRegion (Region range gen) - = Region range (mapGen gen) - - {-# INLINE mapGen #-} - mapGen gen - = case gen of - GenManifest vec - -> GenCursor - P.id - addDim - (\ix -> f $ V.unsafeIndex vec $ S.toIndex sh ix) - - GenCursor makeCursor shiftCursor loadElem - -> GenCursor makeCursor shiftCursor (f . loadElem) +map :: (Shape sh, Repr r a) + => (a -> b) -> Array r sh a -> Array D sh b +{-# INLINE [4] map #-} +map f arr + = case delay arr of + ADelayed sh g -> ADelayed sh (f . g) +-- ZipWith -------------------------------------------------------------------- -- | Combine two arrays, element-wise, with a binary operator. -- If the extent of the two array arguments differ, -- then the resulting array's extent is their intersection. -- -zipWith :: (Shape sh, Elt a, Elt b, Elt c) - => (a -> b -> c) - -> Array sh a - -> Array sh b - -> Array sh c - -{-# INLINE zipWith #-} +zipWith :: (Shape sh, Repr r1 a, Repr r2 b) + => (a -> b -> c) + -> Array r1 sh a -> Array r2 sh b + -> Array D sh c +{-# INLINE [3] zipWith #-} zipWith f arr1 arr2 - | Array sh2 [_] <- arr1 - , Array sh1 [ Region g21 (GenCursor make21 _ load21) - , Region g22 (GenCursor make22 _ load22)] <- arr2 - - = let {-# INLINE load21' #-} - load21' ix = f (arr1 `unsafeIndex` ix) (load21 $ make21 ix) - - {-# INLINE load22' #-} - load22' ix = f (arr1 `unsafeIndex` ix) (load22 $ make22 ix) - - in Array (S.intersectDim sh1 sh2) - [ Region g21 (GenCursor P.id addDim load21') - , Region g22 (GenCursor P.id addDim load22') ] + = arr1 `deepSeqArray` arr2 `deepSeqArray` + let + {-# INLINE get #-} + get ix = f (arr1 `unsafeIndex` ix) (arr2 `unsafeIndex` ix) - | P.otherwise - = let {-# INLINE getElem' #-} - getElem' ix = f (arr1 `unsafeIndex` ix) (arr2 `unsafeIndex` ix) - in fromFunction - (S.intersectDim (extent arr1) (extent arr2)) - getElem' + in fromFunction + (intersectDim (extent arr1) (extent arr2)) + get {-# INLINE (+^) #-} @@ -107,3 +66,109 @@ {-# INLINE (/^) #-} (/^) = zipWith (/) + + +-- Combine -------------------------------------------------------------------- +-- | Combining versions of @map@ and @zipWith@ that preserve the representation +-- of cursored and partitioned arrays. +-- +-- For cursored (@C@) arrays, the cursoring of the source array is preserved. +-- +-- For partitioned (@P@) arrays, the worker function is fused with each array +-- partition separately, instead of treating the whole array as a single +-- bulk object. +-- +-- Preserving the cursored and\/or paritioned representation of an array +-- is will make follow-on computation more efficient than if the array was +-- converted to a vanilla Delayed (@D@) array as with plain `map` and `zipWith`. +-- +-- If the source array is not cursored or partitioned then `cmap` and +-- `czipWith` are identical to the plain functions. +-- +class Combine r1 a r2 b | r1 -> r2 where + + -- | Combining @map@. + cmap :: Shape sh + => (a -> b) + -> Array r1 sh a + -> Array r2 sh b + + -- | Combining @zipWith@. + -- If you have a cursored or partitioned source array then use that as + -- the third argument (corresponding to @r1@ here) + czipWith + :: (Shape sh, Repr r c) + => (c -> a -> b) + -> Array r sh c + -> Array r1 sh a + -> Array r2 sh b + + +-- ByteString ------------------------- +instance Combine B Word8 D b where + cmap = map + czipWith = zipWith + + +-- Cursored --------------------------- +instance Combine C a C b where + {-# INLINE [4] cmap #-} + cmap f (ACursored sh makec shiftc loadc) + = ACursored sh makec shiftc (f . loadc) + + {-# INLINE [3] czipWith #-} + czipWith f arr1 (ACursored sh makec shiftc loadc) + = let {-# INLINE makec' #-} + makec' ix = (ix, makec ix) + + {-# INLINE shiftc' #-} + shiftc' off (ix, cur) = (addDim off ix, shiftc off cur) + + {-# INLINE load' #-} + load' (ix, cur) = f (arr1 `unsafeIndex` ix) (loadc cur) + + in ACursored + (intersectDim (extent arr1) sh) + makec' shiftc' load' + + +-- Delayed ---------------------------- +instance Combine D a D b where + cmap = map + czipWith = zipWith + + +-- ForeignPtr ------------------------- +instance Storable a => Combine F a D b where + cmap = map + czipWith = zipWith + + +-- Partitioned ------------------------ +instance (Combine r11 a r21 b + , Combine r12 a r22 b) + => Combine (P r11 r12) a (P r21 r22) b where + + {-# INLINE [4] cmap #-} + cmap f (APart sh range arr1 arr2) + = APart sh range (cmap f arr1) (cmap f arr2) + + {-# INLINE [3] czipWith #-} + czipWith f arr1 (APart sh range arr21 arr22) + = APart sh range (czipWith f arr1 arr21) + (czipWith f arr1 arr22) + + +-- Unboxed ---------------------------- +instance Unbox a => Combine U a D b where + cmap = map + czipWith = zipWith + + +-- Undefined -------------------------- +instance Combine X a D b where + cmap = map + czipWith = zipWith + + + diff --git a/Data/Array/Repa/Operators/Modify.hs b/Data/Array/Repa/Operators/Modify.hs deleted file mode 100644 --- a/Data/Array/Repa/Operators/Modify.hs +++ /dev/null @@ -1,53 +0,0 @@ -{-# OPTIONS_HADDOCK hide #-} - -module Data.Array.Repa.Operators.Modify - ( -- * Bulk updates - (//)) -where -import Data.Array.Repa.Shape -import Data.Array.Repa.Internals.Elt -import Data.Array.Repa.Internals.Base - -{- -stage :: String -stage = "Data.Array.Repa.Operators.Modify" --} - --- Bulk updates ---------------------------------------------------------------- --- | For each pair @(sh, a)@ from the list of index/value pairs, replace the --- element at position @sh@ by @a@. --- --- > update <5,9,2,7> [(2,1),(0,3),(2,8)] = <3,9,8,7> --- -{-# INLINE (//) #-} -(//) :: (Shape sh, Elt a) => Array sh a -> [(sh,a)] -> Array sh a -(//) arr us - = fromFunction - (extent arr) - (\sh -> case lookup sh us of - Just a -> a - Nothing -> index arr sh) - -{- --- For each pair @(sh, a)@ from the array of index/value pairs, replace the --- element at position @sh@ by @a@. --- --- > update <5,9,2,7> <(2,1),(0,3),(2,8)> = <3,9,8,7> --- -{-# INLINE update #-} -update :: Shape sh - => Array sh a -- ^ initial array - -> Array sh (sh, a) -- ^ array of shape/value pairs - -> Array sh a -update _arr _us = error $ stage ++ ".update: not defined yet" - - --- Same as 'update', but without bounds checks --- -{-# INLINE unsafeUpdate #-} -unsafeUpdate :: Shape sh - => Array sh a - -> Array sh (sh, a) - -> Array sh a -unsafeUpdate _arr _us = error $ stage ++ ".unsafeUpdate: not defined yet" --} diff --git a/Data/Array/Repa/Operators/Reduction.hs b/Data/Array/Repa/Operators/Reduction.hs --- a/Data/Array/Repa/Operators/Reduction.hs +++ b/Data/Array/Repa/Operators/Reduction.hs @@ -1,82 +1,141 @@ -{-# OPTIONS_HADDOCK hide #-} -{-# LANGUAGE BangPatterns, ExplicitForAll, TypeOperators #-} +{-# LANGUAGE BangPatterns, ExplicitForAll, TypeOperators, MagicHash #-} module Data.Array.Repa.Operators.Reduction - ( fold, foldAll - , sum, sumAll) + ( foldS, foldP + , foldAllS, foldAllP + , sumS, sumP + , sumAllS, sumAllP) where +import Data.Array.Repa.Base import Data.Array.Repa.Index -import Data.Array.Repa.Internals.Elt -import Data.Array.Repa.Internals.Base +import Data.Array.Repa.Eval.Elt +import Data.Array.Repa.Repr.Unboxed import Data.Array.Repa.Shape as S import qualified Data.Vector.Unboxed as V import qualified Data.Vector.Unboxed.Mutable as M import Prelude hiding (sum) - -import Data.Array.Repa.Internals.EvalReduction +import qualified Data.Array.Repa.Eval.Reduction as E import System.IO.Unsafe +import GHC.Exts +-- foldS ---------------------------------------------------------------------- +-- | Sequential reduction of the innermost dimension of an arbitrary rank array. +-- +-- Combine this with `transpose` to fold any other dimension. +foldS :: (Shape sh, Elt a, Unbox a, Repr r a) + => (a -> a -> a) + -> a + -> Array r (sh :. Int) a + -> Array U sh a +{-# INLINE [2] foldS #-} +foldS f z arr + = let sh@(sz :. n') = extent arr + !(I# n) = n' + in unsafePerformIO + $ do mvec <- M.unsafeNew (S.size sz) + E.foldS mvec (\ix -> arr `unsafeIndex` fromIndex sh (I# ix)) f z n + !vec <- V.unsafeFreeze mvec + return $ fromUnboxed sz vec --- | Reduction of the innermost dimension of an arbitrary rank array. The first --- argument needs to be an /associative/ operator. The starting element must --- be neutral with respect to the operator, for example @0@ is neutral with --- respect to @(+)@ as @0 + a = a@. These restrictions are required to support --- parallel evaluation, as the starting element may be used multiple --- times depending on the number of threads. --- Combine this with `transpose` to fold any other dimension. -fold :: (Shape sh, Elt a) +-- | Parallel reduction of the innermost dimension of an arbitray rank array. +-- +-- The first argument needs to be an associative sequential operator. +-- The starting element must be neutral with respect to the operator, for +-- example @0@ is neutral with respect to @(+)@ as @0 + a = a@. +-- These restrictions are required to support parallel evaluation, as the +-- starting element may be used multiple times depending on the number of threads. +foldP :: (Shape sh, Elt a, Unbox a, Repr r a) => (a -> a -> a) -> a - -> Array (sh :. Int) a - -> Array sh a -{-# INLINE [1] fold #-} -fold f z arr + -> Array r (sh :. Int) a + -> Array U sh a +{-# INLINE [2] foldP #-} +foldP f z arr = let sh@(sz :. n) = extent arr in case rank sh of - -- specialise rank-1 arrays, else one thread does all the work. We can't - -- match against the shape constructor, otherwise type error: (sz ~ Z) + -- specialise rank-1 arrays, else one thread does all the work. + -- We can't match against the shape constructor, + -- otherwise type error: (sz ~ Z) -- - 1 -> let !x = V.singleton $ foldAll f z arr - in Array sz [Region RangeAll (GenManifest x)] + 1 -> let !vec = V.singleton $ foldAllP f z arr + in fromUnboxed sz vec _ -> unsafePerformIO $ do mvec <- M.unsafeNew (S.size sz) - foldP mvec (\ix -> arr `unsafeIndex` fromIndex sh ix) f z n + E.foldP mvec (\ix -> arr `unsafeIndex` fromIndex sh ix) f z n !vec <- V.unsafeFreeze mvec - return $ Array sz [Region RangeAll (GenManifest vec)] + return $ fromUnboxed sz vec --- | Reduction of an array of arbitrary rank to a single scalar value. The first --- argument needs to be an /associative/ operator. The starting element must --- be neutral with respect to the operator, for example @0@ is neutral with --- respect to @(+)@ as @0 + a = a@. These restrictions are required to support --- parallel evaluation, as the starting element may be used multiple --- times depending on the number of threads. -foldAll :: (Shape sh, Elt a) +-- foldAll -------------------------------------------------------------------- +-- | Sequential reduction of an array of arbitrary rank to a single scalar value. +-- +foldAllS :: (Shape sh, Elt a, Unbox a, Repr r a) => (a -> a -> a) -> a - -> Array sh a + -> Array r sh a -> a -{-# INLINE [1] foldAll #-} -foldAll f z arr +{-# INLINE [2] foldAllS #-} +foldAllS f z arr + = arr `deepSeqArray` + let !ex = extent arr + !(I# n) = size ex + in E.foldAllS + (\ix -> arr `unsafeIndex` fromIndex ex (I# ix)) + f z n + + +-- | Parallel reduction of an array of arbitrary rank to a single scalar value. +-- +-- The first argument needs to be an associative sequential operator. +-- The starting element must be neutral with respect to the operator, +-- for example @0@ is neutral with respect to @(+)@ as @0 + a = a@. +-- These restrictions are required to support parallel evaluation, as the +-- starting element may be used multiple times depending on the number of threads. +foldAllP :: (Shape sh, Elt a, Unbox a, Repr r a) + => (a -> a -> a) + -> a + -> Array r sh a + -> a +{-# INLINE [2] foldAllP #-} +foldAllP f z arr = let sh = extent arr n = size sh - in unsafePerformIO $ foldAllP (\ix -> arr `unsafeIndex` fromIndex sh ix) f z n + in unsafePerformIO + $ E.foldAllP (\ix -> arr `unsafeIndex` fromIndex sh ix) f z n --- | Sum the innermost dimension of an array. -sum :: (Shape sh, Elt a, Num a) - => Array (sh :. Int) a - -> Array sh a -{-# INLINE sum #-} -sum arr = fold (+) 0 arr +-- sum ------------------------------------------------------------------------ +-- | Sequential sum the innermost dimension of an array. +sumS :: (Shape sh, Num a, Elt a, Unbox a, Repr r a) + => Array r (sh :. Int) a + -> Array U sh a +{-# INLINE [4] sumS #-} +sumS = foldS (+) 0 --- | Sum all the elements of an array. -sumAll :: (Shape sh, Elt a, Num a) - => Array sh a +-- | Sequential sum the innermost dimension of an array. +sumP :: (Shape sh, Num a, Elt a, Unbox a, Repr r a) + => Array r (sh :. Int) a + -> Array U sh a +{-# INLINE [4] sumP #-} +sumP = foldP (+) 0 + + +-- sumAll --------------------------------------------------------------------- +-- | Sequential sum of all the elements of an array. +sumAllS :: (Shape sh, Elt a, Unbox a, Num a, Repr r a) + => Array r sh a -> a -{-# INLINE sumAll #-} -sumAll arr = foldAll (+) 0 arr +{-# INLINE [4] sumAllS #-} +sumAllS = foldAllS (+) 0 + + +-- | Parallel sum all the elements of an array. +sumAllP :: (Shape sh, Elt a, Unbox a, Num a, Repr r a) + => Array r sh a + -> a +{-# INLINE [4] sumAllP #-} +sumAllP = foldAllP (+) 0 diff --git a/Data/Array/Repa/Operators/Select.hs b/Data/Array/Repa/Operators/Select.hs deleted file mode 100644 --- a/Data/Array/Repa/Operators/Select.hs +++ /dev/null @@ -1,44 +0,0 @@ -{-# OPTIONS_HADDOCK hide #-} -{-# LANGUAGE BangPatterns #-} - -module Data.Array.Repa.Operators.Select - (select) -where -import Data.Array.Repa.Index -import Data.Array.Repa.Internals.Elt -import Data.Array.Repa.Internals.Base -import Data.Array.Repa.Internals.Select -import qualified Data.Vector.Unboxed as V -import System.IO.Unsafe - - --- | Produce an array by applying a predicate to a range of integers. --- If the predicate matches, then use the second function to generate --- the element. --- --- This is a low-level function helpful for writing filtering operations on arrays. --- Use the integer as the index into the array you're filtering. --- -select :: Elt a - => (Int -> Bool) -- ^ If the Int matches this predicate, - -> (Int -> a) -- ^ ... then pass it to this fn to produce a value - -> Int -- ^ Range between 0 and this maximum. - -> Array DIM1 a -- ^ Array containing produced values. - -{-# INLINE select #-} -select match produce len - = unsafePerformIO - $ do (sh, vec) <- selectIO - return $ sh `seq` vec `seq` - Array sh [Region RangeAll (GenManifest vec)] - - where {-# INLINE selectIO #-} - selectIO - = do vecs <- selectChunkedP match produce len - vecs' <- mapM V.unsafeFreeze vecs - - -- TODO: avoid copy. - let result = V.concat vecs' - - return (Z :. V.length result, result) - diff --git a/Data/Array/Repa/Operators/Selection.hs b/Data/Array/Repa/Operators/Selection.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Operators/Selection.hs @@ -0,0 +1,43 @@ +{-# LANGUAGE BangPatterns #-} +module Data.Array.Repa.Operators.Selection + (select) +where +import Data.Array.Repa.Index +import Data.Array.Repa.Base +import Data.Array.Repa.Eval.Selection +import Data.Array.Repa.Repr.Unboxed as U +import qualified Data.Vector.Unboxed as V +import System.IO.Unsafe + + +-- | Produce an array by applying a predicate to a range of integers. +-- If the predicate matches, then use the second function to generate +-- the element. +-- +-- * This is a low-level function helpful for writing filtering +-- operations on arrays. +-- +-- * Use the integer as the index into the array you're filtering. +-- +select :: Unbox a + => (Int -> Bool) -- ^ If the Int matches this predicate, + -> (Int -> a) -- ^ ... then pass it to this fn to produce a value + -> Int -- ^ Range between 0 and this maximum. + -> Array U DIM1 a -- ^ Array containing produced values. + +{-# INLINE [2] select #-} +select match produce len + = unsafePerformIO + $ do (sh, vec) <- selectIO + return $ sh `seq` vec `seq` + fromUnboxed sh vec + + where {-# INLINE selectIO #-} + selectIO + = do vecs <- selectChunkedP match produce len + vecs' <- mapM V.unsafeFreeze vecs + + -- TODO: avoid copy somehow. + let result = V.concat vecs' + + return (Z :. V.length result, result) diff --git a/Data/Array/Repa/Operators/Traversal.hs b/Data/Array/Repa/Operators/Traversal.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Operators/Traversal.hs @@ -0,0 +1,120 @@ +-- Generic Traversal +module Data.Array.Repa.Operators.Traversal + ( traverse, unsafeTraverse + , traverse2, unsafeTraverse2 + , traverse3, unsafeTraverse3 + , traverse4, unsafeTraverse4) +where +import Data.Array.Repa.Base +import Data.Array.Repa.Shape +import Data.Array.Repa.Repr.Delayed + + +-- | Unstructured traversal. +traverse, unsafeTraverse + :: forall r sh sh' a b + . (Shape sh, Shape sh', Repr r a) + => Array r sh a -- ^ Source array. + -> (sh -> sh') -- ^ Function to produce the extent of the result. + -> ((sh -> a) -> sh' -> b) -- ^ Function to produce elements of the result. + -- It is passed a lookup function to get elements of the source. + -> Array D sh' b + +{-# INLINE [4] traverse #-} +traverse arr transExtent newElem + = arr `deepSeqArray` + fromFunction (transExtent (extent arr)) (newElem (index arr)) + +{-# INLINE [4] unsafeTraverse #-} +unsafeTraverse arr transExtent newElem + = arr `deepSeqArray` + fromFunction (transExtent (extent arr)) (newElem (unsafeIndex arr)) + + +-- | Unstructured traversal over two arrays at once. +traverse2, unsafeTraverse2 + :: forall r1 r2 sh sh' sh'' a b c + . ( Shape sh, Shape sh', Shape sh'' + , Repr r1 a, Repr r2 b) + => Array r1 sh a -- ^ First source array. + -> Array r2 sh' b -- ^ Second source array. + -> (sh -> sh' -> sh'') -- ^ Function to produce the extent of the result. + -> ((sh -> a) -> (sh' -> b) + -> (sh'' -> c)) -- ^ Function to produce elements of the result. + -- It is passed lookup functions to get elements of the + -- source arrays. + -> Array D sh'' c + +{-# INLINE [4] traverse2 #-} +traverse2 arrA arrB transExtent newElem + = arrA `deepSeqArray` arrB `deepSeqArray` + fromFunction (transExtent (extent arrA) (extent arrB)) + (newElem (index arrA) (index arrB)) + +{-# INLINE [4] unsafeTraverse2 #-} +unsafeTraverse2 arrA arrB transExtent newElem + = arrA `deepSeqArray` arrB `deepSeqArray` + fromFunction (transExtent (extent arrA) (extent arrB)) + (newElem (unsafeIndex arrA) (unsafeIndex arrB)) + + +-- | Unstructured traversal over three arrays at once. +traverse3, unsafeTraverse3 + :: forall r1 r2 r3 + sh1 sh2 sh3 sh4 + a b c d + . ( Shape sh1, Shape sh2, Shape sh3, Shape sh4 + , Repr r1 a, Repr r2 b, Repr r3 c) + => Array r1 sh1 a + -> Array r2 sh2 b + -> Array r3 sh3 c + -> (sh1 -> sh2 -> sh3 -> sh4) + -> ( (sh1 -> a) -> (sh2 -> b) + -> (sh3 -> c) + -> sh4 -> d ) + -> Array D sh4 d + +{-# INLINE [4] traverse3 #-} +traverse3 arrA arrB arrC transExtent newElem + = arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray` + fromFunction (transExtent (extent arrA) (extent arrB) (extent arrC)) + (newElem (index arrA) (index arrB) (index arrC)) + +{-# INLINE [4] unsafeTraverse3 #-} +unsafeTraverse3 arrA arrB arrC transExtent newElem + = arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray` + fromFunction (transExtent (extent arrA) (extent arrB) (extent arrC)) + (newElem (unsafeIndex arrA) (unsafeIndex arrB) (unsafeIndex arrC)) + + +-- | Unstructured traversal over four arrays at once. +traverse4, unsafeTraverse4 + :: forall r1 r2 r3 r4 + sh1 sh2 sh3 sh4 sh5 + a b c d e + . ( Shape sh1, Shape sh2, Shape sh3, Shape sh4, Shape sh5 + , Repr r1 a, Repr r2 b, Repr r3 c, Repr r4 d) + => Array r1 sh1 a + -> Array r2 sh2 b + -> Array r3 sh3 c + -> Array r4 sh4 d + -> (sh1 -> sh2 -> sh3 -> sh4 -> sh5 ) + -> ( (sh1 -> a) -> (sh2 -> b) + -> (sh3 -> c) -> (sh4 -> d) + -> sh5 -> e ) + -> Array D sh5 e + +{-# INLINE [4] traverse4 #-} +traverse4 arrA arrB arrC arrD transExtent newElem + = arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray` arrD `deepSeqArray` + fromFunction (transExtent (extent arrA) (extent arrB) (extent arrC) (extent arrD)) + (newElem (index arrA) (index arrB) (index arrC) (index arrD)) + + +{-# INLINE [4] unsafeTraverse4 #-} +unsafeTraverse4 arrA arrB arrC arrD transExtent newElem + = arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray` arrD `deepSeqArray` + fromFunction (transExtent (extent arrA) (extent arrB) (extent arrC) (extent arrD)) + (newElem (unsafeIndex arrA) (unsafeIndex arrB) (unsafeIndex arrC) (unsafeIndex arrD)) + + diff --git a/Data/Array/Repa/Operators/Traverse.hs b/Data/Array/Repa/Operators/Traverse.hs deleted file mode 100644 --- a/Data/Array/Repa/Operators/Traverse.hs +++ /dev/null @@ -1,126 +0,0 @@ -{-# OPTIONS_HADDOCK hide #-} -{-# LANGUAGE ExplicitForAll #-} - -module Data.Array.Repa.Operators.Traverse - ( traverse, unsafeTraverse - , traverse2, unsafeTraverse2 - , traverse3, unsafeTraverse3 - , traverse4, unsafeTraverse4) -where -import Data.Array.Repa.Internals.Elt -import Data.Array.Repa.Internals.Base -import Data.Array.Repa.Shape as S - --- Generic Traversal ----------------------------------------------------------------------------- --- | Unstructured traversal. -traverse - :: forall sh sh' a b - . (Shape sh, Shape sh', Elt a) - => Array sh a -- ^ Source array. - -> (sh -> sh') -- ^ Function to produce the extent of the result. - -> ((sh -> a) -> sh' -> b) -- ^ Function to produce elements of the result. - -- It is passed a lookup function to get elements of the source. - -> Array sh' b - -{-# INLINE traverse #-} -traverse arr transExtent newElem - = arr `deepSeqArray` - fromFunction (transExtent (extent arr)) (newElem (arr !)) - - -{-# INLINE unsafeTraverse #-} -unsafeTraverse arr transExtent newElem - = arr `deepSeqArray` - fromFunction (transExtent (extent arr)) (newElem (unsafeIndex arr)) - - --- | Unstructured traversal over two arrays at once. -traverse2, unsafeTraverse2 - :: forall sh sh' sh'' a b c - . ( Shape sh, Shape sh', Shape sh'' - , Elt a, Elt b) - => Array sh a -- ^ First source array. - -> Array sh' b -- ^ Second source array. - -> (sh -> sh' -> sh'') -- ^ Function to produce the extent of the result. - -> ((sh -> a) -> (sh' -> b) - -> (sh'' -> c)) -- ^ Function to produce elements of the result. - -- It is passed lookup functions to get elements of the - -- source arrays. - -> Array sh'' c - -{-# INLINE traverse2 #-} -traverse2 arrA arrB transExtent newElem - = arrA `deepSeqArray` arrB `deepSeqArray` - fromFunction - (transExtent (extent arrA) (extent arrB)) - (newElem (arrA !) (arrB !)) - -{-# INLINE unsafeTraverse2 #-} -unsafeTraverse2 arrA arrB transExtent newElem - = arrA `deepSeqArray` arrB `deepSeqArray` - fromFunction - (transExtent (extent arrA) (extent arrB)) - (newElem (unsafeIndex arrA) (unsafeIndex arrB)) - - --- | Unstructured traversal over three arrays at once. -traverse3, unsafeTraverse3 - :: forall sh1 sh2 sh3 sh4 - a b c d - . ( Shape sh1, Shape sh2, Shape sh3, Shape sh4 - , Elt a, Elt b, Elt c) - => Array sh1 a - -> Array sh2 b - -> Array sh3 c - -> (sh1 -> sh2 -> sh3 -> sh4) - -> ( (sh1 -> a) -> (sh2 -> b) - -> (sh3 -> c) - -> sh4 -> d ) - -> Array sh4 d - -{-# INLINE traverse3 #-} -traverse3 arrA arrB arrC transExtent newElem - = arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray` - fromFunction - (transExtent (extent arrA) (extent arrB) (extent arrC)) - (newElem (arrA !) (arrB !) (arrC !)) - -{-# INLINE unsafeTraverse3 #-} -unsafeTraverse3 arrA arrB arrC transExtent newElem - = arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray` - fromFunction - (transExtent (extent arrA) (extent arrB) (extent arrC)) - (newElem (unsafeIndex arrA) (unsafeIndex arrB) (unsafeIndex arrC)) - - --- | Unstructured traversal over four arrays at once. -traverse4, unsafeTraverse4 - :: forall sh1 sh2 sh3 sh4 sh5 - a b c d e - . ( Shape sh1, Shape sh2, Shape sh3, Shape sh4, Shape sh5 - , Elt a, Elt b, Elt c, Elt d) - => Array sh1 a - -> Array sh2 b - -> Array sh3 c - -> Array sh4 d - -> (sh1 -> sh2 -> sh3 -> sh4 -> sh5 ) - -> ( (sh1 -> a) -> (sh2 -> b) - -> (sh3 -> c) -> (sh4 -> d) - -> sh5 -> e ) - -> Array sh5 e - -{-# INLINE traverse4 #-} -traverse4 arrA arrB arrC arrD transExtent newElem - = arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray` arrD `deepSeqArray` - fromFunction - (transExtent (extent arrA) (extent arrB) (extent arrC) (extent arrD)) - (newElem (arrA !) (arrB !) (arrC !) (arrD !)) - - -{-# INLINE unsafeTraverse4 #-} -unsafeTraverse4 arrA arrB arrC arrD transExtent newElem - = arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray` arrD `deepSeqArray` - fromFunction - (transExtent (extent arrA) (extent arrB) (extent arrC) (extent arrD)) - (newElem (unsafeIndex arrA) (unsafeIndex arrB) (unsafeIndex arrC) (unsafeIndex arrD)) - diff --git a/Data/Array/Repa/Properties.hs b/Data/Array/Repa/Properties.hs deleted file mode 100644 --- a/Data/Array/Repa/Properties.hs +++ /dev/null @@ -1,115 +0,0 @@ -{-# LANGUAGE ScopedTypeVariables #-} - -module Data.Array.Repa.Properties - ( props_DataArrayRepaIndex - , props_DataArrayRepa) -where -import Data.Array.Repa as R -import Data.Array.Repa.Arbitrary -import qualified Data.Array.Repa.Shape as S -import qualified Data.Vector.Unboxed as V -import Control.Monad -import Test.QuickCheck -import Prelude as P hiding (compare) - -stage :: String -stage = "Data.Array.Repa.Properties" - -compare :: (Eq a, Show a) => a -> a -> Property -compare ans ref = printTestCase message (ref == ans) - where - message = unlines ["*** Expected:", show ref - ,"*** Received:", show ans ] - - --- Data.Array.Repa.Index -------------------------------------------------------------------------- --- | QuickCheck properties for "Data.Array.Repa.Index". -props_DataArrayRepaIndex :: [(String, Property)] -props_DataArrayRepaIndex - = [(stage P.++ "." P.++ name, test) | (name, test) - <- [ ("toIndexFromIndex/DIM1", property prop_toIndexFromIndex_DIM1) - , ("toIndexFromIndex/DIM2", property prop_toIndexFromIndex_DIM2) ]] - -prop_toIndexFromIndex_DIM1 sh ix - = (sizeIsValid sh) - ==> (inShape sh ix) - ==> fromIndex sh (toIndex sh ix) `compare` (ix :: DIM1) - -prop_toIndexFromIndex_DIM2 - = forAll arbitraryShape $ \(sh :: DIM2) -> - forAll (genInShape2 sh) $ \(ix :: DIM2) -> - fromIndex sh (toIndex sh ix) `compare` ix - - - - --- Data.Array.Repa -------------------------------------------------------------------------------- --- | QuickCheck properties for "Data.Array.Repa" and its children. -props_DataArrayRepa :: [(String, Property)] -props_DataArrayRepa - = props_DataArrayRepaIndex - P.++ [(stage P.++ "." P.++ name, test) | (name, test) - <- [ ("id_force/DIM5", property prop_id_force_DIM5) - , ("id_toScalarUnit", property prop_id_toScalarUnit) - , ("id_toListFromList/DIM3", property prop_id_toListFromList_DIM3) - , ("id_transpose/DIM4", property prop_id_transpose_DIM4) - , ("reshapeTransposeSize/DIM3", property prop_reshapeTranspose_DIM3) - , ("appendIsAppend/DIM3", property prop_appendIsAppend_DIM3) - , ("sumIsSum/DIM3", property prop_sumIsSum_DIM3) - , ("sumAllIsSum/DIM3", property prop_sumAllIsSum_DIM3) ]] - - --- The Eq instance uses fold and zipWith. -prop_id_force_DIM5 - = forAll (arbitrarySmallArray 10) $ \(arr :: Array DIM5 Int) -> - force arr `compare` arr - -prop_id_toScalarUnit (x :: Int) - = toScalar (singleton x) `compare` x - --- Conversions ------------------------ -prop_id_toListFromList_DIM3 - = forAll (arbitrarySmallShape 10) $ \(sh :: DIM3) -> - forAll (arbitraryListOfLength (S.size sh)) $ \(xx :: [Int]) -> - toList (fromList sh xx) `compare` xx - --- Index Space Transforms ------------- -prop_id_transpose_DIM4 - = forAll (arbitrarySmallArray 20) $ \(arr :: Array DIM3 Int) -> - transpose (transpose arr) `compare` arr - --- A reshaped array has the same size and sum as the original -prop_reshapeTranspose_DIM3 - = forAll (arbitrarySmallArray 20) $ \(arr :: Array DIM3 Int) -> - let arr' = transpose arr - sh' = extent arr' - in (S.size (extent (reshape sh' arr)) `compare` S.size (extent arr)) - .&&. (sumAll arr' `compare` sumAll arr) - -prop_appendIsAppend_DIM3 - = forAll (arbitrarySmallArray 20) $ \(arr1 :: Array DIM3 Int) -> - sumAll (append arr1 arr1) `compare` (2 * sumAll arr1) - --- Reductions -------------------------- -prop_sumIsSum_DIM3 - = forAll (arbitrarySmallArray 20) $ \(arr :: Array DIM3 Int) -> - let sh :. sz = extent arr - elemFn ix = V.foldl' (+) 0 - $ V.map (\i -> arr ! (ix :. i)) - (V.enumFromTo 0 (sz-1)) - in - R.fold (+) 0 arr `compare` fromFunction sh elemFn - -prop_sumAllIsSum_DIM3 - = forAll (arbitrarySmallShape 20) $ \(sh :: DIM3) -> - forAll (arbitraryListOfLength (S.size sh)) $ \(xx :: [Int]) -> - sumAll (fromList sh xx) `compare` P.sum xx - - --- Utils ------------------------------------------------------------------------------------------ -genInShape2 :: DIM2 -> Gen DIM2 -genInShape2 (Z :. yMax :. xMax) - = do y <- liftM (`mod` yMax) $ arbitrary - x <- liftM (`mod` xMax) $ arbitrary - return $ Z :. y :. x - diff --git a/Data/Array/Repa/Repr/ByteString.hs b/Data/Array/Repa/Repr/ByteString.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Repr/ByteString.hs @@ -0,0 +1,57 @@ + +module Data.Array.Repa.Repr.ByteString + ( B, Array (..) + , fromByteString, toByteString) +where +import Data.Array.Repa.Shape +import Data.Array.Repa.Base +import Data.Array.Repa.Repr.Delayed +import Data.Word +import qualified Data.ByteString as B +import qualified Data.ByteString.Unsafe as BU +import Data.ByteString (ByteString) + + +-- | Strict ByteStrings arrays are represented as ForeignPtr buffers of Word8 +data B +data instance Array B sh Word8 + = AByteString sh !ByteString + +deriving instance Show sh + => Show (Array B sh Word8) + + +-- Repr ----------------------------------------------------------------------- +-- | Read elements from a `ByteString`. +instance Repr B Word8 where + {-# INLINE linearIndex #-} + linearIndex (AByteString _ bs) ix + = bs `B.index` ix + + {-# INLINE unsafeLinearIndex #-} + unsafeLinearIndex (AByteString _ bs) ix + = bs `BU.unsafeIndex` ix + + {-# INLINE extent #-} + extent (AByteString sh _) + = sh + + {-# INLINE deepSeqArray #-} + deepSeqArray (AByteString sh bs) x + = sh `deepSeq` bs `seq` x + + +-- Conversions ---------------------------------------------------------------- +-- | O(1). Wrap a `ByteString` as an array. +fromByteString + :: Shape sh + => sh -> ByteString -> Array B sh Word8 +{-# INLINE fromByteString #-} +fromByteString sh bs + = AByteString sh bs + + +-- | O(1). Unpack a `ByteString` from an array. +toByteString :: Array B sh Word8 -> ByteString +{-# INLINE toByteString #-} +toByteString (AByteString _ bs) = bs diff --git a/Data/Array/Repa/Repr/Cursored.hs b/Data/Array/Repa/Repr/Cursored.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Repr/Cursored.hs @@ -0,0 +1,109 @@ +{-# LANGUAGE MagicHash #-} +module Data.Array.Repa.Repr.Cursored + ( C, Array (..) + , makeCursored) +where +import Data.Array.Repa.Base +import Data.Array.Repa.Shape +import Data.Array.Repa.Index +import Data.Array.Repa.Repr.Delayed +import Data.Array.Repa.Repr.Undefined +import Data.Array.Repa.Eval.Fill +import Data.Array.Repa.Eval.Elt +import Data.Array.Repa.Eval.Cursored +import GHC.Exts + +-- | Cursored Arrays. +-- These are produced by Repa's stencil functions, and help the fusion +-- framework to share index compuations between array elements. +-- +-- The basic idea is described in ``Efficient Parallel Stencil Convolution'', +-- Ben Lippmeier and Gabriele Keller, Haskell 2011 -- though the underlying +-- array representation has changed since this paper was published. +data C + +data instance Array C sh e + = forall cursor. ACursored + { cursoredExtent :: sh + + -- | Make a cursor to a particular element. + , makeCursor :: sh -> cursor + + -- | Shift the cursor by an offset, to get to another element. + , shiftCursor :: sh -> cursor -> cursor + + -- | Load\/compute the element at the given cursor. + , loadCursor :: cursor -> e } + + +-- Repr ----------------------------------------------------------------------- +-- | Compute elements of a cursored array. +instance Repr C a where + {-# INLINE index #-} + index (ACursored _ makec _ loadc) + = loadc . makec + + {-# INLINE unsafeIndex #-} + unsafeIndex = index + + {-# INLINE linearIndex #-} + linearIndex (ACursored sh makec _ loadc) + = loadc . makec . fromIndex sh + + {-# INLINE extent #-} + extent (ACursored sh _ _ _) + = sh + + {-# INLINE deepSeqArray #-} + deepSeqArray (ACursored sh makec shiftc loadc) y + = sh `deepSeq` makec `seq` shiftc `seq` loadc `seq` y + + +-- Fill ----------------------------------------------------------------------- +-- | Compute all elements in an rank-2 array. +instance (Fillable r2 e, Elt e) => Fill C r2 DIM2 e where + {-# INLINE fillP #-} + fillP (ACursored (Z :. h :. w) makec shiftc loadc) marr + = fillCursoredBlock2P + (unsafeWriteMArr marr) + makec shiftc loadc + w 0 0 (w - 1) (h - 1) + + {-# INLINE fillS #-} + fillS (ACursored (Z :. (I# h) :. (I# w)) makec shiftc loadc) marr + = fillCursoredBlock2S + (unsafeWriteMArr marr) + makec shiftc loadc + w 0# 0# (w -# 1#) (h -# 1#) + + +-- | Compute a range of elements in a rank-2 array. +instance (Fillable r2 e, Elt e) => FillRange C r2 DIM2 e where + {-# INLINE fillRangeP #-} + fillRangeP (ACursored (Z :. _h :. w) makec shiftc loadc) marr + (Z :. y0 :. x0) (Z :. y1 :. x1) + = fillCursoredBlock2P + (unsafeWriteMArr marr) + makec shiftc loadc + w x0 y0 x1 y1 + + {-# INLINE fillRangeS #-} + fillRangeS (ACursored (Z :. _h :. (I# w)) makec shiftc loadc) marr + (Z :. (I# y0) :. (I# x0)) + (Z :. (I# y1) :. (I# x1)) + = fillCursoredBlock2S + (unsafeWriteMArr marr) + makec shiftc loadc + w x0 y0 x1 y1 + +-- Conversions ---------------------------------------------------------------- +-- | Define a new cursored array. +makeCursored + :: sh + -> (sh -> cursor) -- ^ Create a cursor for an index. + -> (sh -> cursor -> cursor) -- ^ Shift a cursor by an offset. + -> (cursor -> e) -- ^ Compute the element at the cursor. + -> Array C sh e + +{-# INLINE makeCursored #-} +makeCursored = ACursored diff --git a/Data/Array/Repa/Repr/Delayed.hs b/Data/Array/Repa/Repr/Delayed.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Repr/Delayed.hs @@ -0,0 +1,99 @@ +{-# LANGUAGE MagicHash #-} +module Data.Array.Repa.Repr.Delayed + ( D, Array(..) + , fromFunction, toFunction + , delay) +where +import Data.Array.Repa.Eval.Elt +import Data.Array.Repa.Eval.Cursored +import Data.Array.Repa.Eval.Chunked +import Data.Array.Repa.Eval.Fill +import Data.Array.Repa.Index +import Data.Array.Repa.Shape +import Data.Array.Repa.Base +import GHC.Exts + +-- | Delayed arrays are represented as functions from the index to element value. +data D +data instance Array D sh e + = ADelayed + sh + (sh -> e) + + +-- Repr ----------------------------------------------------------------------- +-- | Compute elements of a delayed array. +instance Repr D a where + {-# INLINE index #-} + index (ADelayed _ f) ix = f ix + + {-# INLINE linearIndex #-} + linearIndex (ADelayed sh f) ix = f (fromIndex sh ix) + + {-# INLINE extent #-} + extent (ADelayed sh _) + = sh + + {-# INLINE deepSeqArray #-} + deepSeqArray (ADelayed sh f) y + = sh `deepSeq` f `seq` y + + +-- Fill ----------------------------------------------------------------------- +-- | Compute all elements in an array. +instance (Fillable r2 e, Shape sh) => Fill D r2 sh e where + {-# INLINE [4] fillP #-} + fillP (ADelayed sh getElem) marr + = fillChunkedP (size sh) (unsafeWriteMArr marr) (getElem . fromIndex sh) + + {-# INLINE [4] fillS #-} + fillS (ADelayed sh getElem) marr + = fillChunkedS (size sh) (unsafeWriteMArr marr) (getElem . fromIndex sh) + + +-- | Compute a range of elements in a rank-2 array. +instance (Fillable r2 e, Elt e) => FillRange D r2 DIM2 e where + {-# INLINE [1] fillRangeP #-} + fillRangeP (ADelayed (Z :. _h :. w) getElem) marr + (Z :. y0 :. x0) (Z :. y1 :. x1) + = fillBlock2P (unsafeWriteMArr marr) + getElem + w x0 y0 x1 y1 + + {-# INLINE [1] fillRangeS #-} + fillRangeS (ADelayed (Z :. _h :. (I# w)) getElem) marr + (Z :. (I# y0) :. (I# x0)) (Z :. (I# y1) :. (I# x1)) + = fillBlock2S (unsafeWriteMArr marr) + getElem + w x0 y0 x1 y1 + + +-- Conversions ---------------------------------------------------------------- +-- | O(1). Wrap a function as a delayed array. +fromFunction :: sh -> (sh -> a) -> Array D sh a +{-# INLINE fromFunction #-} +fromFunction sh f + = ADelayed sh f + + +-- | O(1). Produce the extent of an array and a function to retrieve an arbitrary element. +toFunction + :: (Shape sh, Repr r1 a) + => Array r1 sh a -> (sh, sh -> a) +{-# INLINE toFunction #-} +toFunction arr + = case delay arr of + ADelayed sh f -> (sh, f) + + +-- | O(1). Delay an array. +-- This wraps the internal representation to be a function from +-- indices to elements, so consumers don't need to worry about +-- what the previous representation was. +-- +delay :: (Shape sh, Repr r e) + => Array r sh e -> Array D sh e +{-# INLINE delay #-} +delay arr = ADelayed (extent arr) (index arr) + + diff --git a/Data/Array/Repa/Repr/ForeignPtr.hs b/Data/Array/Repa/Repr/ForeignPtr.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Repr/ForeignPtr.hs @@ -0,0 +1,111 @@ + +module Data.Array.Repa.Repr.ForeignPtr + ( F, Array (..) + , fromForeignPtr, toForeignPtr + , computeIntoS, computeIntoP) +where +import Data.Array.Repa.Shape +import Data.Array.Repa.Base +import Data.Array.Repa.Eval.Fill +import Data.Array.Repa.Repr.Delayed +import Foreign.Storable +import Foreign.ForeignPtr +import Foreign.Marshal.Alloc +import System.IO.Unsafe + +-- | Arrays represented as foreign buffers in the C heap. +data F +data instance Array F sh e + = AForeignPtr !sh !Int !(ForeignPtr e) + +-- Repr ----------------------------------------------------------------------- +-- | Read elements from a foreign buffer. +instance Storable a => Repr F a where + {-# INLINE linearIndex #-} + linearIndex (AForeignPtr _ len fptr) ix + | ix < len + = unsafePerformIO + $ withForeignPtr fptr + $ \ptr -> peekElemOff ptr ix + + | otherwise + = error "Repa: foreign array index out of bounds" + + {-# INLINE unsafeLinearIndex #-} + unsafeLinearIndex (AForeignPtr _ _ fptr) ix + = unsafePerformIO + $ withForeignPtr fptr + $ \ptr -> peekElemOff ptr ix + + {-# INLINE extent #-} + extent (AForeignPtr sh _ _) + = sh + + {-# INLINE deepSeqArray #-} + deepSeqArray (AForeignPtr sh len fptr) x + = sh `deepSeq` len `seq` fptr `seq` x + + +-- Fill ----------------------------------------------------------------------- +-- | Filling of foreign buffers. +instance Storable e => Fillable F e where + data MArr F e + = FPArr !Int !(ForeignPtr e) + + {-# INLINE newMArr #-} + newMArr n + = do let (proxy :: e) = undefined + ptr <- mallocBytes (sizeOf proxy * n) + _ <- peek ptr `asTypeOf` return proxy + + fptr <- newForeignPtr finalizerFree ptr + return $ FPArr n fptr + + {-# INLINE unsafeWriteMArr #-} + unsafeWriteMArr (FPArr _ fptr) !ix !x + = withForeignPtr fptr + $ \ptr -> pokeElemOff ptr ix x + + {-# INLINE unsafeFreezeMArr #-} + unsafeFreezeMArr !sh (FPArr len fptr) + = return $ AForeignPtr sh len fptr + + +-- Conversions ---------------------------------------------------------------- +-- | O(1). Wrap a `ForeignPtr` as an array. +fromForeignPtr + :: Shape sh + => sh -> ForeignPtr e -> Array F sh e +{-# INLINE fromForeignPtr #-} +fromForeignPtr !sh !fptr + = AForeignPtr sh (size sh) fptr + + +-- | O(1). Unpack a `ForeignPtr` from an array. +toForeignPtr :: Array F sh e -> ForeignPtr e +{-# INLINE toForeignPtr #-} +toForeignPtr (AForeignPtr _ _ fptr) + = fptr + + +-- | Compute an array sequentially and write the elements into a foreign +-- buffer without intermediate copying. If you want to copy a +-- pre-existing manifest array to a foreign buffer then `delay` it first. +computeIntoS + :: Fill r1 F sh e + => ForeignPtr e -> Array r1 sh e -> IO () +{-# INLINE computeIntoS #-} +computeIntoS !fptr !arr + = fillS arr (FPArr 0 fptr) + + +-- | Compute an array in parallel and write the elements into a foreign +-- buffer without intermediate copying. If you want to copy a +-- pre-existing manifest array to a foreign buffer then `delay` it first. +computeIntoP + :: Fill r1 F sh e + => ForeignPtr e -> Array r1 sh e -> IO () +{-# INLINE computeIntoP #-} +computeIntoP !fptr !arr + = fillP arr (FPArr 0 fptr) + diff --git a/Data/Array/Repa/Repr/Partitioned.hs b/Data/Array/Repa/Repr/Partitioned.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Repr/Partitioned.hs @@ -0,0 +1,82 @@ + + +module Data.Array.Repa.Repr.Partitioned + ( P, Array (..) + , Range(..) + , inRange) +where +import Data.Array.Repa.Base +import Data.Array.Repa.Shape +import Data.Array.Repa.Eval +import Data.Array.Repa.Repr.Delayed +import Data.Array.Repa.Repr.Undefined + + +-- | Partitioned arrays. +-- The last partition takes priority +-- +-- These are produced by Repa's support functions and allow arrays to be defined +-- using a different element function for each partition. +-- +-- The basic idea is described in ``Efficient Parallel Stencil Convolution'', +-- Ben Lippmeier and Gabriele Keller, Haskell 2011 -- though the underlying +-- array representation has changed since this paper was published. +-- +data P r1 r2 + +data instance Array (P r1 r2) sh e + = APart sh -- size of the whole array + (Range sh) (Array r1 sh e) -- if in range use this array + (Array r2 sh e) -- otherwise use this array + +data Range sh + = Range sh sh -- indices defining the range + (sh -> Bool) -- predicate to check whether were in range + +-- | Check whether an index is within the given range. +{-# INLINE inRange #-} +inRange :: Range sh -> sh -> Bool +inRange (Range _ _ p) ix + = p ix + + +-- Repr ----------------------------------------------------------------------- +-- | Read elements from a partitioned array. +instance (Repr r1 e, Repr r2 e) => Repr (P r1 r2) e where + {-# INLINE index #-} + index (APart _ range arr1 arr2) ix + | inRange range ix = index arr1 ix + | otherwise = index arr2 ix + + {-# INLINE linearIndex #-} + linearIndex arr@(APart sh _ _ _) ix + = index arr $ fromIndex sh ix + + {-# INLINE extent #-} + extent (APart sh _ _ _) + = sh + + {-# INLINE deepSeqArray #-} + deepSeqArray (APart sh range arr1 arr2) y + = sh `deepSeq` range `deepSeqRange` arr1 `deepSeqArray` arr2 `deepSeqArray` y + + +{-# INLINE deepSeqRange #-} +deepSeqRange :: Shape sh => Range sh -> b -> b +deepSeqRange (Range low high f) y + = low `deepSeq` high `deepSeq` f `seq` y + + +-- Fill ----------------------------------------------------------------------- +instance ( FillRange r1 r3 sh e, Fill r2 r3 sh e + , Fillable r3 e) + => Fill (P r1 r2) r3 sh e where + {-# INLINE fillP #-} + fillP (APart _ (Range ix10 ix11 _) arr1 arr2) marr + = do fillRangeP arr1 marr ix10 ix11 + fillP arr2 marr + + {-# INLINE fillS #-} + fillS (APart _ (Range ix10 ix11 _) arr1 arr2) marr + = do fillRangeS arr1 marr ix10 ix11 + fillS arr2 marr diff --git a/Data/Array/Repa/Repr/Unboxed.hs b/Data/Array/Repa/Repr/Unboxed.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Repr/Unboxed.hs @@ -0,0 +1,233 @@ + +module Data.Array.Repa.Repr.Unboxed + ( U, U.Unbox, Array (..) + , computeUnboxedS, computeUnboxedP + , fromListUnboxed + , fromUnboxed, toUnboxed + + , zip, zip3, zip4, zip5, zip6 + , unzip, unzip3, unzip4, unzip5, unzip6) +where +import Data.Array.Repa.Shape as R +import Data.Array.Repa.Base as R +import Data.Array.Repa.Eval as R +import Data.Array.Repa.Repr.Delayed as R +import qualified Data.Vector.Unboxed as U +import qualified Data.Vector.Unboxed.Mutable as UM +import Control.Monad +import Prelude hiding (zip, zip3, unzip, unzip3) + +-- | Unboxed arrays are represented as unboxed vectors. +-- +-- The implementation of `Data.Vector.Unboxed` is based on type families and +-- picks an efficient, specialised representation for every element type. In +-- particular, unboxed vectors of pairs are represented as pairs of unboxed +-- vectors. This is the most efficient representation for numerical data. +-- +data U +data instance U.Unbox e => Array U sh e + = AUnboxed sh !(U.Vector e) + +deriving instance (Show sh, Show e, U.Unbox e) + => Show (Array U sh e) + +-- Repr ----------------------------------------------------------------------- +-- | Read elements from an unboxed vector array. +instance U.Unbox a => Repr U a where + {-# INLINE linearIndex #-} + linearIndex (AUnboxed _ vec) ix + = vec U.! ix + + {-# INLINE unsafeLinearIndex #-} + unsafeLinearIndex (AUnboxed _ vec) ix + = vec `U.unsafeIndex` ix + + {-# INLINE extent #-} + extent (AUnboxed sh _) + = sh + + {-# INLINE deepSeqArray #-} + deepSeqArray (AUnboxed sh vec) x + = sh `deepSeq` vec `seq` x + + +-- Fill ----------------------------------------------------------------------- +-- | Filling of unboxed vector arrays. +instance U.Unbox e => Fillable U e where + data MArr U e + = UMArr (UM.IOVector e) + + {-# INLINE newMArr #-} + newMArr n + = liftM UMArr (UM.new n) + + {-# INLINE unsafeWriteMArr #-} + unsafeWriteMArr (UMArr v) ix + = UM.unsafeWrite v ix + + {-# INLINE unsafeFreezeMArr #-} + unsafeFreezeMArr sh (UMArr mvec) + = do vec <- U.unsafeFreeze mvec + return $ AUnboxed sh vec + + +-- Conversions ---------------------------------------------------------------- +-- | Sequential computation of array elements.. +-- +-- * This is an alias for `computeS` with a more specific type. +-- +computeUnboxedS + :: Fill r1 U sh e + => Array r1 sh e -> Array U sh e +{-# INLINE computeUnboxedS #-} +computeUnboxedS = computeS + + +-- | Parallel computation of array elements. +-- +-- * This is an alias for `computeP` with a more specific type. +-- +computeUnboxedP + :: Fill r1 U sh e + => Array r1 sh e -> Array U sh e +{-# INLINE computeUnboxedP #-} +computeUnboxedP = computeP + + +-- | O(n). Convert a list to an unboxed vector array. +-- +-- * This is an alias for `fromList` with a more specific type. +-- +fromListUnboxed + :: (Shape sh, U.Unbox a) + => sh -> [a] -> Array U sh a +{-# INLINE fromListUnboxed #-} +fromListUnboxed = R.fromList + + +-- | O(1). Wrap an unboxed vector as an array. +fromUnboxed + :: (Shape sh, U.Unbox e) + => sh -> U.Vector e -> Array U sh e +{-# INLINE fromUnboxed #-} +fromUnboxed sh vec + = AUnboxed sh vec + + +-- | O(1). Unpack an unboxed vector from an array. +toUnboxed + :: U.Unbox e + => Array U sh e -> U.Vector e +{-# INLINE toUnboxed #-} +toUnboxed (AUnboxed _ vec) + = vec + +-- Zip ------------------------------------------------------------------------ +-- | O(1). Zip some unboxed arrays. +-- The shapes must be identical else `error`. +zip :: (Shape sh, U.Unbox a, U.Unbox b) + => Array U sh a -> Array U sh b + -> Array U sh (a, b) +{-# INLINE zip #-} +zip (AUnboxed sh1 vec1) (AUnboxed sh2 vec2) + | sh1 /= sh2 = error "Repa: zip array shapes not identical" + | otherwise = AUnboxed sh1 (U.zip vec1 vec2) + + +-- | O(1). Zip some unboxed arrays. +-- The shapes must be identical else `error`. +zip3 :: (Shape sh, U.Unbox a, U.Unbox b, U.Unbox c) + => Array U sh a -> Array U sh b -> Array U sh c + -> Array U sh (a, b, c) +{-# INLINE zip3 #-} +zip3 (AUnboxed sh1 vec1) (AUnboxed sh2 vec2) (AUnboxed sh3 vec3) + | sh1 /= sh2 || sh1 /= sh3 + = error "Repa: zip array shapes not identical" + | otherwise = AUnboxed sh1 (U.zip3 vec1 vec2 vec3) + + +-- | O(1). Zip some unboxed arrays. +-- The shapes must be identical else `error`. +zip4 :: (Shape sh, U.Unbox a, U.Unbox b, U.Unbox c, U.Unbox d) + => Array U sh a -> Array U sh b -> Array U sh c -> Array U sh d + -> Array U sh (a, b, c, d) +{-# INLINE zip4 #-} +zip4 (AUnboxed sh1 vec1) (AUnboxed sh2 vec2) (AUnboxed sh3 vec3) (AUnboxed sh4 vec4) + | sh1 /= sh2 || sh1 /= sh3 || sh1 /= sh4 + = error "Repa: zip array shapes not identical" + | otherwise = AUnboxed sh1 (U.zip4 vec1 vec2 vec3 vec4) + + +-- | O(1). Zip some unboxed arrays. +-- The shapes must be identical else `error`. +zip5 :: (Shape sh, U.Unbox a, U.Unbox b, U.Unbox c, U.Unbox d, U.Unbox e) + => Array U sh a -> Array U sh b -> Array U sh c -> Array U sh d -> Array U sh e + -> Array U sh (a, b, c, d, e) +{-# INLINE zip5 #-} +zip5 (AUnboxed sh1 vec1) (AUnboxed sh2 vec2) (AUnboxed sh3 vec3) (AUnboxed sh4 vec4) (AUnboxed sh5 vec5) + | sh1 /= sh2 || sh1 /= sh3 || sh1 /= sh4 || sh1 /= sh5 + = error "Repa: zip array shapes not identical" + | otherwise = AUnboxed sh1 (U.zip5 vec1 vec2 vec3 vec4 vec5) + + +-- | O(1). Zip some unboxed arrays. +-- The shapes must be identical else `error`. +zip6 :: (Shape sh, U.Unbox a, U.Unbox b, U.Unbox c, U.Unbox d, U.Unbox e, U.Unbox f) + => Array U sh a -> Array U sh b -> Array U sh c -> Array U sh d -> Array U sh e -> Array U sh f + -> Array U sh (a, b, c, d, e, f) +{-# INLINE zip6 #-} +zip6 (AUnboxed sh1 vec1) (AUnboxed sh2 vec2) (AUnboxed sh3 vec3) (AUnboxed sh4 vec4) (AUnboxed sh5 vec5) (AUnboxed sh6 vec6) + | sh1 /= sh2 || sh1 /= sh3 || sh1 /= sh4 || sh1 /= sh5 || sh1 /= sh6 + = error "Repa: zip array shapes not identical" + | otherwise = AUnboxed sh1 (U.zip6 vec1 vec2 vec3 vec4 vec5 vec6) + + +-- Unzip ---------------------------------------------------------------------- +-- | O(1). Unzip an unboxed array. +unzip :: (U.Unbox a, U.Unbox b) + => Array U sh (a, b) + -> (Array U sh a, Array U sh b) +{-# INLINE unzip #-} +unzip (AUnboxed sh vec) + = let (as, bs) = U.unzip vec + in (AUnboxed sh as, AUnboxed sh bs) + + +-- | O(1). Unzip an unboxed array. +unzip3 :: (U.Unbox a, U.Unbox b, U.Unbox c) + => Array U sh (a, b, c) + -> (Array U sh a, Array U sh b, Array U sh c) +{-# INLINE unzip3 #-} +unzip3 (AUnboxed sh vec) + = let (as, bs, cs) = U.unzip3 vec + in (AUnboxed sh as, AUnboxed sh bs, AUnboxed sh cs) + + +-- | O(1). Unzip an unboxed array. +unzip4 :: (U.Unbox a, U.Unbox b, U.Unbox c, U.Unbox d) + => Array U sh (a, b, c, d) + -> (Array U sh a, Array U sh b, Array U sh c, Array U sh d) +{-# INLINE unzip4 #-} +unzip4 (AUnboxed sh vec) + = let (as, bs, cs, ds) = U.unzip4 vec + in (AUnboxed sh as, AUnboxed sh bs, AUnboxed sh cs, AUnboxed sh ds) + + +-- | O(1). Unzip an unboxed array. +unzip5 :: (U.Unbox a, U.Unbox b, U.Unbox c, U.Unbox d, U.Unbox e) + => Array U sh (a, b, c, d, e) + -> (Array U sh a, Array U sh b, Array U sh c, Array U sh d, Array U sh e) +{-# INLINE unzip5 #-} +unzip5 (AUnboxed sh vec) + = let (as, bs, cs, ds, es) = U.unzip5 vec + in (AUnboxed sh as, AUnboxed sh bs, AUnboxed sh cs, AUnboxed sh ds, AUnboxed sh es) + + +-- | O(1). Unzip an unboxed array. +unzip6 :: (U.Unbox a, U.Unbox b, U.Unbox c, U.Unbox d, U.Unbox e, U.Unbox f) + => Array U sh (a, b, c, d, e, f) + -> (Array U sh a, Array U sh b, Array U sh c, Array U sh d, Array U sh e, Array U sh f) +{-# INLINE unzip6 #-} +unzip6 (AUnboxed sh vec) + = let (as, bs, cs, ds, es, fs) = U.unzip6 vec + in (AUnboxed sh as, AUnboxed sh bs, AUnboxed sh cs, AUnboxed sh ds, AUnboxed sh es, AUnboxed sh fs) diff --git a/Data/Array/Repa/Repr/Undefined.hs b/Data/Array/Repa/Repr/Undefined.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Repr/Undefined.hs @@ -0,0 +1,41 @@ + +module Data.Array.Repa.Repr.Undefined + ( X, Array (..)) +where +import Data.Array.Repa.Base +import Data.Array.Repa.Shape +import Data.Array.Repa.Eval.Fill + + +-- | An array with undefined elements. +-- +-- * This is normally used as the last representation in a partitioned array, +-- as the previous partitions are expected to provide full coverage. +data X +data instance Array X sh e + = AUndefined sh + + +-- | Undefined array elements. Inspecting them yields `error`. +-- +instance Repr X e where + {-# INLINE deepSeqArray #-} + deepSeqArray _ x + = x + + {-# INLINE extent #-} + extent (AUndefined sh) + = sh + + {-# INLINE index #-} + index (AUndefined _) _ = error "Repa: array element is undefined." + + {-# INLINE linearIndex #-} + linearIndex (AUndefined _) _ = error "Repa: array element is undefined." + + +instance (Shape sh, Fillable r2 e, Num e) => Fill X r2 sh e where + fillS _ _ = return () + fillP _ _ = return () + + diff --git a/Data/Array/Repa/Repr/Vector.hs b/Data/Array/Repa/Repr/Vector.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Repr/Vector.hs @@ -0,0 +1,112 @@ + +module Data.Array.Repa.Repr.Vector + ( V, Array (..) + , computeVectorS, computeVectorP + , fromListVector + , fromVector + , toVector) +where +import Data.Array.Repa.Shape +import Data.Array.Repa.Base +import Data.Array.Repa.Eval +import qualified Data.Vector as V +import qualified Data.Vector.Mutable as VM +import Control.Monad + +-- | Arrays represented as boxed vectors. +-- +-- This representation should only be used when your element type doesn't +-- have an `Unbox` instsance. If it does, then use the Unboxed `U` +-- representation will be faster. +data V +data instance Array V sh e + = AVector sh !(V.Vector e) + +deriving instance (Show sh, Show e) + => Show (Array V sh e) + +-- Repr ----------------------------------------------------------------------- +-- | Read elements from a boxed vector array. +instance Repr V a where + {-# INLINE linearIndex #-} + linearIndex (AVector _ vec) ix + = vec V.! ix + + {-# INLINE unsafeLinearIndex #-} + unsafeLinearIndex (AVector _ vec) ix + = vec `V.unsafeIndex` ix + + {-# INLINE extent #-} + extent (AVector sh _) + = sh + + {-# INLINE deepSeqArray #-} + deepSeqArray (AVector sh vec) x + = sh `deepSeq` vec `seq` x + + +-- Fill ----------------------------------------------------------------------- +-- | Filling of boxed vector arrays. +instance Fillable V e where + data MArr V e + = MVec (VM.IOVector e) + + {-# INLINE newMArr #-} + newMArr n + = liftM MVec (VM.new n) + + {-# INLINE unsafeWriteMArr #-} + unsafeWriteMArr (MVec v) ix + = VM.unsafeWrite v ix + + {-# INLINE unsafeFreezeMArr #-} + unsafeFreezeMArr sh (MVec mvec) + = do vec <- V.unsafeFreeze mvec + return $ AVector sh vec + + +-- Conversions ---------------------------------------------------------------- +-- | Sequential computation of array elements. +-- +-- * This is an alias for `compute` with a more specific type. +-- +computeVectorS + :: Fill r1 V sh e + => Array r1 sh e -> Array V sh e +{-# INLINE computeVectorS #-} +computeVectorS = computeS + + +-- | Parallel computation of array elements. +computeVectorP + :: Fill r1 V sh e + => Array r1 sh e -> Array V sh e +{-# INLINE computeVectorP #-} +computeVectorP = computeP + + +-- | O(n). Convert a list to a boxed vector array. +-- +-- * This is an alias for `fromList` with a more specific type. +-- +fromListVector :: Shape sh => sh -> [a] -> Array V sh a +{-# INLINE fromListVector #-} +fromListVector = fromList + + +-- | O(1). Wrap a boxed vector as an array. +fromVector + :: Shape sh + => sh -> V.Vector e -> Array V sh e +{-# INLINE fromVector #-} +fromVector sh vec + = AVector sh vec + + +-- | O(1). Unpack a boxed vector from an array. +toVector :: Array V sh e -> V.Vector e +{-# INLINE toVector #-} +toVector (AVector _ vec) + = vec + + diff --git a/Data/Array/Repa/Shape.hs b/Data/Array/Repa/Shape.hs --- a/Data/Array/Repa/Shape.hs +++ b/Data/Array/Repa/Shape.hs @@ -7,7 +7,7 @@ , showShape ) where --- Shape ------------------------------------------------------------------------------------------ +-- Shape ---------------------------------------------------------------------- -- | Class of types that can be used as array shapes and indices. class Eq sh => Shape sh where diff --git a/Data/Array/Repa/Slice.hs b/Data/Array/Repa/Slice.hs --- a/Data/Array/Repa/Slice.hs +++ b/Data/Array/Repa/Slice.hs @@ -10,7 +10,7 @@ , Slice (..)) where import Data.Array.Repa.Index -import Prelude hiding (replicate, drop) +import Prelude hiding (replicate, drop) -- | Select all indices at a certain position. @@ -47,37 +47,37 @@ instance Slice Z where - {-# INLINE sliceOfFull #-} + {-# INLINE [1] sliceOfFull #-} sliceOfFull _ _ = Z - {-# INLINE fullOfSlice #-} + {-# INLINE [1] fullOfSlice #-} fullOfSlice _ _ = Z instance Slice (Any sh) where - {-# INLINE sliceOfFull #-} + {-# INLINE [1] sliceOfFull #-} sliceOfFull _ sh = sh - {-# INLINE fullOfSlice #-} + {-# INLINE [1] fullOfSlice #-} fullOfSlice _ sh = sh instance Slice sl => Slice (sl :. Int) where - {-# INLINE sliceOfFull #-} + {-# INLINE [1] sliceOfFull #-} sliceOfFull (fsl :. _) (ssl :. _) = sliceOfFull fsl ssl - {-# INLINE fullOfSlice #-} + {-# INLINE [1] fullOfSlice #-} fullOfSlice (fsl :. n) ssl = fullOfSlice fsl ssl :. n instance Slice sl => Slice (sl :. All) where - {-# INLINE sliceOfFull #-} + {-# INLINE [1] sliceOfFull #-} sliceOfFull (fsl :. All) (ssl :. s) = sliceOfFull fsl ssl :. s - {-# INLINE fullOfSlice #-} + {-# INLINE [1] fullOfSlice #-} fullOfSlice (fsl :. All) (ssl :. s) = fullOfSlice fsl ssl :. s diff --git a/Data/Array/Repa/Specialised/Dim2.hs b/Data/Array/Repa/Specialised/Dim2.hs --- a/Data/Array/Repa/Specialised/Dim2.hs +++ b/Data/Array/Repa/Specialised/Dim2.hs @@ -7,7 +7,9 @@ , clampToBorder2 , makeBordered2) where -import Data.Array.Repa +import Data.Array.Repa.Index +import Data.Array.Repa.Repr.Partitioned +import Data.Array.Repa.Repr.Undefined -- | Check if an index lies inside the given extent. @@ -63,19 +65,22 @@ | otherwise = sh :. y :. x --- | Make a 2D partitioned array given two generators, one to produce elements in the --- border region, and one to produce values in the internal region. +-- | Make a 2D partitioned array from two others, one to produce the elements +-- in the internal region, and one to produce elements in the border region. +-- The two arrays must have the same extent. -- The border must be the same width on all sides. +-- +-- TODO: Check arrays have same extent. +-- makeBordered2 - :: Elt a - => DIM2 -- ^ Extent of array. + :: DIM2 -- ^ Extent of array. -> Int -- ^ Width of border. - -> Generator DIM2 a -- ^ Generator for border elements. - -> Generator DIM2 a -- ^ Generator for internal elements. - -> Array DIM2 a + -> Array r1 DIM2 a -- ^ Array for internal elements. + -> Array r2 DIM2 a -- ^ Array for border elements. + -> Array (P r1 (P r2 (P r2 (P r2 (P r2 X))))) DIM2 a {-# INLINE makeBordered2 #-} -makeBordered2 sh@(_ :. aHeight :. aWidth) borderWidth genInternal genBorder +makeBordered2 sh@(_ :. aHeight :. aWidth) borderWidth arrInternal arrBorder = let -- minimum and maximum indicies of values in the inner part of the image. !xMin = borderWidth @@ -83,26 +88,22 @@ !xMax = aWidth - borderWidth - 1 !yMax = aHeight - borderWidth - 1 - -- | Range of values where some of the data needed by the stencil is outside the image. - rectsBorder - = [ Rect (Z :. 0 :. 0) (Z :. yMin -1 :. aWidth - 1) -- bot - , Rect (Z :. yMax + 1 :. 0) (Z :. aHeight - 1 :. aWidth - 1) -- top - , Rect (Z :. yMin :. 0) (Z :. yMax :. xMin - 1) -- left - , Rect (Z :. yMin :. xMax + 1) (Z :. yMax :. aWidth - 1) ] -- right - {-# INLINE inBorder #-} - inBorder = not . inInternal - - -- Range of values where we don't need to worry about the border - rectsInternal - = [ Rect (Z :. yMin :. xMin) (Z :. yMax :. xMax ) ] - {-# INLINE inInternal #-} inInternal (Z :. y :. x) = x >= xMin && x <= xMax && y >= yMin && y <= yMax - in Array sh - [ Region (RangeRects inBorder rectsBorder) genInternal - , Region (RangeRects inInternal rectsInternal) genBorder ] + {-# INLINE inBorder #-} + inBorder = not . inInternal + in + -- internal region + APart sh (Range (Z :. yMin :. xMin) (Z :. yMax :. xMax ) inInternal) arrInternal + + -- border regions + $ APart sh (Range (Z :. 0 :. 0) (Z :. yMin -1 :. aWidth - 1) inBorder) arrBorder + $ APart sh (Range (Z :. yMax + 1 :. 0) (Z :. aHeight - 1 :. aWidth - 1) inBorder) arrBorder + $ APart sh (Range (Z :. yMin :. 0) (Z :. yMax :. xMin - 1) inBorder) arrBorder + $ APart sh (Range (Z :. yMin :. xMax + 1) (Z :. yMax :. aWidth - 1) inBorder) arrBorder + $ AUndefined sh diff --git a/Data/Array/Repa/Stencil.hs b/Data/Array/Repa/Stencil.hs --- a/Data/Array/Repa/Stencil.hs +++ b/Data/Array/Repa/Stencil.hs @@ -4,267 +4,16 @@ -- | Efficient computation of stencil based convolutions. -- --- This is specialised for stencils up to 7x7. --- Due to limitations in the GHC optimiser, using larger stencils doesn't work, and will yield `error` --- at runtime. We can probably increase the limit if required -- just ask. --- --- The focus of the stencil is in the center of the 7x7 tile, which has coordinates (0, 0). --- All coefficients in the stencil must fit in the tile, so they can be given X,Y coordinates up to --- +/- 3 positions. The stencil can be any shape, and need not be symmetric -- provided it fits in the 7x7 tile. --- module Data.Array.Repa.Stencil ( Stencil (..) , Boundary (..) -- * Stencil creation. - , makeStencil, makeStencil2 - - -- * Stencil operators. - , mapStencil2, forStencil2 - , mapStencilFrom2, forStencilFrom2 - - -- From Data.Array.Repa.Stencil.Template - , stencil2) + , makeStencil) where -import Data.Array.Repa as R -import Data.Array.Repa.Internals.Base as R +import Data.Array.Repa +import Data.Array.Repa.Base import Data.Array.Repa.Stencil.Base import Data.Array.Repa.Stencil.Template import Data.Array.Repa.Specialised.Dim2 -import qualified Data.Array.Repa.Shape as S -import qualified Data.Vector.Unboxed as V -import Data.List as List -import GHC.Exts -import GHC.Base -import Debug.Trace - --- | A index into the flat array. --- Should be abstract outside the stencil modules. -data Cursor - = Cursor Int - - --- Wrappers --------------------------------------------------------------------------------------- --- | Like `mapStencil2` but with the parameters flipped. -forStencil2 - :: Elt a - => Boundary a - -> Array DIM2 a - -> Stencil DIM2 a - -> Array DIM2 a - -{-# INLINE forStencil2 #-} -forStencil2 boundary arr stencil - = mapStencil2 boundary stencil arr - - --- | Like `mapStencilFrom2` but with the parameters flipped. -forStencilFrom2 - :: (Elt a, Elt b) - => Boundary a - -> Array DIM2 b - -> (b -> a) - -> Stencil DIM2 a - -> Array DIM2 a - -{-# INLINE forStencilFrom2 #-} -forStencilFrom2 boundary arr from stencil - = mapStencilFrom2 boundary stencil arr from - - --- | Apply a stencil to every element of a 2D array. --- The array must be manifest else `error`. -mapStencil2 - :: Elt a - => Boundary a - -> Stencil DIM2 a - -> Array DIM2 a - -> Array DIM2 a - -{-# INLINE mapStencil2 #-} -mapStencil2 boundary stencil arr - = mapStencilFrom2 boundary stencil arr id - - ---------------------------------------------------------------------------------------------------- --- | Apply a stencil to every element of a 2D array. --- The array must be manifest else `error`. -mapStencilFrom2 - :: (Elt a, Elt b) - => Boundary a -- ^ How to handle the boundary of the array. - -> Stencil DIM2 a -- ^ Stencil to apply. - -> Array DIM2 b -- ^ Array to apply stencil to. - -> (b -> a) -- ^ Apply this function to values read from the array before - -- transforming them with the stencil. - -> Array DIM2 a - -{-# INLINE mapStencilFrom2 #-} -mapStencilFrom2 boundary stencil@(StencilStatic sExtent zero load) arr preConvert - = let (_ :. aHeight :. aWidth) = extent arr - (_ :. sHeight :. sWidth) = sExtent - - sHeight2 = sHeight `div` 2 - sWidth2 = sWidth `div` 2 - - -- minimum and maximum indicies of values in the inner part of the image. - !xMin = sWidth2 - !yMin = sHeight2 - !xMax = aWidth - sWidth2 - 1 - !yMax = aHeight - sHeight2 - 1 - - -- Rectangles ----------------------- - -- range of values where we don't need to worry about the border - rectsInternal - = [ Rect (Z :. yMin :. xMin) (Z :. yMax :. xMax ) ] - - {-# INLINE inInternal #-} - inInternal (Z :. y :. x) - = x >= xMin && x <= xMax - && y >= yMin && y <= yMax - - - -- range of values where some of the data needed by the stencil is outside the image. - rectsBorder - = [ Rect (Z :. 0 :. 0) (Z :. yMin -1 :. aWidth - 1) -- bot - , Rect (Z :. yMax + 1 :. 0) (Z :. aHeight - 1 :. aWidth - 1) -- top - , Rect (Z :. yMin :. 0) (Z :. yMax :. xMin - 1) -- left - , Rect (Z :. yMin :. xMax + 1) (Z :. yMax :. aWidth - 1) ] -- right - - {-# INLINE inBorder #-} - inBorder = not . inInternal - - - -- Cursor functions ---------------- - {-# INLINE makeCursor' #-} - makeCursor' (Z :. y :. x) - = Cursor (x + y * aWidth) - - {-# INLINE shiftCursor' #-} - shiftCursor' ix (Cursor off) - = Cursor - $ case ix of - Z :. y :. x -> off + y * aWidth + x - - {-# INLINE getInner' #-} - getInner' cur - = unsafeAppStencilCursor2 shiftCursor' stencil - arr preConvert cur - - {-# INLINE getBorder' #-} - getBorder' cur - = case boundary of - BoundConst c -> c - BoundClamp -> unsafeAppStencilCursor2_clamp addDim stencil - arr preConvert cur - - in Array (extent arr) - [ Region (RangeRects inBorder rectsBorder) - (GenCursor id addDim getBorder') - - , Region (RangeRects inInternal rectsInternal) - (GenCursor makeCursor' shiftCursor' getInner') ] - - -unsafeAppStencilCursor2 - :: (Elt a, Elt b) - => (DIM2 -> Cursor -> Cursor) - -> Stencil DIM2 a - -> Array DIM2 b - -> (b -> a) - -> Cursor - -> a - -{-# INLINE [1] unsafeAppStencilCursor2 #-} -unsafeAppStencilCursor2 shift - stencil@(StencilStatic sExtent zero load) - arr@(Array aExtent [Region RangeAll (GenManifest vec)]) preConvert - cur@(Cursor off) - - | _ :. sHeight :. sWidth <- sExtent - , _ :. aHeight :. aWidth <- aExtent - , sHeight <= 7, sWidth <= 7 - = let - -- Get data from the manifest array. - {-# INLINE [0] getData #-} - getData (Cursor cur) = preConvert $ vec `V.unsafeIndex` cur - - -- Build a function to pass data from the array to our stencil. - {-# INLINE oload #-} - oload oy ox - = let !cur' = shift (Z :. oy :. ox) cur - in load (Z :. oy :. ox) (getData cur') - - in template7x7 oload zero - - --- | Like above, but clamp out of bounds array values to the closest real value. -unsafeAppStencilCursor2_clamp - :: forall a b. (Elt a, Elt b) - => (DIM2 -> DIM2 -> DIM2) - -> Stencil DIM2 a - -> Array DIM2 b - -> (b -> a) - -> DIM2 - -> a - -{-# INLINE [1] unsafeAppStencilCursor2_clamp #-} -unsafeAppStencilCursor2_clamp shift - stencil@(StencilStatic sExtent zero load) - arr@(Array aExtent [Region RangeAll (GenManifest vec)]) preConvert - cur - - | _ :. sHeight :. sWidth <- sExtent - , _ :. aHeight :. aWidth <- aExtent - , sHeight <= 7, sWidth <= 7 - = let - -- Get data from the manifest array. - {-# INLINE [0] getData #-} - getData :: DIM2 -> a - getData (Z :. y :. x) - = wrapLoadX x y - - -- TODO: Inlining this into above makes SpecConstr choke - wrapLoadX :: Int -> Int -> a - wrapLoadX !x !y - | x < 0 = wrapLoadY 0 y - | x >= aWidth = wrapLoadY (aWidth - 1) y - | otherwise = wrapLoadY x y - - {-# INLINE wrapLoadY #-} - wrapLoadY :: Int -> Int -> a - wrapLoadY !x !y - | y < 0 = loadXY x 0 - | y >= aHeight = loadXY x (aHeight - 1) - | otherwise = loadXY x y - - {-# INLINE loadXY #-} - loadXY :: Int -> Int -> a - loadXY !x !y - = preConvert $ vec `V.unsafeIndex` (x + y * aWidth) - - -- Build a function to pass data from the array to our stencil. - {-# INLINE oload #-} - oload oy ox - = let !cur' = shift (Z :. oy :. ox) cur - in load (Z :. oy :. ox) (getData cur') - - in template7x7 oload zero - - - --- | Data template for stencils up to 7x7. -template7x7 - :: (Int -> Int -> a -> a) - -> a -> a - -{-# INLINE [1] template7x7 #-} -template7x7 f zero - = f (-3) (-3) $ f (-3) (-2) $ f (-3) (-1) $ f (-3) 0 $ f (-3) 1 $ f (-3) 2 $ f (-3) 3 - $ f (-2) (-3) $ f (-2) (-2) $ f (-2) (-1) $ f (-2) 0 $ f (-2) 1 $ f (-2) 2 $ f (-2) 3 - $ f (-1) (-3) $ f (-1) (-2) $ f (-1) (-1) $ f (-1) 0 $ f (-1) 1 $ f (-1) 2 $ f (-1) 3 - $ f 0 (-3) $ f 0 (-2) $ f 0 (-1) $ f 0 0 $ f 0 1 $ f 0 2 $ f 0 3 - $ f 1 (-3) $ f 1 (-2) $ f 1 (-1) $ f 1 0 $ f 1 1 $ f 1 2 $ f 1 3 - $ f 2 (-3) $ f 2 (-2) $ f 2 (-1) $ f 2 0 $ f 2 1 $ f 2 2 $ f 2 3 - $ f 3 (-3) $ f 3 (-2) $ f 3 (-1) $ f 3 0 $ f 3 1 $ f 3 2 $ f 3 3 - $ zero diff --git a/Data/Array/Repa/Stencil/Base.hs b/Data/Array/Repa/Stencil/Base.hs --- a/Data/Array/Repa/Stencil/Base.hs +++ b/Data/Array/Repa/Stencil/Base.hs @@ -5,7 +5,6 @@ , Stencil (..) , makeStencil, makeStencil2) where -import Data.Array.Repa.Internals.Elt import Data.Array.Repa.Index -- | How to handle the case when the stencil lies partly outside the array. @@ -18,8 +17,8 @@ deriving (Show) --- | Represents a convolution stencil that we can apply to array. Only statically known stencils --- are supported right now. +-- | Represents a convolution stencil that we can apply to array. +-- Only statically known stencils are supported right now. data Stencil sh a -- | Static stencils are used when the coefficients are fixed, @@ -32,7 +31,7 @@ -- | Make a stencil from a function yielding coefficients at each index. makeStencil - :: (Elt a, Num a) + :: Num a => sh -- ^ Extent of stencil. -> (sh -> Maybe a) -- ^ Get the coefficient at this index. -> Stencil sh a @@ -48,7 +47,7 @@ -- | Wrapper for `makeStencil` that requires a DIM2 stencil. makeStencil2 - :: (Elt a, Num a) + :: Num a => Int -> Int -- ^ extent of stencil -> (DIM2 -> Maybe a) -- ^ Get the coefficient at this index. -> Stencil DIM2 a diff --git a/Data/Array/Repa/Stencil/Dim2.hs b/Data/Array/Repa/Stencil/Dim2.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Repa/Stencil/Dim2.hs @@ -0,0 +1,228 @@ +-- This is specialised for stencils up to 7x7. +-- Due to limitations in the GHC optimiser, using larger stencils doesn't +-- work, and will yield `error` at runtime. We can probably increase the +-- limit if required -- just ask. +-- +-- The focus of the stencil is in the center of the 7x7 tile, which has +-- coordinates (0, 0). All coefficients in the stencil must fit in the tile, +-- so they can be given X,Y coordinates up to +/- 3 positions. +-- The stencil can be any shape, and need not be symmetric -- provided it +-- fits in the 7x7 tile. +-- +module Data.Array.Repa.Stencil.Dim2 + ( + -- * Stencil creation + makeStencil2, stencil2 + + -- * Stencil operators + , PC5, mapStencil2, forStencil2) +where +import Data.Array.Repa +import Data.Array.Repa.Repr.Cursored +import Data.Array.Repa.Repr.Partitioned +import Data.Array.Repa.Repr.Undefined +import Data.Array.Repa.Stencil.Base +import Data.Array.Repa.Stencil.Template + +-- | A index into the flat array. +-- Should be abstract outside the stencil modules. +data Cursor + = Cursor Int + +type PC5 = P C (P D (P D (P D (P D X)))) + + +-- Wrappers ------------------------------------------------------------------- +-- | Like `mapStencil2` but with the parameters flipped. +forStencil2 + :: Repr r a + => Boundary a + -> Array r DIM2 a + -> Stencil DIM2 a + -> Array PC5 DIM2 a + +{-# INLINE forStencil2 #-} +forStencil2 boundary arr stencil + = mapStencil2 boundary stencil arr + + +------------------------------------------------------------------------------- +-- | Apply a stencil to every element of a 2D array. +mapStencil2 + :: Repr r a + => Boundary a -- ^ How to handle the boundary of the array. + -> Stencil DIM2 a -- ^ Stencil to apply. + -> Array r DIM2 a -- ^ Array to apply stencil to. + -> Array PC5 DIM2 a + +{-# INLINE mapStencil2 #-} +mapStencil2 boundary stencil@(StencilStatic sExtent _zero _load) arr + = let sh = extent arr + (_ :. aHeight :. aWidth) = sh + (_ :. sHeight :. sWidth) = sExtent + + sHeight2 = sHeight `div` 2 + sWidth2 = sWidth `div` 2 + + -- minimum and maximum indicies of values in the inner part of the image. + !xMin = sWidth2 + !yMin = sHeight2 + !xMax = aWidth - sWidth2 - 1 + !yMax = aHeight - sHeight2 - 1 + + {-# INLINE inInternal #-} + inInternal (Z :. y :. x) + = x >= xMin && x <= xMax + && y >= yMin && y <= yMax + + {-# INLINE inBorder #-} + inBorder = not . inInternal + + -- Cursor functions ---------------- + {-# INLINE makec #-} + makec (Z :. y :. x) + = Cursor (x + y * aWidth) + + {-# INLINE shiftc #-} + shiftc ix (Cursor off) + = Cursor + $ case ix of + Z :. y :. x -> off + y * aWidth + x + + {-# INLINE getInner' #-} + getInner' cur + = unsafeAppStencilCursor2 shiftc stencil arr cur + + {-# INLINE getBorder' #-} + getBorder' ix + = case boundary of + BoundConst c -> c + BoundClamp -> unsafeAppStencilCursor2_clamp addDim stencil + arr ix + + {-# INLINE arrInternal #-} + arrInternal = makeCursored (extent arr) makec shiftc getInner' + + {-# INLINE arrBorder #-} + arrBorder = fromFunction (extent arr) getBorder' + + in + -- internal region + APart sh (Range (Z :. yMin :. xMin) (Z :. yMax :. xMax ) inInternal) arrInternal + + -- border regions + $ APart sh (Range (Z :. 0 :. 0) (Z :. yMin -1 :. aWidth - 1) inBorder) arrBorder + $ APart sh (Range (Z :. yMax + 1 :. 0) (Z :. aHeight - 1 :. aWidth - 1) inBorder) arrBorder + $ APart sh (Range (Z :. yMin :. 0) (Z :. yMax :. xMin - 1) inBorder) arrBorder + $ APart sh (Range (Z :. yMin :. xMax + 1) (Z :. yMax :. aWidth - 1) inBorder) arrBorder + $ AUndefined sh + + +unsafeAppStencilCursor2 + :: Repr r a + => (DIM2 -> Cursor -> Cursor) + -> Stencil DIM2 a + -> Array r DIM2 a + -> Cursor + -> a + +{-# INLINE unsafeAppStencilCursor2 #-} +unsafeAppStencilCursor2 shift + (StencilStatic sExtent zero loads) + arr cur0 + + | _ :. sHeight :. sWidth <- sExtent + , sHeight <= 7, sWidth <= 7 + = let + -- Get data from the manifest array. + {-# INLINE getData #-} + getData (Cursor cur) = arr `unsafeLinearIndex` cur + + -- Build a function to pass data from the array to our stencil. + {-# INLINE oload #-} + oload oy ox + = let !cur' = shift (Z :. oy :. ox) cur0 + in loads (Z :. oy :. ox) (getData cur') + + in template7x7 oload zero + + | otherwise + = error $ unlines + [ "mapStencil2: Your stencil is too big for this method." + , " It must fit within a 7x7 tile to be compiled statically." ] + + +-- | Like above, but clamp out of bounds array values to the closest real value. +unsafeAppStencilCursor2_clamp + :: forall r a + . Repr r a + => (DIM2 -> DIM2 -> DIM2) + -> Stencil DIM2 a + -> Array r DIM2 a + -> DIM2 + -> a + +{-# INLINE unsafeAppStencilCursor2_clamp #-} +unsafeAppStencilCursor2_clamp shift + (StencilStatic sExtent zero loads) + arr cur + + | _ :. sHeight :. sWidth <- sExtent + , _ :. aHeight :. aWidth <- extent arr + , sHeight <= 7, sWidth <= 7 + = let + -- Get data from the manifest array. + {-# INLINE getData #-} + getData :: DIM2 -> a + getData (Z :. y :. x) + = wrapLoadX x y + + -- TODO: Inlining this into above makes SpecConstr choke + wrapLoadX :: Int -> Int -> a + wrapLoadX !x !y + | x < 0 = wrapLoadY 0 y + | x >= aWidth = wrapLoadY (aWidth - 1) y + | otherwise = wrapLoadY x y + + {-# INLINE wrapLoadY #-} + wrapLoadY :: Int -> Int -> a + wrapLoadY !x !y + | y < 0 = loadXY x 0 + | y >= aHeight = loadXY x (aHeight - 1) + | otherwise = loadXY x y + + {-# INLINE loadXY #-} + loadXY :: Int -> Int -> a + loadXY !x !y + = arr `unsafeIndex` (Z :. y :. x) + + -- Build a function to pass data from the array to our stencil. + {-# INLINE oload #-} + oload oy ox + = let !cur' = shift (Z :. oy :. ox) cur + in loads (Z :. oy :. ox) (getData cur') + + in template7x7 oload zero + + | otherwise + = error $ unlines + [ "mapStencil2: Your stencil is too big for this method." + , " It must fit within a 7x7 tile to be compiled statically." ] + + +-- | Data template for stencils up to 7x7. +template7x7 + :: (Int -> Int -> a -> a) + -> a -> a + +{-# INLINE template7x7 #-} +template7x7 f zero + = f (-3) (-3) $ f (-3) (-2) $ f (-3) (-1) $ f (-3) 0 $ f (-3) 1 $ f (-3) 2 $ f (-3) 3 + $ f (-2) (-3) $ f (-2) (-2) $ f (-2) (-1) $ f (-2) 0 $ f (-2) 1 $ f (-2) 2 $ f (-2) 3 + $ f (-1) (-3) $ f (-1) (-2) $ f (-1) (-1) $ f (-1) 0 $ f (-1) 1 $ f (-1) 2 $ f (-1) 3 + $ f 0 (-3) $ f 0 (-2) $ f 0 (-1) $ f 0 0 $ f 0 1 $ f 0 2 $ f 0 3 + $ f 1 (-3) $ f 1 (-2) $ f 1 (-1) $ f 1 0 $ f 1 1 $ f 1 2 $ f 1 3 + $ f 2 (-3) $ f 2 (-2) $ f 2 (-1) $ f 2 0 $ f 2 1 $ f 2 2 $ f 2 3 + $ f 3 (-3) $ f 3 (-2) $ f 3 (-1) $ f 3 0 $ f 3 1 $ f 3 2 $ f 3 3 + $ zero + diff --git a/Data/Array/Repa/Stencil/Template.hs b/Data/Array/Repa/Stencil/Template.hs --- a/Data/Array/Repa/Stencil/Template.hs +++ b/Data/Array/Repa/Stencil/Template.hs @@ -81,14 +81,17 @@ let fnCoeffs = LamE [VarP ix'] $ CaseE (VarE (mkName "ix")) - $ [ Match (InfixP (InfixP z' (mkName ":.") (LitP (IntegerL oy))) (mkName ":.") (LitP (IntegerL ox))) + $ [ Match (InfixP (InfixP z' (mkName ":.") (LitP (IntegerL oy))) + (mkName ":.") (LitP (IntegerL ox))) (NormalB $ ConE (mkName "Just") `AppE` LitE (IntegerL v)) [] | (oy, ox, v) <- coeffs ] ++ [Match WildP (NormalB $ ConE (mkName "Nothing")) []] return - $ AppE (VarE (mkName "makeStencil2") `AppE` (LitE (IntegerL sizeX)) `AppE` (LitE (IntegerL sizeY))) + $ AppE (VarE (mkName "makeStencil2") + `AppE` (LitE (IntegerL sizeX)) + `AppE` (LitE (IntegerL sizeY))) $ LetE [ PragmaD (InlineP (mkName "coeffs") (InlineSpec True False Nothing)) , ValD (VarP coeffs') (NormalB fnCoeffs) [] ] (VarE (mkName "coeffs")) diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright (c) 2010, University of New South Wales. +Copyright (c) 2010-2012, University of New South Wales. All rights reserved. Redistribution and use in source and binary forms, with or without diff --git a/repa.cabal b/repa.cabal --- a/repa.cabal +++ b/repa.cabal @@ -1,5 +1,5 @@ Name: repa -Version: 2.2.0.1 +Version: 3.0.0.1 License: BSD3 License-file: LICENSE Author: The DPH Team @@ -23,11 +23,12 @@ Library Build-Depends: - base == 4.4.*, + base == 4.5.*, ghc-prim == 0.2.*, vector == 0.9.*, + bytestring == 0.9.*, QuickCheck >= 2.3 && < 2.5, - template-haskell >= 2.5 && < 2.7 + template-haskell >= 2.5 && < 2.8 ghc-options: -Wall -fno-warn-missing-signatures @@ -35,32 +36,52 @@ -funbox-strict-fields -fcpr-off + extensions: + NoMonomorphismRestriction + ExplicitForAll + EmptyDataDecls + BangPatterns + TypeFamilies + MultiParamTypeClasses + FlexibleInstances + FlexibleContexts + StandaloneDeriving + ScopedTypeVariables + PatternGuards + Exposed-modules: - Data.Array.Repa + Data.Array.Repa.Eval.Gang + Data.Array.Repa.Operators.IndexSpace + Data.Array.Repa.Operators.Interleave + Data.Array.Repa.Operators.Mapping + Data.Array.Repa.Operators.Reduction + Data.Array.Repa.Operators.Selection + Data.Array.Repa.Operators.Traversal + Data.Array.Repa.Repr.ByteString + Data.Array.Repa.Repr.Cursored + Data.Array.Repa.Repr.Delayed + Data.Array.Repa.Repr.ForeignPtr + Data.Array.Repa.Repr.Partitioned + Data.Array.Repa.Repr.Unboxed + Data.Array.Repa.Repr.Undefined + Data.Array.Repa.Repr.Vector + Data.Array.Repa.Specialised.Dim2 + Data.Array.Repa.Stencil.Dim2 + Data.Array.Repa.Eval Data.Array.Repa.Index Data.Array.Repa.Shape Data.Array.Repa.Slice Data.Array.Repa.Stencil - Data.Array.Repa.Arbitrary - Data.Array.Repa.Properties - Data.Array.Repa.Specialised.Dim2 + Data.Array.Repa Other-modules: - Data.Array.Repa.Operators.IndexSpace - Data.Array.Repa.Operators.Traverse - Data.Array.Repa.Operators.Interleave - Data.Array.Repa.Operators.Mapping - Data.Array.Repa.Operators.Modify - Data.Array.Repa.Operators.Reduction - Data.Array.Repa.Operators.Select - Data.Array.Repa.Internals.Elt - Data.Array.Repa.Internals.Base - Data.Array.Repa.Internals.Gang - Data.Array.Repa.Internals.EvalChunked - Data.Array.Repa.Internals.EvalBlockwise - Data.Array.Repa.Internals.EvalCursored - Data.Array.Repa.Internals.EvalReduction - Data.Array.Repa.Internals.Forcing - Data.Array.Repa.Internals.Select + Data.Array.Repa.Eval.Chunked + Data.Array.Repa.Eval.Cursored + Data.Array.Repa.Eval.Elt + Data.Array.Repa.Eval.Fill + Data.Array.Repa.Eval.Reduction + Data.Array.Repa.Eval.Selection Data.Array.Repa.Stencil.Base Data.Array.Repa.Stencil.Template + Data.Array.Repa.Base +