diff --git a/Data/Array/Repa/Algorithms/Complex.hs b/Data/Array/Repa/Algorithms/Complex.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Repa/Algorithms/Complex.hs
@@ -0,0 +1,53 @@
+{-# LANGUAGE TypeOperators, TypeSynonymInstances #-}
+
+-- | Strict complex doubles.
+module Data.Array.Repa.Algorithms.Complex
+	( Complex
+	, mag
+	, arg
+	, (:*:)(..))
+where
+import 	Data.Array.Parallel.Base ((:*:)(..))
+
+-- | Strict complex doubles.
+type Complex 
+	= 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
+
+
+instance Fractional Complex where
+  (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
+	
+  fromRational x	= fromRational x :*: 0
+	
+-- | Take the magnitude of a complex number.
+mag :: Complex -> Double
+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)
+ = normaliseAngle $ atan2 im re
+
+ where 	normaliseAngle :: Double -> Double
+	normaliseAngle f
+	 | f < - pi	
+	 = normaliseAngle (f + 2 * pi)
+	
+	 | f > pi
+	 = normaliseAngle (f - 2 * pi)
+
+	 | otherwise
+	 = f
diff --git a/Data/Array/Repa/Algorithms/DFT.hs b/Data/Array/Repa/Algorithms/DFT.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Repa/Algorithms/DFT.hs
@@ -0,0 +1,100 @@
+{-# LANGUAGE TypeOperators, RankNTypes, PatternGuards #-}
+
+-- | Compute the Discrete Fourier Transform (DFT) along the low order dimension
+--   of an array. 
+--
+--   This uses the naive algorithm and takes O(n^2) time. 
+--   However, you can transform an array with an arbitray extent, unlike with FFT which requires
+--   each dimension to be a power of two.
+--
+--   The `dft` and `idft` functions 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 `dftWithRoots` directly.
+--
+--   You can also compute single values of the transform using `dftWithRootsSingle`.
+module Data.Array.Repa.Algorithms.DFT 
+	( dft
+	, idft
+	, dftWithRoots
+	, dftWithRootsSingle)
+where
+import Data.Array.Repa.Algorithms.DFT.Roots
+import Data.Array.Repa.Algorithms.Complex
+import Data.Array.Repa				as A
+
+-- | Compute the DFT along the low order dimension of an array.
+dft 	:: forall sh
+	.  Shape sh
+	=> Array (sh :. Int) Complex
+	-> Array (sh :. Int) Complex
+
+dft v
+ = let	rofu	= calcRootsOfUnity (extent v)
+   in	force $ dftWithRoots rofu v
+
+
+-- | Compute the inverse DFT along the low order dimension of an array.
+idft 	:: forall sh
+	.  Shape sh
+	=> Array (sh :. Int) Complex
+	-> Array (sh :. Int) Complex
+
+idft v
+ = let	_ :. len	= extent v
+	scale		= fromIntegral len :*: 0
+	rofu		= calcInverseRootsOfUnity (extent v)
+   in	force $ A.map (/ scale) $ dftWithRoots rofu v
+
+
+-- | Generic function for computation of forward or inverse DFT.
+--	This function is also useful if you transform many arrays with the same extent, 
+--	and don't want to recompute the roots for each one.
+--	The extent of the given roots must match that of the input array, else `error`.
+dftWithRoots
+	:: forall sh
+	.  Shape sh
+	=> Array (sh :. Int) Complex		-- ^ Roots of unity.
+	-> Array (sh :. Int) Complex		-- ^ Input array.
+	-> Array (sh :. Int) Complex
+
+dftWithRoots rofu arr
+	| _ :. 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 ++ ")"
+
+	| otherwise
+	= traverse arr id (\_ k -> dftWithRootsSingle rofu arr k)
+		
+
+-- | Compute a single value of the DFT.
+--	The extent of the given roots must match that of the input array, else `error`.
+dftWithRootsSingle
+	:: forall sh
+	.  Shape sh
+	=> Array (sh :. Int) Complex 		-- ^ Roots of unity.
+	-> Array (sh :. Int) Complex		-- ^ Input array.
+	-> (sh :. Int)				-- ^ Index of the value we want.
+	-> Complex
+
+{-# INLINE dftWithRootsSingle #-}
+dftWithRootsSingle rofu arrX (_ :. k)
+	| _ :. 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 ++ ")"
+
+	| otherwise
+	= let	sh@(_ :. len)	= extent arrX
+
+		-- All the roots we need to multiply with.
+		wroots		= fromFunction sh elemFn
+		elemFn (sh' :. n) 
+			= 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
new file mode 100644
--- /dev/null
+++ b/Data/Array/Repa/Algorithms/DFT/Center.hs
@@ -0,0 +1,31 @@
+
+-- | 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)
+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
+ = 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
+ = traverse arr id
+	(\get ix@(_ :. y :. x) -> ((-1) ^ (y + x)) * get ix)
diff --git a/Data/Array/Repa/Algorithms/DFT/Roots.hs b/Data/Array/Repa/Algorithms/DFT/Roots.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Repa/Algorithms/DFT/Roots.hs
@@ -0,0 +1,42 @@
+{-# LANGUAGE TypeOperators, RankNTypes #-}
+
+-- | Calculation of roots of unity for the forward and inverse DFT\/FFT.
+module Data.Array.Repa.Algorithms.DFT.Roots
+	( calcRootsOfUnity
+	, calcInverseRootsOfUnity)
+where
+import Data.Array.Repa
+import Data.Array.Repa.Algorithms.Complex
+
+-- | Calculate roots of unity for the forward transform.
+calcRootsOfUnity
+	:: forall sh
+	.  Shape sh
+	=> (sh :. Int) 			-- ^ Length of lowest dimension of result.
+	-> Array (sh :. Int) Complex
+
+calcRootsOfUnity sh@(_ :. n) 
+ = 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))
+
+    len	= fromIntegral n
+
+
+-- | Calculate roots of unity for the inverse transform.
+calcInverseRootsOfUnity
+	:: forall sh
+	.  Shape sh
+	=> (sh :. Int) 			-- ^ Length of lowest dimension of result.
+	-> Array (sh :. Int) Complex
+
+calcInverseRootsOfUnity sh@(_ :. n) 
+ = 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))
+
+    len	= fromIntegral n
diff --git a/Data/Array/Repa/Algorithms/FFT.hs b/Data/Array/Repa/Algorithms/FFT.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Repa/Algorithms/FFT.hs
@@ -0,0 +1,212 @@
+{-# LANGUAGE TypeOperators, PatternGuards, RankNTypes #-}
+
+-- | 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.
+--
+--   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.
+--
+--   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 )
+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
+
+fft v
+ = let	rofu	= calcRootsOfUnity (extent v)
+   in	force $ fftWithRoots rofu v
+
+
+-- | Compute the inverse DFT along the low order dimension of an array.
+ifft	:: Shape sh
+	=> Array (sh :. Int) Complex
+	-> Array (sh :. Int) Complex
+
+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  ++ ")"
+
+	| 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
+	
+
+-- Cube Transform ---------------------------------------------------------------------------------
+-- | Compute the DFT of a 3d cube.
+--   If the array is not a cube then `error`.
+fft3d 	:: 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"
+
+	| 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)) 
+
+		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
+
+ifft3d arrIn
+ 	| Z :. depth :. height :. width	<- extent arrIn
+ 	, (height /= width) || (height /= depth)
+	= error $ "ifft3d: array is not a cube"
+
+	| otherwise
+	= let	rofu		= calcInverseRootsOfUnity (extent arrIn)
+
+		transpose3 arr
+	 	 = traverse arr 
+        		(\(Z :. k :. l :. m)   -> (Z :. l :. m :. k)) 
+            		(\f (Z :. l :. m :. k) -> f (Z :. k :. l :. m)) 
+
+		_ :. 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"
+	
+	| rLen /= vLen
+	= error $  "fft: length of vector (" ++ show vLen ++ ")"
+		++ " does not match the length of the roots (" ++ show rLen ++ ")"
+	
+	| otherwise
+	= fftWithRoots' rofu v
+
+	where	_ :. rLen	= extent rofu
+		_ :. vLen	= extent v
+
+fftWithRoots'
+	:: Shape sh
+	=> Array (sh :. Int) Complex
+	-> Array (sh :. Int) Complex
+        -> Array (sh :. Int) Complex
+
+{-# 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'
+	
+{-# INLINE fft_split #-}
+fft_split rofu v
+ = let 	fft_lr = force $ fftWithRoots' (splitRofu rofu) (splitVector v)
+
+	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))
+
+	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))
+
+   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"
+
+   in	traverse v
+		(\(vSh :. vLen)    -> vSh :. 2 :. (vLen `div` 2)) 
+		vFn'
+        
diff --git a/Data/Array/Repa/Algorithms/Matrix.hs b/Data/Array/Repa/Algorithms/Matrix.hs
new file mode 100644
--- /dev/null
+++ b/Data/Array/Repa/Algorithms/Matrix.hs
@@ -0,0 +1,30 @@
+{-# OPTIONS -fno-warn-incomplete-patterns #-}
+
+-- | Algorithms operating on matrices.
+-- 
+--   These functions should give performance comparable with nested loop C
+--   implementations, but not block-based, cache friendly, SIMD using, vendor
+--   optimised implementions. 
+--   If you care deeply about runtime performance then you may be better off using 
+--   a binding to LAPACK, such as hvector.
+--
+module Data.Array.Repa.Algorithms.Matrix
+	(multiplyMM)
+where
+import Data.Array.Repa
+	
+
+-- | Matrix-matrix multiply.
+--
+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))
diff --git a/LICENSE b/LICENSE
new file mode 100644
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,24 @@
+Copyright (c) 2010, University of New South Wales.
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    * Neither the name of the University of New South Wales nor the
+      names of its contributors may be used to endorse or promote products
+      derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ''AS IS'' AND ANY
+EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL COPYRIGHT HOLDERS BE LIABLE FOR ANY
+DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
diff --git a/Setup.hs b/Setup.hs
new file mode 100644
--- /dev/null
+++ b/Setup.hs
@@ -0,0 +1,2 @@
+import Distribution.Simple
+main = defaultMain
diff --git a/repa-algorithms.cabal b/repa-algorithms.cabal
new file mode 100644
--- /dev/null
+++ b/repa-algorithms.cabal
@@ -0,0 +1,37 @@
+Name:                repa-algorithms
+Version:             1.1.0.0
+License:             BSD3
+License-file:        LICENSE
+Author:              The DPH Team
+Maintainer:          Ben Lippmeier <benl@ouroborus.net>
+Build-Type:          Simple
+Cabal-Version:       >=1.6
+Stability:           experimental
+Category:            Data Structures
+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.*
+
+  ghc-options:
+        -Odph -Wall -fno-warn-missing-signatures
+
+  Exposed-modules:
+        Data.Array.Repa.Algorithms.Complex
+        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