diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,12 +7,22 @@ Policy (PVP)](https://pvp.haskell.org) +## [0.2.0.0] - 2018-04-03 +### Changed + * Update for AoS representation of complex numbers + +### Added + * support for LLVM-6.0 + ## [0.1.0.1] - 2017-09-25 +### Fixed * test-suite: build fix for ghc-8.2 ## [0.1.0.0] - 2017-09-21 * First version. Released on an unsuspecting world. + +[0.2.0.0]: https://github.com/tmcdonell/accelerate-blas/compare/0.1.0.1...0.2.0.0 [0.1.0.1]: https://github.com/tmcdonell/accelerate-blas/compare/0.1.0.0...0.1.0.1 [0.1.0.0]: https://github.com/tmcdonell/accelerate-blas/compare/4c89f4e6c62b8de3f37855ab2e4d27046b2495b2...0.1.0.0 diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/LinearAlgebra.hs +++ /dev/null @@ -1,172 +0,0 @@ -{-# LANGUAGE ConstraintKinds #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE NoImplicitPrelude #-} -{-# LANGUAGE ViewPatterns #-} --- | --- Module : Data.Array.Accelerate.Numeric.LinearAlgebra --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- - -module Data.Array.Accelerate.Numeric.LinearAlgebra ( - - -- * Types - Numeric, Scalar, Vector, Matrix, - - -- * Products - -- ** Vector-vector - (<.>), - (><), - - -- ** Matrix-vector - (#>), (<#), - - -- ** Matrix-matrix - (<>), - - -- * Diagonal - identity, diagonal, - -) where - -import Data.Array.Accelerate as A - -import Data.Array.Accelerate.Numeric.LinearAlgebra.Type -import Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level1 -import Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level2 -import Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level3 - - --- Level 1 --- ------- - --- | An infix synonym for 'dotu'. --- --- >>> let a = fromList (Z:.4) [1..] --- >>> let b = fromList (Z:.4) [-2,0,1,1] --- >>> a <.> b --- Scalar Z [5.0] --- --- >>> let c = fromList (Z:.2) [1:+1, 1:+0] --- >>> let d = fromList (Z:.2) [1:+0, 1:+(-1)] --- >>> c <.> d --- Scalar Z [2.0 :+ 0.0] --- -infixr 8 <.> -(<.>) :: Numeric e => Acc (Vector e) -> Acc (Vector e) -> Acc (Scalar e) -(<.>) = dotu - - --- | Outer product of two vectors --- --- >>> let a = fromList (Z :. 3) [1,2,3] --- >>> let b = fromList (Z :. 3) [5,2,3] --- >>> a >< b --- Matrix (Z :. 3 :. 3) --- [ 5.0, 2.0, 3.0 --- , 10.0, 4.0, 6.0 --- , 15.0, 6.0, 9.0 ] --- -infixr 8 >< -(><) :: Numeric e => Acc (Vector e) -> Acc (Vector e) -> Acc (Matrix e) -(><) x y = xc <> yr - where - xc = reshape (index2 (length x) 1) x - yr = reshape (index2 1 (length y)) y - - --- Level 2 --- ------- - --- | Dense matrix-vector product --- --- >>> let m = fromList (Z :. 2 :. 3) [1..] --- >>> m --- Matrix (Z :. 2 :. 3) --- [ 1.0, 2.0, 3.0 --- , 4.0, 5.0, 6.0 ] --- --- >>> let x = fromList (Z :. 3) [10,20,30] --- --- >>> m #> x --- Vector (Z :. 2) [140.0,320.0] --- --- See 'gemv' for a more general version of this operation. --- -infixr 8 #> -(#>) :: Numeric e => Acc (Matrix e) -> Acc (Vector e) -> Acc (Vector e) -(#>) m x = gemv 1 N m x - - --- | Dense vector-matrix product --- --- >>> let m = fromList (Z :. 2 :. 3) [1..] --- >>> m --- Matrix (Z :. 2 :. 3) --- [1.0,2.0,3.0, --- 4.0,5.0,6.0] --- --- >>> let v = fromList (Z :. 2) [5,10] --- --- >>> v <# m --- Vector (Z :. 3) [45.0,60.0,75.0] --- --- See 'gemv' for a more general version of this operation. --- -infixr 8 <# -(<#) :: Numeric e => Acc (Vector e) -> Acc (Matrix e) -> Acc (Vector e) -(<#) x m = gemv 1 T m x - - --- Level 3 --- ------- - --- | Dense matrix-matrix product --- --- >>> let a = fromList (Z :. 3 :. 5) [1..] --- >>> a --- Matrix (Z:.3:.5) --- [ 1.0, 2.0, 3.0, 4.0, 5.0 --- , 6.0, 7.0, 8.0, 9.0, 10.0 --- , 11.0, 12.0, 13.0, 14.0, 15.0 ] --- --- >>> let b = fromList (Z :. 5 :. 2) [1,3, 0,2, -1,5, 7,7, 6,0] --- >>> b --- Matrix (Z :. 5 :. 2) --- [ 1.0, 3.0 --- , 0.0, 2.0 --- , -1.0, 5.0 --- , 7.0, 7.0 --- , 6.0, 0.0 ] --- --- >>> a <> b --- Matrix (Z :. 3 :. 2) --- [ 56.0, 50.0 --- , 121.0, 135.0 --- , 186.0, 220.0 ] --- --- See 'gemm' for a more general version of this operation. --- -infixr 8 <> -(<>) :: Numeric e => Acc (Matrix e) -> Acc (Matrix e) -> Acc (Matrix e) -(<>) matA matB = gemm 1 N matA N matB - - --- | Create a square identity matrix of the given dimension --- -identity :: Num e => Exp Int -> Acc (Matrix e) -identity n = diagonal (fill (index1 n) 1) - --- | Create a square matrix with the given diagonal --- -diagonal :: Num e => Acc (Vector e) -> Acc (Matrix e) -diagonal v = - let n = length v - zeros = fill (index2 n n) 0 - in - permute const zeros (\(unindex1 -> i) -> index2 i i) v - diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level1.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level1.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level1.hs +++ /dev/null @@ -1,139 +0,0 @@ -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE NoImplicitPrelude #-} -{-# LANGUAGE ScopedTypeVariables #-} --- | --- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level1 --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- --- Level 1 (vector-vector) BLAS operations. --- - -module Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level1 ( - - -- Types - Numeric, Vector, - - -- Level1 operations - sdot, - dotu, - dotc, - asum, - amax, - amin, - -) where - -import Data.Array.Accelerate as A -import Data.Array.Accelerate.Data.Complex as A -import Data.Array.Accelerate.Numeric.LinearAlgebra.Type - - --- | Computes a vector-vector dot product, using double precision accumulation --- of the intermediate result. Includes a scalar (initial) value to be added to --- the inner product. --- --- --- -sdot :: forall e. Numeric e => Exp e -> Acc (Vector e) -> Acc (Vector e) -> Acc (Scalar e) -sdot z xs ys = - case numericR :: NumericR e of - NumericRfloat32 -> map toFloating $ dsdot (toFloating z) (map toFloating xs) (map toFloating ys) - NumericRfloat64 -> dsdot z xs ys - NumericRcomplex32 -> map d2f $ zsdot (f2d z) (map f2d xs) (map f2d ys) - NumericRcomplex64 -> zsdot z xs ys - where - dsdot :: Exp Double -> Acc (Vector Double) -> Acc (Vector Double) -> Acc (Scalar Double) - dsdot z' xs' ys' = fold (+) z' (zipWith (*) xs' ys') - - zsdot :: Exp (Complex Double) -> Acc (Vector (Complex Double)) -> Acc (Vector (Complex Double)) -> Acc (Scalar (Complex Double)) - zsdot z' xs' ys' = fold (+) z' (zipWith (*) xs' ys') - - f2d :: Exp (Complex Float) -> Exp (Complex Double) - f2d c = lift (toFloating (real c) :+ toFloating (imag c)) - - d2f :: Exp (Complex Double) -> Exp (Complex Float) - d2f c = lift (toFloating (real c) :+ toFloating (imag c)) - - --- | Computes a vector-vector dot product --- --- \[ --- res = \sum_i x_i * y_i --- \] --- --- --- -dotu :: Numeric e => Acc (Vector e) -> Acc (Vector e) -> Acc (Scalar e) -dotu xs ys = fold (+) 0 (zipWith (*) xs ys) - - --- | Computes a dot product of a conjugated vector with another vector --- --- \[ --- res = \sum_i \mathrm{conj}(x_i) * y_i --- \] --- --- --- -dotc :: forall e. Numeric (Complex e) - => Acc (Vector (Complex e)) - -> Acc (Vector (Complex e)) - -> Acc (Scalar (Complex e)) -dotc xs ys = - case numericR :: NumericR (Complex e) of - NumericRcomplex32 -> dotu (map conjugate xs) ys - NumericRcomplex64 -> dotu (map conjugate xs) ys - - --- | Computes the sum of magnitudes of the vector elements. For complex values, --- this is given by \(\sum_i \|\mathrm{real}(x_i)\| + \|\mathrm{imag}(x_i)\|\). --- --- --- -asum :: forall e. Numeric e => Acc (Vector e) -> Acc (Scalar (NumericBaseT e)) -asum = - case numericR :: NumericR e of - NumericRfloat32 -> sum . map abs - NumericRfloat64 -> sum . map abs - NumericRcomplex32 -> sum . map mag - NumericRcomplex64 -> sum . map mag - where - mag c = abs (real c) + abs (imag c) - - --- | Return the index of the element with the maximum absolute value. --- --- --- -amax :: forall e. Numeric e => Acc (Vector e) -> Acc (Scalar Int) -amax = - case numericR :: NumericR e of - NumericRfloat32 -> map (indexHead . fst) . fold1 cmp . indexed . map abs - NumericRfloat64 -> map (indexHead . fst) . fold1 cmp . indexed . map abs - NumericRcomplex32 -> map (indexHead . fst) . fold1 cmp . indexed . map mag - NumericRcomplex64 -> map (indexHead . fst) . fold1 cmp . indexed . map mag - where - cmp ix iy = snd ix > snd iy ? ( ix, iy ) - mag c = abs (real c) + abs (imag c) - --- | Return the index of the element with the minimum absolute value. --- --- --- -amin :: forall e. Numeric e => Acc (Vector e) -> Acc (Scalar Int) -amin = - case numericR :: NumericR e of - NumericRfloat32 -> map (indexHead . fst) . fold1 cmp . indexed . map abs - NumericRfloat64 -> map (indexHead . fst) . fold1 cmp . indexed . map abs - NumericRcomplex32 -> map (indexHead . fst) . fold1 cmp . indexed . map mag - NumericRcomplex64 -> map (indexHead . fst) . fold1 cmp . indexed . map mag - where - cmp ix iy = snd ix < snd iy ? ( ix, iy ) - mag c = abs (real c) + abs (imag c) - diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level2.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level2.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level2.hs +++ /dev/null @@ -1,92 +0,0 @@ -{-# LANGUAGE CPP #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE NoImplicitPrelude #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE ViewPatterns #-} --- | --- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level2 --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- --- Level 2 (matrix-vector) BLAS operations. --- - -module Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level2 ( - - -- Types - Numeric, Vector, Matrix, Transpose(..), - - -- Operations - gemv, - -) where - -import Data.Array.Accelerate as A -import Data.Array.Accelerate.Smart as A -import Data.Array.Accelerate.Data.Complex as A -import Data.Array.Accelerate.Numeric.LinearAlgebra.Type - -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND -import qualified Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level2 as CPU -#endif -#ifdef ACCELERATE_LLVM_PTX_BACKEND -import qualified Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level2 as PTX -#endif - - --- | Computes the matrix-vector product of a general matrix. --- --- \[ --- y = \alpha * \mathrm{op}(A) * x --- \] --- --- where: --- --- * 'shape' \(\mathrm{op}(A)\) @= Z :. m :. n@ --- * 'shape' \(x\) @= Z :. n@ --- * 'shape' \(y\) @= Z :. m@ --- --- --- -gemv :: forall e. Numeric e - => Exp e -- ^ \( \alpha \) - -> Transpose -- ^ Operation to apply to A - -> Acc (Matrix e) -- ^ A - -> Acc (Vector e) -- ^ x - -> Acc (Vector e) -- ^ y -gemv alpha opA matA x = go (lift (unit alpha, matA, x)) - where - go = -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND - foreignAcc (CPU.gemv opA) $ -#endif -#ifdef ACCELERATE_LLVM_PTX_BACKEND - foreignAcc (PTX.gemv opA) $ -#endif - (\(unatup3 -> (_, arr, brr)) -> mXv arr brr) - - -- General matrix-vector multiply in pure Accelerate. This is probably not - -- efficient. - -- - mXv :: Acc (Matrix e) -> Acc (Vector e) -> Acc (Vector e) - mXv arr brr - = fold (+) 0 - $ zipWith (\a b -> alpha * a * b) arr' brr' - where - Z :. m :. _ = unlift (shape arr') :: Z :. Exp Int :. Exp Int - - brr' = replicate (lift (Z :. m :. All)) brr - arr' = case opA of - N -> arr - T -> transpose arr - H -> case numericR :: NumericR e of - NumericRcomplex32 -> map conjugate (transpose arr) - NumericRcomplex64 -> map conjugate (transpose arr) - _ -> transpose arr - diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level3.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level3.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level3.hs +++ /dev/null @@ -1,109 +0,0 @@ -{-# LANGUAGE CPP #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE NoImplicitPrelude #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE ViewPatterns #-} --- | --- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level3 --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- --- Level 3 (matrix-matrix) BLAS operations. --- - -module Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level3 ( - - -- Types - Numeric, Matrix, Transpose(..), - - -- Matrix-matrix operations - gemm, - -) where - -import Data.Array.Accelerate as A -import Data.Array.Accelerate.Smart as A -import Data.Array.Accelerate.Data.Complex as A -import Data.Array.Accelerate.Numeric.LinearAlgebra.Type - -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND -import qualified Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level3 as CPU -#endif -#ifdef ACCELERATE_LLVM_PTX_BACKEND -import qualified Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level3 as PTX -#endif - - --- | General matrix-matrix multiply --- --- \[ --- C = \alpha * \mathrm{op}(A) * \mathrm{op}(B) --- \] --- --- where: --- --- * 'shape' \(\mathrm{op}(A)\) @= Z :. m :. k@ --- * 'shape' \(\mathrm{op}(B)\) @= Z :. k :. n@ --- * 'shape' \(C\) @= Z :. m :. n@ --- --- --- -gemm :: forall e. Numeric e - => Exp e -- ^ \( \alpha \) - -> Transpose -- ^ operation to apply to A - -> Acc (Matrix e) -- ^ A - -> Transpose -- ^ operation to apply to B - -> Acc (Matrix e) -- ^ B - -> Acc (Matrix e) -- ^ C -gemm alpha opA matA opB matB = go (lift (unit alpha, matA, matB)) - where - go = -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND - foreignAcc (CPU.gemm opA opB) $ -#endif -#ifdef ACCELERATE_LLVM_PTX_BACKEND - foreignAcc (PTX.gemm opA opB) $ -#endif - (\(unatup3 -> (_, arr, brr)) -> mXm arr brr) - - -- General dense matrix-matrix multiply written in pure Accelerate. This is - -- not efficient due to the memory access patterns. We could probably - -- improve this a little bit with a divide-and-conquer algorithm, for - -- example, but using a foreign implementation will be best. - -- - mXm :: Acc (Matrix e) -> Acc (Matrix e) -> Acc (Matrix e) - mXm arr brr - = fold (+) 0 - $ zipWith (\a b -> alpha * a * b) arrRepl brrRepl - where - Z :. rowsA :. _ = unlift (shape arr') :: Z :. Exp Int :. Exp Int - Z :. colsB :. _ = unlift (shape brr') :: Z :. Exp Int :. Exp Int - -- - arrRepl = replicate (lift $ Z :. All :. colsB :. All) arr' - brrRepl = replicate (lift $ Z :. rowsA :. All :. All) brr' - - -- apply opA - arr' = case opA of - N -> arr - T -> transpose arr - H -> case numericR :: NumericR e of - NumericRcomplex32 -> map conjugate (transpose arr) - NumericRcomplex64 -> map conjugate (transpose arr) - _ -> transpose arr - - -- apply opB and transpose at the same time, which is required for this - -- algorithm - brr' = case opB of - N -> transpose brr - T -> brr - H -> case numericR :: NumericR e of - NumericRcomplex32 -> map conjugate brr - NumericRcomplex64 -> map conjugate brr - _ -> brr - diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Base.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Base.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Base.hs +++ /dev/null @@ -1,106 +0,0 @@ -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE ForeignFunctionInterface #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TypeFamilies #-} --- | --- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Base --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- - -module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Base - where - -import Data.Array.Accelerate.Data.Complex -import Data.Array.Accelerate.Array.Sugar ( Array(..), EltRepr ) -import Data.Array.Accelerate.Array.Data -import Data.Array.Accelerate.Array.Unique -import Data.Array.Accelerate.Numeric.LinearAlgebra.Type - -import Foreign.Marshal.Alloc -import Foreign.Ptr -import Foreign.Storable -import Foreign.Storable.Complex ( ) - -import qualified Blas.Primitive.Types as C - - -encodeTranspose :: Transpose -> C.Transpose -encodeTranspose N = C.NoTrans -encodeTranspose T = C.Trans -encodeTranspose H = C.ConjTrans - - -withArray - :: forall sh e b. Numeric e - => Array sh e - -> (ArrayPtrs (EltRepr e) -> IO b) - -> IO b -withArray (Array _ adata) = withArrayData (numericR::NumericR e) adata - -withArrayData - :: NumericR e - -> ArrayData (EltRepr e) - -> (ArrayPtrs (EltRepr e) -> IO b) - -> IO b -withArrayData NumericRfloat32 (AD_Float ua) f = withUniqueArrayPtr ua f -withArrayData NumericRfloat64 (AD_Double ua) f = withUniqueArrayPtr ua f -withArrayData NumericRcomplex32 (AD_Pair ad1 ad2) f - | AD_Pair AD_Unit (AD_Float ua_re) <- ad1 - , AD_Float ua_im <- ad2 - = withUniqueArrayPtr ua_re $ \p_re -> - withUniqueArrayPtr ua_im $ \p_im -> - f (((),p_re), p_im) - -withArrayData NumericRcomplex64 (AD_Pair ad1 ad2) f - | AD_Pair AD_Unit (AD_Double ua_re) <- ad1 - , AD_Double ua_im <- ad2 - = withUniqueArrayPtr ua_re $ \p_re -> - withUniqueArrayPtr ua_im $ \p_im -> - f (((),p_re), p_im) - - -interleave - :: forall e b. (Storable e, Numeric (Complex e)) - => ArrayPtrs (EltRepr (Complex e)) - -> Int - -> (Ptr (Complex e) -> IO b) - -> IO b -interleave (((), p_re), p_im) n k = do - allocaBytesAligned (n * sizeOf (undefined::Complex e)) 16 $ \p_cplx -> do - () <- case numericR :: NumericR (Complex e) of - NumericRcomplex32 -> c_interleave_f32 0 n p_cplx p_re p_im - NumericRcomplex64 -> c_interleave_f64 0 n p_cplx p_re p_im - -- - k p_cplx - - -deinterleave - :: forall e. (Storable e, Numeric (Complex e)) - => ArrayPtrs (EltRepr (Complex e)) - -> Ptr (Complex e) - -> Int - -> IO () -deinterleave (((), p_re), p_im) p_cplx n = - case numericR :: NumericR (Complex e) of - NumericRcomplex32 -> c_deinterleave_f32 0 n p_re p_im p_cplx - NumericRcomplex64 -> c_deinterleave_f64 0 n p_re p_im p_cplx - - -foreign import ccall unsafe "interleave_f32" - c_interleave_f32 :: Int -> Int -> Ptr (Complex Float) -> Ptr Float -> Ptr Float -> IO () - -foreign import ccall unsafe "interleave_f64" - c_interleave_f64 :: Int -> Int -> Ptr (Complex Double) -> Ptr Double -> Ptr Double -> IO () - -foreign import ccall unsafe "deinterleave_f32" - c_deinterleave_f32 :: Int -> Int -> Ptr Float -> Ptr Float -> Ptr (Complex Float) -> IO () - -foreign import ccall unsafe "deinterleave_f64" - c_deinterleave_f64 :: Int -> Int -> Ptr Double -> Ptr Double -> Ptr (Complex Double) -> IO () - diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level2.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level2.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level2.hs +++ /dev/null @@ -1,72 +0,0 @@ -{-# LANGUAGE GADTs #-} -{-# LANGUAGE ScopedTypeVariables #-} --- | --- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level2 --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- - -module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level2 - where - -import Data.Array.Accelerate as A -import Data.Array.Accelerate.Data.Complex -import Data.Array.Accelerate.LLVM.Native.Foreign -import Data.Array.Accelerate.Numeric.LinearAlgebra.Type -import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Base - -import Foreign.Marshal.Alloc -import Foreign.Storable -import Foreign.Storable.Complex ( ) - -import qualified Blas.Primitive.Types as C -import qualified Blas.Primitive.Unsafe as C - - -gemv :: forall e. Numeric e - => Transpose - -> ForeignAcc ((Scalar e, Matrix e, Vector e) -> Vector e) -gemv opA = ForeignAcc "native.gemv" gemv' - where - gemv' (alpha, matA, vecx) = do - let - Z :. rowsA :. colsA = arrayShape matA - Z :. sizeX = arrayShape vecx - - sizeA = rowsA * colsA - sizeY = case opA of - N -> rowsA - _ -> colsA - - opA' = encodeTranspose opA - alpha' = indexArray alpha Z - -- - vecy <- allocateRemote (Z :. sizeY) :: LLVM Native (Vector e) - () <- liftIO $ do - withArray matA $ \ptr_A -> do - withArray vecx $ \ptr_x -> do - withArray vecy $ \ptr_y -> do - case numericR :: NumericR e of - NumericRfloat32 -> C.sgemv C.RowMajor opA' rowsA colsA alpha' ptr_A colsA ptr_x 1 0 ptr_y 1 - NumericRfloat64 -> C.dgemv C.RowMajor opA' rowsA colsA alpha' ptr_A colsA ptr_x 1 0 ptr_y 1 - -- - NumericRcomplex32 -> do - allocaBytesAligned (sizeY * sizeOf (undefined::Complex e)) 16 $ \ptr_y' -> do - interleave ptr_A sizeA $ \ptr_A' -> do - interleave ptr_x sizeX $ \ptr_x' -> do - C.cgemv C.RowMajor opA' rowsA colsA alpha' ptr_A' colsA ptr_x' 1 0 ptr_y' 1 - deinterleave ptr_y ptr_y' sizeY - -- - NumericRcomplex64 -> do - allocaBytesAligned (sizeY * sizeOf (undefined::Complex e)) 16 $ \ptr_y' -> do - interleave ptr_A sizeA $ \ptr_A' -> do - interleave ptr_x sizeX $ \ptr_x' -> do - C.zgemv C.RowMajor opA' rowsA colsA alpha' ptr_A' colsA ptr_x' 1 0 ptr_y' 1 - deinterleave ptr_y ptr_y' sizeY - -- - return vecy - diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level3.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level3.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level3.hs +++ /dev/null @@ -1,87 +0,0 @@ -{-# LANGUAGE GADTs #-} -{-# LANGUAGE ScopedTypeVariables #-} --- | --- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level3 --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- - -module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level3 - where - -import Data.Array.Accelerate as A -import Data.Array.Accelerate.Data.Complex -import Data.Array.Accelerate.LLVM.Native.Foreign -import Data.Array.Accelerate.Numeric.LinearAlgebra.Type -import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Base - -import Foreign.Marshal.Alloc -import Foreign.Storable -import Foreign.Storable.Complex ( ) - -import qualified Blas.Primitive.Types as C -import qualified Blas.Primitive.Unsafe as C - - --- TODO: check whether it is faster to compute this as column-major order: --- --- https://www.christophlassner.de/using-blas-from-c-with-row-major-data.html --- -gemm :: forall e. Numeric e - => Transpose - -> Transpose - -> ForeignAcc ((Scalar e, Matrix e, Matrix e) -> Matrix e) -gemm opA opB = ForeignAcc "native.gemm" gemm' - where - gemm' (alpha, matA, matB) = do - let - Z :. rowsA :. colsA = arrayShape matA - Z :. rowsB :. colsB = arrayShape matB - - sizeA = rowsA * colsA - sizeB = rowsB * colsB - sizeC = m * n - - (m,k) = case opA of - N -> (rowsA, colsA) - _ -> (colsA, rowsA) - n = case opB of - N -> colsB - _ -> rowsB - - lda = colsA - ldb = colsB - - opA' = encodeTranspose opA - opB' = encodeTranspose opB - alpha' = indexArray alpha Z - -- - matC <- allocateRemote (Z :. m :. n) :: LLVM Native (Matrix e) - () <- liftIO $ do - withArray matA $ \ptr_A -> do - withArray matB $ \ptr_B -> do - withArray matC $ \ptr_C -> do - case numericR :: NumericR e of - NumericRfloat32 -> C.sgemm C.RowMajor opA' opB' m n k alpha' ptr_A lda ptr_B ldb 0 ptr_C n - NumericRfloat64 -> C.dgemm C.RowMajor opA' opB' m n k alpha' ptr_A lda ptr_B ldb 0 ptr_C n - -- - NumericRcomplex32 -> do - allocaBytesAligned (sizeC * sizeOf (undefined::Complex e)) 16 $ \ptr_C' -> do - interleave ptr_A sizeA $ \ptr_A' -> do - interleave ptr_B sizeB $ \ptr_B' -> do - C.cgemm C.RowMajor opA' opB' m n k alpha' ptr_A' lda ptr_B' ldb 0 ptr_C' n - deinterleave ptr_C ptr_C' sizeC - -- - NumericRcomplex64 -> do - allocaBytesAligned (sizeC * sizeOf (undefined::Complex e)) 16 $ \ptr_C' -> do - interleave ptr_A sizeA $ \ptr_A' -> do - interleave ptr_B sizeB $ \ptr_B' -> do - C.zgemm C.RowMajor opA' opB' m n k alpha' ptr_A' lda ptr_B' ldb 0 ptr_C' n - deinterleave ptr_C ptr_C' sizeC - -- - return matC - diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Base.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Base.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Base.hs +++ /dev/null @@ -1,82 +0,0 @@ -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TypeFamilies #-} --- | --- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- - -module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base - where - -import Data.Array.Accelerate.Lifetime -import Data.Array.Accelerate.Array.Sugar ( Array(..), EltRepr ) -import Data.Array.Accelerate.Array.Data -import Data.Array.Accelerate.Numeric.LinearAlgebra.Type - -import Data.Array.Accelerate.LLVM.PTX.Foreign - -import Foreign.CUDA.Ptr ( DevicePtr ) -import qualified Foreign.CUDA.BLAS as C - - -type family DevicePtrs e :: * -type instance DevicePtrs () = () -type instance DevicePtrs Float = DevicePtr Float -type instance DevicePtrs Double = DevicePtr Double -type instance DevicePtrs (a,b) = (DevicePtrs a, DevicePtrs b) - - -encodeTranspose :: Transpose -> C.Operation -encodeTranspose N = C.N -encodeTranspose T = C.T -encodeTranspose H = C.C - - -withArray - :: forall sh e b. Numeric e - => Array sh e - -> Stream - -> (DevicePtrs (EltRepr e) -> LLVM PTX b) - -> LLVM PTX b -withArray (Array _ adata) s k = withArrayData (numericR::NumericR e) adata s k - -withArrayData - :: NumericR e - -> ArrayData (EltRepr e) - -> Stream - -> (DevicePtrs (EltRepr e) -> LLVM PTX b) - -> LLVM PTX b -withArrayData NumericRfloat32 ad s k = - withDevicePtr ad $ \p -> do - r <- k p - e <- checkpoint s - return (Just e,r) -withArrayData NumericRfloat64 ad s k = - withDevicePtr ad $ \p -> do - r <- k p - e <- checkpoint s - return (Just e, r) -withArrayData NumericRcomplex32 (AD_Pair (AD_Pair AD_Unit ad1) ad2) s k = - withDevicePtr ad1 $ \p1 -> - withDevicePtr ad2 $ \p2 -> do - r <- k (((), p1), p2) - e <- checkpoint s - return (Just e, (Just e, r)) -withArrayData NumericRcomplex64 (AD_Pair (AD_Pair AD_Unit ad1) ad2) s k = - withDevicePtr ad1 $ \p1 -> - withDevicePtr ad2 $ \p2 -> do - r <- k (((), p1), p2) - e <- checkpoint s - return (Just e, (Just e, r)) - -withLifetime' :: Lifetime a -> (a -> LLVM PTX b) -> LLVM PTX b -withLifetime' l k = do - r <- k (unsafeGetValue l) - liftIO $ touchLifetime l - return r - diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Context.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Context.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Context.hs +++ /dev/null @@ -1,76 +0,0 @@ -{-# LANGUAGE MagicHash #-} --- | --- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Context --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- - -module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Context ( - - withBLAS - -) where - -import Data.Array.Accelerate.Lifetime -import Data.Array.Accelerate.LLVM.PTX -import Data.Array.Accelerate.LLVM.PTX.Foreign -import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base - -import Control.Monad.State -import Control.Concurrent.MVar -import Data.IntMap.Strict ( IntMap ) -import System.IO.Unsafe -import qualified Data.IntMap.Strict as IM - -import qualified Foreign.CUDA.Driver.Context as CUDA -import qualified Foreign.CUDA.BLAS as BLAS - -import GHC.Ptr -import GHC.Base -import Prelude hiding ( lookup ) - - --- Execute an operation with a cuBLAS handle appropriate for the current --- execution context. --- --- Initial creation of the context is an atomic operation, but subsequently --- multiple threads may use the context concurrently. --- --- --- -withBLAS :: (BLAS.Handle -> LLVM PTX b) -> LLVM PTX b -withBLAS k = do - lc <- gets (deviceContext . ptxContext) - h <- liftIO $ - withLifetime lc $ \ctx -> - modifyMVar handles $ \im -> - let key = toKey ctx in - case IM.lookup key im of - -- handle does not exist yet; create it and add to the global - -- state for reuse - Nothing -> do - h <- BLAS.create - l <- newLifetime h - -- BLAS.setPointerMode h BLAS.Device - BLAS.setAtomicsMode h BLAS.Allowed - addFinalizer lc $ modifyMVar handles (\im' -> return (IM.delete key im', ())) - addFinalizer l $ BLAS.destroy h - return ( IM.insert key l im, l ) - - -- return existing handle - Just h -> return (im, h) - -- - withLifetime' h k - - -toKey :: CUDA.Context -> IM.Key -toKey (CUDA.Context (Ptr addr#)) = I# (addr2Int# addr#) - -{-# NOINLINE handles #-} -handles :: MVar (IntMap (Lifetime BLAS.Handle)) -handles = unsafePerformIO $ newMVar IM.empty - diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level2.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level2.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level2.hs +++ /dev/null @@ -1,131 +0,0 @@ -{-# LANGUAGE GADTs #-} -{-# LANGUAGE ScopedTypeVariables #-} --- | --- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level2 --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- - -module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level2 - where - -import Data.Array.Accelerate as A -import Data.Array.Accelerate.Array.Sugar ( Array(..) ) -import Data.Array.Accelerate.Data.Complex -import Data.Array.Accelerate.LLVM.PTX.Foreign -import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base -import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Context -import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level3 -import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Twine -import Data.Array.Accelerate.Numeric.LinearAlgebra.Type - -import Foreign.Marshal ( with ) -import Foreign.Storable.Complex ( ) - -import qualified Foreign.CUDA.Ptr as CUDA -import qualified Foreign.CUDA.BLAS as BLAS - - --- NOTE: cuBLAS requires matrices to be stored in column-major order --- (Fortran-style), but Accelerate uses C-style arrays in row-major order. --- --- If the operation is N or T, we can just swap the operation. For --- conjugate-transpose (H) operations (on complex valued arguments), since there --- is no conjugate-no-transpose operation, we implement that via 'gemm', which --- I assume is more efficient than ?geam followed by ?gemv. --- -gemv :: Numeric e - => Transpose - -> ForeignAcc ((Scalar e, Matrix e, Vector e) -> Vector e) -gemv opA = ForeignAcc "ptx.gemv" (gemv' numericR opA) - -gemv' :: Numeric e - => NumericR e - -> Transpose - -> Stream - -> (Scalar e, Matrix e, Vector e) - -> LLVM PTX (Vector e) -gemv' NumericRcomplex32 H = as_gemm H -gemv' NumericRcomplex64 H = as_gemm H -gemv' _ t = as_gemv t - - -as_gemm - :: Numeric e - => Transpose - -> Stream - -> (Scalar e, Matrix e, Vector e) - -> LLVM PTX (Vector e) -as_gemm opA stream (alpha, matA, Array sh adata) = do - let matB = Array (sh,1) adata - -- - Array (sh',1) vecy <- gemm' opA N stream (alpha, matA, matB) - return (Array sh' vecy) - -as_gemv - :: forall e. Numeric e - => Transpose - -> Stream - -> (Scalar e, Matrix e, Vector e) - -> LLVM PTX (Vector e) -as_gemv opA stream (alpha, matA, vecx) = do - let - Z :. rowsA :. colsA = arrayShape matA - Z :. sizeX = arrayShape vecx - - sizeA = rowsA * colsA - sizeY = case opA of - N -> rowsA - _ -> colsA - - opA' = encodeTranspose - $ case opA of - N -> T - _ -> N - -- - vecy <- allocateRemote (Z :. sizeY) :: LLVM PTX (Vector e) - alpha' <- indexRemote alpha 0 - () <- do - withArray matA stream $ \ptr_A -> do - withArray vecx stream $ \ptr_x -> do - withArray vecy stream $ \ptr_y -> do - withBLAS $ \hdl -> do - case numericR :: NumericR e of - NumericRfloat32 -> liftIO $ - with alpha' $ \ptr_alpha -> - with 0 $ \ptr_beta -> - BLAS.sgemv hdl opA' colsA rowsA ptr_alpha ptr_A colsA ptr_x 1 ptr_beta ptr_y 1 - - NumericRfloat64 -> liftIO $ - with alpha' $ \ptr_alpha -> - with 0 $ \ptr_beta -> - BLAS.dgemv hdl opA' colsA rowsA ptr_alpha ptr_A colsA ptr_x 1 ptr_beta ptr_y 1 - - NumericRcomplex32 -> do - tmpy <- allocateRemote (Z :. sizeY * 2) :: LLVM PTX (Vector Float) - withArray tmpy stream $ \ptr_y' -> do - interleave ptr_A stream sizeA $ \ptr_A' -> do - interleave ptr_x stream sizeX $ \ptr_x' -> do - liftIO $ do - with alpha' $ \ptr_alpha -> - with 0 $ \ptr_beta -> do - BLAS.cgemv hdl opA' colsA rowsA ptr_alpha ptr_A' colsA ptr_x' 1 ptr_beta (CUDA.castDevPtr ptr_y' :: CUDA.DevicePtr (Complex Float)) 1 - deinterleave ptr_y (CUDA.castDevPtr ptr_y' :: CUDA.DevicePtr (Complex Float)) stream sizeY - - NumericRcomplex64 -> do - tmpy <- allocateRemote (Z :. sizeY * 2) :: LLVM PTX (Vector Double) - withArray tmpy stream $ \ptr_y' -> do - interleave ptr_A stream sizeA $ \ptr_A' -> do - interleave ptr_x stream sizeX $ \ptr_x' -> do - liftIO $ do - with alpha' $ \ptr_alpha -> - with 0 $ \ptr_beta -> do - BLAS.zgemv hdl opA' colsA rowsA ptr_alpha ptr_A' colsA ptr_x' 1 ptr_beta (CUDA.castDevPtr ptr_y' :: CUDA.DevicePtr (Complex Double)) 1 - deinterleave ptr_y (CUDA.castDevPtr ptr_y' :: CUDA.DevicePtr (Complex Double)) stream sizeY - -- - return vecy - diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level3.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level3.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level3.hs +++ /dev/null @@ -1,113 +0,0 @@ -{-# LANGUAGE GADTs #-} -{-# LANGUAGE ScopedTypeVariables #-} --- | --- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level3 --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- - -module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level3 - where - -import Data.Array.Accelerate as A -import Data.Array.Accelerate.Data.Complex -import Data.Array.Accelerate.LLVM.PTX.Foreign -import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base -import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Context -import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Twine -import Data.Array.Accelerate.Numeric.LinearAlgebra.Type - -import Foreign.Marshal ( with ) -import Foreign.Storable.Complex ( ) - -import qualified Foreign.CUDA.Ptr as CUDA -import qualified Foreign.CUDA.BLAS as BLAS - - --- NOTE: cuBLAS requires that matrices are stored in column-major order --- (Fortran-style), but Accelerate uses a C-style convention where matrices are --- stored in row-major order. --- --- At least for matrix-matrix multiply, we can get around this problem by making --- use of the equivalence \( B^T \cdot A^T = (A \cdot B)^T \). --- -gemm :: Numeric e - => Transpose - -> Transpose - -> ForeignAcc ((Scalar e, Matrix e, Matrix e) -> Matrix e) -gemm opA opB = ForeignAcc "ptx.gemm" (gemm' opA opB) - -gemm' - :: forall e. Numeric e - => Transpose - -> Transpose - -> Stream - -> (Scalar e, Matrix e, Matrix e) - -> LLVM PTX (Matrix e) -gemm' opA opB stream (alpha, matA, matB) = do - let - Z :. rowsA :. colsA = arrayShape matA - Z :. rowsB :. colsB = arrayShape matB - - sizeA = rowsA * colsA - sizeB = rowsB * colsB - sizeC = m * n - - (m,k) = case opA of - N -> (rowsA, colsA) - _ -> (colsA, rowsA) - n = case opB of - N -> colsB - _ -> rowsB - - lda = colsA - ldb = colsB - - opA' = encodeTranspose opA - opB' = encodeTranspose opB - -- - matC <- allocateRemote (Z :. m :. n) :: LLVM PTX (Matrix e) - alpha' <- indexRemote alpha 0 - () <- withArray matA stream $ \ptr_A -> do - withArray matB stream $ \ptr_B -> do - withArray matC stream $ \ptr_C -> do - withBLAS $ \hdl -> do - case numericR :: NumericR e of - NumericRfloat32 -> liftIO $ - with alpha' $ \ptr_alpha -> - with 0 $ \ptr_beta -> - BLAS.sgemm hdl opB' opA' n m k ptr_alpha ptr_B ldb ptr_A lda ptr_beta ptr_C n - - NumericRfloat64 -> liftIO $ - with alpha' $ \ptr_alpha -> - with 0 $ \ptr_beta -> - BLAS.dgemm hdl opB' opA' n m k ptr_alpha ptr_B ldb ptr_A lda ptr_beta ptr_C n - - NumericRcomplex32 -> do - tmpC <- allocateRemote (Z :. sizeC * 2) :: LLVM PTX (Vector Float) - withArray tmpC stream $ \ptr_C' -> do - interleave ptr_A stream sizeA $ \ptr_A' -> do - interleave ptr_B stream sizeB $ \ptr_B' -> do - liftIO $ - with alpha' $ \ptr_alpha -> - with 0 $ \ptr_beta -> - BLAS.cgemm hdl opB' opA' n m k ptr_alpha ptr_B' ldb ptr_A' lda ptr_beta (CUDA.castDevPtr ptr_C') n - deinterleave ptr_C (CUDA.castDevPtr ptr_C' :: CUDA.DevicePtr (Complex Float)) stream sizeC - - NumericRcomplex64 -> do - tmpC <- allocateRemote (Z :. sizeC * 2) :: LLVM PTX (Vector Double) - withArray tmpC stream $ \ptr_C' -> do - interleave ptr_A stream sizeA $ \ptr_A' -> do - interleave ptr_B stream sizeB $ \ptr_B' -> do - liftIO $ - with alpha' $ \ptr_alpha -> - with 0 $ \ptr_beta -> - BLAS.zgemm hdl opB' opA' n m k ptr_alpha ptr_B' ldb ptr_A' lda ptr_beta (CUDA.castDevPtr ptr_C') n - deinterleave ptr_C (CUDA.castDevPtr ptr_C' :: CUDA.DevicePtr (Complex Double)) stream sizeC - - return matC - diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Twine.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Twine.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Twine.hs +++ /dev/null @@ -1,176 +0,0 @@ -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE MagicHash #-} -{-# LANGUAGE RecordWildCards #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TemplateHaskell #-} --- | --- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Twine --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- - -module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Twine ( - - interleave, - deinterleave, - -) where - -import Data.Array.Accelerate.Data.Complex -import Data.Array.Accelerate.Array.Sugar ( EltRepr, Vector, Z(..), (:.)(..) ) -import Data.Array.Accelerate.Lifetime -import Data.Array.Accelerate.LLVM.PTX -import Data.Array.Accelerate.LLVM.PTX.Foreign - -import Data.Array.Accelerate.Numeric.LinearAlgebra.Type -import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base - -import Control.Concurrent.MVar -import Control.Monad.State -import Data.ByteString ( ByteString ) -import Data.FileEmbed -import Data.IntMap.Strict ( IntMap ) -import Foreign.Storable.Complex ( ) -import System.IO.Unsafe -import qualified Data.IntMap.Strict as IM - -import Foreign.CUDA.Ptr ( DevicePtr ) -import Foreign.CUDA.Analysis -import qualified Foreign.CUDA.Driver as CUDA -import qualified Foreign.CUDA.Driver.Stream as CUDA - -import GHC.Ptr -import GHC.Base -import Prelude hiding ( lookup ) - - -interleave - :: forall e b. Numeric (Complex e) - => DevicePtrs (EltRepr (Complex e)) - -> Stream - -> Int - -> (DevicePtr (Complex e) -> LLVM PTX b) -- device pointer is in packed representation - -> LLVM PTX b -interleave (((), d_re), d_im) s n k = do - case numericR :: NumericR (Complex e) of - nR@NumericRcomplex32 -> do - cplx <- allocateRemote (Z :. n * 2) :: LLVM PTX (Vector Float) - withTwine nR $ \(_,pack,_) -> do - withArray cplx s $ \d_cplx -> do - withLifetime' s $ \s' -> do - liftIO $ launch pack s' n d_cplx d_re d_im - k (CUDA.castDevPtr d_cplx :: DevicePtr (Complex Float)) - -- - nR@NumericRcomplex64 -> do - cplx <- allocateRemote (Z :. n * 2) :: LLVM PTX (Vector Double) - withTwine nR $ \(_,pack,_) -> do - withArray cplx s $ \d_cplx -> do - withLifetime' s $ \s' -> do - liftIO $ launch pack s' n d_cplx d_re d_im - k (CUDA.castDevPtr d_cplx :: DevicePtr (Complex Double)) - -deinterleave - :: forall e. Numeric (Complex e) - => DevicePtrs (EltRepr (Complex e)) - -> DevicePtr (Complex e) -- in packed representation - -> Stream - -> Int - -> LLVM PTX () -deinterleave (((), d_re), d_im) d_cplx s n = do - case numericR :: NumericR (Complex e) of - nR@NumericRcomplex32 -> do - withTwine nR $ \(_,_,unpack) -> do - withLifetime' s $ \s' -> do - liftIO $ launch unpack s' n d_re d_im (CUDA.castDevPtr d_cplx :: DevicePtr Float) - -- - nR@NumericRcomplex64 -> do - withTwine nR $ \(_,_,unpack) -> do - withLifetime' s $ \s' -> do - liftIO $ launch unpack s' n d_re d_im (CUDA.castDevPtr d_cplx :: DevicePtr Double) - - -withTwine :: NumericR (Complex e) -> ((CUDA.Module, Kernel, Kernel) -> LLVM PTX b) -> LLVM PTX b -withTwine nR k = do - ptx <- gets ptxContext - let lc = deviceContext ptx - prp = deviceProperties ptx - mds = modules nR - -- - mdl <- liftIO $ do - withLifetime lc $ \ctx -> do - modifyMVar mds $ \im -> do - let key = toKey ctx - case IM.lookup key im of - -- Module is not loaded yet; add to the current context and the global - -- state for later reuse - Nothing -> do - mdl <- CUDA.loadData $ case nR of - NumericRcomplex32 -> ptx_twine_f32 - NumericRcomplex64 -> ptx_twine_f64 - pack <- mkKernel "interleave" mdl prp - unpack <- mkKernel "deinterleave" mdl prp - let mkk = (mdl, pack, unpack) - -- - lm <- newLifetime mkk - addFinalizer lc $ modifyMVar mds (\im' -> return (IM.delete key im', ())) - addFinalizer lm $ CUDA.unload mdl - return ( IM.insert key lm im, lm ) - - -- Return existing module - Just lm -> return (im, lm) - -- - withLifetime' mdl k - - -toKey :: CUDA.Context -> IM.Key -toKey (CUDA.Context (Ptr addr#)) = I# (addr2Int# addr#) - - -launch :: Kernel -> CUDA.Stream -> Int -> DevicePtr e -> DevicePtr e -> DevicePtr e -> IO () -launch Kernel{..} s n dx dy dz = - CUDA.launchKernel kernelFun (kernelThreadBlocks n,1,1) (kernelThreadBlockSize,1,1) kernelSharedMemBytes (Just s) - [ CUDA.VArg dx, CUDA.VArg dy, CUDA.VArg dz, CUDA.IArg (fromIntegral n) ] - -mkKernel :: String -> CUDA.Module -> CUDA.DeviceProperties -> IO Kernel -mkKernel name mdl prp = do - fun <- CUDA.getFun mdl name - reg <- CUDA.requires fun CUDA.NumRegs - let - blockSize = 256 - sharedMem = 0 - maxBlocks = maxResidentBlocks prp blockSize reg sharedMem - numBlocks n = maxBlocks `min` ((n + blockSize - 1) `quot` blockSize) - -- - return $ Kernel fun sharedMem blockSize numBlocks name - -data Kernel = Kernel { - kernelFun :: {-# UNPACK #-} !CUDA.Fun - , kernelSharedMemBytes :: {-# UNPACK #-} !Int - , kernelThreadBlockSize :: {-# UNPACK #-} !Int - , kernelThreadBlocks :: (Int -> Int) - , kernelName :: String - } - -modules :: NumericR (Complex e) -> MVar (IntMap (Lifetime (CUDA.Module, Kernel, Kernel))) -modules NumericRcomplex32 = modules_f32 -modules NumericRcomplex64 = modules_f64 - -{-# NOINLINE modules_f32 #-} -modules_f32 :: MVar (IntMap (Lifetime (CUDA.Module, Kernel, Kernel))) -modules_f32 = unsafePerformIO $ newMVar IM.empty - -{-# NOINLINE modules_f64 #-} -modules_f64 :: MVar (IntMap (Lifetime (CUDA.Module, Kernel, Kernel))) -modules_f64 = unsafePerformIO $ newMVar IM.empty - -ptx_twine_f32 :: ByteString -ptx_twine_f32 = $(makeRelativeToProject "cubits/twine_f32.ptx" >>= embedFile) - -ptx_twine_f64 :: ByteString -ptx_twine_f64 = $(makeRelativeToProject "cubits/twine_f64.ptx" >>= embedFile) - diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/Type.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/Type.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/LinearAlgebra/Type.hs +++ /dev/null @@ -1,91 +0,0 @@ -{-# LANGUAGE ConstraintKinds #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE GADTs #-} -{-# LANGUAGE TypeFamilies #-} --- | --- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.Type --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- - -module Data.Array.Accelerate.Numeric.LinearAlgebra.Type - where - -import Data.Array.Accelerate as A -import Data.Array.Accelerate.Data.Complex as A - -import qualified Prelude as P - - --- For explicit dictionary reification, to recover the type the operation should --- be performed at. --- -data NumericR a where - NumericRfloat32 :: NumericR Float - NumericRfloat64 :: NumericR Double - NumericRcomplex32 :: NumericR (Complex Float) - NumericRcomplex64 :: NumericR (Complex Double) - -class (Elt a, Num a) => Numeric a where - numericR :: NumericR a - -instance Numeric Float where - numericR = NumericRfloat32 - -instance Numeric Double where - numericR = NumericRfloat64 - -instance Numeric (Complex Float) where - numericR = NumericRcomplex32 - -instance Numeric (Complex Double) where - numericR = NumericRcomplex64 - --- class Numeric a => RealNumeric a --- --- instance RealNumeric Float --- instance RealNumeric Double - -type family NumericBaseT t where - NumericBaseT Float = Float - NumericBaseT Double = Double - NumericBaseT (Complex Float) = Float - NumericBaseT (Complex Double) = Double - - --- | Matrices as dense two-dimensional arrays in row-major ordering --- -type Matrix e = Array DIM2 e - --- | Orientation of the underlying data. --- --- Accelerate arrays are naturally stored in row-major format. --- -data Orientation - = R -- ^ row major - | C -- ^ column major - deriving (P.Eq, P.Show) - --- | Many operations allow you to implicitly transpose the arguments. For --- a given input matrix @mat@ with dimensions @Z :. m :. n@ (that is; @m@ rows --- and @n@ columns): --- -data Transpose - -- | Leave the matrix as is. - = N - - -- | Treat the matrix as implicitly transposed, with dimensions @Z :. n :. m@. - -- Entry @Z :. j :. i@ is treated as actually being entry @Z :. i :. j@. - | T - - -- | Implicitly transpose and conjugate the input matrix. For complex-valued - -- matrices a given element @mat ! Z:.j:.i == x :+ y@ will be treated as - -- actually being @mat ! Z:.i:.j == x :+ (-y)@. - | H - deriving (P.Eq, P.Show) - diff --git a/Data/Array/Accelerate/Numeric/Sum.hs b/Data/Array/Accelerate/Numeric/Sum.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/Sum.hs +++ /dev/null @@ -1,330 +0,0 @@ -{-# LANGUAGE ConstraintKinds #-} -{-# LANGUAGE DeriveDataTypeable #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE MultiParamTypeClasses #-} -{-# LANGUAGE RebindableSyntax #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE TypeOperators #-} -{-# LANGUAGE ViewPatterns #-} --- | --- Module : Data.Array.Accelerate.Numeric.Sum --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- --- Functions for summing floating point numbers more accurately than the --- straightforward 'Data.Array.Accelerate.sum' operation. --- --- In the worst case, the 'Data.Array.Accelerate.sum' function accumulates error --- at a rate proportional to the number of values being summed. The algorithms --- in this module implement different methods of /compensated summation/, which --- reduce the accumulation of numeric error so that it grows much more slowly --- than the number of inputs (e.g. logarithmically), or remains constant. --- - --- TLM: The standard formulation of the algorithms implemented here are not --- associative; e.g. they would have a type (KBN a -> a -> KBN a). I've --- done what seems like the sensible conversion, but somebody versed in numeric --- analysis should probably look... --- --- See also: --- - -module Data.Array.Accelerate.Numeric.Sum ( - - -- * Summation type class - Summation(..), - sum, - - -- * Kahan-Babuška-Neumaier summation - KBN(..), - kbn, - - -- * Order-2 Kahan-Babuška summation - KB2(..), - kb2, - - -- * Kahan summation - Kahan(..), - kahan, - -) where - -import Data.Array.Accelerate as A hiding ( sum ) -import Data.Array.Accelerate.Type as A -import Data.Array.Accelerate.Smart as A ( Exp(..), PreExp(..) ) -import Data.Array.Accelerate.Product as A -import Data.Array.Accelerate.Array.Sugar as A -import Data.Array.Accelerate.Numeric.Sum.Arithmetic as A - -import Data.Proxy -import Data.Typeable -import Prelude ( Show, fromInteger ) - - --- | Sum an array using a particular compensation scheme. --- --- >>> let xs = [1.0, 1.0e100, 1.0, -1.0e100] :: [Double] --- >>> Prelude.sum xs --- 0.0 --- --- >>> let ys = fromList (Z:.4) [1.0, 1.0e100, 1.0, -1.0e100] :: Vector Double --- >>> sum kbn (use ys) --- Scalar Z [2.0] --- -sum :: (Summation s a, Shape sh) => Proxy s -> Acc (Array (sh:.Int) a) -> Acc (Array sh a) -sum p = A.map (from p) - . A.fold add zero - . A.map (into p) - - --- | A class for the summation of floating-point numbers --- -class (Elt a, Elt (s a)) => Summation s a where - -- | Add a value to the sum - add :: Exp (s a) -> Exp (s a) -> Exp (s a) - - -- | The identity of the summation - zero :: Exp (s a) - - -- | Insert a value into the summation - into :: Proxy s -> Exp a -> Exp (s a) - - -- | Summarise the result of summation - from :: Proxy s -> Exp (s a) -> Exp a - - --- | Kahan-Babuška-Neumaier summation. This is a little more computationally --- costly than plain Kahan summation, but is /always/ at least as accurate. --- -data KBN a = KBN a a - deriving (Show, Typeable) - --- | Return the result of a Kahan-Babuška-Neumaier sum. --- -kbn :: Proxy KBN -kbn = Proxy - -kbnAdd :: (Num a, Ord a, IsFloating a) => Exp (KBN a) -> Exp (KBN a) -> Exp (KBN a) -kbnAdd (unlift -> KBN s1 c1) (unlift -> KBN s2 c2) = lift (KBN s' c') - where - s' = s1 `fadd` s2 - c' = c1 `fadd` c2 `fadd` if abs s1 >= abs s2 - then (s1 `fsub` s') `fadd` s2 - else (s2 `fsub` s') `fadd` s1 - --- instance (Num a, Ord a) => Summation KBN a where --- zero = lift $ KBN (0::Exp a) (0::Exp a) --- add = kbnAdd --- into _ x = lift (KBN x 0) --- from _ x = let KBN s c = unlift x in s + c - -instance Summation KBN Float where - zero = constant (KBN 0 0) - add = kbnAdd - into _ x = lift (KBN x 0) - from _ x = let KBN s c = unlift x in s + c - -instance Summation KBN Double where - zero = constant (KBN 0 0) - add = kbnAdd - into _ x = lift (KBN x 0) - from _ x = let KBN s c = unlift x in s + c - -instance Summation KBN CFloat where - zero = constant (KBN 0 0) - add = kbnAdd - into _ x = lift (KBN x 0) - from _ x = let KBN s c = unlift x in s + c - -instance Summation KBN CDouble where - zero = constant (KBN 0 0) - add = kbnAdd - into _ x = lift (KBN x 0) - from _ x = let KBN s c = unlift x in s + c - -type instance EltRepr (KBN a) = (((), EltRepr a), EltRepr a) - -instance Elt a => Elt (KBN a) where - eltType _ = UnitTuple `PairTuple` eltType (undefined::a) - `PairTuple` eltType (undefined::a) - toElt (((),a),b) = KBN (toElt a) (toElt b) - fromElt (KBN a b) = (((), fromElt a), fromElt b) - -instance Elt a => IsProduct Elt (KBN a) where - type ProdRepr (KBN a) = (((), a), a) - toProd _ (((),a),b) = KBN a b - fromProd _ (KBN a b) = (((),a),b) - prod _ _ = ProdRsnoc $ ProdRsnoc ProdRunit - -instance (Lift Exp a, Elt (Plain a)) => Lift Exp (KBN a) where - type Plain (KBN a) = KBN (Plain a) - lift (KBN a b) = Exp $ Tuple $ NilTup `SnocTup` lift a - `SnocTup` lift b - -instance Elt a => Unlift Exp (KBN (Exp a)) where - unlift t = KBN (Exp $ SuccTupIdx ZeroTupIdx `Prj` t) - (Exp $ ZeroTupIdx `Prj` t) - - --- | Second-order Kahan-Babuška summation. This is more computationally costly --- than Kahan-Babuška-Neumaier summation. Its advantage is that it can lose less --- precision (in admittedly obscure cases). --- --- This method compensates for error in both the sum and the first-order --- compensation term, hence the use of \"second order\" in the name. --- -data KB2 a = KB2 a a a - deriving (Show, Typeable) - --- | Return the result of a second-order Kahan-Babuška sum. --- -kb2 :: Proxy KB2 -kb2 = Proxy - -kb2Add :: (Num a, Ord a, IsFloating a) => Exp (KB2 a) -> Exp (KB2 a) -> Exp (KB2 a) -kb2Add (unlift -> KB2 s1 c1 cc1) (unlift -> KB2 s2 c2 cc2) = lift (KB2 sum' c' cc') - where - sum' = s1 `fadd` s2 - c' = t `fadd` k - cc' = cc1 `fadd` cc2 `fadd` if abs t >= abs k - then (t `fsub` c') `fadd` k - else (k `fsub` c') `fadd` t - t = c1 `fadd` c2 - k = if abs s1 >= abs s2 - then (s1 `fsub` sum') `fadd` s2 - else (s2 `fsub` sum') `fadd` s1 - --- instance (Num a, Ord a) => Summation KB2 a where --- zero = lift $ KB2 (0::Exp a) (0::Exp a) (0::Exp a) --- add = kb2Add --- into _ x = lift (KB2 x 0 0) --- from _ x = let KB2 s c cc = unlift x in s + c + cc - -instance Summation KB2 Float where - zero = constant (KB2 0 0 0) - add = kb2Add - into _ x = lift (KB2 x 0 0) - from _ x = let KB2 s c cc = unlift x in s + c + cc - -instance Summation KB2 Double where - zero = constant (KB2 0 0 0) - add = kb2Add - into _ x = lift (KB2 x 0 0) - from _ x = let KB2 s c cc = unlift x in s + c + cc - -instance Summation KB2 CFloat where - zero = constant (KB2 0 0 0) - add = kb2Add - into _ x = lift (KB2 x 0 0) - from _ x = let KB2 s c cc = unlift x in s + c + cc - -instance Summation KB2 CDouble where - zero = constant (KB2 0 0 0) - add = kb2Add - into _ x = lift (KB2 x 0 0) - from _ x = let KB2 s c cc = unlift x in s + c + cc - -type instance EltRepr (KB2 a) = ((((), EltRepr a), EltRepr a), EltRepr a) - -instance Elt a => Elt (KB2 a) where - eltType _ = UnitTuple `PairTuple` eltType (undefined::a) - `PairTuple` eltType (undefined::a) - `PairTuple` eltType (undefined::a) - toElt ((((),a),b),c) = KB2 (toElt a) (toElt b) (toElt c) - fromElt (KB2 a b c) = ((((), fromElt a), fromElt b), fromElt c) - -instance Elt a => IsProduct Elt (KB2 a) where - type ProdRepr (KB2 a) = ((((), a), a), a) - toProd _ ((((),a),b),c) = KB2 a b c - fromProd _ (KB2 a b c) = ((((),a),b),c) - prod _ _ = ProdRsnoc $ ProdRsnoc $ ProdRsnoc ProdRunit - -instance (Lift Exp a, Elt (Plain a)) => Lift Exp (KB2 a) where - type Plain (KB2 a) = KB2 (Plain a) - lift (KB2 a b c) = Exp $ Tuple $ NilTup `SnocTup` lift a - `SnocTup` lift b - `SnocTup` lift c - -instance Elt a => Unlift Exp (KB2 (Exp a)) where - unlift t = KB2 (Exp $ SuccTupIdx (SuccTupIdx ZeroTupIdx) `Prj` t) - (Exp $ SuccTupIdx ZeroTupIdx `Prj` t) - (Exp $ ZeroTupIdx `Prj` t) - - --- | Kahan summation. This is the least accurate of the compensated summation --- methods. This summation method is included only for completeness. --- -data Kahan a = Kahan a a - deriving (Show, Typeable) - --- | Return the result of a Kahan sum. --- -kahan :: Proxy Kahan -kahan = Proxy - -kahanAdd :: (Num a, IsFloating a) => Exp (Kahan a) -> Exp (Kahan a) -> Exp (Kahan a) -kahanAdd (unlift -> Kahan s1 c1 :: Kahan (Exp a)) (unlift -> Kahan s2 c2) = lift (Kahan s' c') - where - s' = s1 `fadd` y - c' = (s' `fsub` s1) `fsub` y - y = s2 `fsub` c1 `fsub` c2 - --- instance (Num a, Ord a) => Summation Kahan a where --- zero = lift $ Kahan (0::Exp a) (0::Exp a) --- add = kahanAdd --- into _ x = lift (Kahan x 0) --- from _ x = let Kahan s _ = unlift x in s - -instance Summation Kahan Float where - zero = constant (Kahan 0 0) - add = kahanAdd - into _ x = lift (Kahan x 0) - from _ x = let Kahan s _ = unlift x in s - -instance Summation Kahan Double where - zero = constant (Kahan 0 0) - add = kahanAdd - into _ x = lift (Kahan x 0) - from _ x = let Kahan s _ = unlift x in s - -instance Summation Kahan CFloat where - zero = constant (Kahan 0 0) - add = kahanAdd - into _ x = lift (Kahan x 0) - from _ x = let Kahan s _ = unlift x in s - -instance Summation Kahan CDouble where - zero = constant (Kahan 0 0) - add = kahanAdd - into _ x = lift (Kahan x 0) - from _ x = let Kahan s _ = unlift x in s - -type instance EltRepr (Kahan a) = (((), EltRepr a), EltRepr a) - -instance Elt a => Elt (Kahan a) where - eltType _ = UnitTuple `PairTuple` eltType (undefined::a) - `PairTuple` eltType (undefined::a) - toElt (((),a),b) = Kahan (toElt a) (toElt b) - fromElt (Kahan a b) = (((), fromElt a), fromElt b) - -instance Elt a => IsProduct Elt (Kahan a) where - type ProdRepr (Kahan a) = (((), a), a) - toProd _ (((),a),b) = Kahan a b - fromProd _ (Kahan a b) = (((),a),b) - prod _ _ = ProdRsnoc $ ProdRsnoc ProdRunit - -instance (Lift Exp a, Elt (Plain a)) => Lift Exp (Kahan a) where - type Plain (Kahan a) = Kahan (Plain a) - lift (Kahan a b) = Exp $ Tuple $ NilTup `SnocTup` lift a - `SnocTup` lift b - -instance Elt a => Unlift Exp (Kahan (Exp a)) where - unlift t = Kahan (Exp $ SuccTupIdx ZeroTupIdx `Prj` t) - (Exp $ ZeroTupIdx `Prj` t) - diff --git a/Data/Array/Accelerate/Numeric/Sum/Arithmetic.hs b/Data/Array/Accelerate/Numeric/Sum/Arithmetic.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/Sum/Arithmetic.hs +++ /dev/null @@ -1,37 +0,0 @@ -{-# LANGUAGE ConstraintKinds #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE NoImplicitPrelude #-} --- | --- Module : Data.Array.Accelerate.Numeric.Sum.Arithmetic --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- - -module Data.Array.Accelerate.Numeric.Sum.Arithmetic ( - - fadd, fsub, fmul, - -) where - -import Data.Array.Accelerate - -import qualified Data.Array.Accelerate.Numeric.Sum.LLVM.Native as Native -import qualified Data.Array.Accelerate.Numeric.Sum.LLVM.PTX as PTX - - -infixl 6 `fadd` -fadd :: (Num a, IsFloating a) => Exp a -> Exp a -> Exp a -fadd = Native.fadd $ PTX.fadd (+) - -infixl 6 `fsub` -fsub :: (Num a, IsFloating a) => Exp a -> Exp a -> Exp a -fsub = Native.fsub $ PTX.fsub (-) - -infixl 7 `fmul` -fmul :: (Num a, IsFloating a) => Exp a -> Exp a -> Exp a -fmul = Native.fmul $ PTX.fmul (*) - diff --git a/Data/Array/Accelerate/Numeric/Sum/LLVM/Native.hs b/Data/Array/Accelerate/Numeric/Sum/LLVM/Native.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/Sum/LLVM/Native.hs +++ /dev/null @@ -1,58 +0,0 @@ -{-# LANGUAGE CPP #-} --- | --- Module : Data.Array.Accelerate.Numeric.Sum.LLVM.Native --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- - -module Data.Array.Accelerate.Numeric.Sum.LLVM.Native ( - - fadd, fsub, fmul, - -) where - -import Data.Array.Accelerate as A -import Data.Array.Accelerate.Type - -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND -import Data.Array.Accelerate.LLVM.CodeGen.Sugar -import Data.Array.Accelerate.LLVM.Native.Foreign as A -import qualified Data.Array.Accelerate.Numeric.Sum.LLVM.Prim as Prim -#endif - -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND -wrap2 :: (Elt a, Elt b, Elt c) - => String -- name of the operation - -> IRFun1 Native () ((a, b) -> c) -- foreign implementation - -> (Exp a -> Exp b -> Exp c) -- fallback implementation - -> Exp a - -> Exp b - -> Exp c -wrap2 str f g = A.curry (foreignExp (ForeignExp str f) (A.uncurry g)) -#endif - -fadd :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND -fadd = wrap2 "fadd" (Prim.fadd floatingType) -#else -fadd = id -#endif - -fsub :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND -fsub = wrap2 "fsub" (Prim.fsub floatingType) -#else -fsub = id -#endif - -fmul :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND -fmul = wrap2 "fmul" (Prim.fmul floatingType) -#else -fmul = id -#endif - diff --git a/Data/Array/Accelerate/Numeric/Sum/LLVM/PTX.hs b/Data/Array/Accelerate/Numeric/Sum/LLVM/PTX.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/Sum/LLVM/PTX.hs +++ /dev/null @@ -1,58 +0,0 @@ -{-# LANGUAGE CPP #-} --- | --- Module : Data.Array.Accelerate.Numeric.Sum.LLVM.PTX --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- - -module Data.Array.Accelerate.Numeric.Sum.LLVM.PTX ( - - fadd, fsub, fmul, - -) where - -import Data.Array.Accelerate as A -import Data.Array.Accelerate.Type - -#ifdef ACCELERATE_LLVM_PTX_BACKEND -import Data.Array.Accelerate.LLVM.CodeGen.Sugar -import Data.Array.Accelerate.LLVM.PTX.Foreign as A -import qualified Data.Array.Accelerate.Numeric.Sum.LLVM.Prim as Prim -#endif - -#ifdef ACCELERATE_LLVM_PTX_BACKEND -wrap2 :: (Elt a, Elt b, Elt c) - => String -- name of the operation - -> IRFun1 PTX () ((a, b) -> c) -- foreign implementation - -> (Exp a -> Exp b -> Exp c) -- fallback implementation - -> Exp a - -> Exp b - -> Exp c -wrap2 str f g = A.curry (foreignExp (ForeignExp str f) (A.uncurry g)) -#endif - -fadd :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a -#ifdef ACCELERATE_LLVM_PTX_BACKEND -fadd = wrap2 "fadd" (Prim.fadd floatingType) -#else -fadd = id -#endif - -fsub :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a -#ifdef ACCELERATE_LLVM_PTX_BACKEND -fsub = wrap2 "fsub" (Prim.fsub floatingType) -#else -fsub = id -#endif - -fmul :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a -#ifdef ACCELERATE_LLVM_PTX_BACKEND -fmul = wrap2 "fmul" (Prim.fmul floatingType) -#else -fmul = id -#endif - diff --git a/Data/Array/Accelerate/Numeric/Sum/LLVM/Prim.hs b/Data/Array/Accelerate/Numeric/Sum/LLVM/Prim.hs deleted file mode 100644 --- a/Data/Array/Accelerate/Numeric/Sum/LLVM/Prim.hs +++ /dev/null @@ -1,83 +0,0 @@ -{-# LANGUAGE GADTs #-} -{-# LANGUAGE TemplateHaskell #-} -{-# LANGUAGE ViewPatterns #-} --- | --- Module : Data.Array.Accelerate.Numeric.Sum.LLVM.Prim --- Copyright : [2017] Trevor L. McDonell --- License : BSD3 --- --- Maintainer : Trevor L. McDonell --- Stability : experimental --- Portability : non-portable (GHC extensions) --- - -module Data.Array.Accelerate.Numeric.Sum.LLVM.Prim ( - - fadd, fsub, fmul, - -) where - -import Data.Array.Accelerate.Type -import Data.Array.Accelerate.Error - -import Data.Array.Accelerate.LLVM.CodeGen.Downcast ( downcast ) -import Data.Array.Accelerate.LLVM.CodeGen.IR ( IR(..), Operands(..), IROP(..) ) -import Data.Array.Accelerate.LLVM.CodeGen.Monad ( CodeGen, freshName, instr_ ) -import Data.Array.Accelerate.LLVM.CodeGen.Sugar ( IROpenFun1(..) ) -import qualified Data.Array.Accelerate.LLVM.CodeGen.Arithmetic as A -import qualified LLVM.AST.Type.Name as A -import qualified LLVM.AST.Type.Operand as A -import qualified LLVM.AST.Type.Representation as A - -import LLVM.AST.Instruction -import LLVM.AST.Name -import LLVM.AST.Operand -import LLVM.AST.Type - - --- | As (+), but don't allow potentially unsafe floating-point optimisations. --- -fadd :: FloatingType a -> IROpenFun1 arch env aenv ((a,a) -> a) -fadd t = IRFun1 $ A.uncurry (binop FAdd t) - --- | As (-), but don't allow potentially unsafe floating-point optimisations. --- -fsub :: FloatingType a -> IROpenFun1 arch env aenv ((a,a) -> a) -fsub t = IRFun1 $ A.uncurry (binop FSub t) - --- | As (*), but don't allow potentially unsafe floating-point optimisations. --- -fmul :: FloatingType a -> IROpenFun1 arch env aenv ((a,a) -> a) -fmul t = IRFun1 $ A.uncurry (binop FMul t) - -binop :: (FastMathFlags -> Operand -> Operand -> InstructionMetadata -> Instruction) -> FloatingType a -> IR a -> IR a -> CodeGen (IR a) -binop f t (op t -> x) (op t -> y) = do - r <- instr (downcast t) (f fmf (downcast x) (downcast y) md) - return (upcast t r) - - --- Prim --- ---- - -md :: InstructionMetadata -md = [] - -fmf :: FastMathFlags -fmf = NoFastMathFlags - -fresh :: CodeGen Name -fresh = downcast <$> freshName - -instr :: Type -> Instruction -> CodeGen Operand -instr ty ins = do - name <- fresh - instr_ (name := ins) - return (LocalReference ty name) - -upcast :: FloatingType t -> Operand -> IR t -upcast TypeFloat{} (LocalReference (FloatingPointType FloatFP) (UnName x)) = IR $ OP_Float (A.LocalReference A.type' (A.UnName x)) -upcast TypeDouble{} (LocalReference (FloatingPointType DoubleFP) (UnName x)) = IR $ OP_Double (A.LocalReference A.type' (A.UnName x)) -upcast TypeCFloat{} (LocalReference (FloatingPointType FloatFP) (UnName x)) = IR $ OP_CFloat (A.LocalReference A.type' (A.UnName x)) -upcast TypeCDouble{} (LocalReference (FloatingPointType DoubleFP) (UnName x)) = IR $ OP_CDouble (A.LocalReference A.type' (A.UnName x)) -upcast _ _ = $internalError "upcast" "expected local reference" - diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -18,12 +18,6 @@ * **accelerate-llvm-ptx:** FFI bindings to the NVIDIA [cuBLAS] library. -## Complex numbers - -Due to Accelerate's struct-of-array representation of complex numbers, compared -to the C-style array-of-struct representation, calling foreign implementations -of complex-valued operations entails an extra data marshalling step. - [GitHub]: https://github.com/AccelerateHS/accelerate [blas-hs]: http://hackage.haskell.org/package/blas-hs diff --git a/accelerate-blas.cabal b/accelerate-blas.cabal --- a/accelerate-blas.cabal +++ b/accelerate-blas.cabal @@ -1,5 +1,5 @@ name: accelerate-blas -version: 0.1.0.1 +version: 0.2.0.0 synopsis: Numeric Linear Algebra in Accelerate description: Linear systems, matrix decompositions, and other numerical computations for @@ -21,8 +21,6 @@ extra-source-files: README.md CHANGELOG.md - cubits/twine_f32.ptx - cubits/twine_f64.ptx Flag llvm-cpu Description: Enable the LLVM backend for multicore CPUs @@ -48,9 +46,12 @@ Data.Array.Accelerate.Numeric.Sum.LLVM.PTX build-depends: - base >= 4.7 && < 4.11 - , accelerate >= 1.0 && < 1.2 + base >= 4.7 && < 4.12 + , accelerate >= 1.0 && < 1.3 + hs-source-dirs: + src + ghc-options: -O2 -Wall @@ -58,20 +59,10 @@ if flag(llvm-cpu) CPP-options: -DACCELERATE_LLVM_NATIVE_BACKEND build-depends: - accelerate-llvm >= 1.0 && < 1.2 - , accelerate-llvm-native >= 1.0 && < 1.2 + accelerate-llvm >= 1.1 && < 1.3 + , accelerate-llvm-native >= 1.1 && < 1.3 , blas-hs >= 0.1 - , llvm-hs-pure >= 4.0 - , storable-complex >= 0.2 - - cc-options: - -O3 - -Wall - -march=native - - c-sources: - cbits/twine_f32.c - cbits/twine_f64.c + , llvm-hs-pure >= 4.1 && < 6.1 other-modules: Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Base @@ -79,48 +70,105 @@ Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level3 Data.Array.Accelerate.Numeric.Sum.LLVM.Prim - if flag(llvm-ptx) CPP-options: -DACCELERATE_LLVM_PTX_BACKEND build-depends: - accelerate-llvm >= 1.0 && < 1.2 - , accelerate-llvm-ptx >= 1.0 && < 1.2 + accelerate-llvm >= 1.1 && < 1.3 + , accelerate-llvm-ptx >= 1.1 && < 1.3 , bytestring >= 0.9 , containers >= 0.5 , cublas >= 0.3 , cuda >= 0.8 , file-embed >= 0.0.10 - , llvm-hs-pure >= 4.0 + , llvm-hs-pure >= 4.1 && < 6.1 , mtl >= 2.2 - , storable-complex >= 0.2 other-modules: Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Context - Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Twine Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level2 Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level3 Data.Array.Accelerate.Numeric.Sum.LLVM.Prim -test-suite accelerate-blas-test +test-suite test-llvm-native + type: exitcode-stdio-1.0 default-language: Haskell2010 + hs-source-dirs: test + main-is: TestNative.hs + ghc-options: -main-is TestNative + + if !flag(llvm-cpu) + buildable: False + + build-depends: + base >= 4.7 && < 4.12 + , accelerate >= 1.2 + , accelerate-blas + , accelerate-llvm-native + , hedgehog >= 0.5 + , tasty >= 0.11 + , tasty-hedgehog >= 0.1 + + ghc-options: + -Wall + -threaded + -rtsopts + + other-modules: + Test.BLAS + Test.BLAS.Level2 + Test.BLAS.Level3 + Test.Util + +test-suite test-llvm-ptx type: exitcode-stdio-1.0 + default-language: Haskell2010 hs-source-dirs: test - main-is: Main.hs + main-is: TestPTX.hs + ghc-options: -main-is TestPTX + + if !flag(llvm-ptx) + buildable: False + + build-depends: + base >= 4.7 && < 4.12 + , accelerate >= 1.2 + , accelerate-blas + , accelerate-llvm-ptx + , hedgehog >= 0.5 + , tasty >= 0.11 + , tasty-hedgehog >= 0.1 + + ghc-options: + -Wall + -threaded + -rtsopts + other-modules: - Backend - Hedgehog.Gen.Array - Hedgehog.Gen.Shape - Level2 - Level3 - Similar + Test.BLAS + Test.BLAS.Level2 + Test.BLAS.Level3 + Test.Util + +benchmark bench-hmatrix + default-language: Haskell2010 + type: exitcode-stdio-1.0 + hs-source-dirs: bench + main-is: BenchHMatrix.hs + ghc-options: -main-is BenchHMatrix + + -- don't bother if we aren't building one of the real backends + if !flag(llvm-cpu) && !flag(llvm-ptx) + buildable: False + build-depends: - base >= 4.7 && < 4.11 - , accelerate >= 1.0 && < 1.2 - , accelerate-blas - , hedgehog >= 0.5 + base >= 4.7 && < 4.12 + , criterion >= 1.0 + , mwc-random >= 0.8 + , deepseq >= 1.0 + , hmatrix >= 0.17 ghc-options: -O2 @@ -129,53 +177,68 @@ -threaded -with-rtsopts=-N - if flag(llvm-cpu) - CPP-options: -DACCELERATE_LLVM_NATIVE_BACKEND - build-depends: - accelerate-llvm-native >= 1.0 && < 1.2 + other-modules: + Bench.Util + Bench.HMatrix - if flag(llvm-ptx) - CPP-options: -DACCELERATE_LLVM_PTX_BACKEND - build-depends: - accelerate-llvm-ptx >= 1.0 && < 1.2 +benchmark bench-llvm-native + default-language: Haskell2010 + type: exitcode-stdio-1.0 + hs-source-dirs: bench + main-is: BenchNative.hs + ghc-options: -main-is BenchNative + if !flag(llvm-cpu) + buildable: False -benchmark accelerate-blas-bench + build-depends: + base >= 4.7 && < 4.12 + , accelerate + , accelerate-blas + , accelerate-llvm-native + , criterion >= 1.0 + , mwc-random >= 0.8 + , mwc-random-accelerate >= 0.1 + + ghc-options: + -O2 + -Wall + -rtsopts + -threaded + + other-modules: + Bench.Util + Bench.Accelerate + +benchmark bench-llvm-ptx default-language: Haskell2010 type: exitcode-stdio-1.0 hs-source-dirs: bench - main-is: Main.hs - other-modules: - Accelerate - Extra - HMatrix + main-is: BenchPTX.hs + ghc-options: -main-is BenchPTX + if !flag(llvm-ptx) + buildable: False + build-depends: - base >= 4.7 && < 4.11 - , accelerate >= 1.0 && < 1.2 + base >= 4.7 && < 4.12 + , accelerate , accelerate-blas + , accelerate-llvm-ptx , criterion >= 1.0 , mwc-random >= 0.8 , mwc-random-accelerate >= 0.1 - , deepseq >= 1.0 - , hmatrix >= 0.17 ghc-options: -O2 -Wall -rtsopts -threaded - -with-rtsopts=-N - if flag(llvm-cpu) - CPP-options: -DACCELERATE_LLVM_NATIVE_BACKEND - build-depends: - accelerate-llvm-native >= 1.0 && < 1.2 + other-modules: + Bench.Util + Bench.Accelerate - if flag(llvm-ptx) - CPP-options: -DACCELERATE_LLVM_PTX_BACKEND - build-depends: - accelerate-llvm-ptx >= 1.0 && < 1.2 source-repository head type: git @@ -183,7 +246,7 @@ source-repository this type: git - tag: 0.1.0.1 + tag: 0.2.0.0 location: https://github.com/tmcdonell/accelerate-blas -- vim: nospell diff --git a/bench/Accelerate.hs b/bench/Accelerate.hs deleted file mode 100644 --- a/bench/Accelerate.hs +++ /dev/null @@ -1,187 +0,0 @@ -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE CPP #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE PolyKinds #-} -{-# LANGUAGE RankNTypes #-} - -module Accelerate ( - - Backend(..), - benchAcc, - -) where - -import Extra - -import Data.Array.Accelerate ( Acc, Arrays, Elt, Z(..), (:.)(..) ) -import Data.Array.Accelerate.Numeric.LinearAlgebra -import Data.Array.Accelerate.Data.Complex -import Data.Array.Accelerate.System.Random.MWC -import qualified Data.Array.Accelerate as A -import qualified Data.Array.Accelerate.Interpreter as I -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND -import qualified Data.Array.Accelerate.LLVM.Native as CPU -#endif -#ifdef ACCELERATE_LLVM_PTX_BACKEND -import qualified Data.Array.Accelerate.LLVM.PTX as PTX -#endif - -import Criterion.Main -import Data.Proxy -import Text.Printf - - -benchAcc :: Backend -> Benchmark -benchAcc backend = - bgroup (show backend) - [ level2 backend - , level3 backend - ] - - -level2 :: Backend -> Benchmark -level2 backend = - bgroup "matrix-vector" - [ bgroup "(#>)" - [ gemv 200 400 - , gemv 500 1000 - , gemv 1000 2000 - , gemv 2000 3000 - ] - , bgroup "(<#)" - [ gevm 200 400 - , gevm 500 1000 - , gevm 1000 2000 - , gevm 2000 3000 - ] - ] - where - gemv :: Int -> Int -> Benchmark - gemv m n = - let complexity = m * n - - setup :: (Variate e, Elt e) => proxy e -> IO (Matrix e, Vector e) - setup _ = withSystemRandom $ \gen -> do - matA <- randomArrayWith gen uniform (Z :. m :. n) - vecx <- randomArrayWith gen uniform (Z :. n) - return (matA, vecx) - - go :: (Variate e, Numeric e, Show (ArgType e)) => proxy e -> Benchmark - go t = env (setup t) - $ \ ~(matA, vecx) -> bench (showType t) - $ whnf (run2 backend (#>) matA) vecx - in - bgroup (printf "%dx%d" m n) (sdcz go complexity backend) - - gevm :: Int -> Int -> Benchmark - gevm m n = - let complexity = m * n - - setup :: (Variate e, Elt e) => proxy e -> IO (Matrix e, Vector e) - setup _ = withSystemRandom $ \gen -> do - matA <- randomArrayWith gen uniform (Z :. m :. n) - vecx <- randomArrayWith gen uniform (Z :. m) - return (matA, vecx) - - go :: (Variate e, Numeric e, Show (ArgType e)) => proxy e -> Benchmark - go t = env (setup t) - $ \ ~(matA, vecx) -> bench (showType t) - $ whnf (run2 backend (<#) vecx) matA - in - bgroup (printf "%dx%d" m n) (sdcz go complexity backend) - - -level3 :: Backend -> Benchmark -level3 backend = - bgroup "matrix-matrix" - [ bgroup "(<>)" - [ gemm 100 100 100 - , gemm 250 250 250 - , gemm 500 500 500 - , gemm 1000 1000 1000 - ] - ] - where - gemm :: Int -> Int -> Int -> Benchmark - gemm m n k = - let complexity = m * n * k - - setup :: (Variate e, Elt e) => proxy e -> IO (Matrix e, Matrix e) - setup _ = withSystemRandom $ \gen -> do - matA <- randomArrayWith gen uniform (Z :. m :. k) - matB <- randomArrayWith gen uniform (Z :. k :. n) - return (matA, matB) - - go :: (Variate e, Numeric e, Show (ArgType e)) => proxy e -> Benchmark - go t = env (setup t) - $ \ ~(matA, matB) -> bench (showType t) - $ whnf (run2 backend (<>) matA) matB - in - bgroup (printf "%dx%dx%d" m n k) (sdcz go complexity backend) - - -sdcz :: (forall (e :: *). (Variate e, Numeric e, Show (ArgType e)) => Proxy e -> Benchmark) - -> Int - -> Backend - -> [Benchmark] -sdcz go complexity backend = - if maybe True (complexity <=) (complexityLimit backend) - then - [ go (Proxy :: Proxy Float) - , go (Proxy :: Proxy Double) - , go (Proxy :: Proxy (Complex Float)) - , go (Proxy :: Proxy (Complex Double)) - ] - else - [] - -complexityLimit :: Backend -> Maybe Int -complexityLimit Interpreter = Just 50000 -complexityLimit _ = Nothing - - -data Backend = Interpreter -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND - | Native -#endif -#ifdef ACCELERATE_LLVM_PTX_BACKEND - | PTX -#endif - -instance Show Backend where - show Interpreter = "interpreter" -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND - show Native = "llvm-cpu" -#endif -#ifdef ACCELERATE_LLVM_PTX_BACKEND - show PTX = "llvm-ptx" -#endif - -{-# INLINE run #-} -run :: Arrays a => Backend -> Acc a -> a -run Interpreter = I.run -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND -run Native = CPU.run -#endif -#ifdef ACCELERATE_LLVM_PTX_BACKEND -run PTX = PTX.run -#endif - - -{-# INLINE run1 #-} -run1 :: (Arrays a, Arrays b) => Backend -> (Acc a -> Acc b) -> a -> b -run1 Interpreter f = I.run1 f -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND -run1 Native f = CPU.run1 f -#endif -#ifdef ACCELERATE_LLVM_PTX_BACKEND -run1 PTX f = PTX.run1 f -#endif - -{-# INLINE run2 #-} -run2 :: (Arrays a, Arrays b, Arrays c) => Backend -> (Acc a -> Acc b -> Acc c) -> a -> b -> c -run2 b f x y = go (x,y) - where - !go = run1 b (A.uncurry f) - diff --git a/bench/Bench/Accelerate.hs b/bench/Bench/Accelerate.hs new file mode 100644 --- /dev/null +++ b/bench/Bench/Accelerate.hs @@ -0,0 +1,123 @@ +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +-- | +-- Module : Bench.Accelerate +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Bench.Accelerate ( bench_accelerate ) + where + +import Bench.Util + +import Data.Array.Accelerate ( Elt, Z(..), (:.)(..) ) +import Data.Array.Accelerate.Trafo ( Afunction, AfunctionR ) +import Data.Array.Accelerate.Numeric.LinearAlgebra +import Data.Array.Accelerate.Data.Complex +import Data.Array.Accelerate.System.Random.MWC + +import Criterion.Main +import Data.Proxy +import Text.Printf + +type RunN = forall f. Afunction f => f -> AfunctionR f + +bench_accelerate :: RunN -> [Benchmark] +bench_accelerate runN = + [ bench_level2 runN + , bench_level3 runN + ] + + +bench_level2 :: RunN -> Benchmark +bench_level2 runN = + bgroup "matrix-vector" + [ bgroup "(#>)" + [ gemv 200 400 + , gemv 500 1000 + , gemv 1000 2000 + , gemv 2000 3000 + ] + , bgroup "(<#)" + [ gevm 200 400 + , gevm 500 1000 + , gevm 1000 2000 + , gevm 2000 3000 + ] + ] + where + gemv :: Int -> Int -> Benchmark + gemv m n = + let setup :: (Variate e, Elt e) => proxy e -> IO (Matrix e, Vector e) + setup _ = withSystemRandom $ \gen -> do + matA <- randomArrayWith gen uniform (Z :. m :. n) + vecx <- randomArrayWith gen uniform (Z :. n) + return (matA, vecx) + + go :: (Variate e, Numeric e, Show (ArgType e)) => proxy e -> Benchmark + go t = env (setup t) + $ \ ~(matA, vecx) -> bench (showType t) + $ whnf (runN (#>) matA) vecx + in + bgroup (printf "%dx%d" m n) (sdcz go) + + gevm :: Int -> Int -> Benchmark + gevm m n = + let setup :: (Variate e, Elt e) => proxy e -> IO (Matrix e, Vector e) + setup _ = withSystemRandom $ \gen -> do + matA <- randomArrayWith gen uniform (Z :. m :. n) + vecx <- randomArrayWith gen uniform (Z :. m) + return (matA, vecx) + + go :: (Variate e, Numeric e, Show (ArgType e)) => proxy e -> Benchmark + go t = env (setup t) + $ \ ~(matA, vecx) -> bench (showType t) + $ whnf (runN (<#) vecx) matA + in + bgroup (printf "%dx%d" m n) (sdcz go) + + +bench_level3 :: RunN -> Benchmark +bench_level3 runN = + bgroup "matrix-matrix" + [ bgroup "(<>)" + [ gemm 100 100 100 + , gemm 250 250 250 + , gemm 500 500 500 + , gemm 1000 1000 1000 + ] + ] + where + gemm :: Int -> Int -> Int -> Benchmark + gemm m n k = + let setup :: (Variate e, Elt e) => proxy e -> IO (Matrix e, Matrix e) + setup _ = withSystemRandom $ \gen -> do + matA <- randomArrayWith gen uniform (Z :. m :. k) + matB <- randomArrayWith gen uniform (Z :. k :. n) + return (matA, matB) + + go :: (Variate e, Numeric e, Show (ArgType e)) => proxy e -> Benchmark + go t = env (setup t) + $ \ ~(matA, matB) -> bench (showType t) + $ whnf (runN (<>) matA) matB + in + bgroup (printf "%dx%dx%d" m n k) (sdcz go) + +sdcz :: (forall (e :: *). (Variate e, Numeric e, Show (ArgType e)) => Proxy e -> Benchmark) + -> [Benchmark] +sdcz go = + [ go (Proxy :: Proxy Float) + , go (Proxy :: Proxy Double) + , go (Proxy :: Proxy (Complex Float)) + , go (Proxy :: Proxy (Complex Double)) + ] + diff --git a/bench/Bench/HMatrix.hs b/bench/Bench/HMatrix.hs new file mode 100644 --- /dev/null +++ b/bench/Bench/HMatrix.hs @@ -0,0 +1,126 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE RankNTypes #-} +-- | +-- Module : Bench.HMatrix +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Bench.HMatrix ( bench_hmatrix ) + where + +import Bench.Util + +import Control.DeepSeq +import Criterion.Main +import Data.Proxy +import Foreign.Storable +import Numeric.LinearAlgebra hiding ( randomVector ) +import System.Random.MWC +import Text.Printf + + +bench_hmatrix :: [Benchmark] +bench_hmatrix = + [ bench_level2 + , bench_level3 + ] + +bench_level2 :: Benchmark +bench_level2 = + bgroup "matrix-vector" + [ bgroup "(#>)" + [ gemv 200 400 + , gemv 500 1000 + , gemv 1000 2000 + , gemv 2000 3000 + ] + , bgroup "(<#)" + [ gevm 200 400 + , gevm 500 1000 + , gevm 1000 2000 + , gevm 2000 3000 + ] + ] + where + gemv :: Int -> Int -> Benchmark + gemv m n = + let setup :: (Variate e, Storable e) => proxy e -> IO (Matrix e, Vector e) + setup _ = withSystemRandom $ \gen -> do + matA <- randomMatrix gen m n + vecx <- randomVector gen n + return (matA, vecx) + + go :: (Variate e, Numeric e, NFData e, Show (ArgType e)) => proxy e -> Benchmark + go t = env (setup t) + $ \ ~(matA, vecx) -> bench (showType t) + $ whnf (matA #>) vecx + in + bgroup (printf "%dx%d" m n) (sdcz go) + + gevm :: Int -> Int -> Benchmark + gevm m n = + let setup :: (Variate e, Storable e) => proxy e -> IO (Matrix e, Vector e) + setup _ = withSystemRandom $ \gen -> do + matA <- randomMatrix gen m n + vecx <- randomVector gen m + return (matA, vecx) + + go :: (Variate e, Numeric e, NFData e, Show (ArgType e)) => proxy e -> Benchmark + go t = env (setup t) + $ \ ~(matA, vecx) -> bench (showType t) + $ whnf (vecx <#) matA + in + bgroup (printf "%dx%d" m n) (sdcz go) + +bench_level3 :: Benchmark +bench_level3 = + bgroup "matrix-matrix" + [ bgroup "(<>)" + [ gemm 100 100 100 + , gemm 250 250 250 + , gemm 500 500 500 + , gemm 1000 1000 1000 + ] + ] + where + gemm :: Int -> Int -> Int -> Benchmark + gemm m n k = + let + setup :: (Variate e, Storable e) => proxy e -> IO (Matrix e, Matrix e) + setup _ = withSystemRandom $ \gen -> do + matA <- randomMatrix gen m k + matB <- randomMatrix gen k n + return (matA, matB) + + go :: (Variate e, Numeric e, NFData e, Show (ArgType e)) => proxy e -> Benchmark + go t = env (setup t) + $ \ ~(matA, matB) -> bench (showType t) + $ whnf (matA <>) matB + in + bgroup (printf "%dx%dx%d" m n k) (sdcz go) + + +randomVector :: (Variate e, Storable e) => GenIO -> Int -> IO (Vector e) +randomVector = uniformVector + +randomMatrix :: (Variate e, Storable e) => GenIO -> Int -> Int -> IO (Matrix e) +randomMatrix gen m n = do + v <- uniformVector gen (m * n) + return $ reshape n v + +sdcz :: (forall (e :: *). (Variate e, Numeric e, NFData e, Show (ArgType e)) => Proxy e -> Benchmark) + -> [Benchmark] +sdcz go = + [ go (Proxy :: Proxy Float) + , go (Proxy :: Proxy Double) + , go (Proxy :: Proxy (Complex Float)) + , go (Proxy :: Proxy (Complex Double)) + ] + diff --git a/bench/Bench/Util.hs b/bench/Bench/Util.hs new file mode 100644 --- /dev/null +++ b/bench/Bench/Util.hs @@ -0,0 +1,41 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE PolyKinds #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} +-- | +-- Module : Bench.Util +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Bench.Util where + +import Data.Complex +import System.Random.MWC + + +data ArgType (a :: *) = AT + +showType :: forall proxy a. Show (ArgType a) => proxy a -> String +showType _ = show (AT :: ArgType a) + +instance Show (ArgType Float) where show _ = "Float" +instance Show (ArgType Double) where show _ = "Double" +instance Show (ArgType (Complex Float)) where show _ = "Complex Float" +instance Show (ArgType (Complex Double)) where show _ = "Complex Double" + +instance Variate e => Variate (Complex e) where + uniform gen = (:+) <$> uniform gen <*> uniform gen + uniformR r gen = + let (ur:+ui,vr:+vi) = r + in (:+) <$> uniformR (ur,vr) gen <*> uniformR (ui,vi) gen + +infixr 0 $$ +($$) :: (b -> a) -> (c -> d -> b) -> c -> d -> a +(f $$ g) x y = f (g x y) + diff --git a/bench/BenchHMatrix.hs b/bench/BenchHMatrix.hs new file mode 100644 --- /dev/null +++ b/bench/BenchHMatrix.hs @@ -0,0 +1,18 @@ +-- | +-- Module : BenchHMatrix +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module BenchHMatrix where + +import Bench.HMatrix +import Criterion.Main + +main :: IO () +main = defaultMain bench_hmatrix + diff --git a/bench/BenchNative.hs b/bench/BenchNative.hs new file mode 100644 --- /dev/null +++ b/bench/BenchNative.hs @@ -0,0 +1,19 @@ +-- | +-- Module : BenchNative +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module BenchNative where + +import Bench.Accelerate +import Criterion.Main +import Data.Array.Accelerate.LLVM.Native as CPU + +main :: IO () +main = defaultMain (bench_accelerate CPU.runN) + diff --git a/bench/BenchPTX.hs b/bench/BenchPTX.hs new file mode 100644 --- /dev/null +++ b/bench/BenchPTX.hs @@ -0,0 +1,19 @@ +-- | +-- Module : BenchPTX +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module BenchPTX where + +import Bench.Accelerate +import Criterion.Main +import Data.Array.Accelerate.LLVM.PTX as PTX + +main :: IO () +main = defaultMain (bench_accelerate PTX.runN) + diff --git a/bench/Extra.hs b/bench/Extra.hs deleted file mode 100644 --- a/bench/Extra.hs +++ /dev/null @@ -1,33 +0,0 @@ -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE PolyKinds #-} -{-# OPTIONS_GHC -fno-warn-orphans #-} - -module Extra - where - -import Data.Complex -import System.Random.MWC - - -data ArgType (a :: *) = AT - -showType :: forall proxy a. Show (ArgType a) => proxy a -> String -showType _ = show (AT :: ArgType a) - -instance Show (ArgType Float) where show _ = "Float" -instance Show (ArgType Double) where show _ = "Double" -instance Show (ArgType (Complex Float)) where show _ = "ComplexFloat" -instance Show (ArgType (Complex Double)) where show _ = "ComplexDouble" - -instance Variate e => Variate (Complex e) where - uniform gen = (:+) <$> uniform gen <*> uniform gen - uniformR r gen = - let (ur:+ui,vr:+vi) = r - in (:+) <$> uniformR (ur,vr) gen <*> uniformR (ui,vi) gen - -infixr 0 $$ -($$) :: (b -> a) -> (c -> d -> b) -> c -> d -> a -(f $$ g) x y = f (g x y) - diff --git a/bench/HMatrix.hs b/bench/HMatrix.hs deleted file mode 100644 --- a/bench/HMatrix.hs +++ /dev/null @@ -1,122 +0,0 @@ -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE PolyKinds #-} -{-# LANGUAGE RankNTypes #-} - -module HMatrix ( - - benchHMatrix - -) where - -import Extra - -import Control.DeepSeq -import Criterion.Main -import Data.Proxy -import Foreign.Storable -import Numeric.LinearAlgebra hiding ( randomVector ) -import System.Random.MWC -import Text.Printf - - -benchHMatrix :: Benchmark -benchHMatrix = - bgroup "hmatrix" - [ level2 - , level3 - ] - - -level2 :: Benchmark -level2 = - bgroup "matrix-vector" - [ bgroup "(#>)" - [ gemv 200 400 - , gemv 500 1000 - , gemv 1000 2000 - , gemv 2000 3000 - ] - , bgroup "(<#)" - [ gevm 200 400 - , gevm 500 1000 - , gevm 1000 2000 - , gevm 2000 3000 - ] - ] - where - gemv :: Int -> Int -> Benchmark - gemv m n = - let setup :: (Variate e, Storable e) => proxy e -> IO (Matrix e, Vector e) - setup _ = withSystemRandom $ \gen -> do - matA <- randomMatrix gen m n - vecx <- randomVector gen n - return (matA, vecx) - - go :: (Variate e, Numeric e, NFData e, Show (ArgType e)) => proxy e -> Benchmark - go t = env (setup t) - $ \ ~(matA, vecx) -> bench (showType t) - $ whnf (matA #>) vecx - in - bgroup (printf "%dx%d" m n) (sdcz go) - - gevm :: Int -> Int -> Benchmark - gevm m n = - let setup :: (Variate e, Storable e) => proxy e -> IO (Matrix e, Vector e) - setup _ = withSystemRandom $ \gen -> do - matA <- randomMatrix gen m n - vecx <- randomVector gen m - return (matA, vecx) - - go :: (Variate e, Numeric e, NFData e, Show (ArgType e)) => proxy e -> Benchmark - go t = env (setup t) - $ \ ~(matA, vecx) -> bench (showType t) - $ whnf (vecx <#) matA - in - bgroup (printf "%dx%d" m n) (sdcz go) - -level3 :: Benchmark -level3 = - bgroup "matrix-matrix" - [ bgroup "(<>)" - [ gemm 100 100 100 - , gemm 250 250 250 - , gemm 500 500 500 - , gemm 1000 1000 1000 - ] - ] - where - gemm :: Int -> Int -> Int -> Benchmark - gemm m n k = - let - setup :: (Variate e, Storable e) => proxy e -> IO (Matrix e, Matrix e) - setup _ = withSystemRandom $ \gen -> do - matA <- randomMatrix gen m k - matB <- randomMatrix gen k n - return (matA, matB) - - go :: (Variate e, Numeric e, NFData e, Show (ArgType e)) => proxy e -> Benchmark - go t = env (setup t) - $ \ ~(matA, matB) -> bench (showType t) - $ whnf (matA <>) matB - in - bgroup (printf "%dx%dx%d" m n k) (sdcz go) - - -randomVector :: (Variate e, Storable e) => GenIO -> Int -> IO (Vector e) -randomVector = uniformVector - -randomMatrix :: (Variate e, Storable e) => GenIO -> Int -> Int -> IO (Matrix e) -randomMatrix gen m n = do - v <- uniformVector gen (m * n) - return $ reshape n v - -sdcz :: (forall (e :: *). (Variate e, Numeric e, NFData e, Show (ArgType e)) => Proxy e -> Benchmark) - -> [Benchmark] -sdcz go = - [ go (Proxy :: Proxy Float) - , go (Proxy :: Proxy Double) - , go (Proxy :: Proxy (Complex Float)) - , go (Proxy :: Proxy (Complex Double)) - ] - diff --git a/bench/Main.hs b/bench/Main.hs deleted file mode 100644 --- a/bench/Main.hs +++ /dev/null @@ -1,24 +0,0 @@ -{-# LANGUAGE CPP #-} - -module Main where - -import HMatrix -import Accelerate - -import Data.Array.Accelerate.Debug ( accInit ) -import Criterion.Main - - -main :: IO () -main = do - accInit - defaultMain - [ benchHMatrix -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND - , benchAcc Native -#endif -#ifdef ACCELERATE_LLVM_PTX_BACKEND - , benchAcc PTX -#endif - ] - diff --git a/cbits/twine_f32.c b/cbits/twine_f32.c deleted file mode 100644 --- a/cbits/twine_f32.c +++ /dev/null @@ -1,60 +0,0 @@ -/* - * Module : Twine - * Copyright : [2016] Trevor L. McDonell - * License : BSD3 - * - * Maintainer : Trevor L. McDonell - * Stability : experimental - * Portability : non-portable (GHC extensions) - * - * Convert between Accelerate's Struct-of-Array representation of complex - * numbers and the Array-of-Struct representation used by BLAS. - */ - -#include -#include "HsFFI.h" - -#ifdef __cplusplus -extern "C" { -#endif - -void interleave_f32 -( - const StgInt start, - const StgInt end, - complex float * __restrict__ cplx, - const float * __restrict__ real, - const float * __restrict__ imag -) -{ - StgInt i; - for (i = start; i < end; ++i) { - const float re = real[i]; - const float im = imag[i]; - - cplx[i] = re + im * I; - } -} - -void deinterleave_f32 -( - const StgInt start, - const StgInt end, - float * __restrict__ real, - float * __restrict__ imag, - const complex float * __restrict__ cplx -) -{ - StgInt i; - for (i = start; i < end; ++i) { - const complex float c = cplx[i]; - - real[i] = crealf(c); - imag[i] = cimagf(c); - } -} - -#ifdef __cplusplus -} -#endif - diff --git a/cbits/twine_f64.c b/cbits/twine_f64.c deleted file mode 100644 --- a/cbits/twine_f64.c +++ /dev/null @@ -1,60 +0,0 @@ -/* - * Module : Twine - * Copyright : [2016] Trevor L. McDonell - * License : BSD3 - * - * Maintainer : Trevor L. McDonell - * Stability : experimental - * Portability : non-portable (GHC extensions) - * - * Convert between Accelerate's Struct-of-Array representation of complex - * numbers and the Array-of-Struct representation used by BLAS. - */ - -#include -#include "HsFFI.h" - -#ifdef __cplusplus -extern "C" { -#endif - -void interleave_f64 -( - const StgInt start, - const StgInt end, - complex double * __restrict__ cplx, - const double * __restrict__ real, - const double * __restrict__ imag -) -{ - StgInt i; - for (i = start; i < end; ++i) { - const double re = real[i]; - const double im = imag[i]; - - cplx[i] = re + im * I; - } -} - -void deinterleave_f64 -( - const StgInt start, - const StgInt end, - double * __restrict__ real, - double * __restrict__ imag, - const complex double * __restrict__ cplx -) -{ - StgInt i; - for (i = start; i < end; ++i) { - const complex double c = cplx[i]; - - real[i] = creal(c); - imag[i] = cimag(c); - } -} - -#ifdef __cplusplus -} -#endif - diff --git a/cubits/twine_f32.ptx b/cubits/twine_f32.ptx deleted file mode 100644 --- a/cubits/twine_f32.ptx +++ /dev/null @@ -1,108 +0,0 @@ -// -// Generated by NVIDIA NVVM Compiler -// -// Compiler Build ID: CL-21140586 -// Cuda compilation tools, release 8.0, V8.0.44 -// Based on LLVM 3.4svn -// - -.version 5.0 -.target sm_20 -.address_size 64 - - // .globl interleave - -.visible .entry interleave( - .param .u64 interleave_param_0, - .param .u64 interleave_param_1, - .param .u64 interleave_param_2, - .param .u32 interleave_param_3 -) -{ - .reg .pred %p<3>; - .reg .f32 %f<3>; - .reg .b32 %r<11>; - .reg .b64 %rd<12>; - - - ld.param.u64 %rd4, [interleave_param_0]; - ld.param.u64 %rd5, [interleave_param_1]; - ld.param.u64 %rd6, [interleave_param_2]; - ld.param.u32 %r5, [interleave_param_3]; - cvta.to.global.u64 %rd1, %rd4; - cvta.to.global.u64 %rd2, %rd6; - cvta.to.global.u64 %rd3, %rd5; - mov.u32 %r6, %nctaid.x; - mov.u32 %r7, %ntid.x; - mul.lo.s32 %r1, %r6, %r7; - mov.u32 %r8, %ctaid.x; - mov.u32 %r9, %tid.x; - mad.lo.s32 %r10, %r8, %r7, %r9; - setp.ge.s32 %p1, %r10, %r5; - @%p1 bra BB0_2; - -BB0_1: - mul.wide.s32 %rd7, %r10, 4; - add.s64 %rd8, %rd3, %rd7; - add.s64 %rd9, %rd2, %rd7; - mul.wide.s32 %rd10, %r10, 8; - add.s64 %rd11, %rd1, %rd10; - ld.global.f32 %f1, [%rd9]; - ld.global.f32 %f2, [%rd8]; - st.global.v2.f32 [%rd11], {%f2, %f1}; - add.s32 %r10, %r10, %r1; - setp.lt.s32 %p2, %r10, %r5; - @%p2 bra BB0_1; - -BB0_2: - ret; -} - - // .globl deinterleave -.visible .entry deinterleave( - .param .u64 deinterleave_param_0, - .param .u64 deinterleave_param_1, - .param .u64 deinterleave_param_2, - .param .u32 deinterleave_param_3 -) -{ - .reg .pred %p<3>; - .reg .f32 %f<5>; - .reg .b32 %r<11>; - .reg .b64 %rd<12>; - - - ld.param.u64 %rd4, [deinterleave_param_0]; - ld.param.u64 %rd5, [deinterleave_param_1]; - ld.param.u64 %rd6, [deinterleave_param_2]; - ld.param.u32 %r5, [deinterleave_param_3]; - cvta.to.global.u64 %rd1, %rd5; - cvta.to.global.u64 %rd2, %rd4; - cvta.to.global.u64 %rd3, %rd6; - mov.u32 %r6, %nctaid.x; - mov.u32 %r7, %ntid.x; - mul.lo.s32 %r1, %r6, %r7; - mov.u32 %r8, %ctaid.x; - mov.u32 %r9, %tid.x; - mad.lo.s32 %r10, %r8, %r7, %r9; - setp.ge.s32 %p1, %r10, %r5; - @%p1 bra BB1_2; - -BB1_1: - mul.wide.s32 %rd7, %r10, 8; - add.s64 %rd8, %rd3, %rd7; - ld.global.v2.f32 {%f1, %f2}, [%rd8]; - mul.wide.s32 %rd9, %r10, 4; - add.s64 %rd10, %rd2, %rd9; - st.global.f32 [%rd10], %f1; - add.s64 %rd11, %rd1, %rd9; - st.global.f32 [%rd11], %f2; - add.s32 %r10, %r10, %r1; - setp.lt.s32 %p2, %r10, %r5; - @%p2 bra BB1_1; - -BB1_2: - ret; -} - - diff --git a/cubits/twine_f64.ptx b/cubits/twine_f64.ptx deleted file mode 100644 --- a/cubits/twine_f64.ptx +++ /dev/null @@ -1,108 +0,0 @@ -// -// Generated by NVIDIA NVVM Compiler -// -// Compiler Build ID: CL-21140586 -// Cuda compilation tools, release 8.0, V8.0.44 -// Based on LLVM 3.4svn -// - -.version 5.0 -.target sm_20 -.address_size 64 - - // .globl interleave - -.visible .entry interleave( - .param .u64 interleave_param_0, - .param .u64 interleave_param_1, - .param .u64 interleave_param_2, - .param .u32 interleave_param_3 -) -{ - .reg .pred %p<3>; - .reg .b32 %r<11>; - .reg .f64 %fd<3>; - .reg .b64 %rd<12>; - - - ld.param.u64 %rd4, [interleave_param_0]; - ld.param.u64 %rd5, [interleave_param_1]; - ld.param.u64 %rd6, [interleave_param_2]; - ld.param.u32 %r5, [interleave_param_3]; - cvta.to.global.u64 %rd1, %rd4; - cvta.to.global.u64 %rd2, %rd6; - cvta.to.global.u64 %rd3, %rd5; - mov.u32 %r6, %nctaid.x; - mov.u32 %r7, %ntid.x; - mul.lo.s32 %r1, %r6, %r7; - mov.u32 %r8, %ctaid.x; - mov.u32 %r9, %tid.x; - mad.lo.s32 %r10, %r8, %r7, %r9; - setp.ge.s32 %p1, %r10, %r5; - @%p1 bra BB0_2; - -BB0_1: - mul.wide.s32 %rd7, %r10, 8; - add.s64 %rd8, %rd3, %rd7; - add.s64 %rd9, %rd2, %rd7; - mul.wide.s32 %rd10, %r10, 16; - add.s64 %rd11, %rd1, %rd10; - ld.global.f64 %fd1, [%rd9]; - ld.global.f64 %fd2, [%rd8]; - st.global.v2.f64 [%rd11], {%fd2, %fd1}; - add.s32 %r10, %r10, %r1; - setp.lt.s32 %p2, %r10, %r5; - @%p2 bra BB0_1; - -BB0_2: - ret; -} - - // .globl deinterleave -.visible .entry deinterleave( - .param .u64 deinterleave_param_0, - .param .u64 deinterleave_param_1, - .param .u64 deinterleave_param_2, - .param .u32 deinterleave_param_3 -) -{ - .reg .pred %p<3>; - .reg .b32 %r<11>; - .reg .f64 %fd<5>; - .reg .b64 %rd<12>; - - - ld.param.u64 %rd4, [deinterleave_param_0]; - ld.param.u64 %rd5, [deinterleave_param_1]; - ld.param.u64 %rd6, [deinterleave_param_2]; - ld.param.u32 %r5, [deinterleave_param_3]; - cvta.to.global.u64 %rd1, %rd5; - cvta.to.global.u64 %rd2, %rd4; - cvta.to.global.u64 %rd3, %rd6; - mov.u32 %r6, %nctaid.x; - mov.u32 %r7, %ntid.x; - mul.lo.s32 %r1, %r6, %r7; - mov.u32 %r8, %ctaid.x; - mov.u32 %r9, %tid.x; - mad.lo.s32 %r10, %r8, %r7, %r9; - setp.ge.s32 %p1, %r10, %r5; - @%p1 bra BB1_2; - -BB1_1: - mul.wide.s32 %rd7, %r10, 16; - add.s64 %rd8, %rd3, %rd7; - ld.global.v2.f64 {%fd1, %fd2}, [%rd8]; - mul.wide.s32 %rd9, %r10, 8; - add.s64 %rd10, %rd2, %rd9; - st.global.f64 [%rd10], %fd1; - add.s64 %rd11, %rd1, %rd9; - st.global.f64 [%rd11], %fd2; - add.s32 %r10, %r10, %r1; - setp.lt.s32 %p2, %r10, %r5; - @%p2 bra BB1_1; - -BB1_2: - ret; -} - - diff --git a/src/Data/Array/Accelerate/Numeric/LinearAlgebra.hs b/src/Data/Array/Accelerate/Numeric/LinearAlgebra.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/LinearAlgebra.hs @@ -0,0 +1,172 @@ +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE ViewPatterns #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra ( + + -- * Types + Numeric, Scalar, Vector, Matrix, + + -- * Products + -- ** Vector-vector + (<.>), + (><), + + -- ** Matrix-vector + (#>), (<#), + + -- ** Matrix-matrix + (<>), + + -- * Diagonal + identity, diagonal, + +) where + +import Data.Array.Accelerate as A + +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type +import Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level1 +import Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level2 +import Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level3 + + +-- Level 1 +-- ------- + +-- | An infix synonym for 'dotu'. +-- +-- >>> let a = fromList (Z:.4) [1..] +-- >>> let b = fromList (Z:.4) [-2,0,1,1] +-- >>> a <.> b +-- Scalar Z [5.0] +-- +-- >>> let c = fromList (Z:.2) [1:+1, 1:+0] +-- >>> let d = fromList (Z:.2) [1:+0, 1:+(-1)] +-- >>> c <.> d +-- Scalar Z [2.0 :+ 0.0] +-- +infixr 8 <.> +(<.>) :: Numeric e => Acc (Vector e) -> Acc (Vector e) -> Acc (Scalar e) +(<.>) = dotu + + +-- | Outer product of two vectors +-- +-- >>> let a = fromList (Z :. 3) [1,2,3] +-- >>> let b = fromList (Z :. 3) [5,2,3] +-- >>> a >< b +-- Matrix (Z :. 3 :. 3) +-- [ 5.0, 2.0, 3.0 +-- , 10.0, 4.0, 6.0 +-- , 15.0, 6.0, 9.0 ] +-- +infixr 8 >< +(><) :: Numeric e => Acc (Vector e) -> Acc (Vector e) -> Acc (Matrix e) +(><) x y = xc <> yr + where + xc = reshape (index2 (length x) 1) x + yr = reshape (index2 1 (length y)) y + + +-- Level 2 +-- ------- + +-- | Dense matrix-vector product +-- +-- >>> let m = fromList (Z :. 2 :. 3) [1..] +-- >>> m +-- Matrix (Z :. 2 :. 3) +-- [ 1.0, 2.0, 3.0 +-- , 4.0, 5.0, 6.0 ] +-- +-- >>> let x = fromList (Z :. 3) [10,20,30] +-- +-- >>> m #> x +-- Vector (Z :. 2) [140.0,320.0] +-- +-- See 'gemv' for a more general version of this operation. +-- +infixr 8 #> +(#>) :: Numeric e => Acc (Matrix e) -> Acc (Vector e) -> Acc (Vector e) +(#>) m x = gemv 1 N m x + + +-- | Dense vector-matrix product +-- +-- >>> let m = fromList (Z :. 2 :. 3) [1..] +-- >>> m +-- Matrix (Z :. 2 :. 3) +-- [1.0,2.0,3.0, +-- 4.0,5.0,6.0] +-- +-- >>> let v = fromList (Z :. 2) [5,10] +-- +-- >>> v <# m +-- Vector (Z :. 3) [45.0,60.0,75.0] +-- +-- See 'gemv' for a more general version of this operation. +-- +infixr 8 <# +(<#) :: Numeric e => Acc (Vector e) -> Acc (Matrix e) -> Acc (Vector e) +(<#) x m = gemv 1 T m x + + +-- Level 3 +-- ------- + +-- | Dense matrix-matrix product +-- +-- >>> let a = fromList (Z :. 3 :. 5) [1..] +-- >>> a +-- Matrix (Z:.3:.5) +-- [ 1.0, 2.0, 3.0, 4.0, 5.0 +-- , 6.0, 7.0, 8.0, 9.0, 10.0 +-- , 11.0, 12.0, 13.0, 14.0, 15.0 ] +-- +-- >>> let b = fromList (Z :. 5 :. 2) [1,3, 0,2, -1,5, 7,7, 6,0] +-- >>> b +-- Matrix (Z :. 5 :. 2) +-- [ 1.0, 3.0 +-- , 0.0, 2.0 +-- , -1.0, 5.0 +-- , 7.0, 7.0 +-- , 6.0, 0.0 ] +-- +-- >>> a <> b +-- Matrix (Z :. 3 :. 2) +-- [ 56.0, 50.0 +-- , 121.0, 135.0 +-- , 186.0, 220.0 ] +-- +-- See 'gemm' for a more general version of this operation. +-- +infixr 8 <> +(<>) :: Numeric e => Acc (Matrix e) -> Acc (Matrix e) -> Acc (Matrix e) +(<>) matA matB = gemm 1 N matA N matB + + +-- | Create a square identity matrix of the given dimension +-- +identity :: Num e => Exp Int -> Acc (Matrix e) +identity n = diagonal (fill (index1 n) 1) + +-- | Create a square matrix with the given diagonal +-- +diagonal :: Num e => Acc (Vector e) -> Acc (Matrix e) +diagonal v = + let n = length v + zeros = fill (index2 n n) 0 + in + permute const zeros (\(unindex1 -> i) -> index2 i i) v + diff --git a/src/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level1.hs b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level1.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level1.hs @@ -0,0 +1,139 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE ScopedTypeVariables #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level1 +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- +-- Level 1 (vector-vector) BLAS operations. +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level1 ( + + -- Types + Numeric, Vector, + + -- Level1 operations + sdot, + dotu, + dotc, + asum, + amax, + amin, + +) where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Data.Complex as A +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type + + +-- | Computes a vector-vector dot product, using double precision accumulation +-- of the intermediate result. Includes a scalar (initial) value to be added to +-- the inner product. +-- +-- +-- +sdot :: forall e. Numeric e => Exp e -> Acc (Vector e) -> Acc (Vector e) -> Acc (Scalar e) +sdot z xs ys = + case numericR :: NumericR e of + NumericRfloat32 -> map toFloating $ dsdot (toFloating z) (map toFloating xs) (map toFloating ys) + NumericRfloat64 -> dsdot z xs ys + NumericRcomplex32 -> map d2f $ zsdot (f2d z) (map f2d xs) (map f2d ys) + NumericRcomplex64 -> zsdot z xs ys + where + dsdot :: Exp Double -> Acc (Vector Double) -> Acc (Vector Double) -> Acc (Scalar Double) + dsdot z' xs' ys' = fold (+) z' (zipWith (*) xs' ys') + + zsdot :: Exp (Complex Double) -> Acc (Vector (Complex Double)) -> Acc (Vector (Complex Double)) -> Acc (Scalar (Complex Double)) + zsdot z' xs' ys' = fold (+) z' (zipWith (*) xs' ys') + + f2d :: Exp (Complex Float) -> Exp (Complex Double) + f2d c = lift (toFloating (real c) :+ toFloating (imag c)) + + d2f :: Exp (Complex Double) -> Exp (Complex Float) + d2f c = lift (toFloating (real c) :+ toFloating (imag c)) + + +-- | Computes a vector-vector dot product +-- +-- \[ +-- res = \sum_i x_i * y_i +-- \] +-- +-- +-- +dotu :: Numeric e => Acc (Vector e) -> Acc (Vector e) -> Acc (Scalar e) +dotu xs ys = fold (+) 0 (zipWith (*) xs ys) + + +-- | Computes a dot product of a conjugated vector with another vector +-- +-- \[ +-- res = \sum_i \mathrm{conj}(x_i) * y_i +-- \] +-- +-- +-- +dotc :: forall e. Numeric (Complex e) + => Acc (Vector (Complex e)) + -> Acc (Vector (Complex e)) + -> Acc (Scalar (Complex e)) +dotc xs ys = + case numericR :: NumericR (Complex e) of + NumericRcomplex32 -> dotu (map conjugate xs) ys + NumericRcomplex64 -> dotu (map conjugate xs) ys + + +-- | Computes the sum of magnitudes of the vector elements. For complex values, +-- this is given by \(\sum_i \|\mathrm{real}(x_i)\| + \|\mathrm{imag}(x_i)\|\). +-- +-- +-- +asum :: forall e. Numeric e => Acc (Vector e) -> Acc (Scalar (NumericBaseT e)) +asum = + case numericR :: NumericR e of + NumericRfloat32 -> sum . map abs + NumericRfloat64 -> sum . map abs + NumericRcomplex32 -> sum . map mag + NumericRcomplex64 -> sum . map mag + where + mag c = abs (real c) + abs (imag c) + + +-- | Return the index of the element with the maximum absolute value. +-- +-- +-- +amax :: forall e. Numeric e => Acc (Vector e) -> Acc (Scalar Int) +amax = + case numericR :: NumericR e of + NumericRfloat32 -> map (indexHead . fst) . fold1 cmp . indexed . map abs + NumericRfloat64 -> map (indexHead . fst) . fold1 cmp . indexed . map abs + NumericRcomplex32 -> map (indexHead . fst) . fold1 cmp . indexed . map mag + NumericRcomplex64 -> map (indexHead . fst) . fold1 cmp . indexed . map mag + where + cmp ix iy = snd ix > snd iy ? ( ix, iy ) + mag c = abs (real c) + abs (imag c) + +-- | Return the index of the element with the minimum absolute value. +-- +-- +-- +amin :: forall e. Numeric e => Acc (Vector e) -> Acc (Scalar Int) +amin = + case numericR :: NumericR e of + NumericRfloat32 -> map (indexHead . fst) . fold1 cmp . indexed . map abs + NumericRfloat64 -> map (indexHead . fst) . fold1 cmp . indexed . map abs + NumericRcomplex32 -> map (indexHead . fst) . fold1 cmp . indexed . map mag + NumericRcomplex64 -> map (indexHead . fst) . fold1 cmp . indexed . map mag + where + cmp ix iy = snd ix < snd iy ? ( ix, iy ) + mag c = abs (real c) + abs (imag c) + diff --git a/src/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level2.hs b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level2.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level2.hs @@ -0,0 +1,92 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ViewPatterns #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level2 +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- +-- Level 2 (matrix-vector) BLAS operations. +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level2 ( + + -- Types + Numeric, Vector, Matrix, Transpose(..), + + -- Operations + gemv, + +) where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Smart as A +import Data.Array.Accelerate.Data.Complex as A +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type + +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND +import qualified Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level2 as CPU +#endif +#ifdef ACCELERATE_LLVM_PTX_BACKEND +import qualified Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level2 as PTX +#endif + + +-- | Computes the matrix-vector product of a general matrix. +-- +-- \[ +-- y = \alpha * \mathrm{op}(A) * x +-- \] +-- +-- where: +-- +-- * 'shape' \(\mathrm{op}(A)\) @= Z :. m :. n@ +-- * 'shape' \(x\) @= Z :. n@ +-- * 'shape' \(y\) @= Z :. m@ +-- +-- +-- +gemv :: forall e. Numeric e + => Exp e -- ^ \( \alpha \) + -> Transpose -- ^ Operation to apply to A + -> Acc (Matrix e) -- ^ A + -> Acc (Vector e) -- ^ x + -> Acc (Vector e) -- ^ y +gemv alpha opA matA x = go (lift (unit alpha, matA, x)) + where + go = +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND + foreignAcc (CPU.gemv opA) $ +#endif +#ifdef ACCELERATE_LLVM_PTX_BACKEND + foreignAcc (PTX.gemv opA) $ +#endif + (\(unatup3 -> (_, arr, brr)) -> mXv arr brr) + + -- General matrix-vector multiply in pure Accelerate. This is probably not + -- efficient. + -- + mXv :: Acc (Matrix e) -> Acc (Vector e) -> Acc (Vector e) + mXv arr brr + = fold (+) 0 + $ zipWith (\a b -> alpha * a * b) arr' brr' + where + Z :. m :. _ = unlift (shape arr') :: Z :. Exp Int :. Exp Int + + brr' = replicate (lift (Z :. m :. All)) brr + arr' = case opA of + N -> arr + T -> transpose arr + H -> case numericR :: NumericR e of + NumericRcomplex32 -> map conjugate (transpose arr) + NumericRcomplex64 -> map conjugate (transpose arr) + _ -> transpose arr + diff --git a/src/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level3.hs b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level3.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level3.hs @@ -0,0 +1,109 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ViewPatterns #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level3 +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- +-- Level 3 (matrix-matrix) BLAS operations. +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level3 ( + + -- Types + Numeric, Matrix, Transpose(..), + + -- Matrix-matrix operations + gemm, + +) where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Smart as A +import Data.Array.Accelerate.Data.Complex as A +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type + +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND +import qualified Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level3 as CPU +#endif +#ifdef ACCELERATE_LLVM_PTX_BACKEND +import qualified Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level3 as PTX +#endif + + +-- | General matrix-matrix multiply +-- +-- \[ +-- C = \alpha * \mathrm{op}(A) * \mathrm{op}(B) +-- \] +-- +-- where: +-- +-- * 'shape' \(\mathrm{op}(A)\) @= Z :. m :. k@ +-- * 'shape' \(\mathrm{op}(B)\) @= Z :. k :. n@ +-- * 'shape' \(C\) @= Z :. m :. n@ +-- +-- +-- +gemm :: forall e. Numeric e + => Exp e -- ^ \( \alpha \) + -> Transpose -- ^ operation to apply to A + -> Acc (Matrix e) -- ^ A + -> Transpose -- ^ operation to apply to B + -> Acc (Matrix e) -- ^ B + -> Acc (Matrix e) -- ^ C +gemm alpha opA matA opB matB = go (lift (unit alpha, matA, matB)) + where + go = +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND + foreignAcc (CPU.gemm opA opB) $ +#endif +#ifdef ACCELERATE_LLVM_PTX_BACKEND + foreignAcc (PTX.gemm opA opB) $ +#endif + (\(unatup3 -> (_, arr, brr)) -> mXm arr brr) + + -- General dense matrix-matrix multiply written in pure Accelerate. This is + -- not efficient due to the memory access patterns. We could probably + -- improve this a little bit with a divide-and-conquer algorithm, for + -- example, but using a foreign implementation will be best. + -- + mXm :: Acc (Matrix e) -> Acc (Matrix e) -> Acc (Matrix e) + mXm arr brr + = fold (+) 0 + $ zipWith (\a b -> alpha * a * b) arrRepl brrRepl + where + Z :. rowsA :. _ = unlift (shape arr') :: Z :. Exp Int :. Exp Int + Z :. colsB :. _ = unlift (shape brr') :: Z :. Exp Int :. Exp Int + -- + arrRepl = replicate (lift $ Z :. All :. colsB :. All) arr' + brrRepl = replicate (lift $ Z :. rowsA :. All :. All) brr' + + -- apply opA + arr' = case opA of + N -> arr + T -> transpose arr + H -> case numericR :: NumericR e of + NumericRcomplex32 -> map conjugate (transpose arr) + NumericRcomplex64 -> map conjugate (transpose arr) + _ -> transpose arr + + -- apply opB and transpose at the same time, which is required for this + -- algorithm + brr' = case opB of + N -> transpose brr + T -> brr + H -> case numericR :: NumericR e of + NumericRcomplex32 -> map conjugate brr + NumericRcomplex64 -> map conjugate brr + _ -> brr + diff --git a/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Base.hs b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Base.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Base.hs @@ -0,0 +1,48 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Base +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Base + where + +import Data.Array.Accelerate.Array.Sugar ( Array(..), EltRepr ) +import Data.Array.Accelerate.Array.Data +import Data.Array.Accelerate.Array.Unique +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type + +import qualified Blas.Primitive.Types as C + + +encodeTranspose :: Transpose -> C.Transpose +encodeTranspose N = C.NoTrans +encodeTranspose T = C.Trans +encodeTranspose H = C.ConjTrans + + +withArray + :: forall sh e b. Numeric e + => Array sh e + -> (ArrayPtrs (EltRepr e) -> IO b) + -> IO b +withArray (Array _ adata) = withArrayData (numericR::NumericR e) adata + +withArrayData + :: NumericR e + -> ArrayData (EltRepr e) + -> (ArrayPtrs (EltRepr e) -> IO b) + -> IO b +withArrayData NumericRfloat32 (AD_Float ua) f = withUniqueArrayPtr ua f +withArrayData NumericRfloat64 (AD_Double ua) f = withUniqueArrayPtr ua f +withArrayData NumericRcomplex32 (AD_V2 (AD_Float ua)) f = withUniqueArrayPtr ua f +withArrayData NumericRcomplex64 (AD_V2 (AD_Double ua)) f = withUniqueArrayPtr ua f + diff --git a/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level2.hs b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level2.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level2.hs @@ -0,0 +1,54 @@ +{-# LANGUAGE GADTs #-} +{-# LANGUAGE ScopedTypeVariables #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level2 +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level2 + where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.LLVM.Native.Foreign +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Base + +import Foreign.Ptr +import qualified Blas.Primitive.Types as C +import qualified Blas.Primitive.Unsafe as C + + +gemv :: forall e. Numeric e + => Transpose + -> ForeignAcc ((Scalar e, Matrix e, Vector e) -> Vector e) +gemv opA = ForeignAcc "native.gemv" gemv' + where + gemv' (alpha, matA, vecx) = do + let + Z :. rowsA :. colsA = arrayShape matA + + sizeY = case opA of + N -> rowsA + _ -> colsA + + opA' = encodeTranspose opA + alpha' = indexArray alpha Z + -- + vecy <- allocateRemote (Z :. sizeY) :: LLVM Native (Vector e) + () <- liftIO $ do + withArray matA $ \ptr_A -> do + withArray vecx $ \ptr_x -> do + withArray vecy $ \ptr_y -> do + case numericR :: NumericR e of + NumericRfloat32 -> C.sgemv C.RowMajor opA' rowsA colsA alpha' ptr_A colsA ptr_x 1 0 ptr_y 1 + NumericRfloat64 -> C.dgemv C.RowMajor opA' rowsA colsA alpha' ptr_A colsA ptr_x 1 0 ptr_y 1 + NumericRcomplex32 -> C.cgemv C.RowMajor opA' rowsA colsA alpha' (castPtr ptr_A) colsA (castPtr ptr_x) 1 0 (castPtr ptr_y) 1 + NumericRcomplex64 -> C.zgemv C.RowMajor opA' rowsA colsA alpha' (castPtr ptr_A) colsA (castPtr ptr_x) 1 0 (castPtr ptr_y) 1 + -- + return vecy + diff --git a/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level3.hs b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level3.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level3.hs @@ -0,0 +1,67 @@ +{-# LANGUAGE GADTs #-} +{-# LANGUAGE ScopedTypeVariables #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level3 +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level3 + where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.LLVM.Native.Foreign +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Base + +import Foreign.Ptr +import qualified Blas.Primitive.Types as C +import qualified Blas.Primitive.Unsafe as C + + +-- TODO: check whether it is faster to compute this as column-major order: +-- +-- https://www.christophlassner.de/using-blas-from-c-with-row-major-data.html +-- +gemm :: forall e. Numeric e + => Transpose + -> Transpose + -> ForeignAcc ((Scalar e, Matrix e, Matrix e) -> Matrix e) +gemm opA opB = ForeignAcc "native.gemm" gemm' + where + gemm' (alpha, matA, matB) = do + let + Z :. rowsA :. colsA = arrayShape matA + Z :. rowsB :. colsB = arrayShape matB + + (m,k) = case opA of + N -> (rowsA, colsA) + _ -> (colsA, rowsA) + n = case opB of + N -> colsB + _ -> rowsB + + lda = colsA + ldb = colsB + + opA' = encodeTranspose opA + opB' = encodeTranspose opB + alpha' = indexArray alpha Z + -- + matC <- allocateRemote (Z :. m :. n) :: LLVM Native (Matrix e) + () <- liftIO $ do + withArray matA $ \ptr_A -> do + withArray matB $ \ptr_B -> do + withArray matC $ \ptr_C -> do + case numericR :: NumericR e of + NumericRfloat32 -> C.sgemm C.RowMajor opA' opB' m n k alpha' ptr_A lda ptr_B ldb 0 ptr_C n + NumericRfloat64 -> C.dgemm C.RowMajor opA' opB' m n k alpha' ptr_A lda ptr_B ldb 0 ptr_C n + NumericRcomplex32 -> C.cgemm C.RowMajor opA' opB' m n k alpha' (castPtr ptr_A) lda (castPtr ptr_B) ldb 0 (castPtr ptr_C) n + NumericRcomplex64 -> C.zgemm C.RowMajor opA' opB' m n k alpha' (castPtr ptr_A) lda (castPtr ptr_B) ldb 0 (castPtr ptr_C) n + -- + return matC + diff --git a/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Base.hs b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Base.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Base.hs @@ -0,0 +1,82 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base + where + +import Data.Array.Accelerate.Array.Data +import Data.Array.Accelerate.Array.Sugar ( Array(..), EltRepr ) +import Data.Array.Accelerate.Lifetime +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type +import Data.Array.Accelerate.Type + +import Data.Array.Accelerate.LLVM.PTX.Foreign + +import Foreign.CUDA.Ptr ( DevicePtr ) +import qualified Foreign.CUDA.BLAS as C + + +type family DevicePtrs e :: * +type instance DevicePtrs Float = DevicePtr Float +type instance DevicePtrs Double = DevicePtr Double +type instance DevicePtrs (V2 Float) = DevicePtr Float +type instance DevicePtrs (V2 Double) = DevicePtr Double + + +encodeTranspose :: Transpose -> C.Operation +encodeTranspose N = C.N +encodeTranspose T = C.T +encodeTranspose H = C.C + + +withArray + :: forall sh e b. Numeric e + => Array sh e + -> Stream + -> (DevicePtrs (EltRepr e) -> LLVM PTX b) + -> LLVM PTX b +withArray (Array _ adata) s k = withArrayData (numericR::NumericR e) adata s k + +withArrayData + :: NumericR e + -> ArrayData (EltRepr e) + -> Stream + -> (DevicePtrs (EltRepr e) -> LLVM PTX b) + -> LLVM PTX b +withArrayData NumericRfloat32 ad s k = + withDevicePtr ad $ \p -> do + r <- k p + e <- checkpoint s + return (Just e,r) +withArrayData NumericRfloat64 ad s k = + withDevicePtr ad $ \p -> do + r <- k p + e <- checkpoint s + return (Just e, r) +withArrayData NumericRcomplex32 (AD_V2 ad) s k = + withDevicePtr ad $ \p -> do + r <- k p + e <- checkpoint s + return (Just e,r) +withArrayData NumericRcomplex64 (AD_V2 ad) s k = + withDevicePtr ad $ \p -> do + r <- k p + e <- checkpoint s + return (Just e, r) + + +withLifetime' :: Lifetime a -> (a -> LLVM PTX b) -> LLVM PTX b +withLifetime' l k = do + r <- k (unsafeGetValue l) + liftIO $ touchLifetime l + return r + diff --git a/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Context.hs b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Context.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Context.hs @@ -0,0 +1,76 @@ +{-# LANGUAGE MagicHash #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Context +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Context ( + + withBLAS + +) where + +import Data.Array.Accelerate.Lifetime +import Data.Array.Accelerate.LLVM.PTX +import Data.Array.Accelerate.LLVM.PTX.Foreign +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base + +import Control.Monad.State +import Control.Concurrent.MVar +import Data.IntMap.Strict ( IntMap ) +import System.IO.Unsafe +import qualified Data.IntMap.Strict as IM + +import qualified Foreign.CUDA.Driver.Context as CUDA +import qualified Foreign.CUDA.BLAS as BLAS + +import GHC.Ptr +import GHC.Base +import Prelude hiding ( lookup ) + + +-- Execute an operation with a cuBLAS handle appropriate for the current +-- execution context. +-- +-- Initial creation of the context is an atomic operation, but subsequently +-- multiple threads may use the context concurrently. +-- +-- +-- +withBLAS :: (BLAS.Handle -> LLVM PTX b) -> LLVM PTX b +withBLAS k = do + lc <- gets (deviceContext . ptxContext) + h <- liftIO $ + withLifetime lc $ \ctx -> + modifyMVar handles $ \im -> + let key = toKey ctx in + case IM.lookup key im of + -- handle does not exist yet; create it and add to the global + -- state for reuse + Nothing -> do + h <- BLAS.create + l <- newLifetime h + -- BLAS.setPointerMode h BLAS.Device + BLAS.setAtomicsMode h BLAS.Allowed + addFinalizer lc $ modifyMVar handles (\im' -> return (IM.delete key im', ())) + addFinalizer l $ BLAS.destroy h + return ( IM.insert key l im, l ) + + -- return existing handle + Just h -> return (im, h) + -- + withLifetime' h k + + +toKey :: CUDA.Context -> IM.Key +toKey (CUDA.Context (Ptr addr#)) = I# (addr2Int# addr#) + +{-# NOINLINE handles #-} +handles :: MVar (IntMap (Lifetime BLAS.Handle)) +handles = unsafePerformIO $ newMVar IM.empty + diff --git a/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level2.hs b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level2.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level2.hs @@ -0,0 +1,113 @@ +{-# LANGUAGE GADTs #-} +{-# LANGUAGE ScopedTypeVariables #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level2 +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level2 + where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Array.Sugar ( Array(..) ) +import Data.Array.Accelerate.LLVM.PTX.Foreign +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Context +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level3 +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type + +import Foreign.Marshal ( with ) +import qualified Foreign.CUDA.Ptr as CUDA +import qualified Foreign.CUDA.BLAS as BLAS + + +-- NOTE: cuBLAS requires matrices to be stored in column-major order +-- (Fortran-style), but Accelerate uses C-style arrays in row-major order. +-- +-- If the operation is N or T, we can just swap the operation. For +-- conjugate-transpose (H) operations (on complex valued arguments), since there +-- is no conjugate-no-transpose operation, we implement that via 'gemm', which +-- I assume is more efficient than ?geam followed by ?gemv. +-- +gemv :: Numeric e + => Transpose + -> ForeignAcc ((Scalar e, Matrix e, Vector e) -> Vector e) +gemv opA = ForeignAcc "ptx.gemv" (gemv' numericR opA) + +gemv' :: Numeric e + => NumericR e + -> Transpose + -> Stream + -> (Scalar e, Matrix e, Vector e) + -> LLVM PTX (Vector e) +gemv' NumericRcomplex32 H = as_gemm H +gemv' NumericRcomplex64 H = as_gemm H +gemv' _ t = as_gemv t + + +as_gemm + :: Numeric e + => Transpose + -> Stream + -> (Scalar e, Matrix e, Vector e) + -> LLVM PTX (Vector e) +as_gemm opA stream (alpha, matA, Array sh adata) = do + let matB = Array (sh,1) adata + -- + Array (sh',1) vecy <- gemm' opA N stream (alpha, matA, matB) + return (Array sh' vecy) + +as_gemv + :: forall e. Numeric e + => Transpose + -> Stream + -> (Scalar e, Matrix e, Vector e) + -> LLVM PTX (Vector e) +as_gemv opA stream (alpha, matA, vecx) = do + let + Z :. rowsA :. colsA = arrayShape matA + + sizeY = case opA of + N -> rowsA + _ -> colsA + + opA' = encodeTranspose + $ case opA of + N -> T + _ -> N + -- + vecy <- allocateRemote (Z :. sizeY) :: LLVM PTX (Vector e) + alpha' <- indexRemote alpha 0 + () <- do + withArray matA stream $ \ptr_A -> do + withArray vecx stream $ \ptr_x -> do + withArray vecy stream $ \ptr_y -> do + withBLAS $ \hdl -> do + case numericR :: NumericR e of + NumericRfloat32 -> liftIO $ + with alpha' $ \ptr_alpha -> + with 0 $ \ptr_beta -> + BLAS.sgemv hdl opA' colsA rowsA ptr_alpha ptr_A colsA ptr_x 1 ptr_beta ptr_y 1 + + NumericRfloat64 -> liftIO $ + with alpha' $ \ptr_alpha -> + with 0 $ \ptr_beta -> + BLAS.dgemv hdl opA' colsA rowsA ptr_alpha ptr_A colsA ptr_x 1 ptr_beta ptr_y 1 + + NumericRcomplex32 -> liftIO $ + with alpha' $ \ptr_alpha -> + with 0 $ \ptr_beta -> + BLAS.cgemv hdl opA' colsA rowsA ptr_alpha (CUDA.castDevPtr ptr_A) colsA (CUDA.castDevPtr ptr_x) 1 ptr_beta (CUDA.castDevPtr ptr_y) 1 + + NumericRcomplex64 -> liftIO $ + with alpha' $ \ptr_alpha -> + with 0 $ \ptr_beta -> + BLAS.zgemv hdl opA' colsA rowsA ptr_alpha (CUDA.castDevPtr ptr_A) colsA (CUDA.castDevPtr ptr_x) 1 ptr_beta (CUDA.castDevPtr ptr_y) 1 + -- + return vecy + diff --git a/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level3.hs b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level3.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level3.hs @@ -0,0 +1,93 @@ +{-# LANGUAGE GADTs #-} +{-# LANGUAGE ScopedTypeVariables #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level3 +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level3 + where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.LLVM.PTX.Foreign +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Context +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type + +import Foreign.Marshal ( with ) +import qualified Foreign.CUDA.Ptr as CUDA +import qualified Foreign.CUDA.BLAS as BLAS + + +-- NOTE: cuBLAS requires that matrices are stored in column-major order +-- (Fortran-style), but Accelerate uses a C-style convention where matrices are +-- stored in row-major order. +-- +-- At least for matrix-matrix multiply, we can get around this problem by making +-- use of the equivalence \( B^T \cdot A^T = (A \cdot B)^T \). +-- +gemm :: Numeric e + => Transpose + -> Transpose + -> ForeignAcc ((Scalar e, Matrix e, Matrix e) -> Matrix e) +gemm opA opB = ForeignAcc "ptx.gemm" (gemm' opA opB) + +gemm' + :: forall e. Numeric e + => Transpose + -> Transpose + -> Stream + -> (Scalar e, Matrix e, Matrix e) + -> LLVM PTX (Matrix e) +gemm' opA opB stream (alpha, matA, matB) = do + let + Z :. rowsA :. colsA = arrayShape matA + Z :. rowsB :. colsB = arrayShape matB + + (m,k) = case opA of + N -> (rowsA, colsA) + _ -> (colsA, rowsA) + n = case opB of + N -> colsB + _ -> rowsB + + lda = colsA + ldb = colsB + + opA' = encodeTranspose opA + opB' = encodeTranspose opB + -- + matC <- allocateRemote (Z :. m :. n) :: LLVM PTX (Matrix e) + alpha' <- indexRemote alpha 0 + () <- withArray matA stream $ \ptr_A -> do + withArray matB stream $ \ptr_B -> do + withArray matC stream $ \ptr_C -> do + withBLAS $ \hdl -> do + case numericR :: NumericR e of + NumericRfloat32 -> liftIO $ + with alpha' $ \ptr_alpha -> + with 0 $ \ptr_beta -> + BLAS.sgemm hdl opB' opA' n m k ptr_alpha ptr_B ldb ptr_A lda ptr_beta ptr_C n + + NumericRfloat64 -> liftIO $ + with alpha' $ \ptr_alpha -> + with 0 $ \ptr_beta -> + BLAS.dgemm hdl opB' opA' n m k ptr_alpha ptr_B ldb ptr_A lda ptr_beta ptr_C n + + NumericRcomplex32 -> liftIO $ + with alpha' $ \ptr_alpha -> + with 0 $ \ptr_beta -> + BLAS.cgemm hdl opB' opA' n m k ptr_alpha (CUDA.castDevPtr ptr_B) ldb (CUDA.castDevPtr ptr_A) lda ptr_beta (CUDA.castDevPtr ptr_C) n + + NumericRcomplex64 -> liftIO $ + with alpha' $ \ptr_alpha -> + with 0 $ \ptr_beta -> + BLAS.zgemm hdl opB' opA' n m k ptr_alpha (CUDA.castDevPtr ptr_B) ldb (CUDA.castDevPtr ptr_A) lda ptr_beta (CUDA.castDevPtr ptr_C) n + + return matC + diff --git a/src/Data/Array/Accelerate/Numeric/LinearAlgebra/Type.hs b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/Type.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/LinearAlgebra/Type.hs @@ -0,0 +1,100 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE TypeFamilies #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.Type +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.Type ( + +#if MIN_VERSION_accelerate(1,2,0) + Matrix, +#endif + module Data.Array.Accelerate.Numeric.LinearAlgebra.Type, + +) where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Data.Complex as A + +import qualified Prelude as P + + +-- For explicit dictionary reification, to recover the type the operation should +-- be performed at. +-- +data NumericR a where + NumericRfloat32 :: NumericR Float + NumericRfloat64 :: NumericR Double + NumericRcomplex32 :: NumericR (Complex Float) + NumericRcomplex64 :: NumericR (Complex Double) + +class (Elt a, Num a) => Numeric a where + numericR :: NumericR a + +instance Numeric Float where + numericR = NumericRfloat32 + +instance Numeric Double where + numericR = NumericRfloat64 + +instance Numeric (Complex Float) where + numericR = NumericRcomplex32 + +instance Numeric (Complex Double) where + numericR = NumericRcomplex64 + +-- class Numeric a => RealNumeric a +-- +-- instance RealNumeric Float +-- instance RealNumeric Double + +type family NumericBaseT t where + NumericBaseT Float = Float + NumericBaseT Double = Double + NumericBaseT (Complex Float) = Float + NumericBaseT (Complex Double) = Double + + +#if !MIN_VERSION_accelerate(1,2,0) +-- | Matrices as dense two-dimensional arrays in row-major ordering +-- +type Matrix e = Array DIM2 e +#endif + +-- | Orientation of the underlying data. +-- +-- Accelerate arrays are naturally stored in row-major format. +-- +data Orientation + = R -- ^ row major + | C -- ^ column major + deriving (P.Eq, P.Show) + +-- | Many operations allow you to implicitly transpose the arguments. For +-- a given input matrix @mat@ with dimensions @Z :. m :. n@ (that is; @m@ rows +-- and @n@ columns): +-- +data Transpose + -- | Leave the matrix as is. + = N + + -- | Treat the matrix as implicitly transposed, with dimensions @Z :. n :. m@. + -- Entry @Z :. j :. i@ is treated as actually being entry @Z :. i :. j@. + | T + + -- | Implicitly transpose and conjugate the input matrix. For complex-valued + -- matrices a given element @mat ! Z:.j:.i == x :+ y@ will be treated as + -- actually being @mat ! Z:.i:.j == x :+ (-y)@. + | H + deriving (P.Eq, P.Show) + diff --git a/src/Data/Array/Accelerate/Numeric/Sum.hs b/src/Data/Array/Accelerate/Numeric/Sum.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/Sum.hs @@ -0,0 +1,330 @@ +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE RebindableSyntax #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ViewPatterns #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.Sum +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- +-- Functions for summing floating point numbers more accurately than the +-- straightforward 'Data.Array.Accelerate.sum' operation. +-- +-- In the worst case, the 'Data.Array.Accelerate.sum' function accumulates error +-- at a rate proportional to the number of values being summed. The algorithms +-- in this module implement different methods of /compensated summation/, which +-- reduce the accumulation of numeric error so that it grows much more slowly +-- than the number of inputs (e.g. logarithmically), or remains constant. +-- + +-- TLM: The standard formulation of the algorithms implemented here are not +-- associative; e.g. they would have a type (KBN a -> a -> KBN a). I've +-- done what seems like the sensible conversion, but somebody versed in numeric +-- analysis should probably look... +-- +-- See also: +-- + +module Data.Array.Accelerate.Numeric.Sum ( + + -- * Summation type class + Summation(..), + sum, + + -- * Kahan-Babuška-Neumaier summation + KBN(..), + kbn, + + -- * Order-2 Kahan-Babuška summation + KB2(..), + kb2, + + -- * Kahan summation + Kahan(..), + kahan, + +) where + +import Data.Array.Accelerate as A hiding ( sum, fromInteger ) +import Data.Array.Accelerate.Type as A +import Data.Array.Accelerate.Smart as A ( Exp(..), PreExp(..) ) +import Data.Array.Accelerate.Product as A +import Data.Array.Accelerate.Array.Sugar as A +import Data.Array.Accelerate.Numeric.Sum.Arithmetic as A + +import Data.Proxy +import Data.Typeable +import Prelude ( Show, fromInteger ) + + +-- | Sum an array using a particular compensation scheme. +-- +-- >>> let xs = [1.0, 1.0e100, 1.0, -1.0e100] :: [Double] +-- >>> Prelude.sum xs +-- 0.0 +-- +-- >>> let ys = fromList (Z:.4) [1.0, 1.0e100, 1.0, -1.0e100] :: Vector Double +-- >>> sum kbn (use ys) +-- Scalar Z [2.0] +-- +sum :: (Summation s a, Shape sh) => Proxy s -> Acc (Array (sh:.Int) a) -> Acc (Array sh a) +sum p = A.map (from p) + . A.fold add zero + . A.map (into p) + + +-- | A class for the summation of floating-point numbers +-- +class (Elt a, Elt (s a)) => Summation s a where + -- | Add a value to the sum + add :: Exp (s a) -> Exp (s a) -> Exp (s a) + + -- | The identity of the summation + zero :: Exp (s a) + + -- | Insert a value into the summation + into :: Proxy s -> Exp a -> Exp (s a) + + -- | Summarise the result of summation + from :: Proxy s -> Exp (s a) -> Exp a + + +-- | Kahan-Babuška-Neumaier summation. This is a little more computationally +-- costly than plain Kahan summation, but is /always/ at least as accurate. +-- +data KBN a = KBN a a + deriving (Show, Typeable) + +-- | Return the result of a Kahan-Babuška-Neumaier sum. +-- +kbn :: Proxy KBN +kbn = Proxy + +kbnAdd :: (Num a, Ord a, IsFloating a) => Exp (KBN a) -> Exp (KBN a) -> Exp (KBN a) +kbnAdd (unlift -> KBN s1 c1) (unlift -> KBN s2 c2) = lift (KBN s' c') + where + s' = s1 `fadd` s2 + c' = c1 `fadd` c2 `fadd` if abs s1 >= abs s2 + then (s1 `fsub` s') `fadd` s2 + else (s2 `fsub` s') `fadd` s1 + +-- instance (Num a, Ord a) => Summation KBN a where +-- zero = lift $ KBN (0::Exp a) (0::Exp a) +-- add = kbnAdd +-- into _ x = lift (KBN x 0) +-- from _ x = let KBN s c = unlift x in s + c + +instance Summation KBN Float where + zero = constant (KBN 0 0) + add = kbnAdd + into _ x = lift (KBN x 0) + from _ x = let KBN s c = unlift x in s + c + +instance Summation KBN Double where + zero = constant (KBN 0 0) + add = kbnAdd + into _ x = lift (KBN x 0) + from _ x = let KBN s c = unlift x in s + c + +instance Summation KBN CFloat where + zero = constant (KBN 0 0) + add = kbnAdd + into _ x = lift (KBN x 0) + from _ x = let KBN s c = unlift x in s + c + +instance Summation KBN CDouble where + zero = constant (KBN 0 0) + add = kbnAdd + into _ x = lift (KBN x 0) + from _ x = let KBN s c = unlift x in s + c + +type instance EltRepr (KBN a) = (((), EltRepr a), EltRepr a) + +instance Elt a => Elt (KBN a) where + eltType _ = TypeRunit `TypeRpair` eltType (undefined::a) + `TypeRpair` eltType (undefined::a) + toElt (((),a),b) = KBN (toElt a) (toElt b) + fromElt (KBN a b) = (((), fromElt a), fromElt b) + +instance Elt a => IsProduct Elt (KBN a) where + type ProdRepr (KBN a) = (((), a), a) + toProd _ (((),a),b) = KBN a b + fromProd _ (KBN a b) = (((),a),b) + prod _ _ = ProdRsnoc $ ProdRsnoc ProdRunit + +instance (Lift Exp a, Elt (Plain a)) => Lift Exp (KBN a) where + type Plain (KBN a) = KBN (Plain a) + lift (KBN a b) = Exp $ Tuple $ NilTup `SnocTup` lift a + `SnocTup` lift b + +instance Elt a => Unlift Exp (KBN (Exp a)) where + unlift t = KBN (Exp $ SuccTupIdx ZeroTupIdx `Prj` t) + (Exp $ ZeroTupIdx `Prj` t) + + +-- | Second-order Kahan-Babuška summation. This is more computationally costly +-- than Kahan-Babuška-Neumaier summation. Its advantage is that it can lose less +-- precision (in admittedly obscure cases). +-- +-- This method compensates for error in both the sum and the first-order +-- compensation term, hence the use of \"second order\" in the name. +-- +data KB2 a = KB2 a a a + deriving (Show, Typeable) + +-- | Return the result of a second-order Kahan-Babuška sum. +-- +kb2 :: Proxy KB2 +kb2 = Proxy + +kb2Add :: (Num a, Ord a, IsFloating a) => Exp (KB2 a) -> Exp (KB2 a) -> Exp (KB2 a) +kb2Add (unlift -> KB2 s1 c1 cc1) (unlift -> KB2 s2 c2 cc2) = lift (KB2 sum' c' cc') + where + sum' = s1 `fadd` s2 + c' = t `fadd` k + cc' = cc1 `fadd` cc2 `fadd` if abs t >= abs k + then (t `fsub` c') `fadd` k + else (k `fsub` c') `fadd` t + t = c1 `fadd` c2 + k = if abs s1 >= abs s2 + then (s1 `fsub` sum') `fadd` s2 + else (s2 `fsub` sum') `fadd` s1 + +-- instance (Num a, Ord a) => Summation KB2 a where +-- zero = lift $ KB2 (0::Exp a) (0::Exp a) (0::Exp a) +-- add = kb2Add +-- into _ x = lift (KB2 x 0 0) +-- from _ x = let KB2 s c cc = unlift x in s + c + cc + +instance Summation KB2 Float where + zero = constant (KB2 0 0 0) + add = kb2Add + into _ x = lift (KB2 x 0 0) + from _ x = let KB2 s c cc = unlift x in s + c + cc + +instance Summation KB2 Double where + zero = constant (KB2 0 0 0) + add = kb2Add + into _ x = lift (KB2 x 0 0) + from _ x = let KB2 s c cc = unlift x in s + c + cc + +instance Summation KB2 CFloat where + zero = constant (KB2 0 0 0) + add = kb2Add + into _ x = lift (KB2 x 0 0) + from _ x = let KB2 s c cc = unlift x in s + c + cc + +instance Summation KB2 CDouble where + zero = constant (KB2 0 0 0) + add = kb2Add + into _ x = lift (KB2 x 0 0) + from _ x = let KB2 s c cc = unlift x in s + c + cc + +type instance EltRepr (KB2 a) = ((((), EltRepr a), EltRepr a), EltRepr a) + +instance Elt a => Elt (KB2 a) where + eltType _ = TypeRunit `TypeRpair` eltType (undefined::a) + `TypeRpair` eltType (undefined::a) + `TypeRpair` eltType (undefined::a) + toElt ((((),a),b),c) = KB2 (toElt a) (toElt b) (toElt c) + fromElt (KB2 a b c) = ((((), fromElt a), fromElt b), fromElt c) + +instance Elt a => IsProduct Elt (KB2 a) where + type ProdRepr (KB2 a) = ((((), a), a), a) + toProd _ ((((),a),b),c) = KB2 a b c + fromProd _ (KB2 a b c) = ((((),a),b),c) + prod _ _ = ProdRsnoc $ ProdRsnoc $ ProdRsnoc ProdRunit + +instance (Lift Exp a, Elt (Plain a)) => Lift Exp (KB2 a) where + type Plain (KB2 a) = KB2 (Plain a) + lift (KB2 a b c) = Exp $ Tuple $ NilTup `SnocTup` lift a + `SnocTup` lift b + `SnocTup` lift c + +instance Elt a => Unlift Exp (KB2 (Exp a)) where + unlift t = KB2 (Exp $ SuccTupIdx (SuccTupIdx ZeroTupIdx) `Prj` t) + (Exp $ SuccTupIdx ZeroTupIdx `Prj` t) + (Exp $ ZeroTupIdx `Prj` t) + + +-- | Kahan summation. This is the least accurate of the compensated summation +-- methods. This summation method is included only for completeness. +-- +data Kahan a = Kahan a a + deriving (Show, Typeable) + +-- | Return the result of a Kahan sum. +-- +kahan :: Proxy Kahan +kahan = Proxy + +kahanAdd :: (Num a, IsFloating a) => Exp (Kahan a) -> Exp (Kahan a) -> Exp (Kahan a) +kahanAdd (unlift -> Kahan s1 c1 :: Kahan (Exp a)) (unlift -> Kahan s2 c2) = lift (Kahan s' c') + where + s' = s1 `fadd` y + c' = (s' `fsub` s1) `fsub` y + y = s2 `fsub` c1 `fsub` c2 + +-- instance (Num a, Ord a) => Summation Kahan a where +-- zero = lift $ Kahan (0::Exp a) (0::Exp a) +-- add = kahanAdd +-- into _ x = lift (Kahan x 0) +-- from _ x = let Kahan s _ = unlift x in s + +instance Summation Kahan Float where + zero = constant (Kahan 0 0) + add = kahanAdd + into _ x = lift (Kahan x 0) + from _ x = let Kahan s _ = unlift x in s + +instance Summation Kahan Double where + zero = constant (Kahan 0 0) + add = kahanAdd + into _ x = lift (Kahan x 0) + from _ x = let Kahan s _ = unlift x in s + +instance Summation Kahan CFloat where + zero = constant (Kahan 0 0) + add = kahanAdd + into _ x = lift (Kahan x 0) + from _ x = let Kahan s _ = unlift x in s + +instance Summation Kahan CDouble where + zero = constant (Kahan 0 0) + add = kahanAdd + into _ x = lift (Kahan x 0) + from _ x = let Kahan s _ = unlift x in s + +type instance EltRepr (Kahan a) = (((), EltRepr a), EltRepr a) + +instance Elt a => Elt (Kahan a) where + eltType _ = TypeRunit `TypeRpair` eltType (undefined::a) + `TypeRpair` eltType (undefined::a) + toElt (((),a),b) = Kahan (toElt a) (toElt b) + fromElt (Kahan a b) = (((), fromElt a), fromElt b) + +instance Elt a => IsProduct Elt (Kahan a) where + type ProdRepr (Kahan a) = (((), a), a) + toProd _ (((),a),b) = Kahan a b + fromProd _ (Kahan a b) = (((),a),b) + prod _ _ = ProdRsnoc $ ProdRsnoc ProdRunit + +instance (Lift Exp a, Elt (Plain a)) => Lift Exp (Kahan a) where + type Plain (Kahan a) = Kahan (Plain a) + lift (Kahan a b) = Exp $ Tuple $ NilTup `SnocTup` lift a + `SnocTup` lift b + +instance Elt a => Unlift Exp (Kahan (Exp a)) where + unlift t = Kahan (Exp $ SuccTupIdx ZeroTupIdx `Prj` t) + (Exp $ ZeroTupIdx `Prj` t) + diff --git a/src/Data/Array/Accelerate/Numeric/Sum/Arithmetic.hs b/src/Data/Array/Accelerate/Numeric/Sum/Arithmetic.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/Sum/Arithmetic.hs @@ -0,0 +1,37 @@ +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE NoImplicitPrelude #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.Sum.Arithmetic +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.Sum.Arithmetic ( + + fadd, fsub, fmul, + +) where + +import Data.Array.Accelerate + +import qualified Data.Array.Accelerate.Numeric.Sum.LLVM.Native as Native +import qualified Data.Array.Accelerate.Numeric.Sum.LLVM.PTX as PTX + + +infixl 6 `fadd` +fadd :: (Num a, IsFloating a) => Exp a -> Exp a -> Exp a +fadd = Native.fadd $ PTX.fadd (+) + +infixl 6 `fsub` +fsub :: (Num a, IsFloating a) => Exp a -> Exp a -> Exp a +fsub = Native.fsub $ PTX.fsub (-) + +infixl 7 `fmul` +fmul :: (Num a, IsFloating a) => Exp a -> Exp a -> Exp a +fmul = Native.fmul $ PTX.fmul (*) + diff --git a/src/Data/Array/Accelerate/Numeric/Sum/LLVM/Native.hs b/src/Data/Array/Accelerate/Numeric/Sum/LLVM/Native.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/Sum/LLVM/Native.hs @@ -0,0 +1,58 @@ +{-# LANGUAGE CPP #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.Sum.LLVM.Native +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.Sum.LLVM.Native ( + + fadd, fsub, fmul, + +) where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Type + +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND +import Data.Array.Accelerate.LLVM.CodeGen.Sugar +import Data.Array.Accelerate.LLVM.Native.Foreign as A +import qualified Data.Array.Accelerate.Numeric.Sum.LLVM.Prim as Prim +#endif + +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND +wrap2 :: (Elt a, Elt b, Elt c) + => String -- name of the operation + -> IRFun1 Native () ((a, b) -> c) -- foreign implementation + -> (Exp a -> Exp b -> Exp c) -- fallback implementation + -> Exp a + -> Exp b + -> Exp c +wrap2 str f g = A.curry (foreignExp (ForeignExp str f) (A.uncurry g)) +#endif + +fadd :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND +fadd = wrap2 "fadd" (Prim.fadd floatingType) +#else +fadd = id +#endif + +fsub :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND +fsub = wrap2 "fsub" (Prim.fsub floatingType) +#else +fsub = id +#endif + +fmul :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND +fmul = wrap2 "fmul" (Prim.fmul floatingType) +#else +fmul = id +#endif + diff --git a/src/Data/Array/Accelerate/Numeric/Sum/LLVM/PTX.hs b/src/Data/Array/Accelerate/Numeric/Sum/LLVM/PTX.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/Sum/LLVM/PTX.hs @@ -0,0 +1,58 @@ +{-# LANGUAGE CPP #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.Sum.LLVM.PTX +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.Sum.LLVM.PTX ( + + fadd, fsub, fmul, + +) where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Type + +#ifdef ACCELERATE_LLVM_PTX_BACKEND +import Data.Array.Accelerate.LLVM.CodeGen.Sugar +import Data.Array.Accelerate.LLVM.PTX.Foreign as A +import qualified Data.Array.Accelerate.Numeric.Sum.LLVM.Prim as Prim +#endif + +#ifdef ACCELERATE_LLVM_PTX_BACKEND +wrap2 :: (Elt a, Elt b, Elt c) + => String -- name of the operation + -> IRFun1 PTX () ((a, b) -> c) -- foreign implementation + -> (Exp a -> Exp b -> Exp c) -- fallback implementation + -> Exp a + -> Exp b + -> Exp c +wrap2 str f g = A.curry (foreignExp (ForeignExp str f) (A.uncurry g)) +#endif + +fadd :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a +#ifdef ACCELERATE_LLVM_PTX_BACKEND +fadd = wrap2 "fadd" (Prim.fadd floatingType) +#else +fadd = id +#endif + +fsub :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a +#ifdef ACCELERATE_LLVM_PTX_BACKEND +fsub = wrap2 "fsub" (Prim.fsub floatingType) +#else +fsub = id +#endif + +fmul :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a +#ifdef ACCELERATE_LLVM_PTX_BACKEND +fmul = wrap2 "fmul" (Prim.fmul floatingType) +#else +fmul = id +#endif + diff --git a/src/Data/Array/Accelerate/Numeric/Sum/LLVM/Prim.hs b/src/Data/Array/Accelerate/Numeric/Sum/LLVM/Prim.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Numeric/Sum/LLVM/Prim.hs @@ -0,0 +1,88 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE TemplateHaskell #-} +{-# LANGUAGE ViewPatterns #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.Sum.LLVM.Prim +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.Sum.LLVM.Prim ( + + fadd, fsub, fmul, + +) where + +import Data.Array.Accelerate.Type +import Data.Array.Accelerate.Error + +import Data.Array.Accelerate.LLVM.CodeGen.Downcast ( downcast ) +import Data.Array.Accelerate.LLVM.CodeGen.IR ( IR(..), Operands(..), IROP(..) ) +import Data.Array.Accelerate.LLVM.CodeGen.Monad ( CodeGen, freshName, instr_ ) +import Data.Array.Accelerate.LLVM.CodeGen.Sugar ( IROpenFun1(..) ) +import qualified Data.Array.Accelerate.LLVM.CodeGen.Arithmetic as A +import qualified LLVM.AST.Type.Name as A +import qualified LLVM.AST.Type.Operand as A +import qualified LLVM.AST.Type.Representation as A + +import LLVM.AST.Instruction +import LLVM.AST.Name +import LLVM.AST.Operand +import LLVM.AST.Type + + +-- | As (+), but don't allow potentially unsafe floating-point optimisations. +-- +fadd :: FloatingType a -> IROpenFun1 arch env aenv ((a,a) -> a) +fadd t = IRFun1 $ A.uncurry (binop FAdd t) + +-- | As (-), but don't allow potentially unsafe floating-point optimisations. +-- +fsub :: FloatingType a -> IROpenFun1 arch env aenv ((a,a) -> a) +fsub t = IRFun1 $ A.uncurry (binop FSub t) + +-- | As (*), but don't allow potentially unsafe floating-point optimisations. +-- +fmul :: FloatingType a -> IROpenFun1 arch env aenv ((a,a) -> a) +fmul t = IRFun1 $ A.uncurry (binop FMul t) + +binop :: (FastMathFlags -> Operand -> Operand -> InstructionMetadata -> Instruction) -> FloatingType a -> IR a -> IR a -> CodeGen (IR a) +binop f t (op t -> x) (op t -> y) = do + r <- instr (downcast t) (f fmf (downcast x) (downcast y) md) + return (upcast t r) + + +-- Prim +-- ---- + +md :: InstructionMetadata +md = [] + +fmf :: FastMathFlags +#if MIN_VERSION_llvm_hs_pure(6,0,0) +fmf = noFastMathFlags +#else +fmf = NoFastMathFlags +#endif + +fresh :: CodeGen Name +fresh = downcast <$> freshName + +instr :: Type -> Instruction -> CodeGen Operand +instr ty ins = do + name <- fresh + instr_ (name := ins) + return (LocalReference ty name) + +upcast :: FloatingType t -> Operand -> IR t +upcast TypeFloat{} (LocalReference (FloatingPointType FloatFP) (UnName x)) = IR $ OP_Float (A.LocalReference A.type' (A.UnName x)) +upcast TypeDouble{} (LocalReference (FloatingPointType DoubleFP) (UnName x)) = IR $ OP_Double (A.LocalReference A.type' (A.UnName x)) +upcast TypeCFloat{} (LocalReference (FloatingPointType FloatFP) (UnName x)) = IR $ OP_CFloat (A.LocalReference A.type' (A.UnName x)) +upcast TypeCDouble{} (LocalReference (FloatingPointType DoubleFP) (UnName x)) = IR $ OP_CDouble (A.LocalReference A.type' (A.UnName x)) +upcast _ _ = $internalError "upcast" "expected local reference" + diff --git a/test/Backend.hs b/test/Backend.hs deleted file mode 100644 --- a/test/Backend.hs +++ /dev/null @@ -1,58 +0,0 @@ -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE CPP #-} - -module Backend where - -import Data.Array.Accelerate as A -import qualified Data.Array.Accelerate.Interpreter as I -#if ACCELERATE_LLVM_NATIVE_BACKEND -import qualified Data.Array.Accelerate.LLVM.Native as CPU -#endif -#if ACCELERATE_LLVM_PTX_BACKEND -import qualified Data.Array.Accelerate.LLVM.PTX as PTX -#endif - -data Backend = Interpreter -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND - | Native -#endif -#ifdef ACCELERATE_LLVM_PTX_BACKEND - | PTX -#endif - -instance Show Backend where - show Interpreter = "interpreter" -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND - show Native = "llvm-cpu" -#endif -#ifdef ACCELERATE_LLVM_PTX_BACKEND - show PTX = "llvm-ptx" -#endif - -{-# INLINE run #-} -run :: Arrays a => Backend -> Acc a -> a -run Interpreter = I.run -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND -run Native = CPU.run -#endif -#ifdef ACCELERATE_LLVM_PTX_BACKEND -run PTX = PTX.run -#endif - - -{-# INLINE run1 #-} -run1 :: (Arrays a, Arrays b) => Backend -> (Acc a -> Acc b) -> a -> b -run1 Interpreter f = I.run1 f -#ifdef ACCELERATE_LLVM_NATIVE_BACKEND -run1 Native f = CPU.run1 f -#endif -#ifdef ACCELERATE_LLVM_PTX_BACKEND -run1 PTX f = PTX.run1 f -#endif - -{-# INLINE run2 #-} -run2 :: (Arrays a, Arrays b, Arrays c) => Backend -> (Acc a -> Acc b -> Acc c) -> a -> b -> c -run2 b f x y = go (x,y) - where - !go = run1 b (A.uncurry f) - diff --git a/test/Hedgehog/Gen/Array.hs b/test/Hedgehog/Gen/Array.hs deleted file mode 100644 --- a/test/Hedgehog/Gen/Array.hs +++ /dev/null @@ -1,21 +0,0 @@ - -module Hedgehog.Gen.Array where - -import Data.Array.Accelerate as A -import Prelude as P - -import Hedgehog ( Gen ) -import qualified Hedgehog.Gen as Gen -import qualified Hedgehog.Range as Range - - --- Generate an array of the given shape --- -genArray - :: (Shape sh, Elt e) - => sh - -> Gen e - -> Gen (Array sh e) -genArray sh gen = - fromList sh <$> Gen.list (Range.singleton (arraySize sh)) gen - diff --git a/test/Hedgehog/Gen/Shape.hs b/test/Hedgehog/Gen/Shape.hs deleted file mode 100644 --- a/test/Hedgehog/Gen/Shape.hs +++ /dev/null @@ -1,22 +0,0 @@ -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE TypeOperators #-} - -module Hedgehog.Gen.Shape where - -import Data.Array.Accelerate as A - -import Hedgehog ( Gen, Range ) -import Hedgehog.Gen ( int ) - - --- Generate a randomly sized shape of the given dimensionality --- -class GenShape sh where - genShape :: Range Int -> Gen sh - -instance GenShape Z where - genShape _ = return Z - -instance GenShape sh => GenShape (sh :. Int) where - genShape r = (:.) <$> genShape r <*> int r - diff --git a/test/Level2.hs b/test/Level2.hs deleted file mode 100644 --- a/test/Level2.hs +++ /dev/null @@ -1,63 +0,0 @@ -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE RankNTypes #-} -{-# LANGUAGE TemplateHaskell #-} - -module Level2 ( tests ) where - -import Backend -import Similar - -import Data.Array.Accelerate as A -import Data.Array.Accelerate.Data.Complex as A -import Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level2 - -import Hedgehog -import Hedgehog.Gen.Array -import qualified Hedgehog.Gen as Gen -import qualified Hedgehog.Range as Range - -import Data.String -import Text.Printf -import Prelude as P - - -tests :: Backend -> IO Bool -tests backend - = checkParallel - $ Group (fromString $ printf "Tests.Level2.%s" (show backend)) - [ ("gemv.float32", test_gemv backend r f32) - , ("gemv.float64", test_gemv backend r f64) - , ("gemv.complex32", test_gemv backend r c32) - , ("gemv.complex64", test_gemv backend r c64) - ] - where - r = Range.linearFrom 0 1 128 - f32 = Gen.float (Range.linearFracFrom 0 (-1) 1) - f64 = Gen.double (Range.linearFracFrom 0 (-1) 1) - c32 = (:+) <$> f32 <*> f32 - c64 = (:+) <$> f64 <*> f64 - -test_gemv - :: (Numeric e, Similar e) - => Backend - -> Range Int - -> Gen e - -> Property -test_gemv backend r g = - property $ do - alpha <- forAll g - m <- forAll (Gen.int r) - n <- forAll (Gen.int r) - opA <- forAll (Gen.element [N,T,H]) - vecx <- forAll (genArray (Z :. n) g) - matA <- forAll $ case opA of - N -> genArray (Z :. m :. n) g - _ -> genArray (Z :. n :. m) g - -- - let t = gemv (constant alpha) opA (use matA) (use vecx) - -- - run Interpreter t ~~~ run backend t - - - diff --git a/test/Level3.hs b/test/Level3.hs deleted file mode 100644 --- a/test/Level3.hs +++ /dev/null @@ -1,65 +0,0 @@ -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE RankNTypes #-} -{-# LANGUAGE TemplateHaskell #-} - -module Level3 ( tests ) where - -import Backend -import Similar - -import Data.Array.Accelerate as A -import Data.Array.Accelerate.Data.Complex as A -import Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level3 - -import Hedgehog -import Hedgehog.Gen.Array -import qualified Hedgehog.Gen as Gen -import qualified Hedgehog.Range as Range - -import Data.String -import Text.Printf -import Prelude as P - - -tests :: Backend -> IO Bool -tests backend - = checkParallel - $ Group (fromString $ printf "Tests.Level3.%s" (show backend)) - [ ("gemm.float32", test_gemm backend r f32) - , ("gemm.float64", test_gemm backend r f64) - , ("gemm.complex32", test_gemm backend r c32) - , ("gemm.complex64", test_gemm backend r c64) - ] - where - r = Range.linearFrom 0 1 64 - f32 = Gen.float (Range.linearFracFrom 0 (-1) 1) - f64 = Gen.double (Range.linearFracFrom 0 (-1) 1) - c32 = (:+) <$> f32 <*> f32 - c64 = (:+) <$> f64 <*> f64 - -test_gemm - :: (Numeric e, Similar e) - => Backend - -> Range Int - -> Gen e - -> Property -test_gemm backend r g = - property $ do - alpha <- forAll g - m <- forAll (Gen.int r) - n <- forAll (Gen.int r) - k <- forAll (Gen.int r) - opA <- forAll (Gen.element [N,T,H]) - opB <- forAll (Gen.element [N,T,H]) - matA <- forAll $ case opA of - N -> genArray (Z :. m :. k) g - _ -> genArray (Z :. k :. m) g - matB <- forAll $ case opB of - N -> genArray (Z :. k :. n) g - _ -> genArray (Z :. n :. k) g - -- - let t = gemm (constant alpha) opA (use matA) opB (use matB) - -- - run Interpreter t ~~~ run backend t - diff --git a/test/Main.hs b/test/Main.hs deleted file mode 100644 --- a/test/Main.hs +++ /dev/null @@ -1,27 +0,0 @@ -{-# LANGUAGE CPP #-} - -module Main where - -import Backend - -import System.IO -import qualified Level3 -import qualified Level2 - -main :: IO () -main = do - hSetBuffering stdout LineBuffering - hSetBuffering stderr LineBuffering - - sequence_ - [ return True -#if ACCELERATE_LLVM_NATIVE_BACKEND - , Level2.tests Native - , Level3.tests Native -#endif -#if ACCELERATE_LLVM_PTX_BACKEND - , Level2.tests PTX - , Level3.tests PTX -#endif - ] - diff --git a/test/Similar.hs b/test/Similar.hs deleted file mode 100644 --- a/test/Similar.hs +++ /dev/null @@ -1,68 +0,0 @@ -{-# LANGUAGE DefaultSignatures #-} -{-# LANGUAGE FlexibleContexts #-} -{-# LANGUAGE TypeOperators #-} - -module Similar where - -import Data.Complex -import Data.Array.Accelerate ( Array, Shape, Z, (:.)(..), arrayShape, toList ) - -import Hedgehog -import Hedgehog.Internal.Source ( HasCallStack, withFrozenCallStack ) - - -infix 4 ~~~ -(~~~) :: (MonadTest m, Similar a, Show a, HasCallStack) => a -> a -> m () -a ~~~ b = withFrozenCallStack $ Sim a === Sim b - - -data Sim a = Sim a - -instance Similar a => Eq (Sim a) where - Sim a == Sim b = a ~= b - -instance Show a => Show (Sim a) where - show (Sim a) = show a - - --- A class of things that support almost-equality, so that we can disregard --- small amounts of floating-point round-off error. --- -class Similar a where - {-# INLINE (~=) #-} - (~=) :: a -> a -> Bool - default (~=) :: Eq a => a -> a -> Bool - (~=) = (==) - -infix 4 ~= - -instance Similar Float where (~=) = absRelTol 0.00005 0.005 -instance Similar Double where (~=) = absRelTol 0.00005 0.005 - -instance Similar e => Similar (Complex e) where - (r1 :+ i1) ~= (r2 :+ i2) = r1 ~= r2 && i1 ~= i2 - -instance Similar Z -instance (Eq sh, Eq sz) => Similar (sh:.sz) - -instance Similar a => Similar [a] where - [] ~= [] = True - (x:xs) ~= (y:ys) = x ~= y && xs ~= ys - _ ~= _ = False - -instance (Similar e, Eq sh, Shape sh) => Similar (Array sh e) where - a1 ~= a2 = arrayShape a1 == arrayShape a2 - && toList a1 ~= toList a2 - - -{-# INLINEABLE absRelTol #-} -absRelTol :: RealFloat a => a -> a -> a -> a -> Bool -absRelTol epsilonAbs epsilonRel u v - | isInfinite u - && isInfinite v = True - | isNaN u - && isNaN v = True - | abs (u-v) < epsilonAbs = True - | abs u > abs v = abs ((u-v) / u) < epsilonRel - | otherwise = abs ((v-u) / v) < epsilonRel - diff --git a/test/Test/BLAS.hs b/test/Test/BLAS.hs new file mode 100644 --- /dev/null +++ b/test/Test/BLAS.hs @@ -0,0 +1,33 @@ +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ViewPatterns #-} +-- | +-- Module : Test.BLAS +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Test.BLAS ( testBLAS ) + where + +import Test.Util +import Test.BLAS.Level2 +import Test.BLAS.Level3 + +import Test.Tasty + + +testBLAS :: Run -> TestTree +testBLAS run = + testGroup "BLAS" + [ test_level2 run + , test_level3 run + ] + diff --git a/test/Test/BLAS/Level2.hs b/test/Test/BLAS/Level2.hs new file mode 100644 --- /dev/null +++ b/test/Test/BLAS/Level2.hs @@ -0,0 +1,69 @@ +{-# LANGUAGE RankNTypes #-} +-- | +-- Module : Test.BLAS.Level2 +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Test.BLAS.Level2 ( test_level2 ) + where + +import Test.Util as Gen + +import Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level2 + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Test.Similar +import qualified Data.Array.Accelerate.Interpreter as I + +import Hedgehog +import qualified Hedgehog.Gen as Gen +import qualified Hedgehog.Range as Range + +import Test.Tasty +import Test.Tasty.Hedgehog + + +test_gemv + :: (Numeric e, Similar e) + => Run + -> Range Int + -> Gen e + -> Property +test_gemv run r e = + property $ do + alpha <- forAll e + m <- forAll (Gen.int r) + n <- forAll (Gen.int r) + opA <- forAll (Gen.element [N,T,H]) + vecx <- forAll (Gen.array (Z :. n) e) + matA <- forAll $ case opA of + N -> Gen.array (Z :. m :. n) e + _ -> Gen.array (Z :. n :. m) e + -- + let t = gemv (constant alpha) opA (use matA) (use vecx) + -- + I.run t ~~~ run t + + +test_level2 :: Run -> TestTree +test_level2 run = + testGroup "Level2" + [ testGroup "gemv" + [ testProperty "Float" $ test_gemv run r f32 + , testProperty "Double" $ test_gemv run r f64 + , testProperty "Complex Float" $ test_gemv run r c32 + , testProperty "Complex Double" $ test_gemv run r c64 + ] + ] + where + r = Range.linearFrom 0 1 128 + f32 = floating :: Gen Float + f64 = floating :: Gen Double + c32 = complex f32 + c64 = complex f64 + diff --git a/test/Test/BLAS/Level3.hs b/test/Test/BLAS/Level3.hs new file mode 100644 --- /dev/null +++ b/test/Test/BLAS/Level3.hs @@ -0,0 +1,73 @@ +{-# LANGUAGE RankNTypes #-} +-- | +-- Module : Test.BLAS.Level3 +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Test.BLAS.Level3 ( test_level3 ) + where + +import Test.Util as Gen + +import Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level3 + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Test.Similar +import qualified Data.Array.Accelerate.Interpreter as I + +import Hedgehog +import qualified Hedgehog.Gen as Gen +import qualified Hedgehog.Range as Range + +import Test.Tasty +import Test.Tasty.Hedgehog + + +test_gemm + :: (Numeric e, Similar e) + => Run + -> Range Int + -> Gen e + -> Property +test_gemm run r e = + property $ do + alpha <- forAll e + m <- forAll (Gen.int r) + n <- forAll (Gen.int r) + k <- forAll (Gen.int r) + opA <- forAll (Gen.element [N,T,H]) + opB <- forAll (Gen.element [N,T,H]) + matA <- forAll $ case opA of + N -> Gen.array (Z :. m :. k) e + _ -> Gen.array (Z :. k :. m) e + matB <- forAll $ case opB of + N -> Gen.array (Z :. k :. n) e + _ -> Gen.array (Z :. n :. k) e + -- + let t = gemm (constant alpha) opA (use matA) opB (use matB) + -- + I.run t ~~~ run t + + +test_level3 :: Run -> TestTree +test_level3 run = + testGroup "Level3" + [ testGroup "gemm" + [ testProperty "Float" $ test_gemm run r f32 + , testProperty "Double" $ test_gemm run r f64 + , testProperty "Complex Float" $ test_gemm run r c32 + , testProperty "Complex Double" $ test_gemm run r c64 + ] + ] + where + r = Range.linearFrom 0 1 64 + f32 = floating :: Gen Float + f64 = floating :: Gen Double + c32 = complex f32 + c64 = complex f64 + diff --git a/test/Test/Util.hs b/test/Test/Util.hs new file mode 100644 --- /dev/null +++ b/test/Test/Util.hs @@ -0,0 +1,41 @@ +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ViewPatterns #-} +-- | +-- Module : Test.Util +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Test.Util where + +import Data.Array.Accelerate ( Acc, Arrays, Array, Shape, Elt, fromList ) +import Data.Array.Accelerate.Array.Sugar ( size ) +import Data.Array.Accelerate.Trafo ( Afunction, AfunctionR ) +import Data.Array.Accelerate.Data.Complex + +import Hedgehog +import qualified Hedgehog.Gen as Gen +import qualified Hedgehog.Range as Range +import Prelude as P + + +type Run = forall a. Arrays a => Acc a -> a +type RunN = forall f. Afunction f => f -> AfunctionR f + +floating :: P.RealFloat a => Gen a +floating = Gen.realFloat (Range.linearFracFrom 0 (-1) 1) + +complex :: Gen a -> Gen (Complex a) +complex f = (:+) <$> f <*> f + +array :: (Shape sh, Elt e) => sh -> Gen e -> Gen (Array sh e) +array sh gen = fromList sh <$> Gen.list (Range.singleton (size sh)) gen + diff --git a/test/TestNative.hs b/test/TestNative.hs new file mode 100644 --- /dev/null +++ b/test/TestNative.hs @@ -0,0 +1,19 @@ +-- | +-- Module : TestNative +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module TestNative where + +import Test.BLAS +import Test.Tasty +import Data.Array.Accelerate.LLVM.Native as CPU + +main :: IO () +main = defaultMain (testBLAS CPU.run) + diff --git a/test/TestPTX.hs b/test/TestPTX.hs new file mode 100644 --- /dev/null +++ b/test/TestPTX.hs @@ -0,0 +1,19 @@ +-- | +-- Module : TestPTX +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module TestPTX where + +import Test.BLAS +import Test.Tasty +import Data.Array.Accelerate.LLVM.PTX as PTX + +main :: IO () +main = defaultMain (testBLAS PTX.run) +