diff --git a/Data/Array/Repa.hs b/Data/Array/Repa.hs
--- a/Data/Array/Repa.hs
+++ b/Data/Array/Repa.hs
@@ -29,11 +29,8 @@
 --  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:/
+--  /Advice for writing fast code:/
 --
 --  1. Repa does not support nested parallellism. 
 --     This means that you cannot `map` a parallel worker function across
@@ -43,29 +40,59 @@
 --
 --  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`
+--     You need to use `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.
+--  3. Add @INLINE@ pragmas to all leaf-functions in your code, expecially ones
+--     that compute numeric results. Non-inlined lazy function calls can cost
+--     upwards of 50 cycles each, while each numeric operator only costs one (or less).
+--     Inlining leaf functions also ensures they are specialised at the appropriate
+--     numeric types.
+--     
+--  4. Add bang patterns to all function arguments, and all fields of your data
+--     types. In a high-performance Haskell program, the cost of lazy evaluation
+--     can easily dominate the run time if not handled correctly. You don't want
+--     to rely on the strictness analyser in numeric code because if it does not
+--     return a perfect result then the performance of your program will be awful.
+--     This is less of a problem for general Haskell code, and in a different
+--     context relying on strictness analysis is fine.
 --
---  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.
+--  5. Scheduling a parallel computation takes about 200us on an OSX machine. 
+--     You should switch to sequential evaluation functions like `computeS` and
+--     `foldS` for small arrays in inner loops, and at the bottom of a 
+--     divide-and-conquer algorithm. Consider using a `computeP` that evaluates
+--     an array defined using `computeS` or `foldS` for each element.
 --
+--  6. Compile the modules that use Repa with the following flags:
+--     @-Odph -rtsopts -threaded@
+--     @-fno-liberate-case -funfolding-use-threshold1000 -funfolding-keeness-factor1000@
+--     @-fllvm -optlo-O3@
+--     You don't want the liberate-case transform because it tends to duplicate
+--     too much intermediate code, and is not needed if you use bang patterns
+--     as per point 4. The unfolding flags tell the inliner to not to fool around with 
+--     heuristics, and just inline everything. If the binaries become too big then 
+--     split the array part of your program into separate modules and only compile
+--     those with the unfolding flags.
+--
+--  7. Repa writes to the GHC eventlog at the start and end of  each parallel computation.
+--     Use threadscope to see what your program is doing.
+--
+--  8. Follow the advice on program structure in the comment for `deepSeqArrays`
+--
+--  9. When you're sure your program works, switch to the unsafe versions
+--     of functions like `traverse`. These don't do bounds checks.
+-- 
 module Data.Array.Repa
         ( -- * Abstract array representation
-          Array(..)
-        , module Data.Array.Repa.Shape
+          module Data.Array.Repa.Shape
         , module Data.Array.Repa.Index
+        , Array(..)
         , Repr(..), (!), toList
         , deepSeqArrays
 
-        -- * Converting between array representations
+        -- * Computation
         , computeP, computeS
         , copyP,    copyS
-        , now
 
         -- * Concrete array representations
         -- ** Delayed representation
@@ -84,13 +111,15 @@
 	-- ** Index space transformations
 	, reshape
 	, append, (++)
+        , extract
 	, transpose
-	, extend
-	, backpermute,         unsafeBackpermute
+	, backpermute
 	, backpermuteDft
 
+        -- ** Slice transformations
 	, module Data.Array.Repa.Slice
 	, slice
+        , extend
 
 	-- from Data.Array.Repa.Operators.Mapping -------------------
         -- ** Structure preserving operations
@@ -101,10 +130,10 @@
 
 	-- from Data.Array.Repa.Operators.Traversal ------------------
 	-- ** Generic traversal
-	, traverse,            unsafeTraverse
-	, traverse2,           unsafeTraverse2
-	, traverse3,           unsafeTraverse3
-	, traverse4,           unsafeTraverse4
+	, traverse 
+	, traverse2
+	, traverse3
+	, traverse4
 	
 	-- from Data.Array.Repa.Operators.Interleave -----------------
 	-- ** Interleaving
@@ -118,9 +147,11 @@
 	, foldAllP, foldAllS
 	, sumP,     sumS
 	, sumAllP,  sumAllS
+        , equalsP,  equalsS
 	
 	-- from Data.Array.Repa.Operators.Selection ------------------
-	, select)
+        -- ** Selection
+	, selectP)
 where
 import Data.Array.Repa.Base
 import Data.Array.Repa.Shape
diff --git a/Data/Array/Repa/Base.hs b/Data/Array/Repa/Base.hs
--- a/Data/Array/Repa/Base.hs
+++ b/Data/Array/Repa/Base.hs
@@ -12,10 +12,11 @@
 data family Array r sh e
 
 
+-- Repr -----------------------------------------------------------------------
 -- | Class of array representations that we can read elements from.
 --
 class Repr r e where
- -- | O(1). Take the extent of an array.
+ -- | O(1). Take the extent (size) of an array.
  extent       :: Shape sh => Array r sh e -> sh
 
  -- | O(1). Shape polymorphic indexing.
@@ -61,6 +62,32 @@
 --   The implementation of this function has been hand-unwound to work for up to
 --   four arrays. Putting more in the list yields `error`.
 -- 
+--   For functions that are /not/ marked as INLINE, you should apply `deepSeqArrays`
+--   to argument arrays before using them in a @compute@ or @copy@ expression.
+--   For example:
+--
+-- @  processArrays 
+--     :: Monad m 
+--     => Array U DIM2 Int -> Array U DIM2 Int 
+--     -> m (Array U DIM2 Int)
+--  processArrays arr1 arr2
+--   = [arr1, arr2] \`deepSeqArrays\` 
+--     do arr3 <- computeP $ map f arr1
+--        arr4 <- computeP $ zipWith g arr3 arr2
+--        return arr4
+--  @
+--
+--  Applying `deepSeqArrays` tells the GHC simplifier that it's ok to unbox 
+--  size fields and the pointers to the underlying array data at the start
+--  of the function. Without this, they may be unboxed repeatedly when
+--  computing elements in the result arrays, which will make your program slow.
+--
+--  If you INLINE @processArrays@ into the function that computes @arr1@ and @arr2@,
+--  then you don't need to apply `deepSeqArrays`. This is because a pointer
+--  to the underlying data will be passed directly to the consumers and never boxed.
+--
+--  If you're not sure, then just follow the example code above.
+--   
 deepSeqArrays 
         :: (Shape sh, Repr r e)
         => [Array r sh e] -> b -> b
@@ -82,4 +109,6 @@
          -> 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
--- a/Data/Array/Repa/Eval.hs
+++ b/Data/Array/Repa/Eval.hs
@@ -12,8 +12,8 @@
         , fromList
         
         -- * Converting between representations
-        , computeP, computeS
-        , copyP,    copyS
+        , computeS, computeP, suspendedComputeP
+        , copyS,    copyP,    suspendedCopyP
         , now
         
         -- * Chunked filling
@@ -46,8 +46,8 @@
 
 -- | Parallel computation of array elements.
 --
---   * The `Fill` class is defined so that the source array must have a
---     delayed representation (`D` or `C`)
+--   * The source array must have a delayed representation like `D`, `C` or `P`, 
+--     and the result a manifest representation like `U` or `F`.
 --
 --   * If you want to copy data between manifest representations then use
 --    `copyP` instead.
@@ -55,30 +55,47 @@
 --   * 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
+computeP 
+        :: (Shape sh, Fill r1 r2 sh e, Repr r2 e, Monad m)
+        => Array r1 sh e -> m (Array r2 sh e)
+computeP arr = now $ suspendedComputeP arr
 {-# 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
+{-# INLINE [4] computeS #-}
 
 
+-- | Suspended parallel computation of array elements.
+--
+--   This version creates a thunk that will evaluate the array on demand.
+--   If you force it when another parallel computation is already running
+--   then you  will get a runtime warning and evaluation will be sequential. 
+--   Use `deepSeqArray` and `now` to ensure that each array is evaluated
+--   before proceeding to the next one. 
+--  
+--   If unsure then just use the monadic version `computeP`. This one ensures
+--   that each array is fully evaluated before continuing.
+--
+suspendedComputeP 
+        :: Fill r1 r2 sh e
+        => Array r1 sh e -> Array r2 sh e
+suspendedComputeP arr1
+ = arr1 `deepSeqArray` 
+   unsafePerformIO
+ $ do   marr2    <- newMArr (size $ extent arr1) 
+        fillP arr1 marr2
+        unsafeFreezeMArr (extent arr1) marr2
+{-# INLINE [4] suspendedComputeP #-}
 
 
 -- | Parallel copying of arrays.
@@ -87,44 +104,41 @@
 -- 
 --   * 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
+copyP  :: (Shape sh, Fill D r2 sh e, Repr r1 e, Repr r2 e, Monad m)
+        => Array r1 sh e -> m (Array r2 sh e)
+copyP arr = now $ suspendedCopyP arr
 {-# 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
+copyS arr1  = computeS $ delay arr1
 {-# INLINE [4] copyS #-}
-copyS arr1 = computeS $ delay arr1
 
 
-        
+-- | Suspended parallel copy of array elements.
+suspendedCopyP   
+        :: (Repr r1 e, Fill D r2 sh e)
+        => Array r1 sh e -> Array r2 sh e
+suspendedCopyP arr1  = suspendedComputeP $ delay arr1
+{-# INLINE [4] suspendedCopyP #-}
 
--- | 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:
+
+-- | Monadic version of `deepSeqArray`. 
+-- 
+--   Forces an suspended array computation to be completed at this point
+--   in a monadic computation.
 --
---   @ do  arr2 <- now $ computeP $ map f arr1
---     arr3 <- now $ computeP $ zipWith arr2 arr1
---     return arr3
---   @
+-- @ do  let arr2 = suspendedComputeP arr1
+--     ...
+--     arr3 <- now $ arr2
+--     ...
+-- @
 --
 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
+{-# INLINE [4] now #-}
diff --git a/Data/Array/Repa/Eval/Chunked.hs b/Data/Array/Repa/Eval/Chunked.hs
--- a/Data/Array/Repa/Eval/Chunked.hs
+++ b/Data/Array/Repa/Eval/Chunked.hs
@@ -4,7 +4,6 @@
 module Data.Array.Repa.Eval.Chunked
 	( fillChunkedP
 	, fillChunkedS
-        , fillChunkedS'
 	, fillChunkedIOP)
 where
 import Data.Array.Repa.Eval.Gang
@@ -29,21 +28,6 @@
 	 | 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.
diff --git a/Data/Array/Repa/Eval/Cursored.hs b/Data/Array/Repa/Eval/Cursored.hs
--- a/Data/Array/Repa/Eval/Cursored.hs
+++ b/Data/Array/Repa/Eval/Cursored.hs
@@ -11,9 +11,8 @@
 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
+import GHC.Base
+import Prelude                          as P
 
 -- Non-cursored interface -----------------------------------------------------
 -- | Fill a block in a rank-2 array in parallel.
@@ -31,11 +30,11 @@
         :: 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
+	-> 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 #-}
@@ -92,33 +91,33 @@
 	-> (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
+	-> 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
+	!imageWidth !x0 !y0 !x1 !y1
  = 	gangIO theGang fillBlock
- where	!threads	= gangSize theGang
-	!blockWidth	= x1 - x0 + 1
+ where	!(I# threads)  = gangSize theGang
+	!blockWidth    = x1 -# x0 +# 1#
 
 	-- All columns have at least this many pixels.
-	!colChunkLen	= blockWidth `quotInt` threads
+	!colChunkLen   = I# (blockWidth `quotInt#` threads)
 
 	-- Extra pixels that we have to divide between some of the threads.
-	!colChunkSlack	= blockWidth `remInt` threads
+	!colChunkSlack = I# (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
+	 | ix < colChunkSlack = (I# x0) +  ix * (colChunkLen  + 1)
+	 | otherwise	      = (I# x0) + (ix * colChunkLen) + colChunkSlack
 
 	-- Give one column to each thread
 	{-# INLINE fillBlock #-}
@@ -126,8 +125,8 @@
 	fillBlock !ix
 	 = let	!(I# x0')	= colIx ix
 		!(I# x1')	= colIx (ix + 1) - 1
-		!(I# y0')	= y0
-		!(I# y1')	= y1
+		!y0'	  = y0
+		!y1'	  = y1
 	   in	fillCursoredBlock2S
 			write
 			makeCursorFCB shiftCursorFCB getElemFCB
diff --git a/Data/Array/Repa/Eval/Fill.hs b/Data/Array/Repa/Eval/Fill.hs
--- a/Data/Array/Repa/Eval/Fill.hs
+++ b/Data/Array/Repa/Eval/Fill.hs
@@ -26,6 +26,9 @@
  -- | Freeze the mutable array into an immutable Repa array.
  unsafeFreezeMArr :: sh  -> MArr r e -> IO (Array r sh e)
 
+ -- | Ensure the strucure of a mutable array is fully evaluated.
+ deepSeqMArr      :: MArr r e -> a -> a
+
 
 -- | O(n). Construct a manifest array from a list.
 fromList
diff --git a/Data/Array/Repa/Eval/Selection.hs b/Data/Array/Repa/Eval/Selection.hs
--- a/Data/Array/Repa/Eval/Selection.hs
+++ b/Data/Array/Repa/Eval/Selection.hs
@@ -96,31 +96,36 @@
 	-- 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'
+         | ixSrc > ixSrcEnd
+         = do  vecDst   <- VM.new 0
+               writeIORef ref vecDst
 
+         | otherwise
+	 = do  vecDst	<- VM.new (len `div` threads)
+               vecDst'	<- fillChunk ixSrc ixSrcEnd vecDst 0 (VM.length vecDst)
+	       writeIORef ref vecDst'
 
+
 	-- The main filling loop.
 	fillChunk :: Int -> Int -> IOVector a -> Int -> Int -> IO (IOVector a)
-	fillChunk !ixSrc !ixSrcEnd !vecDst !ixDst !ixDstEnd
+	fillChunk !ixSrc !ixSrcEnd !vecDst !ixDst !ixDstLen
          -- 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'
+	 | ixDst >= ixDstLen
+	 = do	let ixDstLen'	= (VM.length vecDst + 1) * 2
+		vecDst' 	<- VM.grow vecDst ixDstLen'
+		fillChunk ixSrc ixSrcEnd vecDst' ixDst ixDstLen'
 
 	 -- 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
+		fillChunk (ixSrc + 1) ixSrcEnd vecDst (ixDst + 1) ixDstLen
 
 	 -- The element doesnt match, so keep going.
 	 | otherwise
-	 =	fillChunk (ixSrc + 1) ixSrcEnd vecDst ixDst ixDstEnd
+	 =	fillChunk (ixSrc + 1) ixSrcEnd vecDst ixDst ixDstLen
 
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
@@ -8,12 +8,8 @@
 	, (:.)	(..)
 
 	-- * Common dimensions.
-	, DIM0
-	, DIM1
-	, DIM2
-	, DIM3
-	, DIM4
-	, DIM5)
+	, DIM0, DIM1, DIM2, DIM3, DIM4, DIM5
+        ,       ix1,  ix2,  ix3,  ix4,  ix5)
 where
 import Data.Array.Repa.Shape
 import GHC.Base 		(quotInt, remInt)
@@ -39,7 +35,33 @@
 type DIM5	= DIM4 :. Int
 
 
--- Shape ------------------------------------------------------------------------------------------
+-- | Helper for index construction.
+--
+--   Use this instead of explicit constructors like @(Z :. (x :: Int))@.
+--   The this is sometimes needed to ensure that 'x' is constrained to 
+--   be in @Int@.
+ix1 :: Int -> DIM1
+ix1 x = Z :. x
+{-# INLINE ix1 #-}
+
+ix2 :: Int -> Int -> DIM2
+ix2 y x = Z :. y :. x
+{-# INLINE ix2 #-}
+
+ix3 :: Int -> Int -> Int -> DIM3
+ix3 z y x = Z :. z :. y :. x
+{-# INLINE ix3 #-}
+
+ix4 :: Int -> Int -> Int -> Int -> DIM4
+ix4 a z y x = Z :. a :. z :. y :. x
+{-# INLINE ix4 #-}
+
+ix5 :: Int -> Int -> Int -> Int -> Int -> DIM5
+ix5 b a z y x = Z :. b :. a :. z :. y :. x
+{-# INLINE ix5 #-}
+
+
+-- Shape ----------------------------------------------------------------------
 instance Shape Z where
 	{-# INLINE [1] rank #-}
 	rank _			= 0
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
@@ -4,10 +4,11 @@
 	( reshape
 	, append, (++)
 	, transpose
-	, extend
-	, slice
+        , extract
 	, backpermute,         unsafeBackpermute
-	, backpermuteDft,      unsafeBackpermuteDft)
+	, backpermuteDft,      unsafeBackpermuteDft
+        , extend,              unsafeExtend 
+        , slice,               unsafeSlice)
 where
 import Data.Array.Repa.Index
 import Data.Array.Repa.Slice
@@ -29,7 +30,6 @@
 	-> Array r1 sh1 e
 	-> Array D  sh2 e
 
-{-# INLINE [3] reshape #-}
 reshape sh2 arr
 	| not $ S.size sh2 == S.size (extent arr)
 	= error 
@@ -38,6 +38,7 @@
 reshape sh2 arr
         = fromFunction sh2 
         $ unsafeIndex arr . fromIndex (extent arr) . toIndex sh2
+{-# INLINE [2] reshape #-}
  
 
 -- | Append two arrays.
@@ -48,7 +49,6 @@
 	-> Array r2 (sh :. Int) e
 	-> Array D  (sh :. Int) e
 
-{-# INLINE [3] append #-}
 append arr1 arr2
  = unsafeTraverse2 arr1 arr2 fnExtent fnElem
  where
@@ -60,9 +60,11 @@
 	fnElem f1 f2 (sh :. i)
       		| i < n		= f1 (sh :. i)
   		| otherwise	= f2 (sh :. (i - n))
+{-# INLINE [2] append #-}
 
-{-# INLINE (++) #-}
+
 (++) arr1 arr2 = append arr1 arr2
+{-# INLINE (++) #-}
 
 
 -- | Transpose the lowest two dimensions of an array.
@@ -73,45 +75,22 @@
 	=> Array r (sh :. Int :. Int) e
 	-> Array D (sh :. Int :. Int) e
 
-{-# INLINE [3] transpose #-}
 transpose arr
  = unsafeTraverse arr
 	(\(sh :. m :. n) 	-> (sh :. n :.m))
 	(\f -> \(sh :. i :. j) 	-> f (sh :. j :. i))
+{-# INLINE [2] transpose #-}
 
 
--- | Extend an array, according to a given slice specification.
-extend
-	:: ( Slice sl
-	   , Shape (FullShape sl)
-	   , Shape (SliceShape sl)
-	   , Repr r e)
-	=> sl
-	-> Array r (SliceShape sl) e
-	-> Array D (FullShape sl)  e
-
-{-# INLINE [3] extend #-}
-extend sl arr
-	= unsafeBackpermute
-		(fullOfSlice sl (extent arr))
-		(sliceOfFull sl)
-		arr
-
--- | Take a slice from an array, according to a given specification.
-slice	:: ( Slice sl
-	   , Shape (FullShape sl)
-	   , Shape (SliceShape sl)
-	   , Repr r e)
-	=> Array r (FullShape sl) e
-	-> sl
-	-> Array D (SliceShape sl) e
-
-{-# INLINE [3] slice #-}
-slice arr sl
-	= unsafeBackpermute
-		(sliceOfFull sl (extent arr))
-		(fullOfSlice sl)
-		arr
+-- | Extract a sub-range of elements from an array.
+extract :: (Shape sh, Repr r e)
+        => sh                   -- ^ Starting index.
+        -> sh                   -- ^ Size of result.
+        -> Array r sh e 
+        -> Array D sh e
+extract start sz arr
+        = fromFunction sz (\ix -> arr `unsafeIndex` (addDim start ix))
+{-# INLINE [2] extract #-}
 
 
 -- | Backwards permutation of an array's elements.
@@ -126,13 +105,13 @@
 	-> Array r  sh1 e 	-- ^ Source array.
 	-> Array D  sh2 e
 
-{-# INLINE [3] backpermute #-}
 backpermute newExtent perm arr
 	= traverse arr (const newExtent) (. perm)
+{-# INLINE [2] backpermute #-}
 
-{-# INLINE [3] unsafeBackpermute #-}
 unsafeBackpermute newExtent perm arr
         = unsafeTraverse arr (const newExtent) (. perm)
+{-# INLINE [2] unsafeBackpermute #-}
 
 
 -- | Default backwards permutation of an array's elements.
@@ -148,19 +127,85 @@
 	-> Array r1 sh1 e	-- ^ Source array.
 	-> Array D  sh2 e
 
-{-# INLINE [3] backpermuteDft #-}
 backpermuteDft arrDft fnIndex arrSrc
 	= fromFunction (extent arrDft) fnElem
 	where	fnElem ix
 		 = case fnIndex ix of
 			Just ix'	-> arrSrc `index` ix'
 			Nothing		-> arrDft `index` ix
+{-# INLINE [2] backpermuteDft #-}
 
-{-# 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
+{-# INLINE [2] unsafeBackpermuteDft #-}
 
+
+
+-- | Extend an array, according to a given slice specification.
+--
+--   For example, to replicate the rows of an array use the following:
+--
+--   @extend arr (Any :. (5::Int) :. All)@
+--
+extend, unsafeExtend
+        :: ( Slice sl
+           , Shape (FullShape sl)
+           , Shape (SliceShape sl)
+           , Repr r e)
+        => sl
+        -> Array r (SliceShape sl) e
+        -> Array D (FullShape sl)  e
+
+extend sl arr
+        = backpermute
+                (fullOfSlice sl (extent arr))
+                (sliceOfFull sl)
+                arr
+{-# INLINE [2] extend #-}
+
+unsafeExtend sl arr
+        = unsafeBackpermute
+                (fullOfSlice sl (extent arr))
+                (sliceOfFull sl)
+                arr
+{-# INLINE [2] unsafeExtend #-}
+
+
+
+-- | Take a slice from an array, according to a given specification.
+--
+--   For example, to take a row from a matrix use the following:
+--
+--   @slice arr (Any :. (5::Int) :. All)@
+--
+--   To take a column use:
+--
+--   @slice arr (Any :. (5::Int))@
+--
+slice, unsafeSlice
+        :: ( Slice sl
+           , Shape (FullShape sl)
+           , Shape (SliceShape sl)
+           , Repr r e)
+        => Array r (FullShape sl) e
+        -> sl
+        -> Array D (SliceShape sl) e
+
+slice arr sl
+        = backpermute
+                (sliceOfFull sl (extent arr))
+                (fullOfSlice sl)
+                arr
+{-# INLINE [2] slice #-}
+
+
+unsafeSlice arr sl
+        = unsafeBackpermute
+                (sliceOfFull sl (extent arr))
+                (fullOfSlice sl)
+                arr
+{-# INLINE [2] unsafeSlice #-}
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
@@ -29,7 +29,7 @@
 	-> Array r2 (sh :. Int) a
 	-> Array D  (sh :. Int) a
 
-{-# INLINE [3] interleave2 #-}
+{-# INLINE [2] interleave2 #-}
 interleave2 arr1 arr2
  = arr1 `deepSeqArray` arr2 `deepSeqArray`
    unsafeTraverse2 arr1 arr2 shapeFn elemFn
@@ -58,7 +58,7 @@
 	-> Array r3 (sh :. Int) a
 	-> Array D  (sh :. Int) a
 
-{-# INLINE [3] interleave3 #-}
+{-# INLINE [2] interleave3 #-}
 interleave3 arr1 arr2 arr3
  = arr1 `deepSeqArray` arr2 `deepSeqArray` arr3 `deepSeqArray`
    unsafeTraverse3 arr1 arr2 arr3 shapeFn elemFn
@@ -90,7 +90,7 @@
 	-> Array r4 (sh :. Int) a
 	-> Array D  (sh :. Int) a
 
-{-# INLINE [3] interleave4 #-}
+{-# INLINE [2] interleave4 #-}
 interleave4 arr1 arr2 arr3 arr4
  = arr1 `deepSeqArray` arr2 `deepSeqArray` arr3 `deepSeqArray` arr4 `deepSeqArray`
    unsafeTraverse4 arr1 arr2 arr3 arr4 shapeFn elemFn
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
@@ -27,7 +27,7 @@
 --
 map     :: (Shape sh, Repr r a)
         => (a -> b) -> Array r sh a -> Array D sh b
-{-# INLINE [4] map #-}
+{-# INLINE [3] map #-}
 map f arr
  = case delay arr of
         ADelayed sh g -> ADelayed sh (f . g)
@@ -42,7 +42,7 @@
         => (a -> b -> c)
         -> Array r1 sh a -> Array r2 sh b
         -> Array D sh c
-{-# INLINE [3] zipWith #-}
+{-# INLINE [2] zipWith #-}
 zipWith f arr1 arr2
  = arr1 `deepSeqArray` arr2 `deepSeqArray`
    let 
@@ -112,11 +112,11 @@
 
 -- Cursored ---------------------------
 instance Combine C a C b where
- {-# INLINE [4] cmap #-}
+ {-# INLINE [3] cmap #-}
  cmap f (ACursored sh makec shiftc loadc)
         = ACursored sh makec shiftc (f . loadc)
 
- {-# INLINE [3] czipWith #-}
+ {-# INLINE [2] czipWith #-}
  czipWith f arr1 (ACursored sh makec shiftc loadc)
   = let {-# INLINE makec' #-}
         makec' ix               = (ix, makec ix)
@@ -149,11 +149,11 @@
         , Combine r12 a r22 b)
        => Combine (P r11 r12) a (P r21 r22) b where
 
- {-# INLINE [4] cmap #-}
+ {-# INLINE [3] cmap #-}
  cmap f (APart sh range arr1 arr2)
         = APart sh range (cmap f arr1) (cmap f arr2)
 
- {-# INLINE [3] czipWith #-}
+ {-# INLINE [2] czipWith #-}
  czipWith f arr1 (APart sh range arr21 arr22)
         = APart sh range (czipWith f arr1 arr21)
                          (czipWith f arr1 arr22)
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,15 +1,17 @@
 {-# LANGUAGE BangPatterns, ExplicitForAll, TypeOperators, MagicHash #-}
-
+{-# OPTIONS -fno-warn-orphans #-}
 module Data.Array.Repa.Operators.Reduction
 	( foldS,        foldP
 	, foldAllS,     foldAllP
 	, sumS,         sumP
-	, sumAllS,      sumAllP)
+	, sumAllS,      sumAllP
+        , equalsS,      equalsP)
 where
 import Data.Array.Repa.Base
 import Data.Array.Repa.Index
-import Data.Array.Repa.Eval.Elt
+import Data.Array.Repa.Eval
 import Data.Array.Repa.Repr.Unboxed
+import Data.Array.Repa.Operators.Mapping        as R
 import Data.Array.Repa.Shape		        as S
 import qualified Data.Vector.Unboxed	        as V
 import qualified Data.Vector.Unboxed.Mutable    as M
@@ -18,24 +20,26 @@
 import System.IO.Unsafe
 import GHC.Exts
 
--- foldS ----------------------------------------------------------------------
+-- fold ----------------------------------------------------------------------
 -- | 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   :: (Shape sh, Elt a, Unbox a, Repr r a)
+        => (a -> a -> a)
+        -> a
+        -> Array r (sh :. Int) a
+        -> Array U sh a
+
 foldS f z arr
- = let  sh@(sz :. n') = extent arr
+ = arr `deepSeqArray`
+   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
+         now $ fromUnboxed sz vec
+{-# INLINE [1] foldS #-}
 
 
 -- | Parallel reduction of the innermost dimension of an arbitray rank array.
@@ -45,27 +49,32 @@
 --   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 r (sh :. Int) a
-	-> Array U sh a
-{-# INLINE [2] foldP #-}
+foldP   
+        :: (Shape sh, Elt a, Unbox a, Repr r a, Monad m)
+        => (a -> a -> a)
+        -> a
+        -> Array r (sh :. Int) a
+        -> m (Array U sh a)
+
 foldP f z arr 
- = let  sh@(sz :. n) = extent arr
+ = arr `deepSeqArray`
+   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 !vec = V.singleton $ foldAllP f z arr
-                in  fromUnboxed sz vec
+           1 -> do
+                x       <- foldAllP f z arr
+                now $ fromUnboxed sz $ V.singleton x
 
-           _ -> unsafePerformIO 
+           _ -> now
+              $ unsafePerformIO 
               $ do mvec   <- M.unsafeNew (S.size sz)
                    E.foldP mvec (\ix -> arr `unsafeIndex` fromIndex sh ix) f z n
                    !vec   <- V.unsafeFreeze mvec
-                   return $ fromUnboxed sz vec
+                   now $ fromUnboxed sz vec
+{-# INLINE [1] foldP #-}
 
 
 -- foldAll --------------------------------------------------------------------
@@ -76,7 +85,7 @@
 	-> a
 	-> Array r sh a
 	-> a
-{-# INLINE [2] foldAllS #-}
+
 foldAllS f z arr 
  = arr `deepSeqArray`
    let  !ex     = extent arr
@@ -84,6 +93,7 @@
    in   E.foldAllS 
                 (\ix -> arr `unsafeIndex` fromIndex ex (I# ix))
                 f z n 
+{-# INLINE [1] foldAllS #-}
 
 
 -- | Parallel reduction of an array of arbitrary rank to a single scalar value.
@@ -93,17 +103,21 @@
 --   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 
+        :: (Shape sh, Elt a, Unbox a, Repr r a, Monad m)
+	=> (a -> a -> a)
+	-> a
+	-> Array r sh a
+	-> m a
+
 foldAllP f z arr 
- = let  sh = extent arr
+ = arr `deepSeqArray`
+   let  sh = extent arr
         n  = size sh
-   in   unsafePerformIO 
+   in   return
+         $ unsafePerformIO 
          $ E.foldAllP (\ix -> arr `unsafeIndex` fromIndex sh ix) f z n
+{-# INLINE [1] foldAllP #-}
 
 
 -- sum ------------------------------------------------------------------------
@@ -111,16 +125,16 @@
 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
+{-# INLINE [3] sumS #-}
 
 
--- | Sequential sum the innermost dimension of an array.
-sumP	:: (Shape sh, Num a, Elt a, Unbox a, Repr r a)
+-- | Parallel sum the innermost dimension of an array.
+sumP	:: (Shape sh, Num a, Elt a, Unbox a, Repr r a, Monad m)
 	=> Array r (sh :. Int) a
-	-> Array U sh a
-{-# INLINE [4] sumP #-}
-sumP = foldP (+) 0
+	-> m (Array U sh a)
+sumP = foldP (+) 0 
+{-# INLINE [3] sumP #-}
 
 
 -- sumAll ---------------------------------------------------------------------
@@ -128,14 +142,43 @@
 sumAllS	:: (Shape sh, Elt a, Unbox a, Num a, Repr r a)
 	=> Array r sh a
 	-> a
-{-# INLINE [4] sumAllS #-}
 sumAllS = foldAllS (+) 0
+{-# INLINE [3] sumAllS #-}
 
 
 -- | Parallel sum all the elements of an array.
-sumAllP	:: (Shape sh, Elt a, Unbox a, Num a, Repr r a)
+sumAllP	:: (Shape sh, Elt a, Unbox a, Num a, Repr r a, Monad m)
 	=> Array r sh a
-	-> a
-{-# INLINE [4] sumAllP #-}
+	-> m a
 sumAllP = foldAllP (+) 0
+{-# INLINE [3] sumAllP #-}
+
+
+-- Equality ------------------------------------------------------------------
+instance (Shape sh, Repr r e, Eq e) => Eq (Array r sh e) where
+ (==) arr1 arr2
+        =   extent arr1 == extent arr2
+        && (foldAllS (&&) True (R.zipWith (==) arr1 arr2))
+
+
+-- | Check whether two arrays have the same shape and contain equal elements,
+--   in parallel.
+equalsP :: (Shape sh, Repr r1 e, Repr r2 e, Eq e, Monad m) 
+        => Array r1 sh e 
+        -> Array r2 sh e
+        -> m Bool
+equalsP arr1 arr2
+ = do   same    <- foldAllP (&&) True (R.zipWith (==) arr1 arr2)
+        return  $ (extent arr1 == extent arr2) && same
+
+
+-- | Check whether two arrays have the same shape and contain equal elements,
+--   sequentially.
+equalsS :: (Shape sh, Repr r1 e, Repr r2 e, Eq e) 
+        => Array r1 sh e 
+        -> Array r2 sh e
+        -> Bool
+equalsS arr1 arr2
+        =   extent arr1 == extent arr2
+        && (foldAllS (&&) True (R.zipWith (==) arr1 arr2))
 
diff --git a/Data/Array/Repa/Operators/Selection.hs b/Data/Array/Repa/Operators/Selection.hs
--- a/Data/Array/Repa/Operators/Selection.hs
+++ b/Data/Array/Repa/Operators/Selection.hs
@@ -1,6 +1,6 @@
 {-# LANGUAGE BangPatterns #-}
 module Data.Array.Repa.Operators.Selection
-	(select)
+	(selectP)
 where
 import Data.Array.Repa.Index
 import Data.Array.Repa.Base
@@ -19,15 +19,15 @@
 --
 --   * Use the integer as the index into the array you're filtering.
 --
-select	:: Unbox a
+selectP	:: (Unbox a, Monad m)
         => (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.
+	-> m (Array U DIM1 a)	-- ^ Array containing produced values.
 
-{-# INLINE [2] select #-}
-select match produce len
- = unsafePerformIO
+selectP match produce len
+ = return
+ $ unsafePerformIO
  $ do   (sh, vec)	<- selectIO
 	return $ sh `seq` vec `seq`
 	         fromUnboxed sh vec
@@ -41,3 +41,4 @@
 		let result	= V.concat vecs'
 
 		return	(Z :. V.length result, result)
+{-# INLINE [1] selectP #-}
diff --git a/Data/Array/Repa/Operators/Traversal.hs b/Data/Array/Repa/Operators/Traversal.hs
--- a/Data/Array/Repa/Operators/Traversal.hs
+++ b/Data/Array/Repa/Operators/Traversal.hs
@@ -20,15 +20,15 @@
 	 				--   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 [3] traverse #-}
 
-{-# INLINE [4] unsafeTraverse #-}
 unsafeTraverse arr transExtent newElem
  = arr `deepSeqArray`
    fromFunction (transExtent (extent arr)) (newElem (unsafeIndex arr))
+{-# INLINE [3] unsafeTraverse #-}
 
 
 -- | Unstructured traversal over two arrays at once.
@@ -45,17 +45,17 @@
 					--   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 [3] traverse2 #-}
 
-{-# INLINE [4] unsafeTraverse2 #-}
 unsafeTraverse2 arrA arrB transExtent newElem
  = arrA `deepSeqArray` arrB `deepSeqArray`
    fromFunction  (transExtent (extent arrA) (extent arrB))
                  (newElem     (unsafeIndex arrA) (unsafeIndex arrB))
+{-# INLINE [3] unsafeTraverse2 #-}
 
 
 -- | Unstructured traversal over three arrays at once.
@@ -74,17 +74,17 @@
            ->  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 [3] traverse3 #-}
 
-{-# 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))
+{-# INLINE [3] unsafeTraverse3 #-}
 
 
 -- | Unstructured traversal over four arrays at once.
@@ -104,17 +104,17 @@
            ->  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 [3] traverse4 #-}
 
 
-{-# 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))
+{-# INLINE [3] unsafeTraverse4 #-}
 
 
diff --git a/Data/Array/Repa/Repr/ByteString.hs b/Data/Array/Repa/Repr/ByteString.hs
--- a/Data/Array/Repa/Repr/ByteString.hs
+++ b/Data/Array/Repa/Repr/ByteString.hs
@@ -24,21 +24,21 @@
 -- Repr -----------------------------------------------------------------------
 -- | Read elements from a `ByteString`.
 instance Repr B Word8 where
- {-# INLINE linearIndex #-}
  linearIndex (AByteString _ bs) ix
         = bs `B.index` ix
+ {-# INLINE linearIndex #-}
 
- {-# INLINE unsafeLinearIndex #-}
  unsafeLinearIndex (AByteString _ bs) ix
         = bs `BU.unsafeIndex` ix
+ {-# INLINE unsafeLinearIndex #-}
 
- {-# INLINE extent #-}
  extent (AByteString sh _)
         = sh
+ {-# INLINE extent #-}
 
- {-# INLINE deepSeqArray #-}
  deepSeqArray (AByteString sh bs) x 
   = sh `deepSeq` bs `seq` x
+ {-# INLINE deepSeqArray #-}
 
 
 -- Conversions ----------------------------------------------------------------
@@ -46,12 +46,12 @@
 fromByteString
         :: Shape sh
         => sh -> ByteString -> Array B sh Word8
-{-# INLINE fromByteString #-}
 fromByteString sh bs
         = AByteString sh bs
+{-# INLINE fromByteString #-}
 
 
 -- | O(1). Unpack a `ByteString` from an array.
 toByteString :: Array B sh Word8 -> ByteString
-{-# INLINE toByteString #-}
 toByteString (AByteString _ bs) = bs
+{-# INLINE toByteString #-}
diff --git a/Data/Array/Repa/Repr/Cursored.hs b/Data/Array/Repa/Repr/Cursored.hs
--- a/Data/Array/Repa/Repr/Cursored.hs
+++ b/Data/Array/Repa/Repr/Cursored.hs
@@ -12,6 +12,7 @@
 import Data.Array.Repa.Eval.Elt
 import Data.Array.Repa.Eval.Cursored
 import GHC.Exts
+import Debug.Trace
 
 -- | Cursored Arrays.
 --   These are produced by Repa's stencil functions, and help the fusion
@@ -39,63 +40,72 @@
 -- Repr -----------------------------------------------------------------------
 -- | Compute elements of a cursored array.
 instance Repr C a where
- {-# INLINE index #-}
  index (ACursored _ makec _ loadc)
         = loadc . makec
+ {-# INLINE index #-}
 
- {-# INLINE unsafeIndex #-}
  unsafeIndex    = index
+ {-# INLINE unsafeIndex #-}
  
- {-# INLINE linearIndex #-}
  linearIndex (ACursored sh makec _ loadc)
         = loadc . makec . fromIndex sh
+ {-# INLINE linearIndex #-}
 
- {-# INLINE extent #-}
  extent (ACursored sh _ _ _)
         = sh
+ {-# INLINE extent #-}
         
- {-# INLINE deepSeqArray #-}
  deepSeqArray (ACursored sh makec shiftc loadc) y
   = sh `deepSeq` makec  `seq` shiftc `seq` loadc `seq` y
+ {-# INLINE deepSeqArray #-}
 
 
 -- 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 
+ fillP (ACursored (Z :. (I# h) :. (I# w)) makec shiftc loadc) marr
+  = do  traceEventIO "Repa.fillP[Cursored]: start"
+        fillCursoredBlock2P 
                 (unsafeWriteMArr marr) 
                 makec shiftc loadc
-                w 0 0 (w - 1) (h - 1) 
-
- {-# INLINE fillS #-}
+                w 0# 0# (w -# 1#) (h -# 1#) 
+        traceEventIO "Repa.fillP[Cursored]: end"
+ {-# INLINE fillP #-}
+        
  fillS (ACursored (Z :. (I# h) :. (I# w)) makec shiftc loadc) marr
-  = fillCursoredBlock2S 
+  = do  traceEventIO "Repa.fillS[Cursored]: start"
+        fillCursoredBlock2S 
                 (unsafeWriteMArr marr) 
                 makec shiftc loadc
                 w 0# 0# (w -# 1#) (h -# 1#) 
-
+        traceEventIO "Repa.fillS[Cursored]: end"
+ {-# INLINE fillS #-}
+        
 
 -- | 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 
+ fillRangeP  (ACursored (Z :. _h :. (I# w)) makec shiftc loadc) marr
+             (Z :. (I# y0) :. (I# x0)) (Z :. (I# y1) :. (I# x1))
+  = do  traceEventIO "Repa.fillRangeP[Cursored]: start"
+        fillCursoredBlock2P 
                 (unsafeWriteMArr marr) 
                 makec shiftc loadc
                 w x0 y0 x1 y1
-
- {-# INLINE fillRangeS #-}
+        traceEventIO "Repa.fillRangeP[Cursored]: end"
+ {-# INLINE fillRangeP #-}
+        
  fillRangeS  (ACursored (Z :. _h :. (I# w)) makec shiftc loadc) marr
              (Z :. (I# y0) :. (I# x0)) 
              (Z :. (I# y1) :. (I# x1))
-  = fillCursoredBlock2S
+  = do  traceEventIO "Repa.fillRangeS[Cursored]: start"
+        fillCursoredBlock2S
                 (unsafeWriteMArr marr) 
                 makec shiftc loadc
                 w x0 y0 x1 y1
- 
+        traceEventIO "Repa.fillRangeS[Cursored]: end"
+ {-# INLINE fillRangeS #-}
+        
+
 -- Conversions ----------------------------------------------------------------
 -- | Define a new cursored array.
 makeCursored 
@@ -105,5 +115,6 @@
         -> (cursor -> e)                -- ^ Compute the element at the cursor.
         -> Array C sh e
 
-{-# INLINE makeCursored #-}
 makeCursored = ACursored
+{-# INLINE makeCursored #-}
+
diff --git a/Data/Array/Repa/Repr/Delayed.hs b/Data/Array/Repa/Repr/Delayed.hs
--- a/Data/Array/Repa/Repr/Delayed.hs
+++ b/Data/Array/Repa/Repr/Delayed.hs
@@ -11,9 +11,13 @@
 import Data.Array.Repa.Index
 import Data.Array.Repa.Shape
 import Data.Array.Repa.Base
+import Debug.Trace
 import GHC.Exts
 
 -- | Delayed arrays are represented as functions from the index to element value.
+--
+--   Every time you index into a delayed array the element at that position 
+--   is recomputed.
 data D
 data instance Array D sh e
         = ADelayed  
@@ -24,66 +28,79 @@
 -- Repr -----------------------------------------------------------------------
 -- | Compute elements of a delayed array.
 instance Repr D a where
- {-# INLINE index #-}
  index       (ADelayed _  f) ix  = f ix
+ {-# INLINE index #-}
 
- {-# INLINE linearIndex #-}
  linearIndex (ADelayed sh f) ix  = f (fromIndex sh ix)
+ {-# INLINE linearIndex #-}
 
- {-# INLINE extent #-}
  extent (ADelayed sh _)
         = sh
+ {-# INLINE extent #-}
 
- {-# INLINE deepSeqArray #-}
  deepSeqArray (ADelayed sh f) y
         = sh `deepSeq` f `seq` y
+ {-# INLINE deepSeqArray #-}
 
 
 -- 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) 
+  = marr `deepSeqMArr` 
+    do  traceEventIO "Repa.fillP[Delayed]: start"
+        fillChunkedP (size sh) (unsafeWriteMArr marr) (getElem . fromIndex sh) 
+        traceEventIO "Repa.fillP[Delayed]: end"
+ {-# INLINE [4] fillP #-}
 
- {-# INLINE [4] fillS #-}
  fillS (ADelayed sh getElem) marr
-  = fillChunkedS (size sh) (unsafeWriteMArr marr) (getElem . fromIndex sh)
+  = marr `deepSeqMArr` 
+    do  traceEventIO "Repa.fillS[Delayed]: start"
+        fillChunkedS (size sh) (unsafeWriteMArr marr) (getElem . fromIndex sh)
+        traceEventIO "Repa.fillS[Delayed]: end"
+ {-# INLINE [4] fillS #-}
 
 
 -- | Compute a range of elements in a rank-2 array.
 instance (Fillable r2 e, Elt e) => FillRange D r2 DIM2 e where
+ fillRangeP  (ADelayed (Z :. _h :. (I# w)) getElem) marr
+             (Z :. (I# y0) :. (I# x0)) (Z :. (I# y1) :. (I# x1))
+  = marr `deepSeqMArr` 
+    do  traceEventIO "Repa.fillRangeP[Delayed]: start"
+        fillBlock2P (unsafeWriteMArr marr) 
+                        getElem
+                        w x0 y0 x1 y1
+        traceEventIO "Repa.fillRangeP[Delayed]: end"
  {-# 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) 
+  = marr `deepSeqMArr`
+    do  traceEventIO "Repa.fillRangeS[Delayed]: start"
+        fillBlock2S (unsafeWriteMArr marr) 
                 getElem
                 w x0 y0 x1 y1
+        traceEventIO "Repa.fillRangeS[Delayed]: end"
+ {-# INLINE [1] fillRangeS #-}
 
 
 -- 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 
+{-# INLINE fromFunction #-}
 
 
--- | O(1). Produce the extent of an array and a function to retrieve an arbitrary element.
+-- | 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)
+{-# INLINE toFunction #-}
 
 
 -- | O(1). Delay an array.
@@ -93,7 +110,7 @@
 --
 delay   :: (Shape sh, Repr r e)
         => Array r sh e -> Array D sh e
+delay arr = ADelayed (extent arr) (unsafeIndex arr)
 {-# 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
--- a/Data/Array/Repa/Repr/ForeignPtr.hs
+++ b/Data/Array/Repa/Repr/ForeignPtr.hs
@@ -12,6 +12,7 @@
 import Foreign.ForeignPtr
 import Foreign.Marshal.Alloc
 import System.IO.Unsafe
+import qualified Foreign.ForeignPtr.Unsafe      as Unsafe
 
 -- | Arrays represented as foreign buffers in the C heap.
 data F
@@ -21,7 +22,6 @@
 -- Repr -----------------------------------------------------------------------
 -- | Read elements from a foreign buffer.
 instance Storable a => Repr F a where
- {-# INLINE linearIndex #-}
  linearIndex (AForeignPtr _ len fptr) ix
   | ix < len  
         = unsafePerformIO 
@@ -30,21 +30,22 @@
   
   | otherwise
   = error "Repa: foreign array index out of bounds"
-
- {-# INLINE unsafeLinearIndex #-}
+ {-# INLINE linearIndex #-}
+ 
  unsafeLinearIndex (AForeignPtr _ _ fptr) ix
         = unsafePerformIO
         $ withForeignPtr fptr 
         $ \ptr -> peekElemOff ptr ix
+ {-# INLINE unsafeLinearIndex #-}
 
- {-# INLINE extent #-}
  extent (AForeignPtr sh _ _)
         = sh
+ {-# INLINE extent #-}
 
- {-# INLINE deepSeqArray #-}
  deepSeqArray (AForeignPtr sh len fptr) x 
   = sh `deepSeq` len `seq` fptr `seq` x
-
+ {-# INLINE deepSeqArray #-}
+ 
 
 -- Fill -----------------------------------------------------------------------
 -- | Filling of foreign buffers.
@@ -52,7 +53,6 @@
  data MArr F e 
   = FPArr !Int !(ForeignPtr e)
 
- {-# INLINE newMArr #-}
  newMArr n
   = do  let (proxy :: e) = undefined
         ptr              <- mallocBytes (sizeOf proxy * n)
@@ -60,32 +60,36 @@
         
         fptr             <- newForeignPtr finalizerFree ptr
         return           $ FPArr n fptr
+ {-# INLINE newMArr #-}
 
- {-# INLINE unsafeWriteMArr #-}
  unsafeWriteMArr (FPArr _ fptr) !ix !x
-  = withForeignPtr fptr
-  $ \ptr -> pokeElemOff ptr ix x
+  = pokeElemOff (Unsafe.unsafeForeignPtrToPtr fptr) ix x
+ {-# INLINE unsafeWriteMArr #-}
 
- {-# INLINE unsafeFreezeMArr #-}
  unsafeFreezeMArr !sh (FPArr len fptr)     
   =     return  $ AForeignPtr sh len fptr
+ {-# INLINE unsafeFreezeMArr #-}
 
+ deepSeqMArr !(FPArr _ ptr) x
+  = Unsafe.unsafeForeignPtrToPtr ptr `seq` x
+ {-# INLINE deepSeqMArr #-}
 
+
 -- 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
+{-# INLINE fromForeignPtr #-}
 
 
 -- | O(1). Unpack a `ForeignPtr` from an array.
 toForeignPtr :: Array F sh e -> ForeignPtr e
-{-# INLINE toForeignPtr #-}
 toForeignPtr (AForeignPtr _ _ fptr)
         = fptr
+{-# INLINE toForeignPtr #-}
 
 
 -- | Compute an array sequentially and write the elements into a foreign
@@ -94,9 +98,9 @@
 computeIntoS
         :: Fill r1 F sh e
         => ForeignPtr e -> Array r1 sh e -> IO ()
-{-# INLINE computeIntoS #-}
 computeIntoS !fptr !arr
  = fillS arr (FPArr 0 fptr)
+{-# INLINE computeIntoS #-}
 
 
 -- | Compute an array in parallel and write the elements into a foreign
@@ -105,7 +109,7 @@
 computeIntoP
         :: Fill r1 F sh e
         => ForeignPtr e -> Array r1 sh e -> IO ()
-{-# INLINE computeIntoP #-}
 computeIntoP !fptr !arr
  = fillP arr (FPArr 0 fptr)
+{-# INLINE computeIntoP #-}
 
diff --git a/Data/Array/Repa/Repr/Partitioned.hs b/Data/Array/Repa/Repr/Partitioned.hs
--- a/Data/Array/Repa/Repr/Partitioned.hs
+++ b/Data/Array/Repa/Repr/Partitioned.hs
@@ -34,49 +34,49 @@
                 (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
+{-# INLINE inRange #-}
 
 
 -- 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 index #-}
 
- {-# INLINE linearIndex #-}
  linearIndex arr@(APart sh _ _ _) ix
         = index arr $ fromIndex sh ix
+ {-# INLINE linearIndex #-}
 
- {-# INLINE extent #-}
  extent (APart sh _ _ _) 
         = sh
+ {-# INLINE extent #-}
 
- {-# INLINE deepSeqArray #-}
  deepSeqArray (APart sh range arr1 arr2) y
   = sh `deepSeq` range `deepSeqRange` arr1 `deepSeqArray` arr2 `deepSeqArray` y
+ {-# INLINE deepSeqArray #-}
 
 
-{-# INLINE deepSeqRange #-}
 deepSeqRange :: Shape sh => Range sh -> b -> b
 deepSeqRange (Range low high f) y
         = low `deepSeq` high `deepSeq` f `seq` y
+{-# INLINE deepSeqRange #-}
 
 
 -- 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 fillP #-}
 
- {-# INLINE fillS #-}
  fillS (APart _ (Range ix10 ix11 _) arr1 arr2) marr
   = do  fillRangeS arr1 marr ix10 ix11
         fillS arr2 marr
+ {-# INLINE fillS #-}
diff --git a/Data/Array/Repa/Repr/Unboxed.hs b/Data/Array/Repa/Repr/Unboxed.hs
--- a/Data/Array/Repa/Repr/Unboxed.hs
+++ b/Data/Array/Repa/Repr/Unboxed.hs
@@ -19,10 +19,11 @@
 
 -- | 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.
+--   The implementation uses @Data.Vector.Unboxed@ which 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
@@ -34,21 +35,21 @@
 -- 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 linearIndex #-}
 
- {-# INLINE unsafeLinearIndex #-}
  unsafeLinearIndex (AUnboxed _ vec) ix
         = vec `U.unsafeIndex` ix
+ {-# INLINE unsafeLinearIndex #-}
 
- {-# INLINE extent #-}
  extent (AUnboxed sh _)
         = sh
+ {-# INLINE extent #-}
 
- {-# INLINE deepSeqArray #-}
  deepSeqArray (AUnboxed sh vec) x 
   = sh `deepSeq` vec `seq` x
+ {-# INLINE deepSeqArray #-}
 
 
 -- Fill -----------------------------------------------------------------------
@@ -57,20 +58,24 @@
  data MArr U e 
   = UMArr (UM.IOVector e)
 
- {-# INLINE newMArr #-}
  newMArr n
   = liftM UMArr (UM.new n)
+ {-# INLINE newMArr #-}
 
- {-# INLINE unsafeWriteMArr #-}
  unsafeWriteMArr (UMArr v) ix
   = UM.unsafeWrite v ix
+ {-# INLINE unsafeWriteMArr #-}
 
- {-# INLINE unsafeFreezeMArr #-}
  unsafeFreezeMArr sh (UMArr mvec)     
   = do  vec     <- U.unsafeFreeze mvec
         return  $  AUnboxed sh vec
+ {-# INLINE unsafeFreezeMArr #-}
 
+ deepSeqMArr (UMArr vec) x
+  = vec `seq` x
+ {-# INLINE deepSeqMArr #-}
 
+
 -- Conversions ----------------------------------------------------------------
 -- | Sequential computation of array elements..
 --
@@ -79,8 +84,8 @@
 computeUnboxedS
         :: Fill r1 U sh e
         => Array r1 sh e -> Array U sh e
-{-# INLINE computeUnboxedS #-}
 computeUnboxedS = computeS
+{-# INLINE computeUnboxedS #-}
 
 
 -- | Parallel computation of array elements.
@@ -88,10 +93,10 @@
 --   * 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 #-}
+        :: (Fill r1 U sh e, Monad m, U.Unbox e)
+        => Array r1 sh e -> m (Array U sh e)
 computeUnboxedP = computeP
+{-# INLINE computeUnboxedP #-}
 
 
 -- | O(n). Convert a list to an unboxed vector array.
@@ -101,37 +106,38 @@
 fromListUnboxed
         :: (Shape sh, U.Unbox a)
         => sh -> [a] -> Array U sh a
-{-# INLINE fromListUnboxed #-}
 fromListUnboxed = R.fromList
+{-# INLINE fromListUnboxed #-}
 
 
 -- | 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
+{-# INLINE fromUnboxed #-}
 
 
 -- | 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
+{-# INLINE toUnboxed #-}
 
+
 -- 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)
+{-# INLINE zip #-}
 
 
 -- | O(1). Zip some unboxed arrays.
@@ -139,11 +145,11 @@
 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)
+{-# INLINE zip3 #-}
 
 
 -- | O(1). Zip some unboxed arrays.
@@ -151,11 +157,11 @@
 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)
+{-# INLINE zip4 #-}
 
 
 -- | O(1). Zip some unboxed arrays.
@@ -163,11 +169,11 @@
 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)
+{-# INLINE zip5 #-}
 
 
 -- | O(1). Zip some unboxed arrays.
@@ -175,11 +181,11 @@
 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)
+{-# INLINE zip6 #-}
  
 
 -- Unzip ----------------------------------------------------------------------
@@ -187,47 +193,47 @@
 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)
+{-# INLINE unzip #-}
 
 
 -- | 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)
+{-# INLINE unzip3 #-}
 
 
 -- | 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)
+{-# INLINE unzip4 #-}
 
 
 -- | 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)
+{-# INLINE unzip5 #-}
 
 
 -- | 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)
+{-# INLINE unzip6 #-}
diff --git a/Data/Array/Repa/Repr/Undefined.hs b/Data/Array/Repa/Repr/Undefined.hs
--- a/Data/Array/Repa/Repr/Undefined.hs
+++ b/Data/Array/Repa/Repr/Undefined.hs
@@ -19,19 +19,19 @@
 -- | Undefined array elements. Inspecting them yields `error`.
 --
 instance Repr X e where
- {-# INLINE deepSeqArray #-}
  deepSeqArray _ x
         = x
+ {-# INLINE deepSeqArray #-}
 
- {-# INLINE extent #-}
  extent (AUndefined sh) 
         = sh
+ {-# INLINE extent #-}
 
- {-# INLINE index #-}
  index (AUndefined _) _        = error "Repa: array element is undefined."
+ {-# INLINE index #-}
         
- {-# INLINE linearIndex #-}
  linearIndex (AUndefined _) _  = error "Repa: array element is undefined."
+ {-# INLINE linearIndex #-}
  
 
 instance (Shape sh, Fillable r2 e, Num e) => Fill X r2 sh e where
diff --git a/Data/Array/Repa/Repr/Vector.hs b/Data/Array/Repa/Repr/Vector.hs
--- a/Data/Array/Repa/Repr/Vector.hs
+++ b/Data/Array/Repa/Repr/Vector.hs
@@ -28,21 +28,21 @@
 -- Repr -----------------------------------------------------------------------
 -- | Read elements from a boxed vector array.
 instance Repr V a where
- {-# INLINE linearIndex #-}
  linearIndex (AVector _ vec) ix
         = vec V.! ix
+ {-# INLINE linearIndex #-}
 
- {-# INLINE unsafeLinearIndex #-}
  unsafeLinearIndex (AVector _ vec) ix
         = vec `V.unsafeIndex` ix
+ {-# INLINE unsafeLinearIndex #-}
 
- {-# INLINE extent #-}
  extent (AVector sh _)
         = sh
+ {-# INLINE extent #-}
 
- {-# INLINE deepSeqArray #-}
  deepSeqArray (AVector sh vec) x 
   = sh `deepSeq` vec `seq` x
+ {-# INLINE deepSeqArray #-}
 
 
 -- Fill -----------------------------------------------------------------------
@@ -51,19 +51,22 @@
  data MArr V e 
   = MVec (VM.IOVector e)
 
- {-# INLINE newMArr #-}
  newMArr n
   = liftM MVec (VM.new n)
+ {-# INLINE newMArr #-}
 
- {-# INLINE unsafeWriteMArr #-}
  unsafeWriteMArr (MVec v) ix
   = VM.unsafeWrite v ix
+ {-# INLINE unsafeWriteMArr #-}
 
- {-# INLINE unsafeFreezeMArr #-}
  unsafeFreezeMArr sh (MVec mvec)     
   = do  vec     <- V.unsafeFreeze mvec
         return  $  AVector sh vec
+ {-# INLINE unsafeFreezeMArr #-}
 
+ deepSeqMArr !_vec x
+  = x
+ {-# INLINE deepSeqMArr #-}
 
 -- Conversions ----------------------------------------------------------------
 -- | Sequential computation of array elements.
@@ -73,16 +76,16 @@
 computeVectorS
         :: Fill r1 V sh e
         => Array r1 sh e -> Array V sh e
-{-# INLINE computeVectorS #-}
 computeVectorS   = computeS
+{-# INLINE computeVectorS #-}
 
 
 -- | Parallel computation of array elements.
 computeVectorP
-        :: Fill r1 V sh e
-        => Array r1 sh e -> Array V sh e
-{-# INLINE computeVectorP #-}
+        :: (Fill r1 V sh e, Monad m)
+        => Array r1 sh e -> m (Array V sh e)
 computeVectorP   = computeP
+{-# INLINE computeVectorP #-}
 
 
 -- | O(n). Convert a list to a boxed vector array.
@@ -90,23 +93,23 @@
 --   * This is an alias for `fromList` with a more specific type.
 --
 fromListVector :: Shape sh => sh -> [a] -> Array V sh a
-{-# INLINE fromListVector #-}
 fromListVector  = fromList
+{-# INLINE fromListVector #-}
 
 
 -- | 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
+{-# INLINE fromVector #-}
 
 
 -- | O(1). Unpack a boxed vector from an array.
 toVector :: Array V sh e -> V.Vector e
-{-# INLINE toVector #-}
 toVector (AVector _ vec)
         = vec
+{-# INLINE toVector #-}
 
 
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
@@ -8,6 +8,7 @@
 	, makeBordered2)
 where
 import Data.Array.Repa.Index
+import Data.Array.Repa.Base
 import Data.Array.Repa.Repr.Partitioned
 import Data.Array.Repa.Repr.Undefined
 
@@ -70,10 +71,9 @@
 --   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
-	:: DIM2			-- ^ Extent of array.
+	:: (Repr r1 a, Repr r2 a)
+        => DIM2			-- ^ Extent of array.
 	-> Int			-- ^ Width of border.
 	-> Array r1 DIM2 a	-- ^ Array for internal elements.
 	-> Array r2 DIM2 a	-- ^ Array for border elements.
@@ -81,21 +81,21 @@
 
 {-# INLINE makeBordered2 #-}
 makeBordered2 sh@(_ :. aHeight :. aWidth) borderWidth arrInternal arrBorder
- = let
+ = checkDims `seq` 
+   let
 	-- minimum and maximum indicies of values in the inner part of the image.
 	!xMin		= borderWidth
 	!yMin		= borderWidth
 	!xMax		= aWidth  - borderWidth  - 1
 	!yMax		= aHeight - borderWidth - 1
 
-
-	{-# INLINE inInternal #-}
 	inInternal (Z :. y :. x)
 		=  x >= xMin && x <= xMax
 		&& y >= yMin && y <= yMax
+        {-# INLINE inInternal #-}
 
-	{-# INLINE inBorder #-}
 	inBorder 	= not . inInternal
+        {-# INLINE inBorder #-}
 
    in	
     --  internal region
@@ -107,3 +107,11 @@
     $   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
+
+ where
+        checkDims
+         = if (extent arrInternal) == (extent arrBorder)
+                then ()
+                else error "makeBordered2: internal and border arrays have different extents"
+        {-# NOINLINE checkDims #-}
+        --  NOINLINE because we don't want the branch in the core code.
diff --git a/Data/Array/Repa/Stencil/Dim2.hs b/Data/Array/Repa/Stencil/Dim2.hs
--- a/Data/Array/Repa/Stencil/Dim2.hs
+++ b/Data/Array/Repa/Stencil/Dim2.hs
@@ -1,3 +1,4 @@
+{-# LANGUAGE MagicHash #-}
 --   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
@@ -23,6 +24,7 @@
 import Data.Array.Repa.Repr.Undefined
 import Data.Array.Repa.Stencil.Base
 import Data.Array.Repa.Stencil.Template
+import GHC.Exts
 
 -- | A index into the flat array.
 --   Should be abstract outside the stencil modules.
@@ -167,34 +169,34 @@
 	   (StencilStatic sExtent zero loads)
 	   arr cur
 
-	| _ :. sHeight :. sWidth	<- sExtent
-	, _ :. aHeight :. aWidth	<- extent arr
+	| _ :. sHeight      :. sWidth	    <- sExtent
+	, _ :. (I# aHeight) :. (I# aWidth)  <- extent arr
 	, sHeight <= 7, sWidth <= 7
 	= let
 		-- Get data from the manifest array.
 		{-# INLINE getData #-}
 		getData :: DIM2 -> a
-		getData (Z :. y :. x)
+		getData (Z :. (I# y) :. (I# x))
 		 = wrapLoadX x y
 
-		-- TODO: Inlining this into above makes SpecConstr choke
-		wrapLoadX :: Int -> Int -> a
+                {-# NOINLINE wrapLoadX #-}
+		wrapLoadX :: Int# -> Int# -> a
 		wrapLoadX !x !y
-		 | x < 0	= wrapLoadY 0      	 y
-		 | x >= aWidth	= wrapLoadY (aWidth - 1) y
+		 | x <# 0#	= wrapLoadY 0#      	   y
+		 | x >=# aWidth	= wrapLoadY (aWidth -# 1#) y
 		 | otherwise    = wrapLoadY x y
 
-		{-# INLINE wrapLoadY #-}
-		wrapLoadY :: Int -> Int -> a
+		{-# NOINLINE wrapLoadY #-}
+		wrapLoadY :: Int# -> Int# -> a
 		wrapLoadY !x !y
-		 | y <  0	= loadXY x 0
-		 | y >= aHeight = loadXY x (aHeight - 1)
-		 | otherwise    = loadXY x y
+		 | y <#  0#	 = loadXY x 0#
+		 | y >=# aHeight = loadXY x (aHeight -# 1#)
+		 | otherwise     = loadXY x y
 
 		{-# INLINE loadXY #-}
-		loadXY :: Int -> Int -> a
+		loadXY :: Int# -> Int# -> a
 		loadXY !x !y
-		 = arr `unsafeIndex` (Z :. y :.  x)
+		 = arr `unsafeIndex` (Z :. (I# y) :.  (I# x))
 
 		-- Build a function to pass data from the array to our stencil.
 		{-# INLINE oload #-}
diff --git a/Data/Array/Repa/Unsafe.hs b/Data/Array/Repa/Unsafe.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Repa/Unsafe.hs
@@ -0,0 +1,14 @@
+
+-- | Functions without sanity or bounds checks.
+module Data.Array.Repa.Unsafe
+        ( unsafeBackpermute
+        , unsafeBackpermuteDft
+        , unsafeSlice
+        , unsafeExtend
+        , unsafeTraverse
+        , unsafeTraverse2
+        , unsafeTraverse3
+        , unsafeTraverse4)
+where
+import Data.Array.Repa.Operators.IndexSpace
+import Data.Array.Repa.Operators.Traversal
diff --git a/repa.cabal b/repa.cabal
--- a/repa.cabal
+++ b/repa.cabal
@@ -1,5 +1,5 @@
 Name:                repa
-Version:             3.0.0.1
+Version:             3.1.0.1
 License:             BSD3
 License-file:        LICENSE
 Author:              The DPH Team
@@ -72,6 +72,7 @@
         Data.Array.Repa.Shape
         Data.Array.Repa.Slice
         Data.Array.Repa.Stencil
+        Data.Array.Repa.Unsafe
         Data.Array.Repa
 
   Other-modules: