diff --git a/Data/Array/Repa.hs b/Data/Array/Repa.hs
--- a/Data/Array/Repa.hs
+++ b/Data/Array/Repa.hs
@@ -3,11 +3,19 @@
 {-# OPTIONS -fno-warn-orphans #-}
 
 -- | See the repa-examples package for examples.
---   
+--
 --   More information at <http://repa.ouroborus.net>.
 --
 --   There is a draft tutorial at <http://www.haskell.org/haskellwiki/Numeric_Haskell:_A_Repa_Tutorial>
---  
+--
+-- @Release Notes:
+--  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.
+-- @
 module Data.Array.Repa
 	( module Data.Array.Repa.Shape
 	, module Data.Array.Repa.Index
@@ -35,10 +43,10 @@
 	, unsafeIndex
 
 	-- * Construction
-	, fromFunction	
+	, fromFunction
 	, fromVector
 	, fromList
-	
+
 	-- from Data.Array.Repa.Interlals.Forcing -------------------
 	-- * Forcing
 	, force, force2
@@ -61,6 +69,10 @@
 	, zipWith
 	, (+^), (-^), (*^), (/^)
 
+        -- from Data.Array.Repa.Operations.Modify -------------------
+        -- * Bulk updates
+        , (//)
+
 	-- from Data.Array.Repa.Operators.Reduction -----------------
 	-- * Reductions
 	, fold,	foldAll
@@ -82,11 +94,11 @@
 	, interleave2
 	, interleave3
 	, interleave4
-	
+
 	-- from Data.Array.Repa.Operators.Select --------------------
 	-- * Selection
 	, select)
-		
+
 where
 import Data.Array.Repa.Index
 import Data.Array.Repa.Slice
@@ -98,11 +110,12 @@
 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 Prelude				hiding (sum, map, zipWith, (++))
 import qualified Prelude		as P
 
 stage	= "Data.Array.Repa"
@@ -118,13 +131,13 @@
 instance (Shape sh, Elt a, Eq a) => Eq (Array sh a) where
 
 	{-# INLINE (==) #-}
-	(==) arr1  arr2 
-		= foldAll (&&) True 
-		$ reshape (Z :. (S.size $ extent arr1)) 
+	(==) arr1  arr2
+		= foldAll (&&) True
+		$ reshape (Z :. (S.size $ extent arr1))
 		$ zipWith (==) arr1 arr2
-		
+
 	{-# INLINE (/=) #-}
-	(/=) a1 a2 
+	(/=) a1 a2
 		= not $ (==) a1 a2
 
 -- Num
@@ -149,12 +162,12 @@
 	signum 		= map signum
 
 	{-# INLINE fromInteger #-}
-	fromInteger n	 = fromFunction failShape (\_ -> fromInteger n) 
+	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 
+withManifest
 	:: (Shape sh, Elt a)
 	=> (Array sh a -> b) -> Array sh a -> b
 
@@ -163,12 +176,12 @@
  = 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' 
+withManifest'
 	:: (Shape sh, Elt a)
 	=> Array sh a -> (Array sh a -> b) -> b
 
@@ -177,8 +190,8 @@
  = case arr of
 	Array sh [Region RangeAll (GenManifest vec)]
 	 -> vec `seq` f (Array sh [Region RangeAll (GenManifest vec)])
-	
+
 	_ -> f (force arr)
 
-	
+
 
diff --git a/Data/Array/Repa/Arbitrary.hs b/Data/Array/Repa/Arbitrary.hs
--- a/Data/Array/Repa/Arbitrary.hs
+++ b/Data/Array/Repa/Arbitrary.hs
@@ -22,24 +22,24 @@
 
 -- | Generate an arbitrary index, which may have 0's for some components.
 instance (Shape sh, Arbitrary sh) => Arbitrary (sh :. Int)  where
-	arbitrary 
+	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 
 
+		return	$ sh1 :. nMaxed
+
 -- | Generate an aribrary shape that does not have 0's for any component.
-arbitraryShape 
-	:: (Shape sh, Arbitrary sh) 
+arbitraryShape
+	:: (Shape sh, Arbitrary sh)
 	=> Gen (sh :. Int)
 
-arbitraryShape 
+arbitraryShape
  = do	sh1		<- arbitrary
 	let sh1Unit	= if size sh1 == 0 then unitDim else sh1
 
@@ -49,13 +49,13 @@
 	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, 
+
+
+-- | Generate an arbitrary shape where each dimension is more than zero,
 --	but less than a specific value.
-arbitrarySmallShape 
+arbitrarySmallShape
 	:: (Shape sh, Arbitrary sh)
 	=> Int
 	-> Gen (sh :. Int)
@@ -68,13 +68,13 @@
 		= case x `mod` maxDim of
 			0	-> 1
 			n	-> n
-						
-	return	$ if True 
+
+	return	$ if True
 			then shapeOfList $ map clamp dims
 			else sh
 
 
-arbitraryListOfLength 
+arbitraryListOfLength
 	:: Arbitrary a
 	=> Int -> Gen [a]
 
@@ -84,10 +84,10 @@
 	= 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 
+arbitrarySmallArray
 	:: (Shape sh, Elt a, Arbitrary sh, Arbitrary a)
 	=> Int
 	-> Gen (Array (sh :. Int) a)
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
@@ -2,7 +2,7 @@
 
 -- | Index types.
 module Data.Array.Repa.Index
-	( 
+	(
 	-- * Index types
 	  Z	(..)
 	, (:.)	(..)
@@ -80,7 +80,7 @@
 	{-# INLINE deepSeq #-}
 	deepSeq Z x		= x
 
-	
+
 instance Shape sh => Shape (sh :. Int) where
 	{-# INLINE rank #-}
 	rank   (sh  :. _)
@@ -93,7 +93,7 @@
 	unitDim = unitDim :. 1
 
 	{-# INLINE intersectDim #-}
-	intersectDim (sh1 :. n1) (sh2 :. n2) 
+	intersectDim (sh1 :. n1) (sh2 :. n2)
 		= (intersectDim sh1 sh2 :. (min n1 n2))
 
 	{-# INLINE addDim #-}
@@ -108,16 +108,16 @@
 	sizeIsValid (sh1 :. n)
 		| size sh1 > 0
 		= n <= maxBound `div` size sh1
-		
+
 		| otherwise
 		= False
-		
+
 	{-# INLINE toIndex #-}
-	toIndex (sh1 :. sh2) (sh1' :. sh2') 
+	toIndex (sh1 :. sh2) (sh1' :. sh2')
 		= toIndex sh1 sh1' * sh2 + sh2'
 
 	{-# INLINE fromIndex #-}
-	fromIndex (ds :. d) n 
+	fromIndex (ds :. d) n
 	 	= fromIndex ds (n `quotInt` d) :. r
 		where
 		-- If we assume that the index is in range, there is no point
@@ -128,7 +128,7 @@
 			| otherwise	= n `remInt` d
 
 	{-# INLINE inShapeRange #-}
-	inShapeRange (zs :. z) (sh1 :. n1) (sh2 :. n2) 
+	inShapeRange (zs :. z) (sh1 :. n1) (sh2 :. n2)
 		= (n2 >= z) && (n2 < n1) && (inShapeRange zs sh1 sh2)
 
 
@@ -138,14 +138,8 @@
 	shapeOfList xx
 	 = case xx of
 		[]	-> error $ stage ++ ".toList: empty list when converting to  (_ :. Int)"
-		x:xs	-> shapeOfList xs :. x			
+		x:xs	-> shapeOfList xs :. x
 
-	{-# INLINE deepSeq #-} 
+	{-# INLINE deepSeq #-}
 	deepSeq (sh :. n) x = deepSeq sh (n `seq` x)
-
-
-
-
-
-
 
diff --git a/Data/Array/Repa/Internals/Base.hs b/Data/Array/Repa/Internals/Base.hs
--- a/Data/Array/Repa/Internals/Base.hs
+++ b/Data/Array/Repa/Internals/Base.hs
@@ -10,7 +10,7 @@
 
 	, singleton, toScalar
 	, extent,    delay
-	
+
 	-- * Predicates
 	, inRange
 
@@ -20,7 +20,7 @@
 	, unsafeIndex
 
 	-- * Construction
-	, fromFunction	
+	, fromFunction
 	, fromVector
 	, fromList)
 where
@@ -33,9 +33,9 @@
 stage	= "Data.Array.Repa.Array"
 
 -- Array ----------------------------------------------------------------------
--- | Repa arrays. 
+-- | Repa arrays.
 data Array sh a
-	= Array 
+	= Array
 	{ -- | The entire extent of the array.
 	  arrayExtent		:: sh
 
@@ -45,7 +45,7 @@
 
 -- | Defines the values in a region of the array.
 data Region sh a
-	= Region 
+	= Region
 	{ -- | The range of elements this region applies to.
 	  regionRange		:: Range sh
 
@@ -54,7 +54,7 @@
 
 
 -- | Represents a range of elements in the array.
-data Range sh 
+data Range sh
 	  -- | Covers the entire array.
 	= RangeAll
 
@@ -71,17 +71,17 @@
 
 -- | Generates array elements for a particular region in the array.
 data Generator sh a
-	-- | Elements are already computed and sitting in this vector. 
+	-- | 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
 
@@ -121,7 +121,7 @@
 	x : xs	-> x `deepSeqArray` xs `deepSeqArrays` y
 
 -- | Ensure the structure for a region is fully evaluated.
-infixr 0 `deepSeqRegion` 
+infixr 0 `deepSeqRegion`
 deepSeqRegion :: Shape sh => Region sh a -> b -> b
 {-# INLINE deepSeqRegion #-}
 deepSeqRegion (Region range gen) x
@@ -192,11 +192,11 @@
 
 
 -- | Unpack an array into delayed form.
-delay 	:: (Shape sh, Elt a) 
-	=> Array sh a 
+delay 	:: (Shape sh, Elt a)
+	=> Array sh a
 	-> (sh, sh -> a)
 
-{-# INLINE delay #-}	
+{-# INLINE delay #-}
 delay arr@(Array sh _)
 	= (sh, (arr !))
 
@@ -208,7 +208,7 @@
 	:: forall sh a
 	.  (Shape sh, Elt a)
 	=> Array sh a
-	-> sh 
+	-> sh
 	-> a
 
 {-# INLINE (!) #-}
@@ -219,14 +219,14 @@
  = 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
 
 
@@ -235,7 +235,7 @@
 	 = case gen of
 		GenManifest vec
 		 -> vec V.! (S.toIndex sh ix')
-		
+
 		GenCursor makeCursor _ loadElem
 		 -> loadElem $ makeCursor ix'
 
@@ -254,7 +254,7 @@
 	:: forall sh a
 	.  (Shape sh, Elt a)
 	=> Array sh a
-	-> sh 
+	-> sh
 	-> Maybe a
 
 {-# INLINE (!?) #-}
@@ -266,15 +266,15 @@
  = 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
 
 
@@ -283,7 +283,7 @@
 	 = case gen of
 		GenManifest vec
 		 -> vec V.!? (S.toIndex sh ix')
-		
+
 		GenCursor makeCursor _ loadElem
 		 -> Just (loadElem $ makeCursor ix')
 
@@ -302,7 +302,7 @@
 	:: forall sh a
 	.  (Shape sh, Elt a)
 	=> Array sh a
-	-> sh 
+	-> sh
 	-> a
 
 {-# INLINE unsafeIndex #-}
@@ -310,14 +310,14 @@
  = 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 #-}
@@ -325,7 +325,7 @@
 	 = case gen of
 		GenManifest vec
 		 -> vec `V.unsafeIndex` (S.toIndex sh ix')
-		
+
 		GenCursor makeCursor _ loadElem
 		 -> loadElem $ makeCursor ix'
 
@@ -335,24 +335,24 @@
 
         unsafeIndex' (Array _ []) _
   	 = zero
-	
 
+
 -- Conversions ------------------------------------------------------------------------------------
 -- | Create a `Delayed` array from a function.
-fromFunction 
+fromFunction
 	:: Shape sh
 	=> sh
 	-> (sh -> a)
 	-> Array sh a
-	
+
 {-# INLINE fromFunction #-}
 fromFunction sh fnElems
 	= sh `S.deepSeq`
-	  Array sh [Region 
-			RangeAll 
+	  Array sh [Region
+			RangeAll
 			(GenCursor id addDim fnElems)]
 
--- | Create a `Manifest` array from an unboxed `Vector`. 
+-- | Create a `Manifest` array from an unboxed `Vector`.
 --	The elements are in row-major order.
 fromVector
 	:: Shape sh
@@ -374,7 +374,7 @@
 	=> sh
 	-> [a]
 	-> Array sh a
-	
+
 {-# INLINE fromList #-}
 fromList sh xx
 	| V.length vec /= S.size sh
@@ -382,7 +382,7 @@
 	 	[ 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)]
 
diff --git a/Data/Array/Repa/Internals/Elt.hs b/Data/Array/Repa/Internals/Elt.hs
--- a/Data/Array/Repa/Internals/Elt.hs
+++ b/Data/Array/Repa/Internals/Elt.hs
@@ -10,7 +10,7 @@
 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 #.
 
@@ -22,7 +22,7 @@
 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 
+	--   Doing a `seq` sometimes isn't enough, because the GHC simplifier can
 	--   erase these, and/or still move around the bindings.
 	touch :: a -> IO ()
 
@@ -50,7 +50,7 @@
 -- Floating -------------------------------------------------------------------
 instance Elt Float where
  {-# INLINE touch #-}
- touch (F# f) 
+ touch (F# f)
   = IO (\state -> case touch# f state of
 			state' -> (# state', () #))
 
@@ -77,7 +77,7 @@
 -- Int ------------------------------------------------------------------------
 instance Elt Int where
  {-# INLINE touch #-}
- touch (I# i) 
+ touch (I# i)
   = IO (\state -> case touch# i state of
 			state' -> (# state', () #))
 
@@ -89,7 +89,7 @@
 
 instance Elt Int8 where
  {-# INLINE touch #-}
- touch (I8# w) 
+ touch (I8# w)
   = IO (\state -> case touch# w state of
 			state' -> (# state', () #))
 
@@ -102,7 +102,7 @@
 
 instance Elt Int16 where
  {-# INLINE touch #-}
- touch (I16# w) 
+ touch (I16# w)
   = IO (\state -> case touch# w state of
 			state' -> (# state', () #))
 
@@ -113,10 +113,36 @@
  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) 
+ touch (W# i)
   = IO (\state -> case touch# i state of
 			state' -> (# state', () #))
 
@@ -129,7 +155,7 @@
 
 instance Elt Word8 where
  {-# INLINE touch #-}
- touch (W8# w) 
+ touch (W8# w)
   = IO (\state -> case touch# w state of
 			state' -> (# state', () #))
 
@@ -142,7 +168,7 @@
 
 instance Elt Word16 where
  {-# INLINE touch #-}
- touch (W16# w) 
+ touch (W16# w)
   = IO (\state -> case touch# w state of
 			state' -> (# state', () #))
 
@@ -153,13 +179,39 @@
  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) 
+ touch (a, b)
   = do	touch a
 	touch b
-	
+
  {-# INLINE zero #-}
  zero = (zero, zero)
 
@@ -169,11 +221,11 @@
 
 instance (Elt a, Elt b, Elt c) => Elt (a, b, c) where
  {-# INLINE touch #-}
- touch (a, b, c) 
+ touch (a, b, c)
   = do	touch a
 	touch b
 	touch c
-	
+
  {-# INLINE zero #-}
  zero = (zero, zero, zero)
 
@@ -183,12 +235,12 @@
 
 instance (Elt a, Elt b, Elt c, Elt d) => Elt (a, b, c, d) where
  {-# INLINE touch #-}
- touch (a, b, c, d) 
+ touch (a, b, c, d)
   = do	touch a
 	touch b
 	touch c
 	touch d
-	
+
  {-# INLINE zero #-}
  zero = (zero, zero, zero, zero)
 
@@ -198,13 +250,13 @@
 
 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) 
+ 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)
 
@@ -214,14 +266,14 @@
 
 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) 
+ 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)
 
diff --git a/Data/Array/Repa/Internals/EvalBlockwise.hs b/Data/Array/Repa/Internals/EvalBlockwise.hs
--- a/Data/Array/Repa/Internals/EvalBlockwise.hs
+++ b/Data/Array/Repa/Internals/EvalBlockwise.hs
@@ -12,35 +12,35 @@
 
 
 -- Blockwise filling ------------------------------------------------------------------------------
-fillVectorBlockwiseP 
+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 
+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        
+	fillBlock !ix
+	 = let	!x0	= colIx ix
 		!x1	= colIx (ix + 1)
 		!y0	= 0
 		!y1	= imageHeight
@@ -67,19 +67,19 @@
  = 	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 ()
@@ -98,7 +98,7 @@
 	=> 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			-- ^ 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)
@@ -108,18 +108,18 @@
 fillVectorBlock !vec !getElemFVB !imageWidth !x0 !y0 !x1 !y1
  = do	-- putStrLn $ "fillVectorBlock: " P.++ show (x0, y0, x1, y1)
 	fillBlock ixStart (ixStart + (x1 - x0))
- where	
+ 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
@@ -129,12 +129,12 @@
 			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
@@ -147,3 +147,4 @@
 	   	 | 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
--- a/Data/Array/Repa/Internals/EvalChunked.hs
+++ b/Data/Array/Repa/Internals/EvalChunked.hs
@@ -24,7 +24,7 @@
 fillChunkedS !vec !getElem
  = fill 0
  where 	!len	= VM.length vec
-	
+
 	fill !ix
 	 | ix >= len	= return ()
 	 | otherwise
@@ -38,13 +38,13 @@
 	=> IOVector a	-- ^ Vector to fill.
 	-> (Int -> a)	-- ^ Fn to get the value at a given index.
 	-> IO ()
-		
+
 {-# INLINE [0] fillChunkedP #-}
 fillChunkedP !vec !getElem
- = 	gangIO theGang 
+ = 	gangIO theGang
 	 $  \thread -> fill (splitIx thread) (splitIx (thread + 1))
 
- where	
+ 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.
@@ -57,10 +57,10 @@
 	splitIx thread
 	 | thread < chunkLeftover = thread * (chunkLen + 1)
 	 | otherwise		  = thread * chunkLen  + chunkLeftover
-	
+
 	-- Evaluate the elements of a single chunk.
 	{-# INLINE fill #-}
-	fill !ix !end 
+	fill !ix !end
 	 | ix >= end		= return ()
 	 | otherwise
 	 = do	VM.unsafeWrite vec ix (getElem ix)
diff --git a/Data/Array/Repa/Internals/EvalCursored.hs b/Data/Array/Repa/Internals/EvalCursored.hs
--- a/Data/Array/Repa/Internals/EvalCursored.hs
+++ b/Data/Array/Repa/Internals/EvalCursored.hs
@@ -30,26 +30,26 @@
 	-> IO ()
 
 {-# INLINE [0] fillCursoredBlock2P #-}
-fillCursoredBlock2P 
+fillCursoredBlock2P
 	!vec
 	!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 ()
@@ -58,8 +58,8 @@
 		!x1'	= colIx (ix + 1) - 1
 		!y0'	= y0
 		!y1'	= y1
-	   in	fillCursoredBlock2 
-			vec 
+	   in	fillCursoredBlock2
+			vec
 			makeCursorFCB shiftCursorFCB getElemFCB
 			imageWidth x0' y0' x1' y1'
 
@@ -73,15 +73,15 @@
 	-> (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				-- ^ 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 
-	!vec 
+fillCursoredBlock2
+	!vec
 	!makeCursor !shiftCursor !getElem
 	!imageWidth !x0 !y0 !x1 !y1
 
@@ -93,7 +93,7 @@
 	 | otherwise
 	 = do	fillLine4 x0
 		fillBlock (y + 1)
-	
+
 	 where	{-# INLINE fillLine4 #-}
 		fillLine4 !x
  	   	 | x + 4 > x1 		= fillLine1 x
@@ -110,7 +110,7 @@
 			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
@@ -121,15 +121,15 @@
 			touch val3
 
 			-- Compute cursor into destination array.
-			let !dstCur0	= x + y * imageWidth				
+			let !dstCur0	= x + y * imageWidth
 			VM.unsafeWrite vec (dstCur0)     val0
 			VM.unsafeWrite vec (dstCur0 + 1) val1
 			VM.unsafeWrite vec (dstCur0 + 2) val2
 			VM.unsafeWrite vec (dstCur0 + 3) val3
 			fillLine4 (x + 4)
-		
+
 		{-# INLINE fillLine1 #-}
-		fillLine1 !x 
+		fillLine1 !x
  	   	 | x > x1		= return ()
 	   	 | otherwise
 	   	 = do	VM.unsafeWrite vec (x + y * imageWidth) (getElem $ makeCursor (Z :. y :. x))
diff --git a/Data/Array/Repa/Internals/EvalReduction.hs b/Data/Array/Repa/Internals/EvalReduction.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Repa/Internals/EvalReduction.hs
@@ -0,0 +1,118 @@
+{-# 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 )
+
+
+-- | 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
+      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 = do
+  mvec <- M.unsafeNew threads
+  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
+
+    {-# 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
--- a/Data/Array/Repa/Internals/Forcing.hs
+++ b/Data/Array/Repa/Internals/Forcing.hs
@@ -23,7 +23,7 @@
 --	The elements come out in row-major order.
 toVector
 	:: (Shape sh, Elt a)
-	=> Array sh a 
+	=> Array sh a
 	-> Vector a
 {-# INLINE toVector #-}
 toVector arr
@@ -48,13 +48,13 @@
 force	:: (Shape sh, Elt a)
 	=> Array sh a -> Array sh a
 
-{-# INLINE [2] force #-}	
+{-# INLINE [2] force #-}
 force arr
  = unsafePerformIO
  $ do	(sh, vec)	<- forceIO arr
-	return $ sh `seq` vec `seq` 
+	return $ sh `seq` vec `seq`
 		 Array sh [Region RangeAll (GenManifest vec)]
-	
+
  where	forceIO arr'
 	 = case arr' of
 		-- Don't force an already forced array.
@@ -66,20 +66,20 @@
 			fillChunkedP mvec (\ix -> arr' `unsafeIndex` fromIndex sh ix)
 			vec	<- V.unsafeFreeze mvec
 			return	(sh, vec)
-		
 
+
 -- | 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 #-}	
+{-# INLINE [2] force2 #-}
 force2 arr
- = unsafePerformIO 
+ = unsafePerformIO
  $ do	(sh, vec)	<- forceIO2 arr
-	return $ sh `seq` vec `seq` 
+	return $ sh `seq` vec `seq`
 		 Array sh [Region RangeAll (GenManifest vec)]
 
  where	forceIO2 arr'
- 	 = arr' `deepSeqArray` 
+ 	 = arr' `deepSeqArray`
 	   case arr' of
 		-- Don't force an already forced array.
 		Array sh [Region RangeAll (GenManifest vec)]
@@ -101,16 +101,16 @@
 			fillRegion2P mvec sh r2
 			vec	<- V.unsafeFreeze mvec
 			return (sh, vec)
-	
+
 		-- Create a vector to hold the new array and load in the regions.
 		Array sh regions
 	 	 -> do	mvec	<- VM.new (S.size sh)
 			mapM_ (fillRegion2P mvec sh) regions
 			vec	<- V.unsafeFreeze mvec
 			return (sh, vec)
-			
 
--- FillRegion2P -----------------------------------------------------------------------------------	
+
+-- 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
@@ -118,24 +118,24 @@
 --   access their source elements from the local neighbourhood, this specialised version
 --   should given better cache performance than plain `fillRegionP`.
 --
-fillRegion2P 
+fillRegion2P
 	:: Elt a
 	=> VM.IOVector a	-- ^ Vector to write elements into.
 	-> DIM2			-- ^ Extent of entire array.
 	-> Region DIM2 a	-- ^ Region to fill.
 	-> IO ()
-	
+
 {-# INLINE [1] fillRegion2P #-}
 fillRegion2P mvec sh@(_ :. height :. width) (Region range gen)
  = mvec `seq` height `seq` width `seq`
-   case range of 
-	RangeAll	
-	 -> fillRect2 mvec sh gen 
-		(Rect 	(Z :. 0          :. 0) 
+   case range of
+	RangeAll
+	 -> fillRect2 mvec sh gen
+		(Rect 	(Z :. 0          :. 0)
 			(Z :. height - 1 :. width - 1))
 
 	RangeRects _ [rect]
-	 -> fillRect2 mvec sh gen rect 
+	 -> fillRect2 mvec sh gen rect
 
 	-- Specialise for the common case of 4 rectangles so we get fusion.
 	-- The following case with mapM_ doesn't fuse because mapM_ isn't completely unrolled.
@@ -148,9 +148,9 @@
 	RangeRects _ rects
 	 -> mapM_ (fillRect2 mvec sh gen) rects
 
-		
+
 -- | Fill a rectangle in a vector.
-fillRect2 
+fillRect2
 	:: Elt a
 	=> VM.IOVector a	-- ^ Vector to write elements into.
 	-> DIM2 		-- ^ Extent of entire array.
@@ -158,15 +158,16 @@
 	-> Rect DIM2		-- ^ Rectangle to fill.
 	-> IO ()
 
-{-# INLINE fillRect2 #-}	
-fillRect2 mvec (_ :. _ :. width) gen (Rect (Z :. y0 :. x0) (Z :. y1 :. x1)) 
- = mvec `seq` width `seq` y0 `seq` x0 `seq` y1 `seq` x1 `seq` 
+{-# INLINE fillRect2 #-}
+fillRect2 mvec (_ :. _ :. width) gen (Rect (Z :. y0 :. x0) (Z :. y1 :. x1))
+ = mvec `seq` width `seq` y0 `seq` x0 `seq` y1 `seq` x1 `seq`
    case gen of
 	GenManifest{}
 	 -> error "fillRegion2P: GenManifest, copy elements."
-	
+
 	-- Cursor based arrays.
 	GenCursor makeCursor shiftCursor loadElem
          -> fillCursoredBlock2P mvec
 		makeCursor shiftCursor loadElem
 		width x0 y0 x1 y1
+
diff --git a/Data/Array/Repa/Internals/Gang.hs b/Data/Array/Repa/Internals/Gang.hs
--- a/Data/Array/Repa/Internals/Gang.hs
+++ b/Data/Array/Repa/Internals/Gang.hs
@@ -2,12 +2,12 @@
 
 -- | Gang Primitives.
 --   Based on DPH code by Roman Leshchinskiy
--- 
+--
 --   Gang primitives.
 --
 #define TRACE_GANG 0
 
-module Data.Array.Repa.Internals.Gang 
+module Data.Array.Repa.Internals.Gang
 	( Gang, seqGang, forkGang, gangSize, gangIO, gangST, traceGang, traceGangST
 	, theGang)
 where
@@ -35,21 +35,21 @@
 
 
 -- Requests ---------------------------------------------------------------------------------------
--- | The 'Req' type encapsulates work requests for individual members of a gang. 
-data Req 
+-- | 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 equest down by writing to the MVar before
-	--   returning to its caller (forkGang) 	
+        --   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 
+newReq p
  = do	mv	<- newEmptyMVar
 	return	$ ReqDo p mv
 
@@ -65,15 +65,15 @@
 
 -- Gang ------------------------------------------------------------------------------------------
 -- | A 'Gang' is a group of threads which execute arbitrary work requests.
---   To get the gang to do work, write Req-uest valuesto its MVars
-data Gang 
+--   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 _ _) 
+  showsPrec p (Gang n _ _)
 	= showString "<<"
         . showsPrec p n
         . showString " threads>>"
@@ -90,7 +90,7 @@
 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"
@@ -98,7 +98,7 @@
 		action threadId
 		end 	<- getGangTime
 		traceGang $ "Worker " ++ show threadId ++ " end (" ++ diffTime start end ++ ")"
-		
+
 		putMVar varDone ()
 		gangWorker threadId varReq
 
@@ -113,18 +113,18 @@
 --     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 
+--
+--   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) 
+	putMVar varReq (ReqShutdown varDone)
 	takeMVar varDone
 	return ()
 
@@ -132,21 +132,21 @@
 -- | Fork a 'Gang' with the given number of threads (at least 1).
 forkGang :: Int -> IO Gang
 forkGang n
- = assert (n > 0) 
- $ do	
+ = 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_ forkOnIO [0..] 
+	zipWithM_ forkOnIO [0..]
 		$ zipWith gangWorker [0 .. n-1] mvs
 
 	-- The gang is currently idle.
 	busy	<- newMVar False
-	
+
 	return $ Gang n mvs busy
 
 
@@ -157,7 +157,7 @@
 
 -- | 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. 
+--   requesting thread.
 --
 --   TODO: We might want to print a configurable warning that this is happening.
 --
@@ -166,7 +166,7 @@
 	-> IO ()
 
 {-# NOINLINE gangIO #-}
-gangIO (Gang n mvs busy) p 
+gangIO (Gang n mvs busy) p
  = do	traceGang   "gangIO: issuing work requests (SEQ_IF_GANG_BUSY)"
 	b <- swapMVar busy True
 
@@ -180,7 +180,7 @@
 				, "  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
@@ -196,7 +196,7 @@
 				--   the particular worker thread it's running on.
 	-> IO ()
 
-parIO n mvs p 
+parIO n mvs p
  = do	traceGang "parIO: begin"
 
 	start 	<- getGangTime
@@ -246,3 +246,4 @@
 
 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
--- a/Data/Array/Repa/Internals/Select.hs
+++ b/Data/Array/Repa/Internals/Select.hs
@@ -12,7 +12,7 @@
 import Data.IORef
 
 -- | Select indices matching a predicate
-selectChunkedS 
+selectChunkedS
 	:: (Shape sh, Unbox a)
 	=> (sh -> Bool)		-- ^ See if this predicate matches.
 	-> (sh -> a)		-- ^  .. and apply fn to the matching index
@@ -25,13 +25,13 @@
  = 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 
+	 , match ixSrc
 	 = do	VM.unsafeWrite vDst nDst (produce ixSrc)
 		fill (nSrc + 1) (nDst + 1)
 
@@ -41,9 +41,9 @@
 
 -- | 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 
+--   The number of elements in the result array depends on how many threads
 --   you're running the program with.
-selectChunkedP 
+selectChunkedP
 	:: forall a
 	.  Unbox a
 	=> (Int -> Bool)	-- ^ See if this predicate matches.
@@ -53,11 +53,11 @@
 
 {-# INLINE selectChunkedP #-}
 selectChunkedP !match !produce !len
- = do	
+ = 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 
+	refs	<- P.replicateM threads
 		$ do	vec	<- VM.new $ len `div` threads
 			newIORef vec
 
@@ -66,12 +66,12 @@
 	 $ \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
@@ -89,29 +89,29 @@
 	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)		
+		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 
+         -- If we've finished selecting elements, then slice the vector down
          -- so it doesn't have any empty space at the end.
-	 | ixSrc >= ixSrcEnd	
+	 | 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
@@ -24,10 +24,10 @@
 -- 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) 
+--
+--   TODO: This only works for arrays with a single region.
+--
+reshape	:: (Shape sh, Shape sh', Elt a)
 	=> sh'
 	-> Array sh a
 	-> Array sh' a
@@ -40,34 +40,34 @@
 reshape sh' (Array sh [Region RangeAll gen])
  = Array sh' [Region RangeAll gen']
  where gen' = case gen of
-		GenManifest vec 
+		GenManifest vec
 	 	 -> GenManifest vec
 
 		GenCursor makeCursor _ loadElem
-	 	 -> GenCursor 
+	 	 -> GenCursor
 			id
 			addDim
 			(loadElem . makeCursor . fromIndex sh . toIndex sh')
 
 reshape _ _
 	= error $ stage P.++ ".reshape: can't reshape a partitioned array"
-	
 
+
 -- | Append two arrays.
 --
-append, (++)	
+append, (++)
 	:: (Shape sh, Elt a)
 	=> Array (sh :. Int) a
 	-> Array (sh :. Int) a
 	-> Array (sh :. Int) a
 
 {-# INLINE append #-}
-append arr1 arr2 
+append arr1 arr2
  = unsafeTraverse2 arr1 arr2 fnExtent fnElem
  where
  	(_ :. n) 	= extent arr1
 
-	fnExtent (sh :. i) (_  :. j) 
+	fnExtent (sh :. i) (_  :. j)
 		= sh :. (i + j)
 
 	fnElem f1 f2 (sh :. i)
@@ -78,15 +78,15 @@
 (++) arr1 arr2 = append arr1 arr2
 
 
--- | Transpose the lowest two dimensions of an array. 
+-- | Transpose the lowest two dimensions of an array.
 --	Transposing an array twice yields the original.
-transpose 
-	:: (Shape sh, Elt a) 
+transpose
+	:: (Shape sh, Elt a)
 	=> Array (sh :. Int :. Int) a
 	-> Array (sh :. Int :. Int) a
 
 {-# INLINE transpose #-}
-transpose arr 
+transpose arr
  = unsafeTraverse arr
 	(\(sh :. m :. n) 	-> (sh :. n :.m))
 	(\f -> \(sh :. i :. j) 	-> f (sh :. j :. i))
@@ -105,8 +105,8 @@
 
 {-# INLINE extend #-}
 extend sl arr
-	= backpermute 
-		(fullOfSlice sl (extent arr)) 
+	= backpermute
+		(fullOfSlice sl (extent arr))
 		(sliceOfFull sl)
 		arr
 
@@ -121,7 +121,7 @@
 
 {-# INLINE slice #-}
 slice arr sl
-	= backpermute 
+	= backpermute
 		(sliceOfFull sl (extent arr))
 		(fullOfSlice sl)
 		arr
@@ -131,7 +131,7 @@
 --	The result array has the same extent as the original.
 backpermute
 	:: forall sh sh' a
-	.  (Shape sh, Shape sh', Elt 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.
@@ -140,15 +140,15 @@
 
 {-# INLINE backpermute #-}
 backpermute newExtent perm arr
-	= traverse arr (const newExtent) (. perm) 
-	
+	= traverse 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) 
+	.  (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.
@@ -158,7 +158,8 @@
 {-# INLINE backpermuteDft #-}
 backpermuteDft arrDft fnIndex arrSrc
 	= fromFunction (extent arrDft) fnElem
-	where	fnElem ix	
+	where	fnElem ix
 		 = case fnIndex ix of
 			Just ix'	-> arrSrc ! ix'
 			Nothing		-> arrDft ! 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
@@ -12,9 +12,9 @@
 import Data.Array.Repa.Operators.Traverse
 import Data.Array.Repa.Shape			as S
 
--- | Interleave the elements of two arrays. 
+-- | Interleave the elements of two arrays.
 --   All the input arrays must have the same extent, else `error`.
---   The lowest dimenion of the result array is twice the size of the inputs.
+--   The lowest dimension of the result array is twice the size of the inputs.
 --
 -- @
 --  interleave2 a1 a2   b1 b2  =>  a1 b1 a2 b2
@@ -26,7 +26,7 @@
 	=> Array (sh :. Int) a
 	-> Array (sh :. Int) a
 	-> Array (sh :. Int) a
-	
+
 {-# INLINE interleave2 #-}
 interleave2 arr1 arr2
  = arr1 `deepSeqArray` arr2 `deepSeqArray`
@@ -36,10 +36,10 @@
 	 | dim1 == dim2
 	 , sh :. len	<- dim1
 	 = sh :. (len * 2)
-	
+
 	 | otherwise
 	 = error "Data.Array.Repa.interleave2: arrays must have same extent"
-		
+
 	elemFn get1 get2 (sh :. ix)
 	 = case ix `mod` 3 of
 		0	-> get1 (sh :. ix `div` 2)
@@ -47,14 +47,14 @@
 		_	-> error "Data.Array.Repa.interleave2: this never happens :-P"
 
 
--- | Interleave the elements of three arrays. 
+-- | 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
-	
+
 {-# INLINE interleave3 #-}
 interleave3 arr1 arr2 arr3
  = arr1 `deepSeqArray` arr2 `deepSeqArray` arr3 `deepSeqArray`
@@ -65,10 +65,10 @@
 	 , dim1 == dim3
 	 , sh :. len	<- dim1
 	 = sh :. (len * 3)
-	
+
 	 | otherwise
 	 = error "Data.Array.Repa.interleave3: arrays must have same extent"
-		
+
 	elemFn get1 get2 get3 (sh :. ix)
 	 = case ix `mod` 3 of
 		0	-> get1 (sh :. ix `div` 3)
@@ -77,7 +77,7 @@
 		_	-> error "Data.Array.Repa.interleave3: this never happens :-P"
 
 
--- | Interleave the elements of four arrays. 
+-- | Interleave the elements of four arrays.
 interleave4
 	:: (Shape sh, Elt a)
 	=> Array (sh :. Int) a
@@ -85,7 +85,7 @@
 	-> Array (sh :. Int) a
 	-> Array (sh :. Int) a
 	-> Array (sh :. Int) a
-	
+
 {-# INLINE interleave4 #-}
 interleave4 arr1 arr2 arr3 arr4
  = arr1 `deepSeqArray` arr2 `deepSeqArray` arr3 `deepSeqArray` arr4 `deepSeqArray`
@@ -97,10 +97,10 @@
 	 , dim1 == dim4
 	 , sh :. len	<- dim1
 	 = sh :. (len * 4)
-	
+
 	 | otherwise
 	 = error "Data.Array.Repa.interleave4: arrays must have same extent"
-		
+
 	elemFn get1 get2 get3 get4 (sh :. ix)
 	 = case ix `mod` 4 of
 		0	-> get1 (sh :. ix `div` 4)
@@ -108,3 +108,4 @@
 		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
@@ -20,7 +20,7 @@
 --
 --   This is specialised for arrays of up to four regions, using more breaks fusion.
 --
-map	:: (Shape sh, Elt a, Elt b) 
+map	:: (Shape sh, Elt a, Elt b)
 	=> (a -> b)
 	-> Array sh a
 	-> Array sh b
@@ -38,16 +38,16 @@
 		[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
@@ -56,17 +56,17 @@
 			P.id
 			addDim
 		 	(\ix -> f $ V.unsafeIndex vec $ S.toIndex sh ix)
-				
+
 		GenCursor makeCursor shiftCursor loadElem
 		 -> GenCursor makeCursor shiftCursor (f . loadElem)
 
 
 -- | Combine two arrays, element-wise, with a binary operator.
---	If the extent of the two array arguments differ, 
+--	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) 
+zipWith :: (Shape sh, Elt a, Elt b, Elt c)
+	=> (a -> b -> c)
 	-> Array sh a
 	-> Array sh b
 	-> Array sh c
@@ -82,7 +82,7 @@
 
 		{-# 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') ]
@@ -106,3 +106,4 @@
 
 {-# INLINE (/^) #-}
 (/^)	= zipWith (/)
+
diff --git a/Data/Array/Repa/Operators/Modify.hs b/Data/Array/Repa/Operators/Modify.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Repa/Operators/Modify.hs
@@ -0,0 +1,53 @@
+{-# 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,5 +1,5 @@
 {-# OPTIONS_HADDOCK hide #-}
-{-# LANGUAGE ExplicitForAll, TypeOperators #-}
+{-# LANGUAGE BangPatterns, ExplicitForAll, TypeOperators #-}
 
 module Data.Array.Repa.Operators.Reduction
 	( fold, foldAll
@@ -8,51 +8,67 @@
 import Data.Array.Repa.Index
 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 Prelude				hiding (sum)
+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 System.IO.Unsafe
 
--- | Sequentially fold the innermost dimension of an array.
---	Combine this with `transpose` to fold any other dimension.
+
+-- | 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)
 	=> (a -> a -> a)
-	-> a 
+	-> a
 	-> Array (sh :. Int) a
 	-> Array sh a
+{-# INLINE [1] fold #-}
+fold 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)
+           --
+           1 -> let !x = V.singleton $ foldAll f z arr
+                in  Array sz [Region RangeAll (GenManifest x)]
 
-{-# INLINE fold #-}
-fold f x arr
- = x `seq` arr `deepSeqArray` 
-   let	sh' :. n	= extent arr
-	elemFn i 	= V.foldl' f x
-			$ V.map	(\ix -> arr ! (i :. ix)) 
-				(V.enumFromTo 0 (n - 1))
-   in	fromFunction sh' elemFn
+           _ -> unsafePerformIO 
+              $ do mvec   <- M.unsafeNew (S.size sz)
+                   foldP mvec (\ix -> arr `unsafeIndex` fromIndex sh ix) f z n
+                   !vec   <- V.unsafeFreeze mvec
+                   return $ Array sz [Region RangeAll (GenManifest vec)]
 
 
--- | Sequentially fold all the elements of an array.
+-- | 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)
 	=> (a -> a -> a)
 	-> a
 	-> Array sh a
 	-> a
-	
-{-# INLINE foldAll #-}
-foldAll f x arr
-	= V.foldl' f x
-	$ V.map ((arr !) . (S.fromIndex (extent arr)))
-	$ V.enumFromTo
-		0
-		((S.size $ extent arr) - 1)
-
+{-# INLINE [1] foldAll #-}
+foldAll f z arr 
+ = let  sh = extent arr
+        n  = size sh
+   in   unsafePerformIO $ 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
 
@@ -61,12 +77,6 @@
 sumAll	:: (Shape sh, Elt a, Num a)
 	=> Array sh a
 	-> a
-
 {-# INLINE sumAll #-}
-sumAll arr
-	= V.foldl' (+) 0
-	$ V.map ((arr !) . (S.fromIndex (extent arr)))
-	$ V.enumFromTo
-		0
-		((S.size $ extent arr) - 1)
+sumAll arr = foldAll (+) 0 arr
 
diff --git a/Data/Array/Repa/Operators/Select.hs b/Data/Array/Repa/Operators/Select.hs
--- a/Data/Array/Repa/Operators/Select.hs
+++ b/Data/Array/Repa/Operators/Select.hs
@@ -16,22 +16,22 @@
 --   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. 
+--   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 -> 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` 
+ = 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
@@ -39,6 +39,6 @@
 
 		-- TODO: avoid copy.
 		let result	= V.concat vecs'
-		
+
 		return	(Z :. V.length result, result)
-		
+
diff --git a/Data/Array/Repa/Operators/Traverse.hs b/Data/Array/Repa/Operators/Traverse.hs
--- a/Data/Array/Repa/Operators/Traverse.hs
+++ b/Data/Array/Repa/Operators/Traverse.hs
@@ -18,19 +18,19 @@
 	.  (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. 
+	-> ((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` 
+ 	= arr `deepSeqArray`
           fromFunction (transExtent (extent arr)) (newElem (arr !))
 
 
 {-# INLINE unsafeTraverse #-}
 unsafeTraverse arr transExtent newElem
- 	= arr `deepSeqArray` 
+ 	= arr `deepSeqArray`
 	  fromFunction (transExtent (extent arr)) (newElem (unsafeIndex arr))
 
 
@@ -38,89 +38,89 @@
 traverse2, unsafeTraverse2
 	:: forall sh sh' sh'' a b c
 	.  ( Shape sh, Shape sh', Shape sh''
-	   , Elt a,    Elt b,     Elt c)
+	   , 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 -> a) -> (sh' -> b)
                       -> (sh'' -> c))		-- ^ Function to produce elements of the result.
-						--   It is passed lookup functions to get elements of the 
+						--   It is passed lookup functions to get elements of the
 						--   source arrays.
-        -> Array sh'' c 
+        -> Array sh'' c
 
 {-# INLINE traverse2 #-}
 traverse2 arrA arrB transExtent newElem
 	= arrA `deepSeqArray` arrB `deepSeqArray`
    	  fromFunction
-		(transExtent (extent arrA) (extent arrB)) 
+		(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)) 
+		(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 
+	          a   b   c   d
 	.  ( Shape sh1, Shape sh2, Shape sh3, Shape sh4
-	   , Elt a,     Elt b,     Elt c,     Elt d)
-        => Array sh1 a 		
-	-> Array sh2 b			
-	-> Array sh3 c			
-        -> (sh1 -> sh2 -> sh3 -> sh4)	
-        -> (  (sh1 -> a) -> (sh2 -> b) 
+	   , 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 )		
+           ->  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)) 
+		(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)) 
+		(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 
+	:: 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,     Elt e)
-        => Array sh1 a 			
-	-> Array sh2 b			
-	-> Array sh3 c			
-	-> Array sh4 d				
-        -> (sh1 -> sh2 -> sh3 -> sh4 -> sh5 )	
-        -> (  (sh1 -> a) -> (sh2 -> b) 
+	   , 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 
+           ->  sh5 -> e )
+        -> Array sh5 e
 
 {-# INLINE traverse4 #-}
 traverse4 arrA arrB arrC arrD transExtent newElem
-	= arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray` arrD `deepSeqArray` 
+	= arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray` arrD `deepSeqArray`
    	  fromFunction
-		(transExtent (extent arrA) (extent arrB) (extent arrC) (extent arrD)) 
+		(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` 
+	= arrA `deepSeqArray` arrB `deepSeqArray` arrC `deepSeqArray` arrD `deepSeqArray`
    	  fromFunction
-		(transExtent (extent arrA) (extent arrB) (extent arrC) (extent arrD)) 
+		(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
--- a/Data/Array/Repa/Properties.hs
+++ b/Data/Array/Repa/Properties.hs
@@ -5,8 +5,9 @@
 	, props_DataArrayRepa)
 where
 import Data.Array.Repa			as R
-import qualified Data.Array.Repa.Shape	as S
 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
@@ -19,7 +20,7 @@
 props_DataArrayRepaIndex :: [(String, Property)]
 props_DataArrayRepaIndex
   = [(stage P.++ "." P.++ name, test) | (name, test)
-     <-	[ ("toIndexFromIndex/DIM1", 	property prop_toIndexFromIndex_DIM1) 
+     <-	[ ("toIndexFromIndex/DIM1", 	property prop_toIndexFromIndex_DIM1)
 	, ("toIndexFromIndex/DIM2", 	property prop_toIndexFromIndex_DIM2) ]]
 
 prop_toIndexFromIndex_DIM1 sh ix
@@ -34,9 +35,9 @@
    	forAll (genInShape2 sh) $ \(ix :: DIM2) ->
 	fromIndex sh (toIndex sh ix) == ix
 
-	
-	
 
+
+
 -- Data.Array.Repa --------------------------------------------------------------------------------
 -- | QuickCheck properties for "Data.Array.Repa" and its children.
 props_DataArrayRepa :: [(String, Property)]
@@ -45,18 +46,19 @@
  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_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) ->
 	arr == force arr
-	
+
 prop_id_toScalarUnit (x :: Int)
  =	toScalar (singleton x) == x
 
@@ -84,9 +86,18 @@
 	sumAll (append arr1 arr1) == (2 * sumAll arr1)
 
 -- Reductions --------------------------
+prop_sumIsSum_DIM3
+  = forAll (arbitrarySmallArray 20)                     $ \(arr :: Array DIM3 Float) ->
+    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 == fromFunction sh elemFn
+
 prop_sumAllIsSum_DIM3
- = 	forAll (arbitrarySmallShape 100)		$ \(sh :: DIM2) ->
-	forAll (arbitraryListOfLength (S.size sh))	$ \(xx :: [Int]) -> 
+ = 	forAll (arbitrarySmallShape 20)		        $ \(sh :: DIM3) ->
+	forAll (arbitraryListOfLength (S.size sh))	$ \(xx :: [Int]) ->
 	sumAll (fromList sh xx) == P.sum xx
 
 
@@ -96,3 +107,4 @@
  = do	y	<- liftM (`mod` yMax) $ arbitrary
 	x	<- liftM (`mod` xMax) $ arbitrary
 	return	$ Z :. y :. x
+
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
@@ -5,13 +5,13 @@
 	( Shape(..)
 	, inShape )
 where
-	
--- Shape ------------------------------------------------------------------------------------------	
+
+-- Shape ------------------------------------------------------------------------------------------
 -- | Class of types that can be used as array shapes and indices.
 class Eq sh => Shape sh where
 
 	-- | Get the number of dimensions in a shape.
-	rank	:: sh -> Int           
+	rank	:: sh -> Int
 
 	-- | The shape of an array of size zero, with a particular dimensionality.
 	zeroDim	:: sh
@@ -26,7 +26,7 @@
 	addDim  :: sh -> sh -> sh
 
 	-- | Get the total number of elements in an array with this shape.
-	size	:: sh -> Int           
+	size	:: sh -> Int
 
 	-- | Check whether this shape is small enough so that its flat
 	--	indices an be represented as `Int`. If this returns `False` then your
@@ -37,13 +37,13 @@
 	-- | Convert an index into its equivalent flat, linear, row-major version.
 	toIndex :: sh	-- ^ Shape of the array.
 		-> sh 	-- ^ Index into the array.
-		-> Int     
+		-> Int
 
 	-- | Inverse of `toIndex`.
-	fromIndex 
+	fromIndex
 		:: sh 	-- ^ Shape of the array.
 		-> Int 	-- ^ Index into linear representation.
-		-> sh   
+		-> sh
 
 	-- | Check whether an index is within a given shape.
 	inShapeRange
@@ -54,7 +54,7 @@
 
 	-- | Convert a shape into its list of dimensions.
 	listOfShape	:: sh -> [Int]
-	
+
 	-- | Convert a list of dimensions to a shape
 	shapeOfList	:: [Int] -> sh
 
@@ -65,7 +65,7 @@
 
 -- | Check whether an index is a part of a given shape.
 inShape :: forall sh
-	.  Shape sh 
+	.  Shape sh
 	=> sh 		-- ^ Shape of the array.
 	-> sh		-- ^ Index.
 	-> Bool
@@ -73,3 +73,4 @@
 {-# INLINE inShape #-}
 inShape sh ix
 	= inShapeRange zeroDim sh ix
+
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
@@ -44,35 +44,35 @@
 
 	-- | Map an index of a slice onto an index of the full shape.
 	fullOfSlice	:: ss -> SliceShape ss -> FullShape  ss
-		
 
+
 instance Slice Z  where
 	{-# INLINE sliceOfFull #-}
 	sliceOfFull _ _		= Z
 
 	{-# INLINE fullOfSlice #-}
 	fullOfSlice _ _		= Z
-	
-	
+
+
 instance Slice (Any sh) where
 	{-# INLINE sliceOfFull #-}
 	sliceOfFull _ sh	= sh
 
 	{-# INLINE fullOfSlice #-}
 	fullOfSlice _ sh	= sh
-	
 
+
 instance Slice sl => Slice (sl :. Int) where
 	{-# INLINE sliceOfFull #-}
-	sliceOfFull (fsl :. _) (ssl :. _)	
+	sliceOfFull (fsl :. _) (ssl :. _)
 		= sliceOfFull fsl ssl
 
 	{-# INLINE fullOfSlice #-}
-	fullOfSlice (fsl :. n) ssl		
+	fullOfSlice (fsl :. n) ssl
 		= fullOfSlice fsl ssl :. n
-	
-	
-instance Slice sl => Slice (sl :. All) where	
+
+
+instance Slice sl => Slice (sl :. All) where
 	{-# INLINE sliceOfFull #-}
 	sliceOfFull (fsl :. All) (ssl :. s)
 		= sliceOfFull fsl ssl :. s
@@ -80,3 +80,4 @@
 	{-# INLINE 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
@@ -11,12 +11,12 @@
 
 
 -- | Check if an index lies inside the given extent.
---   As opposed to `inRange` from "Data.Array.Repa.Index", 
+--   As opposed to `inRange` from "Data.Array.Repa.Index",
 --   this is a short-circuited test that checks that lowest dimension first.
-isInside2 
+isInside2
 	:: DIM2 	-- ^ Extent of array.
 	-> DIM2 	-- ^ Index to check.
-	-> Bool	
+	-> Bool
 
 {-# INLINE isInside2 #-}
 isInside2 ex 	= not . isOutside2 ex
@@ -25,13 +25,13 @@
 -- | Check if an index lies outside the given extent.
 --   As opposed to `inRange` from "Data.Array.Repa.Index",
 --   this is a short-circuited test that checks the lowest dimension first.
-isOutside2 
-	:: DIM2		-- ^ Extent of array. 
+isOutside2
+	:: DIM2		-- ^ Extent of array.
 	-> DIM2		-- ^ Index to check.
 	-> Bool
-	
+
 {-# INLINE isOutside2 #-}
-isOutside2 (_ :. yLen :. xLen) (_ :. yy :. xx) 
+isOutside2 (_ :. yLen :. xLen) (_ :. yy :. xx)
 	| xx < 0	= True
 	| xx >= xLen	= True
 	| yy < 0	= True
@@ -42,7 +42,7 @@
 -- | Given the extent of an array, clamp the components of an index so they
 --   lie within the given array. Outlying indices are clamped to the index
 --   of the nearest border element.
-clampToBorder2 
+clampToBorder2
 	:: DIM2 	-- ^ Extent of array.
 	-> DIM2		-- ^ Index to clamp.
 	-> DIM2
@@ -55,7 +55,7 @@
 	  | x < 0	= clampY y 0
 	  | x >= xLen	= clampY y (xLen - 1)
 	  | otherwise	= clampY y x
-		
+
 	{-# INLINE clampY #-}
 	clampY !y !x
 	  | y < 0	= sh :. 0	   :. x
@@ -63,7 +63,7 @@
 	  | otherwise	= sh :. y	   :. x
 
 
--- | Make a 2D partitioned array given two generators, one to produce elements in the 
+-- | 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.
 --   The border must be the same width on all sides.
 makeBordered2
@@ -85,7 +85,7 @@
 
 	-- | 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 :. 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
@@ -94,16 +94,15 @@
 	inBorder 	= not . inInternal
 
 	-- Range of values where we don't need to worry about the border
-	rectsInternal	
+	rectsInternal
 	 = 	[ Rect (Z :. yMin :. xMin)	   (Z :. yMax :. xMax ) ]
 
 	{-# INLINE inInternal #-}
 	inInternal (Z :. y :. x)
-		=  x >= xMin && x <= xMax 
+		=  x >= xMin && x <= xMax
 		&& y >= yMin && y <= yMax
 
    in	Array sh
 		[ Region (RangeRects inBorder   rectsBorder)    genInternal
 		, Region (RangeRects inInternal rectsInternal)  genBorder ]
-
 
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
@@ -1,4 +1,4 @@
-{-# LANGUAGE 	MagicHash, PatternGuards, BangPatterns, TemplateHaskell, QuasiQuotes, 
+{-# LANGUAGE 	MagicHash, PatternGuards, BangPatterns, TemplateHaskell, QuasiQuotes,
 		ParallelListComp, TypeOperators, ExplicitForAll, ScopedTypeVariables #-}
 {-# OPTIONS -Wnot #-}
 
@@ -22,7 +22,7 @@
 	-- * Stencil operators.
 	, mapStencil2,     forStencil2
 	, mapStencilFrom2, forStencilFrom2
-	
+
 	--  From Data.Array.Repa.Stencil.Template
 	, stencil2)
 where
@@ -40,7 +40,7 @@
 
 -- | A index into the flat array.
 --   Should be abstract outside the stencil modules.
-data Cursor 
+data Cursor
 	= Cursor Int
 
 
@@ -89,7 +89,7 @@
 ---------------------------------------------------------------------------------------------------
 -- | Apply a stencil to every element of a 2D array.
 --   The array must be manifest else `error`.
-mapStencilFrom2 
+mapStencilFrom2
 	:: (Elt a, Elt b)
 	=> Boundary a		-- ^ How to handle the boundary of the array.
 	-> Stencil DIM2 a	-- ^ Stencil to apply.
@@ -114,18 +114,18 @@
 
 	-- Rectangles -----------------------
 	-- range of values where we don't need to worry about the border
-	rectsInternal	
+	rectsInternal
 	 = 	[ Rect (Z :. yMin :. xMin)	   (Z :. yMax :. xMax ) ]
 
 	{-# INLINE inInternal #-}
 	inInternal (Z :. y :. x)
-		=  x >= xMin && x <= xMax 
+		=  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 :. 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
@@ -136,31 +136,31 @@
 
 	-- Cursor functions ----------------
 	{-# INLINE makeCursor' #-}
-	makeCursor' (Z :. y :. x)	
+	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	
+	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') ]
 
@@ -171,7 +171,7 @@
 	-> Stencil DIM2 a
 	-> Array DIM2 b
 	-> (b -> a)
-	-> Cursor 
+	-> Cursor
 	-> a
 
 {-# INLINE [1] unsafeAppStencilCursor2 #-}
@@ -183,17 +183,17 @@
 	| _ :. sHeight :. sWidth	<- sExtent
 	, _ :. aHeight :. aWidth	<- aExtent
 	, sHeight <= 7, sWidth <= 7
-	= let	
+	= 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	
+		oload oy ox
 		 = let	!cur' = shift (Z :. oy :. ox) cur
 		   in	load (Z :. oy :. ox) (getData cur')
-	
+
 	   in	template7x7 oload zero
 
 
@@ -202,7 +202,7 @@
 	:: forall a b. (Elt a, Elt b)
 	=> (DIM2 -> DIM2 -> DIM2)
 	-> Stencil DIM2 a
-	-> Array DIM2 b 
+	-> Array DIM2 b
 	-> (b -> a)
 	-> DIM2
 	-> a
@@ -216,7 +216,7 @@
 	| _ :. sHeight :. sWidth	<- sExtent
 	, _ :. aHeight :. aWidth	<- aExtent
 	, sHeight <= 7, sWidth <= 7
-	= let	
+	= let
 		-- Get data from the manifest array.
 		{-# INLINE [0] getData #-}
 		getData :: DIM2 -> a
@@ -229,25 +229,25 @@
 		 | 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	
+		oload oy ox
 		 = let	!cur' = shift (Z :. oy :. ox) cur
 		   in	load (Z :. oy :. ox) (getData cur')
-	
+
 	   in	template7x7 oload zero
 
 
@@ -267,5 +267,4 @@
 	$ 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
@@ -11,7 +11,7 @@
 -- | How to handle the case when the stencil lies partly outside the array.
 data Boundary a
 	-- | Treat points outside as having a constant value.
-	= BoundConst a	
+	= BoundConst a
 
 	-- | Clamp points outside to the same value as the edge pixel.
 	| BoundClamp
@@ -28,18 +28,18 @@
 	{ stencilExtent	:: !sh
 	, stencilZero	:: !a
 	, stencilAcc	:: !(sh -> a -> a -> a) }
-	
-	
+
+
 -- | Make a stencil from a function yielding coefficients at each index.
 makeStencil
-	:: (Elt a, Num a) 
+	:: (Elt a, Num a)
 	=> sh			-- ^ Extent of stencil.
 	-> (sh -> Maybe a) 	-- ^ Get the coefficient at this index.
 	-> Stencil sh a
 
 {-# INLINE makeStencil #-}
 makeStencil ex getCoeff
- = StencilStatic ex 0 
+ = StencilStatic ex 0
  $ \ix val acc
 	-> case getCoeff ix of
 		Nothing		-> acc
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
@@ -1,6 +1,6 @@
 {-# LANGUAGE TemplateHaskell, QuasiQuotes, ParallelListComp #-}
 
--- | Template 
+-- | Template
 module Data.Array.Repa.Stencil.Template
 	(stencil2)
 where
@@ -10,9 +10,9 @@
 import qualified Data.List	as List
 
 -- | QuasiQuoter for producing a static stencil defintion.
---   
---   A definition like 
---  
+--
+--   A definition like
+--
 --   @
 --     [stencil2|  0 1 0
 --                 1 0 1
@@ -20,7 +20,7 @@
 --   @
 --
 --   Is converted to:
---   
+--
 --   @
 --     makeStencil2 (Z:.3:.3)
 --        (\\ix -> case ix of
@@ -32,7 +32,7 @@
 --   @
 --
 stencil2 :: QuasiQuoter
-stencil2 = QuasiQuoter 
+stencil2 = QuasiQuoter
 		{ quoteExp	= parseStencil2
 		, quotePat	= undefined
 		, quoteType	= undefined
@@ -42,15 +42,15 @@
 -- | Parse a stencil definition.
 --   TODO: make this more robust.
 parseStencil2 :: String -> Q Exp
-parseStencil2 str 
- = let	
+parseStencil2 str
+ = let
 	-- Determine the extent of the stencil based on the layout.
 	-- TODO: make this more robust. In particular, handle blank
 	--       lines at the start of the definition.
 	line1 : _	= lines str
 	sizeX		= fromIntegral $ length $ lines str
 	sizeY		= fromIntegral $ length $ words line1
-	
+
 	-- TODO: this probably doesn't work for stencils who's extents are even.
 	minX		= negate (sizeX `div` 2)
 	minY		= negate (sizeY `div` 2)
@@ -59,7 +59,7 @@
 
 	-- List of coefficients for the stencil.
 	coeffs		= (List.map read $ words str) :: [Integer]
-	
+
    in	makeStencil2' sizeX sizeY
 	 $ filter (\(_, _, v) -> v /= 0)
 	 $ [ (fromIntegral y, fromIntegral x, fromIntegral v)
@@ -80,20 +80,19 @@
 	ix'		<- newName "ix"
 	z'		<- [p| Z |]
 	coeffs'		<- newName "coeffs"
-	
-	let fnCoeffs	
+
+	let fnCoeffs
 		= LamE  [VarP ix']
-	 	$ CaseE (VarE ix') 
+	 	$ CaseE (VarE ix')
 	 	$   [ Match	(InfixP (InfixP z' dot' (LitP (IntegerL oy))) dot' (LitP (IntegerL ox)))
 				(NormalB $ ConE just' `AppE` LitE (IntegerL v))
 				[] | (oy, ox, v) <- coeffs ]
-	  	    ++ [Match WildP 
+	  	    ++ [Match WildP
 				(NormalB $ ConE (mkName "Nothing")) []]
-	
-	return 
+
+	return
 	 $ AppE (VarE makeStencil' `AppE` (LitE (IntegerL sizeX)) `AppE` (LitE (IntegerL sizeY)))
          $ LetE [ PragmaD (InlineP coeffs' (InlineSpec True False Nothing))
 		, ValD 	(VarP coeffs') (NormalB fnCoeffs) [] ]
 		(VarE coeffs')
-			
 
diff --git a/repa.cabal b/repa.cabal
--- a/repa.cabal
+++ b/repa.cabal
@@ -1,5 +1,5 @@
 Name:                repa
-Version:             2.0.2.1
+Version:             2.1.0.1
 License:             BSD3
 License-file:        LICENSE
 Author:              The DPH Team
@@ -50,6 +50,7 @@
         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
@@ -58,6 +59,7 @@
         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.Stencil.Base