accelerate-blas (empty) → 0.1.0.0
raw patch · 38 files changed
+3301/−0 lines, 38 filesdep +acceleratedep +accelerate-blasdep +accelerate-llvmsetup-changed
Dependencies added: accelerate, accelerate-blas, accelerate-llvm, accelerate-llvm-native, accelerate-llvm-ptx, base, blas-hs, bytestring, containers, criterion, cublas, cuda, deepseq, file-embed, hedgehog, hmatrix, llvm-hs-pure, mtl, mwc-random, mwc-random-accelerate, storable-complex
Files
- CHANGELOG.md +11/−0
- Data/Array/Accelerate/Numeric/LinearAlgebra.hs +172/−0
- Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level1.hs +139/−0
- Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level2.hs +92/−0
- Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level3.hs +109/−0
- Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Base.hs +106/−0
- Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level2.hs +72/−0
- Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level3.hs +87/−0
- Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Base.hs +82/−0
- Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Context.hs +76/−0
- Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level2.hs +131/−0
- Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level3.hs +113/−0
- Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Twine.hs +176/−0
- Data/Array/Accelerate/Numeric/LinearAlgebra/Type.hs +91/−0
- Data/Array/Accelerate/Numeric/Sum.hs +330/−0
- Data/Array/Accelerate/Numeric/Sum/Arithmetic.hs +37/−0
- Data/Array/Accelerate/Numeric/Sum/LLVM/Native.hs +58/−0
- Data/Array/Accelerate/Numeric/Sum/LLVM/PTX.hs +58/−0
- Data/Array/Accelerate/Numeric/Sum/LLVM/Prim.hs +83/−0
- LICENSE +30/−0
- README.md +31/−0
- Setup.hs +2/−0
- accelerate-blas.cabal +189/−0
- bench/Accelerate.hs +187/−0
- bench/Extra.hs +33/−0
- bench/HMatrix.hs +122/−0
- bench/Main.hs +24/−0
- cbits/twine_f32.c +60/−0
- cbits/twine_f64.c +60/−0
- cubits/twine_f32.ptx +108/−0
- cubits/twine_f64.ptx +108/−0
- test/Backend.hs +58/−0
- test/Hedgehog/Gen/Array.hs +21/−0
- test/Hedgehog/Gen/Shape.hs +22/−0
- test/Level2.hs +63/−0
- test/Level3.hs +65/−0
- test/Main.hs +27/−0
- test/Similar.hs +68/−0
+ CHANGELOG.md view
@@ -0,0 +1,11 @@+# Revision history for accelerate-blas++Notable changes to the project will be documented in this file.++The format is based on [Keep a Changelog](http://keepachangelog.com/) and the+project adheres to the [Haskell Package Versioning+Policy (PVP)](https://pvp.haskell.org)++## 0.1.0.0 - 2017-09-21+ * First version. Released on an unsuspecting world.+
+ 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+
+ 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)+
+ 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+
+ 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+
+ Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Base.hs view
@@ -0,0 +1,106 @@+{-# 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 view
@@ -0,0 +1,72 @@+{-# 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 view
@@ -0,0 +1,87 @@+{-# 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 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.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 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+
+ Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level2.hs view
@@ -0,0 +1,131 @@+{-# 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 view
@@ -0,0 +1,113 @@+{-# 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 view
@@ -0,0 +1,176 @@+{-# 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 view
@@ -0,0 +1,91 @@+{-# 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 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 )+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 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 (*)+
+ 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+
+ 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+
+ Data/Array/Accelerate/Numeric/Sum/LLVM/Prim.hs view
@@ -0,0 +1,83 @@+{-# 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"+
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2016, Trevor L. McDonell++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Trevor L. McDonell nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ README.md view
@@ -0,0 +1,31 @@+# Numeric linear algebra in Accelerate++[](https://travis-ci.org/tmcdonell/accelerate-blas)++Linear systems, matrix decompositions, and other numerical computations for use+in Accelerate. Most operations are implemented efficiently via FFI calls to BLAS+and LAPACK. For details on Accelerate, refer to the [main repository][GitHub].++Please get in touch to let me know which missing operations you would like see+added to the library. Contributions are also welcome!+++## FFI bindings++ * **accelerate-llvm-native:** FFI bindings are provided by the [blas-hs] package,+ which has several options for which underlying BLAS library to link against;+ see that package for setup details.++ * **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+ [cuBLAS]: http://docs.nvidia.com/cuda/cublas/index.html+
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ accelerate-blas.cabal view
@@ -0,0 +1,189 @@+name: accelerate-blas+version: 0.1.0.0+synopsis: Numeric Linear Algebra in Accelerate+description:+ Linear systems, matrix decompositions, and other numerical computations for+ use in Accelerate. Most operations are implemented efficiently via FFI calls+ to BLAS and LAPACK+ .+ For further information refer to the main /Accelerate/ package:+ <http://hackage.haskell.org/package/accelerate>++license: BSD3+license-file: LICENSE+author: Trevor L. McDonell+maintainer: tmcdonell@cse.unsw.edu.au+category: Math+build-type: Simple+extra-source-files: CHANGELOG.md+cabal-version: >=1.10++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+ Default: True++Flag llvm-ptx+ Description: Enable the LLVM PTX backend for NVIDIA GPUs+ Default: True++library+ default-language: Haskell2010+ exposed-modules:+ Data.Array.Accelerate.Numeric.Sum+ Data.Array.Accelerate.Numeric.LinearAlgebra+ Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level1+ Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level2+ Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level3++ other-modules:+ Data.Array.Accelerate.Numeric.LinearAlgebra.Type+ Data.Array.Accelerate.Numeric.Sum.Arithmetic+ Data.Array.Accelerate.Numeric.Sum.LLVM.Native+ Data.Array.Accelerate.Numeric.Sum.LLVM.PTX++ build-depends:+ base >= 4.7 && < 4.11+ , accelerate >= 1.0 && < 1.2++ ghc-options:+ -O2+ -Wall++ 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+ , 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++ other-modules:+ Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Base+ Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level2+ 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+ , bytestring >= 0.9+ , containers >= 0.5+ , cublas >= 0.3+ , cuda >= 0.8+ , file-embed >= 0.0.10+ , llvm-hs-pure >= 4.0+ , 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+ default-language: Haskell2010+ type: exitcode-stdio-1.0+ hs-source-dirs: test+ main-is: Main.hs+ other-modules:+ Backend+ Hedgehog.Gen.Array+ Hedgehog.Gen.Shape+ Level2+ Level3+ Similar++ build-depends:+ base >= 4.7 && < 4.11+ , accelerate >= 1.0 && < 1.2+ , accelerate-blas+ , hedgehog >= 0.5++ 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++ if flag(llvm-ptx)+ CPP-options: -DACCELERATE_LLVM_PTX_BACKEND+ build-depends:+ accelerate-llvm-ptx >= 1.0 && < 1.2+++benchmark accelerate-blas-bench+ default-language: Haskell2010+ type: exitcode-stdio-1.0+ hs-source-dirs: bench+ main-is: Main.hs+ other-modules:+ Accelerate+ Extra+ HMatrix++ build-depends:+ base >= 4.7 && < 4.11+ , accelerate >= 1.0 && < 1.2+ , accelerate-blas+ , 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++ if flag(llvm-ptx)+ CPP-options: -DACCELERATE_LLVM_PTX_BACKEND+ build-depends:+ accelerate-llvm-ptx >= 1.0 && < 1.2++source-repository head+ type: git+ location: https://github.com/tmcdonell/accelerate-blas++source-repository this+ type: git+ tag: 0.1.0.0+ location: https://github.com/tmcdonell/accelerate-blas++-- vim: nospell
+ bench/Accelerate.hs view
@@ -0,0 +1,187 @@+{-# 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/Extra.hs view
@@ -0,0 +1,33 @@+{-# 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 view
@@ -0,0 +1,122 @@+{-# 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 view
@@ -0,0 +1,24 @@+{-# 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 view
@@ -0,0 +1,60 @@+/*+ * 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 view
@@ -0,0 +1,60 @@+/*+ * 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 view
@@ -0,0 +1,108 @@+//+// 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 view
@@ -0,0 +1,108 @@+//+// 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;+}++
+ test/Backend.hs view
@@ -0,0 +1,58 @@+{-# 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 view
@@ -0,0 +1,21 @@++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 view
@@ -0,0 +1,22 @@+{-# 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 :: Monad m => 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 view
@@ -0,0 +1,63 @@+{-# 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 view
@@ -0,0 +1,65 @@+{-# 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 view
@@ -0,0 +1,27 @@+{-# 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 view
@@ -0,0 +1,68 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE DefaultSignatures #-}+{-# 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 (Sim 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+