diff --git a/Data/Array/Accelerate/Math/Complex.hs b/Data/Array/Accelerate/Math/Complex.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Math/Complex.hs +++ /dev/null @@ -1,142 +0,0 @@ -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE IncoherentInstances #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE TypeSynonymInstances #-} -{-# OPTIONS -fno-warn-orphans #-} - -module Data.Array.Accelerate.Math.Complex ( - - Complex(..), - magnitude, phase, real, imag, conj, - -) where - -import Prelude -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 instance EltRepr (Complex a) = (EltRepr a, EltRepr' a) -type instance EltRepr' (Complex a) = (EltRepr a, EltRepr' a) - -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) - - 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 - 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) - - fromRational x - = 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 - in sqrt (r*r + i*i) - --- | The phase of a complex number, in the range (-pi, pi]. If the magnitude is --- zero, then so is the phase. --- -phase :: (Elt a, IsFloating a) => Exp (Complex a) -> Exp a -phase 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 c = - let r :+ _ = unlift c - in r - --- | Return the imaginary part of a complex number --- -imag :: Elt a => Exp (Complex a) -> Exp a -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) - diff --git a/Data/Array/Accelerate/Math/DFT.hs b/Data/Array/Accelerate/Math/DFT.hs --- a/Data/Array/Accelerate/Math/DFT.hs +++ b/Data/Array/Accelerate/Math/DFT.hs @@ -32,7 +32,7 @@ import Prelude as P hiding ((!!)) import Data.Array.Accelerate as A import Data.Array.Accelerate.Math.DFT.Roots -import Data.Array.Accelerate.Math.Complex +import Data.Array.Accelerate.Data.Complex -- | Compute the DFT along the low order dimension of an array diff --git a/Data/Array/Accelerate/Math/DFT/Centre.hs b/Data/Array/Accelerate/Math/DFT/Centre.hs --- a/Data/Array/Accelerate/Math/DFT/Centre.hs +++ b/Data/Array/Accelerate/Math/DFT/Centre.hs @@ -25,7 +25,7 @@ import Prelude as P import Data.Array.Accelerate as A -import Data.Array.Accelerate.Math.Complex +import Data.Array.Accelerate.Data.Complex -- | Apply the centring transform to a vector diff --git a/Data/Array/Accelerate/Math/DFT/Roots.hs b/Data/Array/Accelerate/Math/DFT/Roots.hs --- a/Data/Array/Accelerate/Math/DFT/Roots.hs +++ b/Data/Array/Accelerate/Math/DFT/Roots.hs @@ -16,7 +16,7 @@ import Prelude as P import Data.Array.Accelerate as A -import Data.Array.Accelerate.Math.Complex +import Data.Array.Accelerate.Data.Complex -- | Calculate the roots of unity for the forward transform diff --git a/Data/Array/Accelerate/Math/FFT.hs b/Data/Array/Accelerate/Math/FFT.hs --- a/Data/Array/Accelerate/Math/FFT.hs +++ b/Data/Array/Accelerate/Math/FFT.hs @@ -34,7 +34,7 @@ import Prelude as P import Data.Array.Accelerate as A import Data.Array.Accelerate.Array.Sugar ( showShape ) -import Data.Array.Accelerate.Math.Complex +import Data.Array.Accelerate.Data.Complex #ifdef ACCELERATE_CUDA_BACKEND import Data.Array.Accelerate.CUDA.Foreign @@ -42,6 +42,7 @@ import Data.Array.Accelerate.Type import Data.Functor +import System.Mem.Weak import System.IO.Unsafe import Foreign.CUDA.FFT import qualified Foreign.CUDA.Driver as CUDA hiding (free) @@ -130,7 +131,7 @@ = let sign = signOfMode mode :: e scale = P.fromIntegral (width * height) #ifdef ACCELERATE_CUDA_BACKEND - sh = (Z:.width:.height) + sh = (Z:.height:.width) arr' = cudaFFT mode sh fft' arr #else arr' = fft' arr @@ -176,7 +177,7 @@ = let sign = signOfMode mode :: e scale = P.fromIntegral (width * height) #ifdef ACCELERATE_CUDA_BACKEND - sh = (Z:.width:.height:.depth) + sh = (Z:.depth:.height:.width) arr' = cudaFFT mode sh fft' arr #else arr' = fft' arr @@ -270,13 +271,14 @@ -- Doing this in unsafePerformIO so it is not reperformed every time the -- AST is evaluated. -- - -- 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" + plan <- case shapeToList sh of + [width] -> plan1D width types 1 + [width, height] -> plan2D height width types + [width, height, depth] -> plan3D depth height width types + _ -> error "Accelerate-fft cannot use CUFFT for arrays of dimensions higher than 3" + addFinalizer plan (destroy plan) + return plan types = case (floatingType :: FloatingType e) of TypeFloat{} -> C2C @@ -356,6 +358,7 @@ -- flattened vector. This allows us to mimic the float2 structure used by CUFFT -- to store complex numbers. -- +{-# NOINLINE interleave #-} interleave :: (Shape sh, Elt e) => Acc (Array sh (Complex e)) -> Acc (Vector e) interleave arr = generate sh swizzle where @@ -368,9 +371,9 @@ -- Deinterleave a vector into a complex array. Assumes the array is even in length. -- +{-# NOINLINE deinterleave #-} 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))) #endif - diff --git a/accelerate-fft.cabal b/accelerate-fft.cabal --- a/accelerate-fft.cabal +++ b/accelerate-fft.cabal @@ -1,7 +1,7 @@ Name: accelerate-fft -Version: 0.14.0.0 +Version: 0.15.0.0 Cabal-version: >= 1.6 -Tested-with: GHC == 7.6.* +Tested-with: GHC >= 7.6 Build-type: Simple Synopsis: FFT using the Accelerate library @@ -31,11 +31,10 @@ Default: True Library - Build-depends: accelerate == 0.14.*, - base == 4.6.* + Build-depends: accelerate == 0.15.*, + base >= 4.6 && < 4.8 - Exposed-modules: Data.Array.Accelerate.Math.Complex - Data.Array.Accelerate.Math.FFT + Exposed-modules: Data.Array.Accelerate.Math.FFT Data.Array.Accelerate.Math.DFT Data.Array.Accelerate.Math.DFT.Centre Data.Array.Accelerate.Math.DFT.Roots @@ -44,9 +43,9 @@ if flag(cuda) CPP-options: -DACCELERATE_CUDA_BACKEND - Build-depends: accelerate-cuda == 0.14.*, - cuda >= 0.5 && < 0.6, - cufft >= 0.1 && < 0.2 + Build-depends: accelerate-cuda == 0.15.*, + cuda >= 0.5, + cufft >= 0.1 -- Don't add the extensions list here. Instead, place individual LANGUAGE -- pragmas in the files that require a specific extension. This means the