accelerate-fft 0.13.0.0 → 0.14.0.0

raw patch · 6 files changed

+173/−108 lines, 6 filesdep ~acceleratedep ~accelerate-cudadep ~basePVP ok

version bump matches the API change (PVP)

Dependency ranges changed: accelerate, accelerate-cuda, base

API changes (from Hackage documentation)

Files

• Data/Array/Accelerate/Math/Complex.hs +98/−35
• Data/Array/Accelerate/Math/DFT.hs +3/−3
• Data/Array/Accelerate/Math/DFT/Centre.hs +3/−3
• Data/Array/Accelerate/Math/DFT/Roots.hs +2/−2
• Data/Array/Accelerate/Math/FFT.hs +60/−58
• accelerate-fft.cabal +7/−7

@@ -1,54 +1,113 @@-{-# LANGUAGE FlexibleInstances    #-}-{-# LANGUAGE IncoherentInstances  #-}-{-# LANGUAGE ScopedTypeVariables  #-}-{-# LANGUAGE TypeSynonymInstances #-}+{-# LANGUAGE FlexibleContexts      #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE IncoherentInstances   #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE ScopedTypeVariables   #-}+{-# LANGUAGE TypeFamilies          #-}+{-# LANGUAGE TypeSynonymInstances  #-} {-# OPTIONS -fno-warn-orphans #-} -module Data.Array.Accelerate.Math.Complex-  where+module Data.Array.Accelerate.Math.Complex ( +  Complex(..),+  magnitude, phase, real, imag, conj,++) where+ import Prelude-import Data.Function-import Data.Array.Accelerate                    as A+import Data.Complex                             ( Complex(..) )+import Data.Array.Accelerate+import Data.Array.Accelerate.Smart+import Data.Array.Accelerate.Tuple+import Data.Array.Accelerate.Array.Sugar -type Complex a = (a, a) -instance (Elt a, IsFloating a) => Num (Exp (Complex a)) where-  c1 + c2       = lift ( on (+) real c1 c2, on (+) imag c1 c2 )-  c1 - c2       = lift ( on (-) real c1 c2, on (-) imag c1 c2 )-  c1 * c2       = let (x,  y)   = unlift c1-                      (x', y')  = unlift c2     :: Complex (Exp a)-                  in lift (x*x'-y*y', x*y'+y*x')+type instance EltRepr  (Complex a) = (EltRepr a, EltRepr' a)+type instance EltRepr' (Complex a) = (EltRepr a, EltRepr' a) -  negate c      = lift ( negate (real c), negate (imag c) )-  abs z         = lift ( magnitude z, constant 0 )-  signum z      = let r         = magnitude z-                      (x, y)    = unlift z-                  in r ==* 0 ? (constant (0,0), lift (x/r, y/r))+instance Elt a => Elt (Complex a) where+  eltType (_::Complex a)        = eltType (undefined :: (a,a))+  toElt (a,b)                   = toElt a :+ toElt' b+  fromElt (a :+ b)              = (fromElt a, fromElt' b) -  fromInteger n-    = lift (constant (fromInteger n), constant 0)+  eltType' (_::Complex a)       = eltType' (undefined :: (a,a))+  toElt' (a,b)                  = toElt a :+ toElt' b+  fromElt' (a :+ b)             = (fromElt a, fromElt' b) +instance IsTuple (Complex a) where+  type TupleRepr (Complex a) = (((), a), a)+  fromTuple (x :+ y)    = (((), x), y)+  toTuple (((), x), y)  = (x :+ y) +instance (Lift Exp a, Elt (Plain a)) => Lift Exp (Complex a) where+  type Plain (Complex a) = Complex (Plain a)+  lift (x1 :+ x2)       = Exp $ Tuple (NilTup `SnocTup` lift x1 `SnocTup` lift x2)++instance Elt a => Unlift Exp (Complex (Exp a)) where+  unlift e+    = let x     = Exp $ SuccTupIdx ZeroTupIdx `Prj` e+          y     = Exp $ ZeroTupIdx `Prj` e+      in+      x :+ y++instance (Elt a, IsFloating a) => Num (Exp (Complex a)) where+  (+)           = lift2 ((+) :: Complex (Exp a) -> Complex (Exp a) -> Complex (Exp a))+  (-)           = lift2 ((-) :: Complex (Exp a) -> Complex (Exp a) -> Complex (Exp a))+  (*)           = lift2 ((*) :: Complex (Exp a) -> Complex (Exp a) -> Complex (Exp a))+  negate        = lift1 (negate :: Complex (Exp a) -> Complex (Exp a))+  signum        = lift1 (signum :: Complex (Exp a) -> Complex (Exp a))+  abs           = lift1 (abs :: Complex (Exp a) -> Complex (Exp a))+  fromInteger n = lift (constant (fromInteger n) :+ 0)++ instance (Elt a, IsFloating a) => Fractional (Exp (Complex a)) where-  c1 / c2-    = let (a,b) = unlift c1-          (c,d) = unlift c2     :: Complex (Exp a)-          den   = c^(2 :: Int) + d^(2 :: Int)-          re    = (a * c + b * d) / den-          im    = (b * c - a * d) / den+  c / c'+    = let x  :+ y       = unlift c+          x' :+ y'      = unlift c'     :: Complex (Exp a)+          den           = x'^(2 :: Int) + y'^(2 :: Int)+          re            = (x * x' + y * y') / den+          im            = (y * x' - x * y') / den       in-      lift (re, im)+      lift (re :+ im)    fromRational x-    = lift (constant (fromRational x), constant 0)+    = lift (constant (fromRational x) :+ constant 0)  +instance (Elt a, IsFloating a, RealFloat a) => Floating (Exp (Complex a)) where+  sqrt z+    = let+          x :+ y        = unlift z+          v'            = abs y / (u'*2)+          u'            = sqrt ((magnitude z + abs x) / 2)+          (u, v)        = unlift ( x <* 0 ? ( lift (v',u'), lift (u',v') ) )+      in+      x ==* 0 &&* y ==* 0 ?+        {- then -} ( 0+        {- else -} , lift (u :+ (y <* 0 ? (-v,v))) )++  pi            = lift (pi :+ constant 0)+  log z         = lift (log (magnitude z) :+ phase z)+  exp           = lift1 (exp :: Complex (Exp a) -> Complex (Exp a))+  sin           = lift1 (sin :: Complex (Exp a) -> Complex (Exp a))+  cos           = lift1 (cos :: Complex (Exp a) -> Complex (Exp a))+  tan           = lift1 (tan :: Complex (Exp a) -> Complex (Exp a))+  sinh          = lift1 (sinh :: Complex (Exp a) -> Complex (Exp a))+  cosh          = lift1 (cosh :: Complex (Exp a) -> Complex (Exp a))+  tanh          = lift1 (tanh :: Complex (Exp a) -> Complex (Exp a))+  asin          = lift1 (asin :: Complex (Exp a) -> Complex (Exp a))+  acos          = lift1 (acos :: Complex (Exp a) -> Complex (Exp a))+  atan          = lift1 (atan :: Complex (Exp a) -> Complex (Exp a))+  asinh         = lift1 (asinh :: Complex (Exp a) -> Complex (Exp a))+  acosh         = lift1 (acosh :: Complex (Exp a) -> Complex (Exp a))+  atanh         = lift1 (atanh :: Complex (Exp a) -> Complex (Exp a))++ -- | Non-negative magnitude of a complex number -- magnitude :: (Elt a, IsFloating a) => Exp (Complex a) -> Exp a magnitude c =-  let (r, i) = unlift c+  let r :+ i    = unlift c   in sqrt (r*r + i*i)  -- | The phase of a complex number, in the range (-pi, pi]. If the magnitude is@@ -56,24 +115,28 @@ -- phase :: (Elt a, IsFloating a) => Exp (Complex a) -> Exp a phase c =-  let (x, y) = unlift c+  let x :+ y    = unlift c   in atan2 y x   -- | Return the real part of a complex number -- real :: Elt a => Exp (Complex a) -> Exp a-real = A.fst+real c =+  let r :+ _    = unlift c+  in  r  -- | Return the imaginary part of a complex number -- imag :: Elt a => Exp (Complex a) -> Exp a-imag = A.snd+imag c =+  let _ :+ i    = unlift c+  in  i  -- | Return the complex conjugate of a complex number, defined as -- -- > conj(Z) = X - iY -- conj :: (Elt a, IsNum a) => Exp (Complex a) -> Exp (Complex a)-conj z = lift (real z, - imag z)+conj z = lift $ real z :+ (- imag z) 

@@ -52,7 +52,7 @@   = let sh      = shape v         n       = indexHead sh         roots   = inverseRootsOfUnity sh-        scale   = lift (A.fromIntegral n, constant 0)+        scale   = lift (A.fromIntegral n :+ constant 0)     in     A.map (/scale) $ dftG roots v @@ -68,7 +68,7 @@      -> Acc (Array (sh:.Int) (Complex e))       -- ^ input array      -> Acc (Array (sh:.Int) (Complex e)) dftG roots arr-  = A.fold (+) (constant (0,0))+  = A.fold (+) (constant (0 :+ 0))   $ A.zipWith (*) arr' roots'   where     base        = shape arr@@ -106,5 +106,5 @@                              (\ix' -> let sh :. n = unlift ix'  :: Exp sh :. Exp Int                                       in  roots ! lift (sh :. (k*n) `mod` l))     in-    A.foldAll (+) (constant (0,0)) $ A.zipWith (*) arr roots'+    A.foldAll (+) (constant (0 :+ 0)) $ A.zipWith (*) arr roots' 

@@ -36,7 +36,7 @@ centre1D arr   = A.generate (shape arr)                (\ix -> let Z :. x = unlift ix           :: Z :. Exp Int-                       in  lift (-1 ** A.fromIntegral x, A.constant 0) * arr!ix)+                       in  lift (((-1) ** A.fromIntegral x) :+ A.constant 0) * arr!ix)  -- | Apply the centring transform to a matrix --@@ -46,7 +46,7 @@ centre2D arr   = A.generate (shape arr)                (\ix -> let Z :. y :. x = unlift ix      :: Z :. Exp Int :. Exp Int-                       in  lift (-1 ** A.fromIntegral (y + x), A.constant 0) * arr!ix)+                       in  lift (((-1) ** A.fromIntegral (y + x)) :+ A.constant 0) * arr!ix)  -- | Apply the centring transform to a 3D array --@@ -56,7 +56,7 @@ centre3D arr   = A.generate (shape arr)                (\ix -> let Z :. z :. y :. x = unlift ix :: Z :. Exp Int :. Exp Int :. Exp Int-                       in  lift (-1 ** A.fromIntegral (z + y + x), A.constant 0) * arr!ix)+                       in  lift (((-1) ** A.fromIntegral (z + y + x)) :+ A.constant 0) * arr!ix)   -- | Apply the shifting transform to a vector

@@ -31,7 +31,7 @@   A.generate sh (\ix -> let i = A.fromIntegral (A.indexHead ix)                             k = 2 * pi * i / n                         in-                        A.lift ( cos k, -sin k ))+                        A.lift ( cos k :+ (-sin k) ))   -- | Calculate the roots of unity for an inverse transform@@ -46,5 +46,5 @@   A.generate sh (\ix -> let i = A.fromIntegral (A.indexHead ix)                             k = 2 * pi * i / n                         in-                        A.lift ( cos k, sin k ))+                        A.lift ( cos k :+ sin k )) 

@@ -42,6 +42,7 @@ import Data.Array.Accelerate.Type  import Data.Functor+import System.IO.Unsafe import Foreign.CUDA.FFT import qualified Foreign.CUDA.Driver            as CUDA hiding (free) #endif@@ -96,7 +97,7 @@     if P.not (isPow2 len)        then error $ unlines               [ "Data.Array.Accelerate.FFT: fft1D"-              , "  Array dimensions must be powers of two, but are: " ++ showShape (Z:.len) ]+              , "  Array dimensions must be powers of two, but are: " P.++ showShape (Z:.len) ]         else case mode of                  Inverse -> A.map (/scale) vec'@@ -141,7 +142,7 @@     if P.not (isPow2 width && isPow2 height)        then error $ unlines               [ "Data.Array.Accelerate.FFT: fft2D"-              , "  Array dimensions must be powers of two, but are: " ++ showShape (Z:.height:.width) ]+              , "  Array dimensions must be powers of two, but are: " P.++ showShape (Z:.height:.width) ]         else case mode of                  Inverse -> A.map (/scale) arr'@@ -188,7 +189,7 @@     if P.not (isPow2 width && isPow2 height && isPow2 depth)        then error $ unlines               [ "Data.Array.Accelerate.FFT: fft3D"-              , "  Array dimensions must be powers of two, but are: " ++ showShape (Z:.depth:.height:.width) ]+              , "  Array dimensions must be powers of two, but are: " P.++ showShape (Z:.depth:.height:.width) ]         else case mode of                  Inverse -> A.map (/scale) arr'@@ -250,7 +251,7 @@               i = A.fromIntegral i'               k = 2*pi*i/n           in-          lift ( cos k, A.constant sign * sin k )+          lift ( cos k :+ A.constant sign * sin k )  #ifdef ACCELERATE_CUDA_BACKEND -- FFT using the CUFFT library to enable high performance for the CUDA backend of@@ -263,69 +264,70 @@         -> (Acc (Array sh (Complex e)) -> Acc (Array sh (Complex e)))         -> Acc (Array sh (Complex e))         -> Acc (Array sh (Complex e))-cudaFFT mode sh p arr = deinterleave sh (foreignAcc ff pure (interleave arr))+cudaFFT mode sh = cudaFFT'   where-    ff          = cudaAcc foreignFFT-    -- Unfortunately the pure version of the function needs to be wrapped in-    -- interleave and deinterleave to match how the foreign version works.-    ---    -- RCE: Do the interleaving and deinterleaving in foreignFFT-    ---    -- TLM: The interleaving might get fused into other parts of the-    --      computation and thus be okay. We should really support multi types-    --      such as float2 instead.+    -- Plan the FFT.+    -- Doing this in unsafePerformIO so it is not reperformed every time the+    -- AST is evaluated.     ---    pure        = interleave . p . deinterleave sh-    sign        = signOfMode mode :: Int--    foreignFFT :: Array DIM1 e -> CIO (Array DIM1 e)-    foreignFFT arr' = do-      -- Create the plan-      -- TODO: Cache this.-      ---      hndl <- liftIO $-        case shapeToList sh of-          [width]                -> plan1D              width types 1-          [height, width]        -> plan2D       height width types-          [depth, height, width] -> plan3D depth height width types-          _                      -> error "Accelerate-fft cannot use CUFFT for arrays of dimensions higher than 3"+    -- RCE: This is currently not destroyed properly. Need to attach a finaliser.+    hndl = unsafePerformIO $ do+            case shapeToList sh of+              [width]                -> plan1D              width types 1+              [height, width]        -> plan2D       height width types+              [depth, height, width] -> plan3D depth height width types+              _                      -> error "Accelerate-fft cannot use CUFFT for arrays of dimensions higher than 3" -      output <- allocateArray (S.shape arr')-      iptr   <- floatingDevicePtr arr'-      optr   <- floatingDevicePtr output+    types = case (floatingType :: FloatingType e) of+              TypeFloat{}   -> C2C+              TypeDouble{}  -> Z2Z+              TypeCFloat{}  -> C2C+              TypeCDouble{} -> Z2Z -      --Execute-      liftIO $ execute hndl iptr optr+    cudaFFT' p arr = deinterleave sh (foreignAcc ff pure (interleave arr))+      where+        ff          = CUDAForeignAcc "foreignFFT" foreignFFT+        -- Unfortunately the pure version of the function needs to be wrapped in+        -- interleave and deinterleave to match how the foreign version works.+        --+        -- RCE: Do the interleaving and deinterleaving in foreignFFT+        --+        -- TLM: The interleaving might get fused into other parts of the+        --      computation and thus be okay. We should really support multi types+        --      such as float2 instead.+        --+        pure        = interleave . p . deinterleave sh+        sign        = signOfMode mode :: Int -      liftIO $ destroy hndl+        foreignFFT :: Array DIM1 e -> CIO (Array DIM1 e)+        foreignFFT arr' = do+          output <- allocateArray (S.shape arr')+          iptr   <- floatingDevicePtr arr'+          optr   <- floatingDevicePtr output -      return output+          --Execute+          liftIO $ execute iptr optr -    types-      = case (floatingType :: FloatingType e) of-          TypeFloat{}   -> C2C-          TypeDouble{}  -> Z2Z-          TypeCFloat{}  -> C2C-          TypeCDouble{} -> Z2Z+          return output -    execute :: Handle -> CUDA.DevicePtr e -> CUDA.DevicePtr e -> IO ()-    execute hndl iptr optr-      = case (floatingType :: FloatingType e) of-          TypeFloat{}   -> execC2C hndl iptr optr sign-          TypeDouble{}  -> execZ2Z hndl iptr optr sign-          TypeCFloat{}  -> execC2C hndl (CUDA.castDevPtr iptr) (CUDA.castDevPtr optr) sign-          TypeCDouble{} -> execZ2Z hndl (CUDA.castDevPtr iptr) (CUDA.castDevPtr optr) sign+        execute :: CUDA.DevicePtr e -> CUDA.DevicePtr e -> IO ()+        execute iptr optr+          = case (floatingType :: FloatingType e) of+              TypeFloat{}   -> execC2C hndl iptr optr sign+              TypeDouble{}  -> execZ2Z hndl iptr optr sign+              TypeCFloat{}  -> execC2C hndl (CUDA.castDevPtr iptr) (CUDA.castDevPtr optr) sign+              TypeCDouble{} -> execZ2Z hndl (CUDA.castDevPtr iptr) (CUDA.castDevPtr optr) sign -    floatingDevicePtr :: Vector e -> CIO (CUDA.DevicePtr e)-    floatingDevicePtr v-      = case (floatingType :: FloatingType e) of-          TypeFloat{}   -> singleDevicePtr v-          TypeDouble{}  -> singleDevicePtr v-          TypeCFloat{}  -> CUDA.castDevPtr <$> singleDevicePtr v-          TypeCDouble{} -> CUDA.castDevPtr <$> singleDevicePtr v+        floatingDevicePtr :: Vector e -> CIO (CUDA.DevicePtr e)+        floatingDevicePtr v+          = case (floatingType :: FloatingType e) of+              TypeFloat{}   -> singleDevicePtr v+              TypeDouble{}  -> singleDevicePtr v+              TypeCFloat{}  -> CUDA.castDevPtr <$> singleDevicePtr v+              TypeCDouble{} -> CUDA.castDevPtr <$> singleDevicePtr v -    singleDevicePtr :: DevicePtrs (EltRepr e) ~ ((),CUDA.DevicePtr b) => Vector e -> CIO (CUDA.DevicePtr b)-    singleDevicePtr v = P.snd <$> devicePtrsOfArray v+        singleDevicePtr :: DevicePtrs (EltRepr e) ~ ((),CUDA.DevicePtr b) => Vector e -> CIO (CUDA.DevicePtr b)+        singleDevicePtr v = P.snd <$> devicePtrsOfArray v #endif  -- Append two arrays. Doesn't do proper bounds checking or intersection...@@ -369,6 +371,6 @@ deinterleave :: (Shape sh, Elt e) => sh -> Acc (Vector e) -> Acc (Array sh (Complex e)) deinterleave (constant -> sh) arr =   generate sh (\ix -> let i = toIndex sh ix * 2-                      in  lift (arr A.!! i, arr A.!! (i+1)))+                      in  lift (arr A.!! i :+ arr A.!! (i+1))) #endif 

@@ -1,7 +1,7 @@ Name:                   accelerate-fft-Version:                0.13.0.0+Version:                0.14.0.0 Cabal-version:          >= 1.6-Tested-with:            GHC >= 7.4+Tested-with:            GHC == 7.6.* Build-type:             Simple  Synopsis:               FFT using the Accelerate library@@ -31,8 +31,8 @@   Default:              True  Library-  Build-depends:        accelerate              == 0.13.*,-                        base                    == 4.*+  Build-depends:        accelerate              == 0.14.*,+                        base                    == 4.6.*    Exposed-modules:      Data.Array.Accelerate.Math.Complex                         Data.Array.Accelerate.Math.FFT@@ -44,9 +44,9 @@    if flag(cuda)     CPP-options:        -DACCELERATE_CUDA_BACKEND-    Build-depends:      accelerate-cuda         == 0.13.*,-                        cuda                    >= 0.5  && < 0.6,-                        cufft                   >= 0.1  && < 0.2+    Build-depends:      accelerate-cuda         == 0.14.*,+                        cuda                    >= 0.5          && < 0.6,+                        cufft                   >= 0.1          && < 0.2    -- Don't add the extensions list here. Instead, place individual LANGUAGE   -- pragmas in the files that require a specific extension. This means the