accelerate-blas 0.1.0.1 → 0.2.0.0
raw patch · 65 files changed
+2453/−3108 lines, 65 filesdep +tastydep +tasty-hedgehogdep −storable-complexdep ~acceleratedep ~accelerate-llvmdep ~accelerate-llvm-native
Dependencies added: tasty, tasty-hedgehog
Dependencies removed: storable-complex
Dependency ranges changed: accelerate, accelerate-llvm, accelerate-llvm-native, accelerate-llvm-ptx, base, llvm-hs-pure
Files
- CHANGELOG.md +10/−0
- Data/Array/Accelerate/Numeric/LinearAlgebra.hs +0/−172
- Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level1.hs +0/−139
- Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level2.hs +0/−92
- Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level3.hs +0/−109
- Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Base.hs +0/−106
- Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level2.hs +0/−72
- Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level3.hs +0/−87
- Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Base.hs +0/−82
- Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Context.hs +0/−76
- Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level2.hs +0/−131
- Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level3.hs +0/−113
- Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Twine.hs +0/−176
- Data/Array/Accelerate/Numeric/LinearAlgebra/Type.hs +0/−91
- Data/Array/Accelerate/Numeric/Sum.hs +0/−330
- Data/Array/Accelerate/Numeric/Sum/Arithmetic.hs +0/−37
- Data/Array/Accelerate/Numeric/Sum/LLVM/Native.hs +0/−58
- Data/Array/Accelerate/Numeric/Sum/LLVM/PTX.hs +0/−58
- Data/Array/Accelerate/Numeric/Sum/LLVM/Prim.hs +0/−83
- README.md +0/−6
- accelerate-blas.cabal +127/−64
- bench/Accelerate.hs +0/−187
- bench/Bench/Accelerate.hs +123/−0
- bench/Bench/HMatrix.hs +126/−0
- bench/Bench/Util.hs +41/−0
- bench/BenchHMatrix.hs +18/−0
- bench/BenchNative.hs +19/−0
- bench/BenchPTX.hs +19/−0
- bench/Extra.hs +0/−33
- bench/HMatrix.hs +0/−122
- bench/Main.hs +0/−24
- cbits/twine_f32.c +0/−60
- cbits/twine_f64.c +0/−60
- cubits/twine_f32.ptx +0/−108
- cubits/twine_f64.ptx +0/−108
- src/Data/Array/Accelerate/Numeric/LinearAlgebra.hs +172/−0
- src/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level1.hs +139/−0
- src/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level2.hs +92/−0
- src/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level3.hs +109/−0
- src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Base.hs +48/−0
- src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level2.hs +54/−0
- src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level3.hs +67/−0
- src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Base.hs +82/−0
- src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Context.hs +76/−0
- src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level2.hs +113/−0
- src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level3.hs +93/−0
- src/Data/Array/Accelerate/Numeric/LinearAlgebra/Type.hs +100/−0
- src/Data/Array/Accelerate/Numeric/Sum.hs +330/−0
- src/Data/Array/Accelerate/Numeric/Sum/Arithmetic.hs +37/−0
- src/Data/Array/Accelerate/Numeric/Sum/LLVM/Native.hs +58/−0
- src/Data/Array/Accelerate/Numeric/Sum/LLVM/PTX.hs +58/−0
- src/Data/Array/Accelerate/Numeric/Sum/LLVM/Prim.hs +88/−0
- test/Backend.hs +0/−58
- test/Hedgehog/Gen/Array.hs +0/−21
- test/Hedgehog/Gen/Shape.hs +0/−22
- test/Level2.hs +0/−63
- test/Level3.hs +0/−65
- test/Main.hs +0/−27
- test/Similar.hs +0/−68
- test/Test/BLAS.hs +33/−0
- test/Test/BLAS/Level2.hs +69/−0
- test/Test/BLAS/Level3.hs +73/−0
- test/Test/Util.hs +41/−0
- test/TestNative.hs +19/−0
- test/TestPTX.hs +19/−0
CHANGELOG.md view
@@ -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
− Data/Array/Accelerate/Numeric/LinearAlgebra.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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-
− Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level1.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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.------ <https://software.intel.com/en-us/mkl-developer-reference-c-cblas-sdot>----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--- \]------ <https://software.intel.com/en-us/mkl-developer-reference-c-cblas-dotu>----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--- \]------ <https://software.intel.com/en-us/mkl-developer-reference-c-cblas-dotc>----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)\|\).------ <https://software.intel.com/en-us/mkl-developer-reference-c-cblas-asum>----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.------ <https://software.intel.com/en-us/mkl-developer-reference-c-cblas-i-amax>----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.------ <https://software.intel.com/en-us/mkl-developer-reference-c-cblas-i-amin>----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)-
− Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level2.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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@------ <https://software.intel.com/en-us/mkl-developer-reference-c-cblas-gemv>----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-
− Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level3.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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@------ <https://software.intel.com/en-us/mkl-developer-reference-c-cblas-gemm>----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-
− Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Base.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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 ()-
− Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level2.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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-
− Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level3.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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-
− Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Base.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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-
− Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Context.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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.------ <http://docs.nvidia.com/cuda/cublas/index.html#thread-safety2>----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-
− Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level2.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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-
− Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level3.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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-
− Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Twine.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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)-
− Data/Array/Accelerate/Numeric/LinearAlgebra/Type.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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)-
− Data/Array/Accelerate/Numeric/Sum.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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: <https://hackage.haskell.org/package/math-functions>-----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)-
− Data/Array/Accelerate/Numeric/Sum/Arithmetic.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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 (*)-
− Data/Array/Accelerate/Numeric/Sum/LLVM/Native.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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-
− Data/Array/Accelerate/Numeric/Sum/LLVM/PTX.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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-
− Data/Array/Accelerate/Numeric/Sum/LLVM/Prim.hs
@@ -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 <tmcdonell@cse.unsw.edu.au>--- 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"-
README.md view
@@ -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
accelerate-blas.cabal view
@@ -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
− bench/Accelerate.hs
@@ -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)-
+ bench/Bench/Accelerate.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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))+ ]+
+ bench/Bench/HMatrix.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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))+ ]+
+ bench/Bench/Util.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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)+
+ bench/BenchHMatrix.hs view
@@ -0,0 +1,18 @@+-- |+-- Module : BenchHMatrix+-- Copyright : [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- Stability : experimental+-- Portability : non-portable (GHC extensions)+--++module BenchHMatrix where++import Bench.HMatrix+import Criterion.Main++main :: IO ()+main = defaultMain bench_hmatrix+
+ bench/BenchNative.hs view
@@ -0,0 +1,19 @@+-- |+-- Module : BenchNative+-- Copyright : [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- 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)+
+ bench/BenchPTX.hs view
@@ -0,0 +1,19 @@+-- |+-- Module : BenchPTX+-- Copyright : [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- 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)+
− bench/Extra.hs
@@ -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)-
− bench/HMatrix.hs
@@ -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))- ]-
− bench/Main.hs
@@ -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- ]-
− cbits/twine_f32.c
@@ -1,60 +0,0 @@-/*- * Module : Twine- * Copyright : [2016] Trevor L. McDonell- * License : BSD3- *- * Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>- * 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 <complex.h>-#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-
− cbits/twine_f64.c
@@ -1,60 +0,0 @@-/*- * Module : Twine- * Copyright : [2016] Trevor L. McDonell- * License : BSD3- *- * Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>- * 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 <complex.h>-#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-
− cubits/twine_f32.ptx
@@ -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;-}--
− cubits/twine_f64.ptx
@@ -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;-}--
+ src/Data/Array/Accelerate/Numeric/LinearAlgebra.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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+
+ src/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level1.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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.+--+-- <https://software.intel.com/en-us/mkl-developer-reference-c-cblas-sdot>+--+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+-- \]+--+-- <https://software.intel.com/en-us/mkl-developer-reference-c-cblas-dotu>+--+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+-- \]+--+-- <https://software.intel.com/en-us/mkl-developer-reference-c-cblas-dotc>+--+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)\|\).+--+-- <https://software.intel.com/en-us/mkl-developer-reference-c-cblas-asum>+--+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.+--+-- <https://software.intel.com/en-us/mkl-developer-reference-c-cblas-i-amax>+--+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.+--+-- <https://software.intel.com/en-us/mkl-developer-reference-c-cblas-i-amin>+--+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)+
+ src/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level2.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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@+--+-- <https://software.intel.com/en-us/mkl-developer-reference-c-cblas-gemv>+--+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+
+ src/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level3.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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@+--+-- <https://software.intel.com/en-us/mkl-developer-reference-c-cblas-gemm>+--+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+
+ src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Base.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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+
+ src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level2.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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+
+ src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level3.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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+
+ src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Base.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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+
+ src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Context.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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.+--+-- <http://docs.nvidia.com/cuda/cublas/index.html#thread-safety2>+--+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+
+ src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level2.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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+
+ src/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level3.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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+
+ src/Data/Array/Accelerate/Numeric/LinearAlgebra/Type.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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)+
+ src/Data/Array/Accelerate/Numeric/Sum.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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: <https://hackage.haskell.org/package/math-functions>+--++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)+
+ src/Data/Array/Accelerate/Numeric/Sum/Arithmetic.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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 (*)+
+ src/Data/Array/Accelerate/Numeric/Sum/LLVM/Native.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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+
+ src/Data/Array/Accelerate/Numeric/Sum/LLVM/PTX.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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+
+ src/Data/Array/Accelerate/Numeric/Sum/LLVM/Prim.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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"+
− test/Backend.hs
@@ -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)-
− test/Hedgehog/Gen/Array.hs
@@ -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-
− test/Hedgehog/Gen/Shape.hs
@@ -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-
− test/Level2.hs
@@ -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---
− test/Level3.hs
@@ -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-
− test/Main.hs
@@ -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- ]-
− test/Similar.hs
@@ -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-
+ test/Test/BLAS.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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+ ]+
+ test/Test/BLAS/Level2.hs view
@@ -0,0 +1,69 @@+{-# LANGUAGE RankNTypes #-}+-- |+-- Module : Test.BLAS.Level2+-- Copyright : [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- 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+
+ test/Test/BLAS/Level3.hs view
@@ -0,0 +1,73 @@+{-# LANGUAGE RankNTypes #-}+-- |+-- Module : Test.BLAS.Level3+-- Copyright : [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- 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+
+ test/Test/Util.hs view
@@ -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 <tmcdonell@cse.unsw.edu.au>+-- 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+
+ test/TestNative.hs view
@@ -0,0 +1,19 @@+-- |+-- Module : TestNative+-- Copyright : [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- 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)+
+ test/TestPTX.hs view
@@ -0,0 +1,19 @@+-- |+-- Module : TestPTX+-- Copyright : [2017] Trevor L. McDonell+-- License : BSD3+--+-- Maintainer : Trevor L. McDonell <tmcdonell@cse.unsw.edu.au>+-- 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)+