diff --git a/Data/Array/Accelerate/Math/Complex.hs b/Data/Array/Accelerate/Math/Complex.hs --- a/Data/Array/Accelerate/Math/Complex.hs +++ b/Data/Array/Accelerate/Math/Complex.hs @@ -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) 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 @@ -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' 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 @@ -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 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 @@ -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 )) 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 @@ -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 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.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