diff --git a/Data/Array/Repa/Algorithms/Complex.hs b/Data/Array/Repa/Algorithms/Complex.hs
--- a/Data/Array/Repa/Algorithms/Complex.hs
+++ b/Data/Array/Repa/Algorithms/Complex.hs
@@ -1,44 +1,58 @@
-{-# LANGUAGE TypeOperators, TypeSynonymInstances #-}
+{-# LANGUAGE TypeOperators, TypeSynonymInstances, FlexibleInstances #-}
 
 -- | Strict complex doubles.
 module Data.Array.Repa.Algorithms.Complex
 	( Complex
 	, mag
-	, arg
-	, (:*:)(..))
+	, arg)
 where
-import 	Data.Array.Parallel.Base ((:*:)(..))
 
--- | Strict complex doubles.
+
+-- | Complex doubles.
 type Complex 
-	= Double :*: Double
+	= (Double, Double)
 
 instance Num Complex where
-  abs x				= (mag x) :*: 0
-  signum (re :*: _)		= signum re :*: 0
-  fromInteger n			= fromInteger n :*: 0.0
-  (r :*: i) + (r' :*: i')	= r+r' :*: i+i'
-  (r :*: i) - (r' :*: i')	= r-r' :*: i-i'
-  (r :*: i) * (r' :*: i')	= r*r' - i*i' :*: r*i' + r'*i
 
+  {-# INLINE abs #-}
+  abs x			= (mag x, 0)
 
+  {-# INLINE signum #-}
+  signum (re, _)	= (signum re, 0)
+
+  {-# INLINE fromInteger #-}
+  fromInteger n		= (fromInteger n, 0.0)
+
+  {-# INLINE (+) #-}
+  (r, i) + (r', i')	= (r+r', i+i')
+
+  {-# INLINE (-) #-}
+  (r, i) - (r', i')	= (r-r', i-i')
+
+  {-# INLINE (*) #-}
+  (r, i) * (r', i')	= (r*r' - i*i', r*i' + r'*i)
+
+
 instance Fractional Complex where
-  (a :*: b) / (c :*: d)		
+  {-# INLINE (/) #-}
+  (a, b) / (c, d)		
  	= let	den	= c^(2 :: Int) + d^(2 :: Int)
 		re	= (a * c + b * d) / den
 		im	= (b * c - a * d) / den
-	  in	re :*: im
+	  in	(re, im)
 	
-  fromRational x	= fromRational x :*: 0
+  fromRational x	= (fromRational x, 0)
 	
 -- | Take the magnitude of a complex number.
 mag :: Complex -> Double
-mag (r :*: i)	= sqrt (r * r + i * i)
+{-# INLINE mag #-}
+mag (r, i)	= sqrt (r * r + i * i)
 
 
 -- | Take the argument (phase) of a complex number, in the range [-pi .. pi].
 arg :: Complex -> Double
-arg (re :*: im)
+{-# INLINE arg #-}
+arg (re, im)
  = normaliseAngle $ atan2 im re
 
  where 	normaliseAngle :: Double -> Double
diff --git a/Data/Array/Repa/Algorithms/Convolve.hs b/Data/Array/Repa/Algorithms/Convolve.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Repa/Algorithms/Convolve.hs
@@ -0,0 +1,167 @@
+{-# LANGUAGE BangPatterns, PackageImports #-}
+{-# OPTIONS -Wall -fno-warn-missing-signatures -fno-warn-incomplete-patterns #-}
+
+-- | Old support for stencil based convolutions. 
+--
+--   NOTE: This is slated to be merged with the new Stencil support in the next version
+--         of Repa. We'll still expose the `convolve` function though.
+--
+module Data.Array.Repa.Algorithms.Convolve
+	( convolve
+
+	, GetOut
+	, outAs
+	, outClamp
+	, convolveOut )
+where
+import Data.Array.Repa 					as A
+import qualified Data.Vector.Unboxed			as V
+import qualified Data.Array.Repa.Shape			as S
+import Prelude						as P
+
+
+-- Plain Convolve ---------------------------------------------------------------------------------
+-- | Image-kernel convolution,
+--   which takes a function specifying what value to return when the kernel doesn't apply.
+convolve
+	:: (Elt a, Num a)
+	=> (DIM2 -> a) 		-- ^ Use this function to get border elements where the kernel apply.
+	-> Array DIM2 a		-- ^ Kernel to use in the convolution.
+	-> Array DIM2 a		-- ^ Input image.
+	-> Array DIM2 a
+
+{-# INLINE convolve #-}
+convolve makeOut
+ 	kernel@(Array       (_ :. krnHeight :. krnWidth) [Region RangeAll (GenManifest krnVec)])
+  	 image@(Array imgSh@(_ :. imgHeight :. imgWidth) [Region RangeAll (GenManifest imgVec)])
+
+ = kernel `deepSeqArray` image `deepSeqArray` 
+   force $ unsafeTraverse image id update
+ where	
+	!krnHeight2	= krnHeight `div` 2
+	!krnWidth2	= krnWidth  `div` 2
+
+	-- If we're too close to the edge of the input image then
+	-- we can't apply the stencil because we don't have enough data.
+	!borderLeft	= krnWidth2
+	!borderRight	= imgWidth   - krnWidth2  - 1
+	!borderUp	= krnHeight2
+	!borderDown	= imgHeight  - krnHeight2 - 1
+
+	{-# INLINE update #-}
+	update _ ix@(_ :. j :. i)
+ 	 | i < borderLeft	= makeOut ix
+ 	 | i > borderRight	= makeOut ix
+  	 | j < borderUp		= makeOut ix
+ 	 | j > borderDown	= makeOut ix
+	 | otherwise		= stencil j i
+
+	-- The actual stencil function.
+	{-# INLINE stencil #-}
+	stencil j i
+	 = let	imgStart = S.toIndex imgSh (Z :. j - krnHeight2 :. i - krnWidth2)
+	   in	integrate 0 0 0 imgStart 0
+
+	{-# INLINE integrate #-}
+	integrate !acc !x !y !imgCur !krnCur  
+	 | y >= krnHeight
+	 = acc
+
+	 | x >= krnWidth
+	 = integrate acc 0 (y + 1) (imgCur + imgWidth - krnWidth) krnCur 
+	
+	 | otherwise
+	 = let	imgZ	= imgVec `V.unsafeIndex` imgCur 
+		krnZ	= krnVec `V.unsafeIndex` krnCur 
+		here	= imgZ * krnZ 
+	   in	integrate (acc + here) (x + 1) y (imgCur + 1) (krnCur + 1)
+
+
+-- Convolve Out -----------------------------------------------------------------------------------
+-- | A function that gets out of range elements from an image.
+type GetOut a
+	= (DIM2 -> a) 	-- ^ The original get function.
+	-> DIM2 	-- ^ The shape of the image.
+	-> DIM2 	-- ^ Index of element we were trying to get.
+	-> a
+
+
+-- | Use the provided value for every out-of-range element.
+outAs :: a -> GetOut a
+{-# INLINE outAs #-}
+outAs x _ _ _ = x
+
+
+-- | If the requested element is out of range use
+--   the closest one from the real image.
+outClamp :: GetOut a
+{-# INLINE outClamp #-}
+outClamp get (_ :. yLen :. xLen) (sh :. j :. i)
+ = clampX j i
+ where 	{-# INLINE clampX #-}
+	clampX !y !x
+	  | x < 0	= clampY y 0
+	  | x >= xLen	= clampY y (xLen - 1)
+	  | otherwise	= clampY y x
+		
+	{-# INLINE clampY #-}
+	clampY !y !x
+	  | y < 0	= get (sh :. 0 		:. x)
+	  | y >= yLen	= get (sh :. (yLen - 1) :. x)
+	  | otherwise	= get (sh :. y 		:. x)
+
+
+-- | Image-kernel convolution, 
+--   which takes a function specifying what value to use for out-of-range elements.
+convolveOut
+	:: (Elt a, Num a)
+	=> GetOut a		-- ^ Use this fn to get out of range elements.
+	-> Array DIM2 a		-- ^ Kernel
+	-> Array DIM2 a		-- ^ Image
+	-> Array DIM2 a
+
+{-# INLINE convolveOut #-}
+convolveOut getOut
+ 	kernel@(Array krnSh@(_ :. krnHeight :. krnWidth) _)
+  	 image@(Array imgSh@(_ :. imgHeight :. imgWidth) _)
+
+ = kernel `deepSeqArray` image `deepSeqArray` 
+   force $ unsafeTraverse image id stencil
+ where	
+	!krnHeight2	= krnHeight `div` 2
+	!krnWidth2	= krnWidth  `div` 2
+        !krnSize	= S.size krnSh
+
+	-- If we're too close to the edge of the input image then
+	-- we can't apply the stencil because we don't have enough data.
+	!borderLeft	= krnWidth2
+	!borderRight	= imgWidth   - krnWidth2  - 1
+	!borderUp	= krnHeight2
+	!borderDown	= imgHeight  - krnHeight2 - 1
+
+	-- The actual stencil function.
+	{-# INLINE stencil #-}
+	stencil get (_ :. j :. i)
+	 = let
+		{-# INLINE get' #-}
+		get' ix@(_ :. y :. x)
+		 | x < borderLeft	= getOut get imgSh ix
+		 | x > borderRight	= getOut get imgSh ix
+		 | y < borderUp		= getOut get imgSh ix
+		 | y > borderDown	= getOut get imgSh ix
+		 | otherwise		= get ix
+
+		!ikrnWidth'	= i - krnWidth2
+		!jkrnHeight'	= j - krnHeight2
+
+		{-# INLINE integrate #-}
+		integrate !count !acc
+		 | count == krnSize		= acc
+		 | otherwise
+		 = let	!ix@(sh :. y :. x)	= S.fromIndex krnSh count
+			!ix'			= sh :. y + jkrnHeight' :. x + ikrnWidth'
+			!here			= kernel `unsafeIndex` ix * (get' ix')
+		   in	integrate (count + 1) (acc + here)
+
+	   in	integrate 0 0
+
diff --git a/Data/Array/Repa/Algorithms/DFT.hs b/Data/Array/Repa/Algorithms/DFT.hs
--- a/Data/Array/Repa/Algorithms/DFT.hs
+++ b/Data/Array/Repa/Algorithms/DFT.hs
@@ -22,6 +22,7 @@
 import Data.Array.Repa.Algorithms.DFT.Roots
 import Data.Array.Repa.Algorithms.Complex
 import Data.Array.Repa				as A
+import Prelude					as P
 
 -- | Compute the DFT along the low order dimension of an array.
 dft 	:: forall sh
@@ -42,7 +43,7 @@
 
 idft v
  = let	_ :. len	= extent v
-	scale		= fromIntegral len :*: 0
+	scale		= (fromIntegral len, 0)
 	rofu		= calcInverseRootsOfUnity (extent v)
    in	force $ A.map (/ scale) $ dftWithRoots rofu v
 
@@ -62,8 +63,8 @@
 	| _ :. rLen 	<- extent rofu
 	, _ :. vLen 	<- extent arr
 	, rLen /= vLen
-	= error $  "dftWithRoots: length of vector (" ++ show vLen ++ ")"
-		++ " does not match the length of the roots (" ++ show rLen ++ ")"
+	= error $    "dftWithRoots: length of vector (" P.++ show vLen P.++ ")"
+		P.++ " does not match the length of the roots (" P.++ show rLen P.++ ")"
 
 	| otherwise
 	= traverse arr id (\_ k -> dftWithRootsSingle rofu arr k)
@@ -84,8 +85,8 @@
 	| _ :. rLen 	<- extent rofu
 	, _ :. vLen 	<- extent arrX
 	, rLen /= vLen
-	= error $  "dftWithRootsSingle: length of vector (" ++ show vLen ++ ")"
-		++ " does not match the length of the roots (" ++ show rLen ++ ")"
+	= error $    "dftWithRootsSingle: length of vector (" P.++ show vLen P.++ ")"
+		P.++ " does not match the length of the roots (" P.++ show rLen P.++ ")"
 
 	| otherwise
 	= let	sh@(_ :. len)	= extent arrX
@@ -93,7 +94,7 @@
 		-- All the roots we need to multiply with.
 		wroots		= fromFunction sh elemFn
 		elemFn (sh' :. n) 
-			= rofu !: (sh' :. (k * n) `mod` len)
+			= rofu ! (sh' :. (k * n) `mod` len)
 
 	  in  A.sumAll $ A.zipWith (*) arrX wroots
 
diff --git a/Data/Array/Repa/Algorithms/DFT/Center.hs b/Data/Array/Repa/Algorithms/DFT/Center.hs
--- a/Data/Array/Repa/Algorithms/DFT/Center.hs
+++ b/Data/Array/Repa/Algorithms/DFT/Center.hs
@@ -2,30 +2,32 @@
 -- | Applying these transforms to the input of a DFT causes the output 
 --   to be centered so that the zero frequency is in the middle. 
 module Data.Array.Repa.Algorithms.DFT.Center
-	( centerVector
-	, centerMatrix)
+	( center1d
+	, center2d
+	, center3d)
 where
 import Data.Array.Repa
 import Data.Array.Repa.Algorithms.Complex
 
-
 -- | Apply the centering transform to a vector.
-centerVector
-	:: Array DIM1 Complex
-	-> Array DIM1 Complex
-
-{-# INLINE centerVector #-}
-centerVector arr
+center1d :: Array DIM1 Complex -> Array DIM1 Complex
+{-# INLINE center1d #-}
+center1d arr
  = traverse arr id
 	(\get ix@(_ :. x) -> ((-1) ^ x) * get ix)
 
 
 -- | Apply the centering transform to a matrix.
-centerMatrix
-	:: Array DIM2 Complex
-	-> Array DIM2 Complex
-
-{-# INLINE centerMatrix #-}
-centerMatrix arr
+center2d :: Array DIM2 Complex -> Array DIM2 Complex
+{-# INLINE center2d #-}
+center2d arr
  = traverse arr id
 	(\get ix@(_ :. y :. x) -> ((-1) ^ (y + x)) * get ix)
+
+
+-- | Apply the centering transform to a 3d array.
+center3d :: Array DIM3 Complex -> Array DIM3 Complex
+{-# INLINE center3d #-}
+center3d arr
+ = traverse arr id
+	(\get ix@(_ :. z :. y :. x) -> ((-1) ^ (z + y + x)) * get ix)
diff --git a/Data/Array/Repa/Algorithms/DFT/Roots.hs b/Data/Array/Repa/Algorithms/DFT/Roots.hs
--- a/Data/Array/Repa/Algorithms/DFT/Roots.hs
+++ b/Data/Array/Repa/Algorithms/DFT/Roots.hs
@@ -19,8 +19,9 @@
  = force $ fromFunction sh f
  where
     f :: Shape sh => (sh :. Int) -> Complex
-    f (_ :. i) =      (cos  (2 * pi * (fromIntegral i) / len))
-		:*: (- sin  (2 * pi * (fromIntegral i) / len))
+    f (_ :. i) 
+	= ( cos  (2 * pi * (fromIntegral i) / len)
+	  , - sin  (2 * pi * (fromIntegral i) / len))
 
     len	= fromIntegral n
 
@@ -36,7 +37,8 @@
  = force $ fromFunction sh f
  where
     f :: Shape sh => (sh :. Int) -> Complex
-    f (_ :. i) =      (cos  (2 * pi * (fromIntegral i) / len))
-		:*:   (sin  (2 * pi * (fromIntegral i) / len))
+    f (_ :. i) 
+	= ( cos  (2 * pi * (fromIntegral i) / len)
+	  , sin  (2 * pi * (fromIntegral i) / len))
 
     len	= fromIntegral n
diff --git a/Data/Array/Repa/Algorithms/FFT.hs b/Data/Array/Repa/Algorithms/FFT.hs
--- a/Data/Array/Repa/Algorithms/FFT.hs
+++ b/Data/Array/Repa/Algorithms/FFT.hs
@@ -1,212 +1,185 @@
-{-# LANGUAGE TypeOperators, PatternGuards, RankNTypes #-}
+{-# LANGUAGE TypeOperators, PatternGuards, RankNTypes, ScopedTypeVariables, BangPatterns, FlexibleContexts #-}
+{-# OPTIONS -fno-warn-incomplete-patterns #-}
 
--- | Fast computation of Discrete Fourier Transforms using the Cooley-Tuckey algorithm.
---
---   Time complexity is O(n log n) in the size of the input.
---
---   Input dimensions must be powers of two, else `error`.
---
---   The `fft` and `ifft` functions (and friends) also compute the roots of unity needed.
---   If you need to transform several arrays with the same extent then it is faster to
---   compute the roots once using `calcRootsOfUnity` or `calcInverseRootsOfUnity`, 
---   then call `fftWithRoots` directly.
+-- | Fast computation of Discrete Fourier Transforms using the Cooley-Tuckey algorithm. 
+--   Time complexity is O(n log n) in the size of the input. 
 --
---   The inverse transforms provided also perform post-scaling so that `ifft` is the true inverse of `fft`. 
---   If you don't want that then call `fftWithRoots` directly.
+--   This uses a naive divide-and-conquer algorithm, the absolute performance is about
+--   50x slower than FFTW in estimate mode.
 --
---   The functions `fft2d` and `fft3d` require their inputs to be squares (and cubes) respectively. 
---   This allows them to reuse the same roots-of-unity when transforming along each axis. If you 
---   need to transform rectanglular arrays then call `fftWithRoots` directly.
 module Data.Array.Repa.Algorithms.FFT
-	( fft,   ifft
-	, fft2d, ifft2d
-	, fft3d, ifft3d
-	, fftWithRoots )
+	( Mode(..)
+	, isPowerOfTwo
+	, fft3d
+	, fft2d
+	, fft1d)
 where
-import Data.Array.Repa.Algorithms.DFT.Roots
 import Data.Array.Repa.Algorithms.Complex
 import Data.Array.Repa				as A
-import Data.Ratio
 
--- Vector Transform -------------------------------------------------------------------------------
--- | Compute the DFT along the low order dimension of an array.
-fft	:: Shape sh
-	=> Array (sh :. Int) Complex
-	-> Array (sh :. Int) Complex
+data Mode
+	= Forward
+	| Reverse
+	| Inverse
+	deriving (Show, Eq)
 
-fft v
- = let	rofu	= calcRootsOfUnity (extent v)
-   in	force $ fftWithRoots rofu v
+{-# INLINE signOfMode #-}
+signOfMode :: Mode -> Double
+signOfMode mode
+ = case mode of
+	Forward		-> (-1)
+	Reverse		->   1
+	Inverse		->   1
 
 
--- | Compute the inverse DFT along the low order dimension of an array.
-ifft	:: Shape sh
-	=> Array (sh :. Int) Complex
-	-> Array (sh :. Int) Complex
+{-# INLINE isPowerOfTwo #-}
+-- | Check if an `Int` is a power of two.
+isPowerOfTwo :: Int -> Bool
+isPowerOfTwo x
+ = let	r	= (log (fromIntegral x) / log 2) :: Double
+   in	ceiling r == (floor r :: Int)
 
-ifft v
- = let	_ :. len	= extent v
-	scale		= fromIntegral len :*: 0
-	rofu		= calcInverseRootsOfUnity (extent v)
-   in	force $ A.map (/ scale) $ fftWithRoots rofu v
 
 
--- Matrix Transform -------------------------------------------------------------------------------
--- | Compute the DFT of a square matrix.
---   If the matrix is not square then `error`.
-fft2d 	:: Array DIM2 Complex
-	-> Array DIM2 Complex
 
-fft2d arr
- 	| Z :. height :. width	<- extent arr
- 	, height /= width	
-	= error $ "fft2d: height of matrix (" ++ show height ++ ")"
-		++  " does not match width (" ++ show width  ++ ")"
-
-	| otherwise
-	= let	rofu		= calcRootsOfUnity (extent arr)
-  		fftTrans 	= transpose . fftWithRoots rofu
-   	  in	force $ fftTrans $ fftTrans arr
-
-
--- | Compute the inverse DFT of a square matrix. 
-ifft2d	:: Array DIM2 Complex
-	-> Array DIM2 Complex
-	
-ifft2d arr
- 	| Z :. height :. width	<- extent arr
- 	, height /= width	
-	= error $ "fft2d: height of matrix (" ++ show height ++ ")"
-		++  " does not match width (" ++ show width  ++ ")"
+-- 3D Transform -----------------------------------------------------------------------------------
+-- | Compute the DFT of a 3d array. Array dimensions must be powers of two else `error`.
+fft3d 	:: Mode
+	-> Array DIM3 Complex
+	-> Array DIM3 Complex
 
-	| otherwise
-	= let	_ :. height :. width = extent arr
-		scale		= fromIntegral (height * width) :*: 0
-		rofu		= calcInverseRootsOfUnity (extent arr)
-		fftTrans	= transpose . fftWithRoots rofu
-	  in	force $ A.map (/ scale) $ fftTrans $ fftTrans arr
-	
+fft3d mode arr
+ = let	_ :. depth :. height :. width	= extent arr
+	!sign	= signOfMode mode
+	!scale 	= fromIntegral (depth * width * height) 
+		
+   in	if not (isPowerOfTwo depth && isPowerOfTwo height && isPowerOfTwo width)
+	 then error "Data.Array.Repa.Algorithms.FFT: fft3d -- array dimensions must be powers of two."
+	 else arr `deepSeqArray` 
+		case mode of
+			Forward	-> fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr
+			Reverse	-> fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr
+			Inverse	-> force $ A.map (/ scale) 
+					$ fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr
 
--- Cube Transform ---------------------------------------------------------------------------------
--- | Compute the DFT of a 3d cube.
---   If the array is not a cube then `error`.
-fft3d 	:: Array DIM3 Complex
+fftTrans3d 
+	:: Double
+	-> Array DIM3 Complex 
 	-> Array DIM3 Complex
 
-fft3d arrIn
- 	| Z :. depth :. height :. width	<- extent arrIn
- 	, (height /= width) || (height /= depth)
-	= error $ "fft3d: array is not a cube"
+{-# NOINLINE fftTrans3d #-}
+fftTrans3d sign arr'
+ = let 	arr		= force arr'
+	(sh :. len)	= extent arr
+   in	force $ rotate3d $ fft sign sh len arr
 
-	| otherwise
-	= let	rofu		= calcRootsOfUnity (extent arrIn)
 
-		transpose3 arr
-	 	 = traverse arr 
-        		(\(Z :. k :. l :. m)   -> (Z :. l :. m :. k)) 
-            		(\f (Z :. l :. m :. k) -> f (Z :. k :. l :. m)) 
+rotate3d :: Array DIM3 Complex -> Array DIM3 Complex
+{-# INLINE rotate3d #-}
+rotate3d arr
+ = backpermute (sh :. m :. k :. l) f arr
+ where	(sh :. k :. l :. m)		= extent arr
+	f (sh' :. m' :. k' :. l')	= sh' :. k' :. l' :. m'
 
-		fftTrans	= transpose3 . fftWithRoots rofu
-	
-  	  in	force $ fftTrans $ fftTrans $ fftTrans arrIn
 
 
--- | Compute the inverse DFT of a 3d cube.
---   If the array is not a cube then `error`.
-ifft3d 	:: Array DIM3 Complex
-	-> Array DIM3 Complex
+-- Matrix Transform -------------------------------------------------------------------------------
+-- | Compute the DFT of a matrix. Array dimensions must be powers of two else `error`.
+fft2d 	:: Mode
+	-> Array DIM2 Complex
+	-> Array DIM2 Complex
 
-ifft3d arrIn
- 	| Z :. depth :. height :. width	<- extent arrIn
- 	, (height /= width) || (height /= depth)
-	= error $ "ifft3d: array is not a cube"
+fft2d mode arr
+ = let	_ :. height :. width	= extent arr
+	sign	= signOfMode mode
+	scale 	= fromIntegral (width * height) 
+		
+   in	if not (isPowerOfTwo height && isPowerOfTwo width)
+	 then error "Data.Array.Repa.Algorithms.FFT: fft2d -- array dimensions must be powers of two."
+	 else arr `deepSeqArray` 
+		case mode of
+			Forward	-> fftTrans2d sign $ fftTrans2d sign arr
+			Reverse	-> fftTrans2d sign $ fftTrans2d sign arr
+			Inverse	-> force $ A.map (/ scale) $ fftTrans2d sign $ fftTrans2d sign arr
 
-	| otherwise
-	= let	rofu		= calcInverseRootsOfUnity (extent arrIn)
+fftTrans2d 
+	:: Double
+	-> Array DIM2 Complex 
+	-> Array DIM2 Complex
 
-		transpose3 arr
-	 	 = traverse arr 
-        		(\(Z :. k :. l :. m)   -> (Z :. l :. m :. k)) 
-            		(\f (Z :. l :. m :. k) -> f (Z :. k :. l :. m)) 
+{-# NOINLINE fftTrans2d #-}
+fftTrans2d sign arr'
+ = let 	arr		= force arr'
+	(sh :. len)	= extent arr
+   in	force $ transpose $ fft sign sh len arr
 
-		_ :. depth :. height :. width 
-				= extent arrIn
-		scale		= fromIntegral (height * width * depth) :*: 0
 
-		fftTrans	= transpose3 . fftWithRoots rofu
-	  in	force $ A.map (/ scale) $ fftTrans $ fftTrans $ fftTrans arrIn
-
-	
--- Worker -----------------------------------------------------------------------------------------
--- | Generic function for computation of forward or inverse Discrete Fourier Transforms.
---	Computation is along the low order dimension of the array.
-fftWithRoots	
-	:: forall sh
-	.  Shape sh
-	=> Array (sh :. Int) Complex		-- ^ Roots of unity.
-	-> Array (sh :. Int) Complex		-- ^ Input values.
-        -> Array (sh :. Int) Complex
-
-fftWithRoots rofu v
-	| not $ (denominator $ toRational (logBase (2 :: Double) $ fromIntegral vLen)) == 1
-	= error $ "fft: vector length of " ++ show vLen ++ " is not a power of 2"
+-- Vector Transform -------------------------------------------------------------------------------
+-- | Compute the DFT of a vector. Array dimensions must be powers of two else `error`.
+fft1d	:: Mode 
+	-> Array DIM1 Complex 
+	-> Array DIM1 Complex
 	
-	| rLen /= vLen
-	= error $  "fft: length of vector (" ++ show vLen ++ ")"
-		++ " does not match the length of the roots (" ++ show rLen ++ ")"
+fft1d mode arr
+ = let	_ :. len	= extent arr
+	sign	= signOfMode mode
+	scale	= fromIntegral len
 	
-	| otherwise
-	= fftWithRoots' rofu v
+   in	if not $ isPowerOfTwo len
+	 then error "Data.Array.Repa.Algorithms.FFT: fft1d -- array dimensions must be powers of two."
+	 else arr `deepSeqArray`
+		case mode of
+			Forward	-> fftTrans1d sign arr
+			Reverse	-> fftTrans1d sign arr
+			Inverse -> force $ A.map (/ scale) $ fftTrans1d sign arr
 
-	where	_ :. rLen	= extent rofu
-		_ :. vLen	= extent v
+fftTrans1d
+	:: Double 
+	-> Array DIM1 Complex
+	-> Array DIM1 Complex
 
-fftWithRoots'
-	:: Shape sh
-	=> Array (sh :. Int) Complex
-	-> Array (sh :. Int) Complex
-        -> Array (sh :. Int) Complex
+{-# NOINLINE fftTrans1d #-}
+fftTrans1d sign arr'
+ = let	arr		= force arr'
+	(sh :. len)	= extent arr
+   in	fft sign sh len arr
 
-{-# INLINE fftWithRoots' #-}
-fftWithRoots' rofu v
- = case extent v of
-	_ :. 2	-> fft_two   v
-	_	-> fft_split rofu v
 
-{-# INLINE fft_two #-}
-fft_two v
- = let	vFn' vFn (sh :. 0)  = vFn (sh :. 0) + vFn (sh :. 1)
-	vFn' vFn (sh :. 1)  = vFn (sh :. 0) - vFn (sh :. 1)
-	vFn' _   _          = error "Data.Array.Repa.Algorithms.FFT fft_two fail"
-   in	traverse v id vFn'
+-- Rank Generalised Worker ------------------------------------------------------------------------
+{-# INLINE fft #-}
+fft !sign !sh !lenVec !vec
+ = go lenVec 0 1
+ where	go !len !offset !stride
+	 | len == 2
+	 = force $ fromFunction (sh :. 2) swivel
 	
-{-# INLINE fft_split #-}
-fft_split rofu v
- = let 	fft_lr = force $ fftWithRoots' (splitRofu rofu) (splitVector v)
+	 | otherwise
+	 = combine len 
+		(go (len `div` 2) offset            (stride * 2))
+		(go (len `div` 2) (offset + stride) (stride * 2))
 
-	fft_l  = traverse2 fft_lr rofu 
- 		   (\(sh :. 2 :. n) _ -> sh :. n)
-		   (\f r (sh :. i)    -> f (sh :. 0 :. i) + r (sh :. i) * f (sh :. 1 :. i))
+	 where	swivel (sh' :. ix)
+		 = case ix of
+			0	-> (vec `unsafeIndex` (sh' :. offset)) + (vec `unsafeIndex` (sh' :. (offset + stride)))
+			1	-> (vec `unsafeIndex` (sh' :. offset)) - (vec `unsafeIndex` (sh' :. (offset + stride)))
 
-	fft_r  = traverse2 fft_lr rofu 
-		   (\(sh :. 2 :. n) _ -> sh :. n)
-		   (\f r (sh :. i)    -> f (sh :. 0 :. i) - r (sh :. i) * f (sh :. 1 :. i))
+		{-# INLINE combine #-}
+		combine !len' 	evens@(Array _ [Region RangeAll GenManifest{}]) 
+				 odds@(Array _ [Region RangeAll GenManifest{}])
+ 	 	 = evens `deepSeqArray` odds `deepSeqArray`
+   	   	   let	odds'	= unsafeTraverse odds id (\get ix@(_ :. k) -> twiddle sign k len' * get ix) 
+   	   	   in	force 	$ (evens +^ odds') A.++ (evens -^ odds')
 
-   in	fft_l +:+ fft_r
 
-{-# INLINE splitRofu #-}
-splitRofu rofu
- = traverse rofu
-	(\(rSh :. rLen) 	-> rSh :. (2::Int) :. (rLen `div` 2))
-	(\rFn (sh :. _ :. i) 	-> rFn (sh :. 2*i))
-
-{-# INLINE splitVector #-}
-splitVector v 
- = let	vFn' vFn (sh :. 0 :. i) = vFn (sh :. 2*i)
-	vFn' vFn (sh :. 1 :. i) = vFn (sh :. 2*i+1)
-	vFn' _   _              = error "Data.Array.Repa.Algorithms.FFT splitVector fail"
+-- Compute a twiddle factor.
+twiddle :: Double
+	-> Int 			-- index
+	-> Int 			-- length
+	-> Complex
 
-   in	traverse v
-		(\(vSh :. vLen)    -> vSh :. 2 :. (vLen `div` 2)) 
-		vFn'
-        
+{-# INLINE twiddle #-}
+twiddle sign k' n'
+ 	=  (cos (2 * pi * k / n), sign * sin  (2 * pi * k / n))
+	where 	k	= fromIntegral k'
+		n	= fromIntegral n'
+      
diff --git a/Data/Array/Repa/Algorithms/Iterate.hs b/Data/Array/Repa/Algorithms/Iterate.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Repa/Algorithms/Iterate.hs
@@ -0,0 +1,45 @@
+{-# LANGUAGE BangPatterns, PackageImports #-}
+{-# OPTIONS -Wall -fno-warn-missing-signatures -fno-warn-incomplete-patterns #-}
+module Data.Array.Repa.Algorithms.Iterate
+	(iterateBlockwise, iterateBlockwise')
+where
+import Data.Array.Repa
+
+-- | Iterate array transformation function a fixed number of times, applying `force2` between
+--   each iteration.
+iterateBlockwise
+	:: (Elt a, Num a)
+	=> Int					-- ^ Number of iterations to run for.
+	-> (Array DIM2 a -> Array DIM2 a)	-- ^ Fn to step the array.
+	-> Array DIM2 a				-- ^ Initial array value.
+	-> Array DIM2 a
+	
+{-# INLINE iterateBlockwise #-}
+iterateBlockwise steps f arrInit@(Array shInit [Region RangeAll (GenManifest vecInit)])
+ = arrInit `deepSeqArray`
+  goSolve steps shInit vecInit
+
+ where	-- NOTE: We manually unpack the current array into its shape and vector to
+	--	 stop GHC from unboxing the vector again for every loop. deepSeqing
+	--	 the arrays at the start of solveLaplace makes the unboxings happen
+	--	 at that point in the corresponding core code.
+	goSolve !i !shCurrent !vecCurrent
+	 = let	!arrCurrent	= fromVector shCurrent vecCurrent
+	   in   if i == 0 
+		 then arrCurrent
+		 else let arrNew@(Array _ [Region RangeAll (GenManifest _)]) 
+				= force2 $ f arrCurrent
+		      in  goSolve (i - 1) (extent arrNew) (toVector arrNew)
+
+
+-- | As above, but with the parameters flipped.
+iterateBlockwise'
+	:: (Elt a, Num a)
+	=> Int					-- ^ Number of iterations to run for.
+	-> Array DIM2 a				-- ^ Initial array value.
+	-> (Array DIM2 a -> Array DIM2 a)	-- ^ Fn to step the array.
+	-> Array DIM2 a
+
+{-# INLINE iterateBlockwise' #-}
+iterateBlockwise' steps arr fn
+	= iterateBlockwise steps fn arr
diff --git a/Data/Array/Repa/Algorithms/Matrix.hs b/Data/Array/Repa/Algorithms/Matrix.hs
--- a/Data/Array/Repa/Algorithms/Matrix.hs
+++ b/Data/Array/Repa/Algorithms/Matrix.hs
@@ -1,4 +1,5 @@
 {-# OPTIONS -fno-warn-incomplete-patterns #-}
+{-# LANGUAGE PackageImports #-}
 
 -- | Algorithms operating on matrices.
 -- 
@@ -11,20 +12,30 @@
 module Data.Array.Repa.Algorithms.Matrix
 	(multiplyMM)
 where
-import Data.Array.Repa
-	
+import Data.Array.Repa	as A
 
 -- | Matrix-matrix multiply.
---
-multiplyMM	
+multiplyMM
 	:: Array DIM2 Double
 	-> Array DIM2 Double
 	-> Array DIM2 Double
 
-multiplyMM  arr1 arr2
- = multiplyMM' (force arr1) (force arr2)
- where multiplyMM' arr1'@Manifest{} arr2'@Manifest{}
-	= fold (+) 0 
-  	$ traverse2 arr1' (force $ transpose arr2')
-      		(\(sh :. m1 :. n1) -> \(_ :. n2 :. _m2) -> (sh :. m1 :. n2 :. n1))
-		(\f1 -> \f2 -> \(sh :. i :. j :. k) -> f1 (sh :. i :. k) * f2 (sh :. j :. k))
+{-# NOINLINE multiplyMM #-}
+multiplyMM arr@(Array _ [Region RangeAll (GenManifest _)])
+	   brr@(Array _ [Region RangeAll (GenManifest _)])
+ = [arr, brr] `deepSeqArrays`
+   A.force $ A.sum (A.zipWith (*) arrRepl brrRepl)
+ where	trr@(Array _ [Region RangeAll (GenManifest _)])
+			= force $ transpose2D brr
+	arrRepl		= trr `deepSeqArray` A.extend (Z :. All   :. colsB :. All) arr
+	brrRepl		= trr `deepSeqArray` A.extend (Z :. rowsA :. All   :. All) trr
+	(Z :. _     :. rowsA) = extent arr
+	(Z :. colsB :. _    ) = extent brr
+	
+
+transpose2D :: Elt e => Array DIM2 e -> Array DIM2 e
+{-# INLINE transpose2D #-}
+transpose2D arr
+ = backpermute new_extent swap arr
+ where	swap (Z :. i :. j)	= Z :. j :. i
+	new_extent		= swap (extent arr)
diff --git a/Data/Array/Repa/Algorithms/Randomish.hs b/Data/Array/Repa/Algorithms/Randomish.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Repa/Algorithms/Randomish.hs
@@ -0,0 +1,115 @@
+{-# LANGUAGE BangPatterns #-}
+
+module Data.Array.Repa.Algorithms.Randomish
+ 	( randomishIntArray
+	, randomishIntVector
+	, randomishDoubleArray
+	, randomishDoubleVector)
+where
+import Data.Word
+import Data.Vector.Unboxed			(Vector)
+import Data.Array.Repa				as R
+import qualified Data.Vector.Unboxed.Mutable	as MV
+import qualified Data.Vector.Unboxed		as V
+import qualified Data.Vector.Generic		as G
+
+
+-- | Use the ''minimal standard'' Lehmer generator to quickly generate some random
+--   numbers with reasonable statistical properties. By ''reasonable'' we mean good
+--   enough for games and test data, but not cryptography or anything where the
+--   quality of the randomness really matters. 
+--
+--   By nature of the algorithm, the maximum value in the output is clipped
+--   to (valMin + 2^31 - 1)
+-- 
+--   From ''Random Number Generators: Good ones are hard to find''
+--   Stephen K. Park and Keith W. Miller.
+--   Communications of the ACM, Oct 1988, Volume 31, Number 10.
+--
+randomishIntArray
+	:: Shape sh
+	=> sh 			-- ^ Shape of array
+	-> Int 			-- ^ Minumum value in output.
+	-> Int 			-- ^ Maximum value in output.
+	-> Int 			-- ^ Random seed.	
+	-> Array sh Int		-- ^ Array of randomish numbers.
+
+randomishIntArray !sh !valMin !valMax !seed
+	= fromVector sh $ randomishIntVector (R.size sh) valMin valMax seed
+
+
+randomishIntVector 
+	:: Int 			-- ^ Length of vector.
+	-> Int 			-- ^ Minumum value in output.
+	-> Int 			-- ^ Maximum value in output.
+	-> Int 			-- ^ Random seed.	
+	-> Vector Int		-- ^ Vector of randomish numbers.
+
+randomishIntVector !len !valMin' !valMax' !seed'
+ = let	-- a magic number
+	-- (don't change it, the randomness depends on this specific number).
+	multiplier :: Word64
+	multiplier = 16807
+
+	-- a merzenne prime
+	-- (don't change it, the randomness depends on this specific number).
+	modulus	:: Word64
+	modulus	= 2^(31 :: Integer) - 1
+
+	-- if the seed is 0 all the numbers in the sequence are also 0.
+	seed	
+	 | seed' == 0	= 1
+	 | otherwise	= seed'
+
+	!valMin	= fromIntegral valMin'
+	!valMax	= fromIntegral valMax' + 1
+	!range	= valMax - valMin
+
+	{-# INLINE f #-}
+	f x		= multiplier * x `mod` modulus
+ in G.create 
+     $ do	
+	vec		<- MV.new len
+
+	let go !ix !x 
+	  	| ix == len	= return ()
+		| otherwise
+		= do	let x'	= f x
+			MV.write vec ix $ fromIntegral $ (x `mod` range) + valMin
+			go (ix + 1) x'
+
+	go 0 (f $ f $ f $ fromIntegral seed)
+	return vec
+
+
+-- | Generate some randomish doubles with terrible statistical properties.
+--   This just takes randomish ints then scales them, so there's not much randomness in low-order bits.
+randomishDoubleArray
+	:: Shape sh
+	=> sh 			-- ^ Shape of array
+	-> Double		-- ^ Minumum value in output.
+	-> Double		-- ^ Maximum value in output.
+	-> Int 			-- ^ Random seed.	
+	-> Array sh Double	-- ^ Array of randomish numbers.
+
+randomishDoubleArray !sh !valMin !valMax !seed
+	= fromVector sh $ randomishDoubleVector (R.size sh) valMin valMax seed
+
+
+-- | Generate some randomish doubles with terrible statistical properties.
+--   This just takes randmish ints then scales them, so there's not much randomness in low-order bits.
+randomishDoubleVector
+	:: Int			-- ^ Length of vector
+	-> Double		-- ^ Minimum value in output
+	-> Double		-- ^ Maximum value in output
+	-> Int			-- ^ Random seed.
+	-> Vector Double	-- ^ Vector of randomish doubles.
+
+randomishDoubleVector !len !valMin !valMax !seed
+ = let	range	= valMax - valMin
+
+	mx	= 2^(30 :: Integer) - 1
+	mxf	= fromIntegral mx
+	ints	= randomishIntVector len 0 mx seed
+	
+   in	V.map (\n -> valMin + (fromIntegral n / mxf) * range) ints
diff --git a/repa-algorithms.cabal b/repa-algorithms.cabal
--- a/repa-algorithms.cabal
+++ b/repa-algorithms.cabal
@@ -1,5 +1,5 @@
 Name:                repa-algorithms
-Version:             1.1.0.0
+Version:             2.0.0.1
 License:             BSD3
 License-file:        LICENSE
 Author:              The DPH Team
@@ -11,27 +11,33 @@
 Homepage:            http://trac.haskell.org/repa
 Bug-reports:         http://trac.haskell.org/repa/newticket
 Description:
-        NOTE: You must use the GHC head branch > 6.13.20100309 to get decent performance.
         Reusable algorithms using the Repa array library.
 
 Synopsis:
         Algorithms using the Repa array library.
 
-Tested-with: GHC == 6.13.20100309, GHC == 6.12.1
-
 Library
   Build-Depends: 
         base                 == 4.*,
-        dph-base             == 0.4.*,
-        repa                 == 1.1.*
+        vector               == 0.7.*,
+        repa                 == 2.0.*
 
   ghc-options:
-        -Odph -Wall -fno-warn-missing-signatures
+        -Wall -fno-warn-missing-signatures
+        -Odph
+        -fsimplifier-phases=4
+        -fstrictness-before=5
+        -funfolding-use-threshold=30
+        -funbox-strict-fields
+        -fcpr-off
 
   Exposed-modules:
         Data.Array.Repa.Algorithms.Complex
+        Data.Array.Repa.Algorithms.Randomish
         Data.Array.Repa.Algorithms.DFT
         Data.Array.Repa.Algorithms.DFT.Roots
         Data.Array.Repa.Algorithms.DFT.Center
         Data.Array.Repa.Algorithms.FFT
         Data.Array.Repa.Algorithms.Matrix
+        Data.Array.Repa.Algorithms.Convolve
+        Data.Array.Repa.Algorithms.Iterate