repa-algorithms 3.4.1.1 → 3.4.1.2
raw patch · 3 files changed
+42/−37 lines, 3 filesdep ~basenew-uploader
Dependency ranges changed: base
Files
- Data/Array/Repa/Algorithms/FFT.hs +36/−33
- Data/Array/Repa/Algorithms/Matrix.hs +1/−1
- repa-algorithms.cabal +5/−3
Data/Array/Repa/Algorithms/FFT.hs view
@@ -1,8 +1,8 @@ {-# 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. +-- | Fast computation of Discrete Fourier Transforms using the Cooley-Tuckey algorithm.+-- Time complexity is O(n log n) in the size of the input. -- -- This uses a naive divide-and-conquer algorithm, the absolute performance is about -- 50x slower than FFTW in estimate mode.@@ -14,10 +14,13 @@ , fft2dP , fft1dP) where+ import Data.Array.Repa.Algorithms.Complex import Data.Array.Repa as R import Data.Array.Repa.Eval as R import Data.Array.Repa.Unsafe as R++import Data.Bits ((.&.)) import Prelude as P @@ -37,13 +40,14 @@ {-# INLINE signOfMode #-} --- | Check if an `Int` is a power of two.+-- | Check if an `Int` is a power of two. Assumes `n` is a natural number.++-- The implementation can be found in Henry S. Warren, Jr.'s book+-- Hacker's delight, Chapter 2. isPowerOfTwo :: Int -> Bool-isPowerOfTwo n- | 0 <- n = True- | 2 <- n = True- | n `mod` 2 == 0 = isPowerOfTwo (n `div` 2)- | otherwise = False+isPowerOfTwo 0 = True+isPowerOfTwo 1 = False+isPowerOfTwo n = (n .&. (n-1)) == 0 {-# INLINE isPowerOfTwo #-} @@ -56,29 +60,29 @@ fft3dP mode arr = let _ :. depth :. height :. width = extent arr !sign = signOfMode mode- !scale = fromIntegral (depth * width * height) - + !scale = fromIntegral (depth * width * height)+ in if not (isPowerOfTwo depth && isPowerOfTwo height && isPowerOfTwo width) then error $ unlines [ "Data.Array.Repa.Algorithms.FFT: fft3d" , " Array dimensions must be powers of two,"- , " but the provided array is " + , " but the provided array is " P.++ show height P.++ "x" P.++ show width P.++ "x" P.++ show depth ]- - else arr `deepSeqArray` ++ else arr `deepSeqArray` case mode of Forward -> now $ fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr Reverse -> now $ fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr Inverse -> computeP- $ R.map (/ scale) + $ R.map (/ scale) $ fftTrans3d sign $ fftTrans3d sign $ fftTrans3d sign arr {-# INLINE fft3dP #-} -fftTrans3d +fftTrans3d :: Source r Complex => Double- -> Array r DIM3 Complex + -> Array r DIM3 Complex -> Array U DIM3 Complex fftTrans3d sign arr@@ -87,7 +91,7 @@ {-# INLINE fftTrans3d #-} -rotate3d +rotate3d :: Source r Complex => Array r DIM3 Complex -> Array D DIM3 Complex rotate3d arr@@ -107,15 +111,15 @@ fft2dP mode arr = let _ :. height :. width = extent arr sign = signOfMode mode- scale = fromIntegral (width * height) - + scale = fromIntegral (width * height)+ in if not (isPowerOfTwo height && isPowerOfTwo width) then error $ unlines [ "Data.Array.Repa.Algorithms.FFT: fft2d" , " Array dimensions must be powers of two," , " but the provided array is " P.++ show height P.++ "x" P.++ show width ]- - else arr `deepSeqArray` ++ else arr `deepSeqArray` case mode of Forward -> now $ fftTrans2d sign $ fftTrans2d sign arr Reverse -> now $ fftTrans2d sign $ fftTrans2d sign arr@@ -126,7 +130,7 @@ fftTrans2d :: Source r Complex => Double- -> Array r DIM2 Complex + -> Array r DIM2 Complex -> Array U DIM2 Complex fftTrans2d sign arr@@ -138,20 +142,20 @@ -- Vector Transform ------------------------------------------------------------------------------- -- | Compute the DFT of a vector. Array dimensions must be powers of two else `error`. fft1dP :: (Source r Complex, Monad m)- => Mode - -> Array r DIM1 Complex + => Mode+ -> Array r DIM1 Complex -> m (Array U DIM1 Complex) fft1dP mode arr = let _ :. len = extent arr sign = signOfMode mode scale = fromIntegral len- + in if not $ isPowerOfTwo len- then error $ unlines + then error $ unlines [ "Data.Array.Repa.Algorithms.FFT: fft1d" , " Array dimensions must be powers of two, " , " but the provided array is " P.++ show len ]- + else arr `deepSeqArray` case mode of Forward -> now $ fftTrans1d sign arr@@ -162,7 +166,7 @@ fftTrans1d :: Source r Complex- => Double + => Double -> Array r DIM1 Complex -> Array U DIM1 Complex @@ -174,7 +178,7 @@ -- Rank Generalised Worker ------------------------------------------------------------------------ fft :: (Shape sh, Source r Complex)- => Double -> sh -> Int + => Double -> sh -> Int -> Array r (sh :. Int) Complex -> Array U (sh :. Int) Complex @@ -183,9 +187,9 @@ where go !len !offset !stride | len == 2 = suspendedComputeP $ fromFunction (sh :. 2) swivel- + | otherwise- = combine len + = combine len (go (len `div` 2) offset (stride * 2)) (go (len `div` 2) (offset + stride) (stride * 2)) @@ -197,7 +201,7 @@ {-# INLINE combine #-} combine !len' evens odds = evens `deepSeqArray` odds `deepSeqArray`- let odds' = unsafeTraverse odds id (\get ix@(_ :. k) -> twiddle sign k len' * get ix) + let odds' = unsafeTraverse odds id (\get ix@(_ :. k) -> twiddle sign k len' * get ix) in suspendedComputeP $ (evens +^ odds') R.++ (evens -^ odds') {-# INLINE fft #-} @@ -213,4 +217,3 @@ where k = fromIntegral k' n = fromIntegral n' {-# INLINE twiddle #-}-
Data/Array/Repa/Algorithms/Matrix.hs view
@@ -56,7 +56,7 @@ do trr <- transpose2P brr let (Z :. h1 :. _) = extent arr let (Z :. _ :. w2) = extent brr- computeP + trr `deepSeqArray` computeP $ fromFunction (Z :. h1 :. w2) $ \ix -> R.sumAllS $ R.zipWith (*)
repa-algorithms.cabal view
@@ -1,5 +1,5 @@ Name: repa-algorithms-Version: 3.4.1.1+Version: 3.4.1.2 License: BSD3 License-file: LICENSE Author: The DPH Team@@ -18,8 +18,8 @@ Library Build-Depends:- base >= 4.8 && < 4.10- , vector == 0.11.*+ base >= 4.8 && < 4.11+ , vector >= 0.11 && < 0.13 , repa == 3.4.1.* ghc-options:@@ -61,3 +61,5 @@ Data.Array.Repa.Algorithms.Matrix Data.Array.Repa.Algorithms.Pixel Data.Array.Repa.Algorithms.Randomish++-- vim: nospell