diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 --- /dev/null +++ b/CHANGELOG.md @@ -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. + diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/LinearAlgebra.hs @@ -0,0 +1,172 @@ +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE ViewPatterns #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra ( + + -- * Types + Numeric, Scalar, Vector, Matrix, + + -- * Products + -- ** Vector-vector + (<.>), + (><), + + -- ** Matrix-vector + (#>), (<#), + + -- ** Matrix-matrix + (<>), + + -- * Diagonal + identity, diagonal, + +) where + +import Data.Array.Accelerate as A + +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type +import Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level1 +import Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level2 +import Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level3 + + +-- Level 1 +-- ------- + +-- | An infix synonym for 'dotu'. +-- +-- >>> let a = fromList (Z:.4) [1..] +-- >>> let b = fromList (Z:.4) [-2,0,1,1] +-- >>> a <.> b +-- Scalar Z [5.0] +-- +-- >>> let c = fromList (Z:.2) [1:+1, 1:+0] +-- >>> let d = fromList (Z:.2) [1:+0, 1:+(-1)] +-- >>> c <.> d +-- Scalar Z [2.0 :+ 0.0] +-- +infixr 8 <.> +(<.>) :: Numeric e => Acc (Vector e) -> Acc (Vector e) -> Acc (Scalar e) +(<.>) = dotu + + +-- | Outer product of two vectors +-- +-- >>> let a = fromList (Z :. 3) [1,2,3] +-- >>> let b = fromList (Z :. 3) [5,2,3] +-- >>> a >< b +-- Matrix (Z :. 3 :. 3) +-- [ 5.0, 2.0, 3.0 +-- , 10.0, 4.0, 6.0 +-- , 15.0, 6.0, 9.0 ] +-- +infixr 8 >< +(><) :: Numeric e => Acc (Vector e) -> Acc (Vector e) -> Acc (Matrix e) +(><) x y = xc <> yr + where + xc = reshape (index2 (length x) 1) x + yr = reshape (index2 1 (length y)) y + + +-- Level 2 +-- ------- + +-- | Dense matrix-vector product +-- +-- >>> let m = fromList (Z :. 2 :. 3) [1..] +-- >>> m +-- Matrix (Z :. 2 :. 3) +-- [ 1.0, 2.0, 3.0 +-- , 4.0, 5.0, 6.0 ] +-- +-- >>> let x = fromList (Z :. 3) [10,20,30] +-- +-- >>> m #> x +-- Vector (Z :. 2) [140.0,320.0] +-- +-- See 'gemv' for a more general version of this operation. +-- +infixr 8 #> +(#>) :: Numeric e => Acc (Matrix e) -> Acc (Vector e) -> Acc (Vector e) +(#>) m x = gemv 1 N m x + + +-- | Dense vector-matrix product +-- +-- >>> let m = fromList (Z :. 2 :. 3) [1..] +-- >>> m +-- Matrix (Z :. 2 :. 3) +-- [1.0,2.0,3.0, +-- 4.0,5.0,6.0] +-- +-- >>> let v = fromList (Z :. 2) [5,10] +-- +-- >>> v <# m +-- Vector (Z :. 3) [45.0,60.0,75.0] +-- +-- See 'gemv' for a more general version of this operation. +-- +infixr 8 <# +(<#) :: Numeric e => Acc (Vector e) -> Acc (Matrix e) -> Acc (Vector e) +(<#) x m = gemv 1 T m x + + +-- Level 3 +-- ------- + +-- | Dense matrix-matrix product +-- +-- >>> let a = fromList (Z :. 3 :. 5) [1..] +-- >>> a +-- Matrix (Z:.3:.5) +-- [ 1.0, 2.0, 3.0, 4.0, 5.0 +-- , 6.0, 7.0, 8.0, 9.0, 10.0 +-- , 11.0, 12.0, 13.0, 14.0, 15.0 ] +-- +-- >>> let b = fromList (Z :. 5 :. 2) [1,3, 0,2, -1,5, 7,7, 6,0] +-- >>> b +-- Matrix (Z :. 5 :. 2) +-- [ 1.0, 3.0 +-- , 0.0, 2.0 +-- , -1.0, 5.0 +-- , 7.0, 7.0 +-- , 6.0, 0.0 ] +-- +-- >>> a <> b +-- Matrix (Z :. 3 :. 2) +-- [ 56.0, 50.0 +-- , 121.0, 135.0 +-- , 186.0, 220.0 ] +-- +-- See 'gemm' for a more general version of this operation. +-- +infixr 8 <> +(<>) :: Numeric e => Acc (Matrix e) -> Acc (Matrix e) -> Acc (Matrix e) +(<>) matA matB = gemm 1 N matA N matB + + +-- | Create a square identity matrix of the given dimension +-- +identity :: Num e => Exp Int -> Acc (Matrix e) +identity n = diagonal (fill (index1 n) 1) + +-- | Create a square matrix with the given diagonal +-- +diagonal :: Num e => Acc (Vector e) -> Acc (Matrix e) +diagonal v = + let n = length v + zeros = fill (index2 n n) 0 + in + permute const zeros (\(unindex1 -> i) -> index2 i i) v + diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level1.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level1.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level1.hs @@ -0,0 +1,139 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE ScopedTypeVariables #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level1 +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- +-- Level 1 (vector-vector) BLAS operations. +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level1 ( + + -- Types + Numeric, Vector, + + -- Level1 operations + sdot, + dotu, + dotc, + asum, + amax, + amin, + +) where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Data.Complex as A +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type + + +-- | Computes a vector-vector dot product, using double precision accumulation +-- of the intermediate result. Includes a scalar (initial) value to be added to +-- the inner product. +-- +-- +-- +sdot :: forall e. Numeric e => Exp e -> Acc (Vector e) -> Acc (Vector e) -> Acc (Scalar e) +sdot z xs ys = + case numericR :: NumericR e of + NumericRfloat32 -> map toFloating $ dsdot (toFloating z) (map toFloating xs) (map toFloating ys) + NumericRfloat64 -> dsdot z xs ys + NumericRcomplex32 -> map d2f $ zsdot (f2d z) (map f2d xs) (map f2d ys) + NumericRcomplex64 -> zsdot z xs ys + where + dsdot :: Exp Double -> Acc (Vector Double) -> Acc (Vector Double) -> Acc (Scalar Double) + dsdot z' xs' ys' = fold (+) z' (zipWith (*) xs' ys') + + zsdot :: Exp (Complex Double) -> Acc (Vector (Complex Double)) -> Acc (Vector (Complex Double)) -> Acc (Scalar (Complex Double)) + zsdot z' xs' ys' = fold (+) z' (zipWith (*) xs' ys') + + f2d :: Exp (Complex Float) -> Exp (Complex Double) + f2d c = lift (toFloating (real c) :+ toFloating (imag c)) + + d2f :: Exp (Complex Double) -> Exp (Complex Float) + d2f c = lift (toFloating (real c) :+ toFloating (imag c)) + + +-- | Computes a vector-vector dot product +-- +-- \[ +-- res = \sum_i x_i * y_i +-- \] +-- +-- +-- +dotu :: Numeric e => Acc (Vector e) -> Acc (Vector e) -> Acc (Scalar e) +dotu xs ys = fold (+) 0 (zipWith (*) xs ys) + + +-- | Computes a dot product of a conjugated vector with another vector +-- +-- \[ +-- res = \sum_i \mathrm{conj}(x_i) * y_i +-- \] +-- +-- +-- +dotc :: forall e. Numeric (Complex e) + => Acc (Vector (Complex e)) + -> Acc (Vector (Complex e)) + -> Acc (Scalar (Complex e)) +dotc xs ys = + case numericR :: NumericR (Complex e) of + NumericRcomplex32 -> dotu (map conjugate xs) ys + NumericRcomplex64 -> dotu (map conjugate xs) ys + + +-- | Computes the sum of magnitudes of the vector elements. For complex values, +-- this is given by \(\sum_i \|\mathrm{real}(x_i)\| + \|\mathrm{imag}(x_i)\|\). +-- +-- +-- +asum :: forall e. Numeric e => Acc (Vector e) -> Acc (Scalar (NumericBaseT e)) +asum = + case numericR :: NumericR e of + NumericRfloat32 -> sum . map abs + NumericRfloat64 -> sum . map abs + NumericRcomplex32 -> sum . map mag + NumericRcomplex64 -> sum . map mag + where + mag c = abs (real c) + abs (imag c) + + +-- | Return the index of the element with the maximum absolute value. +-- +-- +-- +amax :: forall e. Numeric e => Acc (Vector e) -> Acc (Scalar Int) +amax = + case numericR :: NumericR e of + NumericRfloat32 -> map (indexHead . fst) . fold1 cmp . indexed . map abs + NumericRfloat64 -> map (indexHead . fst) . fold1 cmp . indexed . map abs + NumericRcomplex32 -> map (indexHead . fst) . fold1 cmp . indexed . map mag + NumericRcomplex64 -> map (indexHead . fst) . fold1 cmp . indexed . map mag + where + cmp ix iy = snd ix > snd iy ? ( ix, iy ) + mag c = abs (real c) + abs (imag c) + +-- | Return the index of the element with the minimum absolute value. +-- +-- +-- +amin :: forall e. Numeric e => Acc (Vector e) -> Acc (Scalar Int) +amin = + case numericR :: NumericR e of + NumericRfloat32 -> map (indexHead . fst) . fold1 cmp . indexed . map abs + NumericRfloat64 -> map (indexHead . fst) . fold1 cmp . indexed . map abs + NumericRcomplex32 -> map (indexHead . fst) . fold1 cmp . indexed . map mag + NumericRcomplex64 -> map (indexHead . fst) . fold1 cmp . indexed . map mag + where + cmp ix iy = snd ix < snd iy ? ( ix, iy ) + mag c = abs (real c) + abs (imag c) + diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level2.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level2.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level2.hs @@ -0,0 +1,92 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ViewPatterns #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level2 +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- +-- Level 2 (matrix-vector) BLAS operations. +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level2 ( + + -- Types + Numeric, Vector, Matrix, Transpose(..), + + -- Operations + gemv, + +) where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Smart as A +import Data.Array.Accelerate.Data.Complex as A +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type + +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND +import qualified Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level2 as CPU +#endif +#ifdef ACCELERATE_LLVM_PTX_BACKEND +import qualified Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level2 as PTX +#endif + + +-- | Computes the matrix-vector product of a general matrix. +-- +-- \[ +-- y = \alpha * \mathrm{op}(A) * x +-- \] +-- +-- where: +-- +-- * 'shape' \(\mathrm{op}(A)\) @= Z :. m :. n@ +-- * 'shape' \(x\) @= Z :. n@ +-- * 'shape' \(y\) @= Z :. m@ +-- +-- +-- +gemv :: forall e. Numeric e + => Exp e -- ^ \( \alpha \) + -> Transpose -- ^ Operation to apply to A + -> Acc (Matrix e) -- ^ A + -> Acc (Vector e) -- ^ x + -> Acc (Vector e) -- ^ y +gemv alpha opA matA x = go (lift (unit alpha, matA, x)) + where + go = +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND + foreignAcc (CPU.gemv opA) $ +#endif +#ifdef ACCELERATE_LLVM_PTX_BACKEND + foreignAcc (PTX.gemv opA) $ +#endif + (\(unatup3 -> (_, arr, brr)) -> mXv arr brr) + + -- General matrix-vector multiply in pure Accelerate. This is probably not + -- efficient. + -- + mXv :: Acc (Matrix e) -> Acc (Vector e) -> Acc (Vector e) + mXv arr brr + = fold (+) 0 + $ zipWith (\a b -> alpha * a * b) arr' brr' + where + Z :. m :. _ = unlift (shape arr') :: Z :. Exp Int :. Exp Int + + brr' = replicate (lift (Z :. m :. All)) brr + arr' = case opA of + N -> arr + T -> transpose arr + H -> case numericR :: NumericR e of + NumericRcomplex32 -> map conjugate (transpose arr) + NumericRcomplex64 -> map conjugate (transpose arr) + _ -> transpose arr + diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level3.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level3.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/LinearAlgebra/BLAS/Level3.hs @@ -0,0 +1,109 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE GADTs #-} +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ViewPatterns #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level3 +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- +-- Level 3 (matrix-matrix) BLAS operations. +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.BLAS.Level3 ( + + -- Types + Numeric, Matrix, Transpose(..), + + -- Matrix-matrix operations + gemm, + +) where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Smart as A +import Data.Array.Accelerate.Data.Complex as A +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type + +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND +import qualified Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level3 as CPU +#endif +#ifdef ACCELERATE_LLVM_PTX_BACKEND +import qualified Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level3 as PTX +#endif + + +-- | General matrix-matrix multiply +-- +-- \[ +-- C = \alpha * \mathrm{op}(A) * \mathrm{op}(B) +-- \] +-- +-- where: +-- +-- * 'shape' \(\mathrm{op}(A)\) @= Z :. m :. k@ +-- * 'shape' \(\mathrm{op}(B)\) @= Z :. k :. n@ +-- * 'shape' \(C\) @= Z :. m :. n@ +-- +-- +-- +gemm :: forall e. Numeric e + => Exp e -- ^ \( \alpha \) + -> Transpose -- ^ operation to apply to A + -> Acc (Matrix e) -- ^ A + -> Transpose -- ^ operation to apply to B + -> Acc (Matrix e) -- ^ B + -> Acc (Matrix e) -- ^ C +gemm alpha opA matA opB matB = go (lift (unit alpha, matA, matB)) + where + go = +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND + foreignAcc (CPU.gemm opA opB) $ +#endif +#ifdef ACCELERATE_LLVM_PTX_BACKEND + foreignAcc (PTX.gemm opA opB) $ +#endif + (\(unatup3 -> (_, arr, brr)) -> mXm arr brr) + + -- General dense matrix-matrix multiply written in pure Accelerate. This is + -- not efficient due to the memory access patterns. We could probably + -- improve this a little bit with a divide-and-conquer algorithm, for + -- example, but using a foreign implementation will be best. + -- + mXm :: Acc (Matrix e) -> Acc (Matrix e) -> Acc (Matrix e) + mXm arr brr + = fold (+) 0 + $ zipWith (\a b -> alpha * a * b) arrRepl brrRepl + where + Z :. rowsA :. _ = unlift (shape arr') :: Z :. Exp Int :. Exp Int + Z :. colsB :. _ = unlift (shape brr') :: Z :. Exp Int :. Exp Int + -- + arrRepl = replicate (lift $ Z :. All :. colsB :. All) arr' + brrRepl = replicate (lift $ Z :. rowsA :. All :. All) brr' + + -- apply opA + arr' = case opA of + N -> arr + T -> transpose arr + H -> case numericR :: NumericR e of + NumericRcomplex32 -> map conjugate (transpose arr) + NumericRcomplex64 -> map conjugate (transpose arr) + _ -> transpose arr + + -- apply opB and transpose at the same time, which is required for this + -- algorithm + brr' = case opB of + N -> transpose brr + T -> brr + H -> case numericR :: NumericR e of + NumericRcomplex32 -> map conjugate brr + NumericRcomplex64 -> map conjugate brr + _ -> brr + diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Base.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Base.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Base.hs @@ -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 +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Base + where + +import Data.Array.Accelerate.Data.Complex +import Data.Array.Accelerate.Array.Sugar ( Array(..), EltRepr ) +import Data.Array.Accelerate.Array.Data +import Data.Array.Accelerate.Array.Unique +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type + +import Foreign.Marshal.Alloc +import Foreign.Ptr +import Foreign.Storable +import Foreign.Storable.Complex ( ) + +import qualified Blas.Primitive.Types as C + + +encodeTranspose :: Transpose -> C.Transpose +encodeTranspose N = C.NoTrans +encodeTranspose T = C.Trans +encodeTranspose H = C.ConjTrans + + +withArray + :: forall sh e b. Numeric e + => Array sh e + -> (ArrayPtrs (EltRepr e) -> IO b) + -> IO b +withArray (Array _ adata) = withArrayData (numericR::NumericR e) adata + +withArrayData + :: NumericR e + -> ArrayData (EltRepr e) + -> (ArrayPtrs (EltRepr e) -> IO b) + -> IO b +withArrayData NumericRfloat32 (AD_Float ua) f = withUniqueArrayPtr ua f +withArrayData NumericRfloat64 (AD_Double ua) f = withUniqueArrayPtr ua f +withArrayData NumericRcomplex32 (AD_Pair ad1 ad2) f + | AD_Pair AD_Unit (AD_Float ua_re) <- ad1 + , AD_Float ua_im <- ad2 + = withUniqueArrayPtr ua_re $ \p_re -> + withUniqueArrayPtr ua_im $ \p_im -> + f (((),p_re), p_im) + +withArrayData NumericRcomplex64 (AD_Pair ad1 ad2) f + | AD_Pair AD_Unit (AD_Double ua_re) <- ad1 + , AD_Double ua_im <- ad2 + = withUniqueArrayPtr ua_re $ \p_re -> + withUniqueArrayPtr ua_im $ \p_im -> + f (((),p_re), p_im) + + +interleave + :: forall e b. (Storable e, Numeric (Complex e)) + => ArrayPtrs (EltRepr (Complex e)) + -> Int + -> (Ptr (Complex e) -> IO b) + -> IO b +interleave (((), p_re), p_im) n k = do + allocaBytesAligned (n * sizeOf (undefined::Complex e)) 16 $ \p_cplx -> do + () <- case numericR :: NumericR (Complex e) of + NumericRcomplex32 -> c_interleave_f32 0 n p_cplx p_re p_im + NumericRcomplex64 -> c_interleave_f64 0 n p_cplx p_re p_im + -- + k p_cplx + + +deinterleave + :: forall e. (Storable e, Numeric (Complex e)) + => ArrayPtrs (EltRepr (Complex e)) + -> Ptr (Complex e) + -> Int + -> IO () +deinterleave (((), p_re), p_im) p_cplx n = + case numericR :: NumericR (Complex e) of + NumericRcomplex32 -> c_deinterleave_f32 0 n p_re p_im p_cplx + NumericRcomplex64 -> c_deinterleave_f64 0 n p_re p_im p_cplx + + +foreign import ccall unsafe "interleave_f32" + c_interleave_f32 :: Int -> Int -> Ptr (Complex Float) -> Ptr Float -> Ptr Float -> IO () + +foreign import ccall unsafe "interleave_f64" + c_interleave_f64 :: Int -> Int -> Ptr (Complex Double) -> Ptr Double -> Ptr Double -> IO () + +foreign import ccall unsafe "deinterleave_f32" + c_deinterleave_f32 :: Int -> Int -> Ptr Float -> Ptr Float -> Ptr (Complex Float) -> IO () + +foreign import ccall unsafe "deinterleave_f64" + c_deinterleave_f64 :: Int -> Int -> Ptr Double -> Ptr Double -> Ptr (Complex Double) -> IO () + diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level2.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level2.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level2.hs @@ -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 +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level2 + where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Data.Complex +import Data.Array.Accelerate.LLVM.Native.Foreign +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Base + +import Foreign.Marshal.Alloc +import Foreign.Storable +import Foreign.Storable.Complex ( ) + +import qualified Blas.Primitive.Types as C +import qualified Blas.Primitive.Unsafe as C + + +gemv :: forall e. Numeric e + => Transpose + -> ForeignAcc ((Scalar e, Matrix e, Vector e) -> Vector e) +gemv opA = ForeignAcc "native.gemv" gemv' + where + gemv' (alpha, matA, vecx) = do + let + Z :. rowsA :. colsA = arrayShape matA + Z :. sizeX = arrayShape vecx + + sizeA = rowsA * colsA + sizeY = case opA of + N -> rowsA + _ -> colsA + + opA' = encodeTranspose opA + alpha' = indexArray alpha Z + -- + vecy <- allocateRemote (Z :. sizeY) :: LLVM Native (Vector e) + () <- liftIO $ do + withArray matA $ \ptr_A -> do + withArray vecx $ \ptr_x -> do + withArray vecy $ \ptr_y -> do + case numericR :: NumericR e of + NumericRfloat32 -> C.sgemv C.RowMajor opA' rowsA colsA alpha' ptr_A colsA ptr_x 1 0 ptr_y 1 + NumericRfloat64 -> C.dgemv C.RowMajor opA' rowsA colsA alpha' ptr_A colsA ptr_x 1 0 ptr_y 1 + -- + NumericRcomplex32 -> do + allocaBytesAligned (sizeY * sizeOf (undefined::Complex e)) 16 $ \ptr_y' -> do + interleave ptr_A sizeA $ \ptr_A' -> do + interleave ptr_x sizeX $ \ptr_x' -> do + C.cgemv C.RowMajor opA' rowsA colsA alpha' ptr_A' colsA ptr_x' 1 0 ptr_y' 1 + deinterleave ptr_y ptr_y' sizeY + -- + NumericRcomplex64 -> do + allocaBytesAligned (sizeY * sizeOf (undefined::Complex e)) 16 $ \ptr_y' -> do + interleave ptr_A sizeA $ \ptr_A' -> do + interleave ptr_x sizeX $ \ptr_x' -> do + C.zgemv C.RowMajor opA' rowsA colsA alpha' ptr_A' colsA ptr_x' 1 0 ptr_y' 1 + deinterleave ptr_y ptr_y' sizeY + -- + return vecy + diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level3.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level3.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/Native/Level3.hs @@ -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 +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Level3 + where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Data.Complex +import Data.Array.Accelerate.LLVM.Native.Foreign +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.Native.Base + +import Foreign.Marshal.Alloc +import Foreign.Storable +import Foreign.Storable.Complex ( ) + +import qualified Blas.Primitive.Types as C +import qualified Blas.Primitive.Unsafe as C + + +-- TODO: check whether it is faster to compute this as column-major order: +-- +-- https://www.christophlassner.de/using-blas-from-c-with-row-major-data.html +-- +gemm :: forall e. Numeric e + => Transpose + -> Transpose + -> ForeignAcc ((Scalar e, Matrix e, Matrix e) -> Matrix e) +gemm opA opB = ForeignAcc "native.gemm" gemm' + where + gemm' (alpha, matA, matB) = do + let + Z :. rowsA :. colsA = arrayShape matA + Z :. rowsB :. colsB = arrayShape matB + + sizeA = rowsA * colsA + sizeB = rowsB * colsB + sizeC = m * n + + (m,k) = case opA of + N -> (rowsA, colsA) + _ -> (colsA, rowsA) + n = case opB of + N -> colsB + _ -> rowsB + + lda = colsA + ldb = colsB + + opA' = encodeTranspose opA + opB' = encodeTranspose opB + alpha' = indexArray alpha Z + -- + matC <- allocateRemote (Z :. m :. n) :: LLVM Native (Matrix e) + () <- liftIO $ do + withArray matA $ \ptr_A -> do + withArray matB $ \ptr_B -> do + withArray matC $ \ptr_C -> do + case numericR :: NumericR e of + NumericRfloat32 -> C.sgemm C.RowMajor opA' opB' m n k alpha' ptr_A lda ptr_B ldb 0 ptr_C n + NumericRfloat64 -> C.dgemm C.RowMajor opA' opB' m n k alpha' ptr_A lda ptr_B ldb 0 ptr_C n + -- + NumericRcomplex32 -> do + allocaBytesAligned (sizeC * sizeOf (undefined::Complex e)) 16 $ \ptr_C' -> do + interleave ptr_A sizeA $ \ptr_A' -> do + interleave ptr_B sizeB $ \ptr_B' -> do + C.cgemm C.RowMajor opA' opB' m n k alpha' ptr_A' lda ptr_B' ldb 0 ptr_C' n + deinterleave ptr_C ptr_C' sizeC + -- + NumericRcomplex64 -> do + allocaBytesAligned (sizeC * sizeOf (undefined::Complex e)) 16 $ \ptr_C' -> do + interleave ptr_A sizeA $ \ptr_A' -> do + interleave ptr_B sizeB $ \ptr_B' -> do + C.zgemm C.RowMajor opA' opB' m n k alpha' ptr_A' lda ptr_B' ldb 0 ptr_C' n + deinterleave ptr_C ptr_C' sizeC + -- + return matC + diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Base.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Base.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Base.hs @@ -0,0 +1,82 @@ +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base + where + +import Data.Array.Accelerate.Lifetime +import Data.Array.Accelerate.Array.Sugar ( Array(..), EltRepr ) +import Data.Array.Accelerate.Array.Data +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type + +import Data.Array.Accelerate.LLVM.PTX.Foreign + +import Foreign.CUDA.Ptr ( DevicePtr ) +import qualified Foreign.CUDA.BLAS as C + + +type family DevicePtrs e :: * +type instance DevicePtrs () = () +type instance DevicePtrs Float = DevicePtr Float +type instance DevicePtrs Double = DevicePtr Double +type instance DevicePtrs (a,b) = (DevicePtrs a, DevicePtrs b) + + +encodeTranspose :: Transpose -> C.Operation +encodeTranspose N = C.N +encodeTranspose T = C.T +encodeTranspose H = C.C + + +withArray + :: forall sh e b. Numeric e + => Array sh e + -> Stream + -> (DevicePtrs (EltRepr e) -> LLVM PTX b) + -> LLVM PTX b +withArray (Array _ adata) s k = withArrayData (numericR::NumericR e) adata s k + +withArrayData + :: NumericR e + -> ArrayData (EltRepr e) + -> Stream + -> (DevicePtrs (EltRepr e) -> LLVM PTX b) + -> LLVM PTX b +withArrayData NumericRfloat32 ad s k = + withDevicePtr ad $ \p -> do + r <- k p + e <- checkpoint s + return (Just e,r) +withArrayData NumericRfloat64 ad s k = + withDevicePtr ad $ \p -> do + r <- k p + e <- checkpoint s + return (Just e, r) +withArrayData NumericRcomplex32 (AD_Pair (AD_Pair AD_Unit ad1) ad2) s k = + withDevicePtr ad1 $ \p1 -> + withDevicePtr ad2 $ \p2 -> do + r <- k (((), p1), p2) + e <- checkpoint s + return (Just e, (Just e, r)) +withArrayData NumericRcomplex64 (AD_Pair (AD_Pair AD_Unit ad1) ad2) s k = + withDevicePtr ad1 $ \p1 -> + withDevicePtr ad2 $ \p2 -> do + r <- k (((), p1), p2) + e <- checkpoint s + return (Just e, (Just e, r)) + +withLifetime' :: Lifetime a -> (a -> LLVM PTX b) -> LLVM PTX b +withLifetime' l k = do + r <- k (unsafeGetValue l) + liftIO $ touchLifetime l + return r + diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Context.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Context.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Context.hs @@ -0,0 +1,76 @@ +{-# LANGUAGE MagicHash #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Context +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Context ( + + withBLAS + +) where + +import Data.Array.Accelerate.Lifetime +import Data.Array.Accelerate.LLVM.PTX +import Data.Array.Accelerate.LLVM.PTX.Foreign +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base + +import Control.Monad.State +import Control.Concurrent.MVar +import Data.IntMap.Strict ( IntMap ) +import System.IO.Unsafe +import qualified Data.IntMap.Strict as IM + +import qualified Foreign.CUDA.Driver.Context as CUDA +import qualified Foreign.CUDA.BLAS as BLAS + +import GHC.Ptr +import GHC.Base +import Prelude hiding ( lookup ) + + +-- Execute an operation with a cuBLAS handle appropriate for the current +-- execution context. +-- +-- Initial creation of the context is an atomic operation, but subsequently +-- multiple threads may use the context concurrently. +-- +-- +-- +withBLAS :: (BLAS.Handle -> LLVM PTX b) -> LLVM PTX b +withBLAS k = do + lc <- gets (deviceContext . ptxContext) + h <- liftIO $ + withLifetime lc $ \ctx -> + modifyMVar handles $ \im -> + let key = toKey ctx in + case IM.lookup key im of + -- handle does not exist yet; create it and add to the global + -- state for reuse + Nothing -> do + h <- BLAS.create + l <- newLifetime h + -- BLAS.setPointerMode h BLAS.Device + BLAS.setAtomicsMode h BLAS.Allowed + addFinalizer lc $ modifyMVar handles (\im' -> return (IM.delete key im', ())) + addFinalizer l $ BLAS.destroy h + return ( IM.insert key l im, l ) + + -- return existing handle + Just h -> return (im, h) + -- + withLifetime' h k + + +toKey :: CUDA.Context -> IM.Key +toKey (CUDA.Context (Ptr addr#)) = I# (addr2Int# addr#) + +{-# NOINLINE handles #-} +handles :: MVar (IntMap (Lifetime BLAS.Handle)) +handles = unsafePerformIO $ newMVar IM.empty + diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level2.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level2.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level2.hs @@ -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 +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level2 + where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Array.Sugar ( Array(..) ) +import Data.Array.Accelerate.Data.Complex +import Data.Array.Accelerate.LLVM.PTX.Foreign +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Context +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level3 +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Twine +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type + +import Foreign.Marshal ( with ) +import Foreign.Storable.Complex ( ) + +import qualified Foreign.CUDA.Ptr as CUDA +import qualified Foreign.CUDA.BLAS as BLAS + + +-- NOTE: cuBLAS requires matrices to be stored in column-major order +-- (Fortran-style), but Accelerate uses C-style arrays in row-major order. +-- +-- If the operation is N or T, we can just swap the operation. For +-- conjugate-transpose (H) operations (on complex valued arguments), since there +-- is no conjugate-no-transpose operation, we implement that via 'gemm', which +-- I assume is more efficient than ?geam followed by ?gemv. +-- +gemv :: Numeric e + => Transpose + -> ForeignAcc ((Scalar e, Matrix e, Vector e) -> Vector e) +gemv opA = ForeignAcc "ptx.gemv" (gemv' numericR opA) + +gemv' :: Numeric e + => NumericR e + -> Transpose + -> Stream + -> (Scalar e, Matrix e, Vector e) + -> LLVM PTX (Vector e) +gemv' NumericRcomplex32 H = as_gemm H +gemv' NumericRcomplex64 H = as_gemm H +gemv' _ t = as_gemv t + + +as_gemm + :: Numeric e + => Transpose + -> Stream + -> (Scalar e, Matrix e, Vector e) + -> LLVM PTX (Vector e) +as_gemm opA stream (alpha, matA, Array sh adata) = do + let matB = Array (sh,1) adata + -- + Array (sh',1) vecy <- gemm' opA N stream (alpha, matA, matB) + return (Array sh' vecy) + +as_gemv + :: forall e. Numeric e + => Transpose + -> Stream + -> (Scalar e, Matrix e, Vector e) + -> LLVM PTX (Vector e) +as_gemv opA stream (alpha, matA, vecx) = do + let + Z :. rowsA :. colsA = arrayShape matA + Z :. sizeX = arrayShape vecx + + sizeA = rowsA * colsA + sizeY = case opA of + N -> rowsA + _ -> colsA + + opA' = encodeTranspose + $ case opA of + N -> T + _ -> N + -- + vecy <- allocateRemote (Z :. sizeY) :: LLVM PTX (Vector e) + alpha' <- indexRemote alpha 0 + () <- do + withArray matA stream $ \ptr_A -> do + withArray vecx stream $ \ptr_x -> do + withArray vecy stream $ \ptr_y -> do + withBLAS $ \hdl -> do + case numericR :: NumericR e of + NumericRfloat32 -> liftIO $ + with alpha' $ \ptr_alpha -> + with 0 $ \ptr_beta -> + BLAS.sgemv hdl opA' colsA rowsA ptr_alpha ptr_A colsA ptr_x 1 ptr_beta ptr_y 1 + + NumericRfloat64 -> liftIO $ + with alpha' $ \ptr_alpha -> + with 0 $ \ptr_beta -> + BLAS.dgemv hdl opA' colsA rowsA ptr_alpha ptr_A colsA ptr_x 1 ptr_beta ptr_y 1 + + NumericRcomplex32 -> do + tmpy <- allocateRemote (Z :. sizeY * 2) :: LLVM PTX (Vector Float) + withArray tmpy stream $ \ptr_y' -> do + interleave ptr_A stream sizeA $ \ptr_A' -> do + interleave ptr_x stream sizeX $ \ptr_x' -> do + liftIO $ do + with alpha' $ \ptr_alpha -> + with 0 $ \ptr_beta -> do + BLAS.cgemv hdl opA' colsA rowsA ptr_alpha ptr_A' colsA ptr_x' 1 ptr_beta (CUDA.castDevPtr ptr_y' :: CUDA.DevicePtr (Complex Float)) 1 + deinterleave ptr_y (CUDA.castDevPtr ptr_y' :: CUDA.DevicePtr (Complex Float)) stream sizeY + + NumericRcomplex64 -> do + tmpy <- allocateRemote (Z :. sizeY * 2) :: LLVM PTX (Vector Double) + withArray tmpy stream $ \ptr_y' -> do + interleave ptr_A stream sizeA $ \ptr_A' -> do + interleave ptr_x stream sizeX $ \ptr_x' -> do + liftIO $ do + with alpha' $ \ptr_alpha -> + with 0 $ \ptr_beta -> do + BLAS.zgemv hdl opA' colsA rowsA ptr_alpha ptr_A' colsA ptr_x' 1 ptr_beta (CUDA.castDevPtr ptr_y' :: CUDA.DevicePtr (Complex Double)) 1 + deinterleave ptr_y (CUDA.castDevPtr ptr_y' :: CUDA.DevicePtr (Complex Double)) stream sizeY + -- + return vecy + diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level3.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level3.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Level3.hs @@ -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 +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Level3 + where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Data.Complex +import Data.Array.Accelerate.LLVM.PTX.Foreign +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Context +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Twine +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type + +import Foreign.Marshal ( with ) +import Foreign.Storable.Complex ( ) + +import qualified Foreign.CUDA.Ptr as CUDA +import qualified Foreign.CUDA.BLAS as BLAS + + +-- NOTE: cuBLAS requires that matrices are stored in column-major order +-- (Fortran-style), but Accelerate uses a C-style convention where matrices are +-- stored in row-major order. +-- +-- At least for matrix-matrix multiply, we can get around this problem by making +-- use of the equivalence \( B^T \cdot A^T = (A \cdot B)^T \). +-- +gemm :: Numeric e + => Transpose + -> Transpose + -> ForeignAcc ((Scalar e, Matrix e, Matrix e) -> Matrix e) +gemm opA opB = ForeignAcc "ptx.gemm" (gemm' opA opB) + +gemm' + :: forall e. Numeric e + => Transpose + -> Transpose + -> Stream + -> (Scalar e, Matrix e, Matrix e) + -> LLVM PTX (Matrix e) +gemm' opA opB stream (alpha, matA, matB) = do + let + Z :. rowsA :. colsA = arrayShape matA + Z :. rowsB :. colsB = arrayShape matB + + sizeA = rowsA * colsA + sizeB = rowsB * colsB + sizeC = m * n + + (m,k) = case opA of + N -> (rowsA, colsA) + _ -> (colsA, rowsA) + n = case opB of + N -> colsB + _ -> rowsB + + lda = colsA + ldb = colsB + + opA' = encodeTranspose opA + opB' = encodeTranspose opB + -- + matC <- allocateRemote (Z :. m :. n) :: LLVM PTX (Matrix e) + alpha' <- indexRemote alpha 0 + () <- withArray matA stream $ \ptr_A -> do + withArray matB stream $ \ptr_B -> do + withArray matC stream $ \ptr_C -> do + withBLAS $ \hdl -> do + case numericR :: NumericR e of + NumericRfloat32 -> liftIO $ + with alpha' $ \ptr_alpha -> + with 0 $ \ptr_beta -> + BLAS.sgemm hdl opB' opA' n m k ptr_alpha ptr_B ldb ptr_A lda ptr_beta ptr_C n + + NumericRfloat64 -> liftIO $ + with alpha' $ \ptr_alpha -> + with 0 $ \ptr_beta -> + BLAS.dgemm hdl opB' opA' n m k ptr_alpha ptr_B ldb ptr_A lda ptr_beta ptr_C n + + NumericRcomplex32 -> do + tmpC <- allocateRemote (Z :. sizeC * 2) :: LLVM PTX (Vector Float) + withArray tmpC stream $ \ptr_C' -> do + interleave ptr_A stream sizeA $ \ptr_A' -> do + interleave ptr_B stream sizeB $ \ptr_B' -> do + liftIO $ + with alpha' $ \ptr_alpha -> + with 0 $ \ptr_beta -> + BLAS.cgemm hdl opB' opA' n m k ptr_alpha ptr_B' ldb ptr_A' lda ptr_beta (CUDA.castDevPtr ptr_C') n + deinterleave ptr_C (CUDA.castDevPtr ptr_C' :: CUDA.DevicePtr (Complex Float)) stream sizeC + + NumericRcomplex64 -> do + tmpC <- allocateRemote (Z :. sizeC * 2) :: LLVM PTX (Vector Double) + withArray tmpC stream $ \ptr_C' -> do + interleave ptr_A stream sizeA $ \ptr_A' -> do + interleave ptr_B stream sizeB $ \ptr_B' -> do + liftIO $ + with alpha' $ \ptr_alpha -> + with 0 $ \ptr_beta -> + BLAS.zgemm hdl opB' opA' n m k ptr_alpha ptr_B' ldb ptr_A' lda ptr_beta (CUDA.castDevPtr ptr_C') n + deinterleave ptr_C (CUDA.castDevPtr ptr_C' :: CUDA.DevicePtr (Complex Double)) stream sizeC + + return matC + diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Twine.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Twine.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/LinearAlgebra/LLVM/PTX/Twine.hs @@ -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 +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Twine ( + + interleave, + deinterleave, + +) where + +import Data.Array.Accelerate.Data.Complex +import Data.Array.Accelerate.Array.Sugar ( EltRepr, Vector, Z(..), (:.)(..) ) +import Data.Array.Accelerate.Lifetime +import Data.Array.Accelerate.LLVM.PTX +import Data.Array.Accelerate.LLVM.PTX.Foreign + +import Data.Array.Accelerate.Numeric.LinearAlgebra.Type +import Data.Array.Accelerate.Numeric.LinearAlgebra.LLVM.PTX.Base + +import Control.Concurrent.MVar +import Control.Monad.State +import Data.ByteString ( ByteString ) +import Data.FileEmbed +import Data.IntMap.Strict ( IntMap ) +import Foreign.Storable.Complex ( ) +import System.IO.Unsafe +import qualified Data.IntMap.Strict as IM + +import Foreign.CUDA.Ptr ( DevicePtr ) +import Foreign.CUDA.Analysis +import qualified Foreign.CUDA.Driver as CUDA +import qualified Foreign.CUDA.Driver.Stream as CUDA + +import GHC.Ptr +import GHC.Base +import Prelude hiding ( lookup ) + + +interleave + :: forall e b. Numeric (Complex e) + => DevicePtrs (EltRepr (Complex e)) + -> Stream + -> Int + -> (DevicePtr (Complex e) -> LLVM PTX b) -- device pointer is in packed representation + -> LLVM PTX b +interleave (((), d_re), d_im) s n k = do + case numericR :: NumericR (Complex e) of + nR@NumericRcomplex32 -> do + cplx <- allocateRemote (Z :. n * 2) :: LLVM PTX (Vector Float) + withTwine nR $ \(_,pack,_) -> do + withArray cplx s $ \d_cplx -> do + withLifetime' s $ \s' -> do + liftIO $ launch pack s' n d_cplx d_re d_im + k (CUDA.castDevPtr d_cplx :: DevicePtr (Complex Float)) + -- + nR@NumericRcomplex64 -> do + cplx <- allocateRemote (Z :. n * 2) :: LLVM PTX (Vector Double) + withTwine nR $ \(_,pack,_) -> do + withArray cplx s $ \d_cplx -> do + withLifetime' s $ \s' -> do + liftIO $ launch pack s' n d_cplx d_re d_im + k (CUDA.castDevPtr d_cplx :: DevicePtr (Complex Double)) + +deinterleave + :: forall e. Numeric (Complex e) + => DevicePtrs (EltRepr (Complex e)) + -> DevicePtr (Complex e) -- in packed representation + -> Stream + -> Int + -> LLVM PTX () +deinterleave (((), d_re), d_im) d_cplx s n = do + case numericR :: NumericR (Complex e) of + nR@NumericRcomplex32 -> do + withTwine nR $ \(_,_,unpack) -> do + withLifetime' s $ \s' -> do + liftIO $ launch unpack s' n d_re d_im (CUDA.castDevPtr d_cplx :: DevicePtr Float) + -- + nR@NumericRcomplex64 -> do + withTwine nR $ \(_,_,unpack) -> do + withLifetime' s $ \s' -> do + liftIO $ launch unpack s' n d_re d_im (CUDA.castDevPtr d_cplx :: DevicePtr Double) + + +withTwine :: NumericR (Complex e) -> ((CUDA.Module, Kernel, Kernel) -> LLVM PTX b) -> LLVM PTX b +withTwine nR k = do + ptx <- gets ptxContext + let lc = deviceContext ptx + prp = deviceProperties ptx + mds = modules nR + -- + mdl <- liftIO $ do + withLifetime lc $ \ctx -> do + modifyMVar mds $ \im -> do + let key = toKey ctx + case IM.lookup key im of + -- Module is not loaded yet; add to the current context and the global + -- state for later reuse + Nothing -> do + mdl <- CUDA.loadData $ case nR of + NumericRcomplex32 -> ptx_twine_f32 + NumericRcomplex64 -> ptx_twine_f64 + pack <- mkKernel "interleave" mdl prp + unpack <- mkKernel "deinterleave" mdl prp + let mkk = (mdl, pack, unpack) + -- + lm <- newLifetime mkk + addFinalizer lc $ modifyMVar mds (\im' -> return (IM.delete key im', ())) + addFinalizer lm $ CUDA.unload mdl + return ( IM.insert key lm im, lm ) + + -- Return existing module + Just lm -> return (im, lm) + -- + withLifetime' mdl k + + +toKey :: CUDA.Context -> IM.Key +toKey (CUDA.Context (Ptr addr#)) = I# (addr2Int# addr#) + + +launch :: Kernel -> CUDA.Stream -> Int -> DevicePtr e -> DevicePtr e -> DevicePtr e -> IO () +launch Kernel{..} s n dx dy dz = + CUDA.launchKernel kernelFun (kernelThreadBlocks n,1,1) (kernelThreadBlockSize,1,1) kernelSharedMemBytes (Just s) + [ CUDA.VArg dx, CUDA.VArg dy, CUDA.VArg dz, CUDA.IArg (fromIntegral n) ] + +mkKernel :: String -> CUDA.Module -> CUDA.DeviceProperties -> IO Kernel +mkKernel name mdl prp = do + fun <- CUDA.getFun mdl name + reg <- CUDA.requires fun CUDA.NumRegs + let + blockSize = 256 + sharedMem = 0 + maxBlocks = maxResidentBlocks prp blockSize reg sharedMem + numBlocks n = maxBlocks `min` ((n + blockSize - 1) `quot` blockSize) + -- + return $ Kernel fun sharedMem blockSize numBlocks name + +data Kernel = Kernel { + kernelFun :: {-# UNPACK #-} !CUDA.Fun + , kernelSharedMemBytes :: {-# UNPACK #-} !Int + , kernelThreadBlockSize :: {-# UNPACK #-} !Int + , kernelThreadBlocks :: (Int -> Int) + , kernelName :: String + } + +modules :: NumericR (Complex e) -> MVar (IntMap (Lifetime (CUDA.Module, Kernel, Kernel))) +modules NumericRcomplex32 = modules_f32 +modules NumericRcomplex64 = modules_f64 + +{-# NOINLINE modules_f32 #-} +modules_f32 :: MVar (IntMap (Lifetime (CUDA.Module, Kernel, Kernel))) +modules_f32 = unsafePerformIO $ newMVar IM.empty + +{-# NOINLINE modules_f64 #-} +modules_f64 :: MVar (IntMap (Lifetime (CUDA.Module, Kernel, Kernel))) +modules_f64 = unsafePerformIO $ newMVar IM.empty + +ptx_twine_f32 :: ByteString +ptx_twine_f32 = $(makeRelativeToProject "cubits/twine_f32.ptx" >>= embedFile) + +ptx_twine_f64 :: ByteString +ptx_twine_f64 = $(makeRelativeToProject "cubits/twine_f64.ptx" >>= embedFile) + diff --git a/Data/Array/Accelerate/Numeric/LinearAlgebra/Type.hs b/Data/Array/Accelerate/Numeric/LinearAlgebra/Type.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/LinearAlgebra/Type.hs @@ -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 +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.LinearAlgebra.Type + where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Data.Complex as A + +import qualified Prelude as P + + +-- For explicit dictionary reification, to recover the type the operation should +-- be performed at. +-- +data NumericR a where + NumericRfloat32 :: NumericR Float + NumericRfloat64 :: NumericR Double + NumericRcomplex32 :: NumericR (Complex Float) + NumericRcomplex64 :: NumericR (Complex Double) + +class (Elt a, Num a) => Numeric a where + numericR :: NumericR a + +instance Numeric Float where + numericR = NumericRfloat32 + +instance Numeric Double where + numericR = NumericRfloat64 + +instance Numeric (Complex Float) where + numericR = NumericRcomplex32 + +instance Numeric (Complex Double) where + numericR = NumericRcomplex64 + +-- class Numeric a => RealNumeric a +-- +-- instance RealNumeric Float +-- instance RealNumeric Double + +type family NumericBaseT t where + NumericBaseT Float = Float + NumericBaseT Double = Double + NumericBaseT (Complex Float) = Float + NumericBaseT (Complex Double) = Double + + +-- | Matrices as dense two-dimensional arrays in row-major ordering +-- +type Matrix e = Array DIM2 e + +-- | Orientation of the underlying data. +-- +-- Accelerate arrays are naturally stored in row-major format. +-- +data Orientation + = R -- ^ row major + | C -- ^ column major + deriving (P.Eq, P.Show) + +-- | Many operations allow you to implicitly transpose the arguments. For +-- a given input matrix @mat@ with dimensions @Z :. m :. n@ (that is; @m@ rows +-- and @n@ columns): +-- +data Transpose + -- | Leave the matrix as is. + = N + + -- | Treat the matrix as implicitly transposed, with dimensions @Z :. n :. m@. + -- Entry @Z :. j :. i@ is treated as actually being entry @Z :. i :. j@. + | T + + -- | Implicitly transpose and conjugate the input matrix. For complex-valued + -- matrices a given element @mat ! Z:.j:.i == x :+ y@ will be treated as + -- actually being @mat ! Z:.i:.j == x :+ (-y)@. + | H + deriving (P.Eq, P.Show) + diff --git a/Data/Array/Accelerate/Numeric/Sum.hs b/Data/Array/Accelerate/Numeric/Sum.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/Sum.hs @@ -0,0 +1,330 @@ +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE RebindableSyntax #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE TypeOperators #-} +{-# LANGUAGE ViewPatterns #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.Sum +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- +-- Functions for summing floating point numbers more accurately than the +-- straightforward 'Data.Array.Accelerate.sum' operation. +-- +-- In the worst case, the 'Data.Array.Accelerate.sum' function accumulates error +-- at a rate proportional to the number of values being summed. The algorithms +-- in this module implement different methods of /compensated summation/, which +-- reduce the accumulation of numeric error so that it grows much more slowly +-- than the number of inputs (e.g. logarithmically), or remains constant. +-- + +-- TLM: The standard formulation of the algorithms implemented here are not +-- associative; e.g. they would have a type (KBN a -> a -> KBN a). I've +-- done what seems like the sensible conversion, but somebody versed in numeric +-- analysis should probably look... +-- +-- See also: +-- + +module Data.Array.Accelerate.Numeric.Sum ( + + -- * Summation type class + Summation(..), + sum, + + -- * Kahan-Babuška-Neumaier summation + KBN(..), + kbn, + + -- * Order-2 Kahan-Babuška summation + KB2(..), + kb2, + + -- * Kahan summation + Kahan(..), + kahan, + +) where + +import Data.Array.Accelerate as A hiding ( sum ) +import Data.Array.Accelerate.Type as A +import Data.Array.Accelerate.Smart as A ( Exp(..), PreExp(..) ) +import Data.Array.Accelerate.Product as A +import Data.Array.Accelerate.Array.Sugar as A +import Data.Array.Accelerate.Numeric.Sum.Arithmetic as A + +import Data.Proxy +import Data.Typeable +import Prelude ( Show, fromInteger ) + + +-- | Sum an array using a particular compensation scheme. +-- +-- >>> let xs = [1.0, 1.0e100, 1.0, -1.0e100] :: [Double] +-- >>> Prelude.sum xs +-- 0.0 +-- +-- >>> let ys = fromList (Z:.4) [1.0, 1.0e100, 1.0, -1.0e100] :: Vector Double +-- >>> sum kbn (use ys) +-- Scalar Z [2.0] +-- +sum :: (Summation s a, Shape sh) => Proxy s -> Acc (Array (sh:.Int) a) -> Acc (Array sh a) +sum p = A.map (from p) + . A.fold add zero + . A.map (into p) + + +-- | A class for the summation of floating-point numbers +-- +class (Elt a, Elt (s a)) => Summation s a where + -- | Add a value to the sum + add :: Exp (s a) -> Exp (s a) -> Exp (s a) + + -- | The identity of the summation + zero :: Exp (s a) + + -- | Insert a value into the summation + into :: Proxy s -> Exp a -> Exp (s a) + + -- | Summarise the result of summation + from :: Proxy s -> Exp (s a) -> Exp a + + +-- | Kahan-Babuška-Neumaier summation. This is a little more computationally +-- costly than plain Kahan summation, but is /always/ at least as accurate. +-- +data KBN a = KBN a a + deriving (Show, Typeable) + +-- | Return the result of a Kahan-Babuška-Neumaier sum. +-- +kbn :: Proxy KBN +kbn = Proxy + +kbnAdd :: (Num a, Ord a, IsFloating a) => Exp (KBN a) -> Exp (KBN a) -> Exp (KBN a) +kbnAdd (unlift -> KBN s1 c1) (unlift -> KBN s2 c2) = lift (KBN s' c') + where + s' = s1 `fadd` s2 + c' = c1 `fadd` c2 `fadd` if abs s1 >= abs s2 + then (s1 `fsub` s') `fadd` s2 + else (s2 `fsub` s') `fadd` s1 + +-- instance (Num a, Ord a) => Summation KBN a where +-- zero = lift $ KBN (0::Exp a) (0::Exp a) +-- add = kbnAdd +-- into _ x = lift (KBN x 0) +-- from _ x = let KBN s c = unlift x in s + c + +instance Summation KBN Float where + zero = constant (KBN 0 0) + add = kbnAdd + into _ x = lift (KBN x 0) + from _ x = let KBN s c = unlift x in s + c + +instance Summation KBN Double where + zero = constant (KBN 0 0) + add = kbnAdd + into _ x = lift (KBN x 0) + from _ x = let KBN s c = unlift x in s + c + +instance Summation KBN CFloat where + zero = constant (KBN 0 0) + add = kbnAdd + into _ x = lift (KBN x 0) + from _ x = let KBN s c = unlift x in s + c + +instance Summation KBN CDouble where + zero = constant (KBN 0 0) + add = kbnAdd + into _ x = lift (KBN x 0) + from _ x = let KBN s c = unlift x in s + c + +type instance EltRepr (KBN a) = (((), EltRepr a), EltRepr a) + +instance Elt a => Elt (KBN a) where + eltType _ = UnitTuple `PairTuple` eltType (undefined::a) + `PairTuple` eltType (undefined::a) + toElt (((),a),b) = KBN (toElt a) (toElt b) + fromElt (KBN a b) = (((), fromElt a), fromElt b) + +instance Elt a => IsProduct Elt (KBN a) where + type ProdRepr (KBN a) = (((), a), a) + toProd _ (((),a),b) = KBN a b + fromProd _ (KBN a b) = (((),a),b) + prod _ _ = ProdRsnoc $ ProdRsnoc ProdRunit + +instance (Lift Exp a, Elt (Plain a)) => Lift Exp (KBN a) where + type Plain (KBN a) = KBN (Plain a) + lift (KBN a b) = Exp $ Tuple $ NilTup `SnocTup` lift a + `SnocTup` lift b + +instance Elt a => Unlift Exp (KBN (Exp a)) where + unlift t = KBN (Exp $ SuccTupIdx ZeroTupIdx `Prj` t) + (Exp $ ZeroTupIdx `Prj` t) + + +-- | Second-order Kahan-Babuška summation. This is more computationally costly +-- than Kahan-Babuška-Neumaier summation. Its advantage is that it can lose less +-- precision (in admittedly obscure cases). +-- +-- This method compensates for error in both the sum and the first-order +-- compensation term, hence the use of \"second order\" in the name. +-- +data KB2 a = KB2 a a a + deriving (Show, Typeable) + +-- | Return the result of a second-order Kahan-Babuška sum. +-- +kb2 :: Proxy KB2 +kb2 = Proxy + +kb2Add :: (Num a, Ord a, IsFloating a) => Exp (KB2 a) -> Exp (KB2 a) -> Exp (KB2 a) +kb2Add (unlift -> KB2 s1 c1 cc1) (unlift -> KB2 s2 c2 cc2) = lift (KB2 sum' c' cc') + where + sum' = s1 `fadd` s2 + c' = t `fadd` k + cc' = cc1 `fadd` cc2 `fadd` if abs t >= abs k + then (t `fsub` c') `fadd` k + else (k `fsub` c') `fadd` t + t = c1 `fadd` c2 + k = if abs s1 >= abs s2 + then (s1 `fsub` sum') `fadd` s2 + else (s2 `fsub` sum') `fadd` s1 + +-- instance (Num a, Ord a) => Summation KB2 a where +-- zero = lift $ KB2 (0::Exp a) (0::Exp a) (0::Exp a) +-- add = kb2Add +-- into _ x = lift (KB2 x 0 0) +-- from _ x = let KB2 s c cc = unlift x in s + c + cc + +instance Summation KB2 Float where + zero = constant (KB2 0 0 0) + add = kb2Add + into _ x = lift (KB2 x 0 0) + from _ x = let KB2 s c cc = unlift x in s + c + cc + +instance Summation KB2 Double where + zero = constant (KB2 0 0 0) + add = kb2Add + into _ x = lift (KB2 x 0 0) + from _ x = let KB2 s c cc = unlift x in s + c + cc + +instance Summation KB2 CFloat where + zero = constant (KB2 0 0 0) + add = kb2Add + into _ x = lift (KB2 x 0 0) + from _ x = let KB2 s c cc = unlift x in s + c + cc + +instance Summation KB2 CDouble where + zero = constant (KB2 0 0 0) + add = kb2Add + into _ x = lift (KB2 x 0 0) + from _ x = let KB2 s c cc = unlift x in s + c + cc + +type instance EltRepr (KB2 a) = ((((), EltRepr a), EltRepr a), EltRepr a) + +instance Elt a => Elt (KB2 a) where + eltType _ = UnitTuple `PairTuple` eltType (undefined::a) + `PairTuple` eltType (undefined::a) + `PairTuple` eltType (undefined::a) + toElt ((((),a),b),c) = KB2 (toElt a) (toElt b) (toElt c) + fromElt (KB2 a b c) = ((((), fromElt a), fromElt b), fromElt c) + +instance Elt a => IsProduct Elt (KB2 a) where + type ProdRepr (KB2 a) = ((((), a), a), a) + toProd _ ((((),a),b),c) = KB2 a b c + fromProd _ (KB2 a b c) = ((((),a),b),c) + prod _ _ = ProdRsnoc $ ProdRsnoc $ ProdRsnoc ProdRunit + +instance (Lift Exp a, Elt (Plain a)) => Lift Exp (KB2 a) where + type Plain (KB2 a) = KB2 (Plain a) + lift (KB2 a b c) = Exp $ Tuple $ NilTup `SnocTup` lift a + `SnocTup` lift b + `SnocTup` lift c + +instance Elt a => Unlift Exp (KB2 (Exp a)) where + unlift t = KB2 (Exp $ SuccTupIdx (SuccTupIdx ZeroTupIdx) `Prj` t) + (Exp $ SuccTupIdx ZeroTupIdx `Prj` t) + (Exp $ ZeroTupIdx `Prj` t) + + +-- | Kahan summation. This is the least accurate of the compensated summation +-- methods. This summation method is included only for completeness. +-- +data Kahan a = Kahan a a + deriving (Show, Typeable) + +-- | Return the result of a Kahan sum. +-- +kahan :: Proxy Kahan +kahan = Proxy + +kahanAdd :: (Num a, IsFloating a) => Exp (Kahan a) -> Exp (Kahan a) -> Exp (Kahan a) +kahanAdd (unlift -> Kahan s1 c1 :: Kahan (Exp a)) (unlift -> Kahan s2 c2) = lift (Kahan s' c') + where + s' = s1 `fadd` y + c' = (s' `fsub` s1) `fsub` y + y = s2 `fsub` c1 `fsub` c2 + +-- instance (Num a, Ord a) => Summation Kahan a where +-- zero = lift $ Kahan (0::Exp a) (0::Exp a) +-- add = kahanAdd +-- into _ x = lift (Kahan x 0) +-- from _ x = let Kahan s _ = unlift x in s + +instance Summation Kahan Float where + zero = constant (Kahan 0 0) + add = kahanAdd + into _ x = lift (Kahan x 0) + from _ x = let Kahan s _ = unlift x in s + +instance Summation Kahan Double where + zero = constant (Kahan 0 0) + add = kahanAdd + into _ x = lift (Kahan x 0) + from _ x = let Kahan s _ = unlift x in s + +instance Summation Kahan CFloat where + zero = constant (Kahan 0 0) + add = kahanAdd + into _ x = lift (Kahan x 0) + from _ x = let Kahan s _ = unlift x in s + +instance Summation Kahan CDouble where + zero = constant (Kahan 0 0) + add = kahanAdd + into _ x = lift (Kahan x 0) + from _ x = let Kahan s _ = unlift x in s + +type instance EltRepr (Kahan a) = (((), EltRepr a), EltRepr a) + +instance Elt a => Elt (Kahan a) where + eltType _ = UnitTuple `PairTuple` eltType (undefined::a) + `PairTuple` eltType (undefined::a) + toElt (((),a),b) = Kahan (toElt a) (toElt b) + fromElt (Kahan a b) = (((), fromElt a), fromElt b) + +instance Elt a => IsProduct Elt (Kahan a) where + type ProdRepr (Kahan a) = (((), a), a) + toProd _ (((),a),b) = Kahan a b + fromProd _ (Kahan a b) = (((),a),b) + prod _ _ = ProdRsnoc $ ProdRsnoc ProdRunit + +instance (Lift Exp a, Elt (Plain a)) => Lift Exp (Kahan a) where + type Plain (Kahan a) = Kahan (Plain a) + lift (Kahan a b) = Exp $ Tuple $ NilTup `SnocTup` lift a + `SnocTup` lift b + +instance Elt a => Unlift Exp (Kahan (Exp a)) where + unlift t = Kahan (Exp $ SuccTupIdx ZeroTupIdx `Prj` t) + (Exp $ ZeroTupIdx `Prj` t) + diff --git a/Data/Array/Accelerate/Numeric/Sum/Arithmetic.hs b/Data/Array/Accelerate/Numeric/Sum/Arithmetic.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/Sum/Arithmetic.hs @@ -0,0 +1,37 @@ +{-# LANGUAGE ConstraintKinds #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE NoImplicitPrelude #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.Sum.Arithmetic +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.Sum.Arithmetic ( + + fadd, fsub, fmul, + +) where + +import Data.Array.Accelerate + +import qualified Data.Array.Accelerate.Numeric.Sum.LLVM.Native as Native +import qualified Data.Array.Accelerate.Numeric.Sum.LLVM.PTX as PTX + + +infixl 6 `fadd` +fadd :: (Num a, IsFloating a) => Exp a -> Exp a -> Exp a +fadd = Native.fadd $ PTX.fadd (+) + +infixl 6 `fsub` +fsub :: (Num a, IsFloating a) => Exp a -> Exp a -> Exp a +fsub = Native.fsub $ PTX.fsub (-) + +infixl 7 `fmul` +fmul :: (Num a, IsFloating a) => Exp a -> Exp a -> Exp a +fmul = Native.fmul $ PTX.fmul (*) + diff --git a/Data/Array/Accelerate/Numeric/Sum/LLVM/Native.hs b/Data/Array/Accelerate/Numeric/Sum/LLVM/Native.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/Sum/LLVM/Native.hs @@ -0,0 +1,58 @@ +{-# LANGUAGE CPP #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.Sum.LLVM.Native +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.Sum.LLVM.Native ( + + fadd, fsub, fmul, + +) where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Type + +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND +import Data.Array.Accelerate.LLVM.CodeGen.Sugar +import Data.Array.Accelerate.LLVM.Native.Foreign as A +import qualified Data.Array.Accelerate.Numeric.Sum.LLVM.Prim as Prim +#endif + +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND +wrap2 :: (Elt a, Elt b, Elt c) + => String -- name of the operation + -> IRFun1 Native () ((a, b) -> c) -- foreign implementation + -> (Exp a -> Exp b -> Exp c) -- fallback implementation + -> Exp a + -> Exp b + -> Exp c +wrap2 str f g = A.curry (foreignExp (ForeignExp str f) (A.uncurry g)) +#endif + +fadd :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND +fadd = wrap2 "fadd" (Prim.fadd floatingType) +#else +fadd = id +#endif + +fsub :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND +fsub = wrap2 "fsub" (Prim.fsub floatingType) +#else +fsub = id +#endif + +fmul :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a +#ifdef ACCELERATE_LLVM_NATIVE_BACKEND +fmul = wrap2 "fmul" (Prim.fmul floatingType) +#else +fmul = id +#endif + diff --git a/Data/Array/Accelerate/Numeric/Sum/LLVM/PTX.hs b/Data/Array/Accelerate/Numeric/Sum/LLVM/PTX.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/Sum/LLVM/PTX.hs @@ -0,0 +1,58 @@ +{-# LANGUAGE CPP #-} +-- | +-- Module : Data.Array.Accelerate.Numeric.Sum.LLVM.PTX +-- Copyright : [2017] Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Trevor L. McDonell +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.Sum.LLVM.PTX ( + + fadd, fsub, fmul, + +) where + +import Data.Array.Accelerate as A +import Data.Array.Accelerate.Type + +#ifdef ACCELERATE_LLVM_PTX_BACKEND +import Data.Array.Accelerate.LLVM.CodeGen.Sugar +import Data.Array.Accelerate.LLVM.PTX.Foreign as A +import qualified Data.Array.Accelerate.Numeric.Sum.LLVM.Prim as Prim +#endif + +#ifdef ACCELERATE_LLVM_PTX_BACKEND +wrap2 :: (Elt a, Elt b, Elt c) + => String -- name of the operation + -> IRFun1 PTX () ((a, b) -> c) -- foreign implementation + -> (Exp a -> Exp b -> Exp c) -- fallback implementation + -> Exp a + -> Exp b + -> Exp c +wrap2 str f g = A.curry (foreignExp (ForeignExp str f) (A.uncurry g)) +#endif + +fadd :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a +#ifdef ACCELERATE_LLVM_PTX_BACKEND +fadd = wrap2 "fadd" (Prim.fadd floatingType) +#else +fadd = id +#endif + +fsub :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a +#ifdef ACCELERATE_LLVM_PTX_BACKEND +fsub = wrap2 "fsub" (Prim.fsub floatingType) +#else +fsub = id +#endif + +fmul :: (IsFloating a, Elt a) => (Exp a -> Exp a -> Exp a) -> Exp a -> Exp a -> Exp a +#ifdef ACCELERATE_LLVM_PTX_BACKEND +fmul = wrap2 "fmul" (Prim.fmul floatingType) +#else +fmul = id +#endif + diff --git a/Data/Array/Accelerate/Numeric/Sum/LLVM/Prim.hs b/Data/Array/Accelerate/Numeric/Sum/LLVM/Prim.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Numeric/Sum/LLVM/Prim.hs @@ -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 +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Numeric.Sum.LLVM.Prim ( + + fadd, fsub, fmul, + +) where + +import Data.Array.Accelerate.Type +import Data.Array.Accelerate.Error + +import Data.Array.Accelerate.LLVM.CodeGen.Downcast ( downcast ) +import Data.Array.Accelerate.LLVM.CodeGen.IR ( IR(..), Operands(..), IROP(..) ) +import Data.Array.Accelerate.LLVM.CodeGen.Monad ( CodeGen, freshName, instr_ ) +import Data.Array.Accelerate.LLVM.CodeGen.Sugar ( IROpenFun1(..) ) +import qualified Data.Array.Accelerate.LLVM.CodeGen.Arithmetic as A +import qualified LLVM.AST.Type.Name as A +import qualified LLVM.AST.Type.Operand as A +import qualified LLVM.AST.Type.Representation as A + +import LLVM.AST.Instruction +import LLVM.AST.Name +import LLVM.AST.Operand +import LLVM.AST.Type + + +-- | As (+), but don't allow potentially unsafe floating-point optimisations. +-- +fadd :: FloatingType a -> IROpenFun1 arch env aenv ((a,a) -> a) +fadd t = IRFun1 $ A.uncurry (binop FAdd t) + +-- | As (-), but don't allow potentially unsafe floating-point optimisations. +-- +fsub :: FloatingType a -> IROpenFun1 arch env aenv ((a,a) -> a) +fsub t = IRFun1 $ A.uncurry (binop FSub t) + +-- | As (*), but don't allow potentially unsafe floating-point optimisations. +-- +fmul :: FloatingType a -> IROpenFun1 arch env aenv ((a,a) -> a) +fmul t = IRFun1 $ A.uncurry (binop FMul t) + +binop :: (FastMathFlags -> Operand -> Operand -> InstructionMetadata -> Instruction) -> FloatingType a -> IR a -> IR a -> CodeGen (IR a) +binop f t (op t -> x) (op t -> y) = do + r <- instr (downcast t) (f fmf (downcast x) (downcast y) md) + return (upcast t r) + + +-- Prim +-- ---- + +md :: InstructionMetadata +md = [] + +fmf :: FastMathFlags +fmf = NoFastMathFlags + +fresh :: CodeGen Name +fresh = downcast <$> freshName + +instr :: Type -> Instruction -> CodeGen Operand +instr ty ins = do + name <- fresh + instr_ (name := ins) + return (LocalReference ty name) + +upcast :: FloatingType t -> Operand -> IR t +upcast TypeFloat{} (LocalReference (FloatingPointType FloatFP) (UnName x)) = IR $ OP_Float (A.LocalReference A.type' (A.UnName x)) +upcast TypeDouble{} (LocalReference (FloatingPointType DoubleFP) (UnName x)) = IR $ OP_Double (A.LocalReference A.type' (A.UnName x)) +upcast TypeCFloat{} (LocalReference (FloatingPointType FloatFP) (UnName x)) = IR $ OP_CFloat (A.LocalReference A.type' (A.UnName x)) +upcast TypeCDouble{} (LocalReference (FloatingPointType DoubleFP) (UnName x)) = IR $ OP_CDouble (A.LocalReference A.type' (A.UnName x)) +upcast _ _ = $internalError "upcast" "expected local reference" + diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,31 @@ +# Numeric linear algebra in Accelerate + +[![Build Status](https://travis-ci.org/tmcdonell/accelerate-blas.svg?branch=master)](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 + diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/accelerate-blas.cabal b/accelerate-blas.cabal new file mode 100644 --- /dev/null +++ b/accelerate-blas.cabal @@ -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: + + +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 diff --git a/bench/Accelerate.hs b/bench/Accelerate.hs new file mode 100644 --- /dev/null +++ b/bench/Accelerate.hs @@ -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) + diff --git a/bench/Extra.hs b/bench/Extra.hs new file mode 100644 --- /dev/null +++ b/bench/Extra.hs @@ -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) + diff --git a/bench/HMatrix.hs b/bench/HMatrix.hs new file mode 100644 --- /dev/null +++ b/bench/HMatrix.hs @@ -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)) + ] + diff --git a/bench/Main.hs b/bench/Main.hs new file mode 100644 --- /dev/null +++ b/bench/Main.hs @@ -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 + ] + diff --git a/cbits/twine_f32.c b/cbits/twine_f32.c new file mode 100644 --- /dev/null +++ b/cbits/twine_f32.c @@ -0,0 +1,60 @@ +/* + * Module : Twine + * Copyright : [2016] Trevor L. McDonell + * License : BSD3 + * + * Maintainer : Trevor L. McDonell + * Stability : experimental + * Portability : non-portable (GHC extensions) + * + * Convert between Accelerate's Struct-of-Array representation of complex + * numbers and the Array-of-Struct representation used by BLAS. + */ + +#include +#include "HsFFI.h" + +#ifdef __cplusplus +extern "C" { +#endif + +void interleave_f32 +( + const StgInt start, + const StgInt end, + complex float * __restrict__ cplx, + const float * __restrict__ real, + const float * __restrict__ imag +) +{ + StgInt i; + for (i = start; i < end; ++i) { + const float re = real[i]; + const float im = imag[i]; + + cplx[i] = re + im * I; + } +} + +void deinterleave_f32 +( + const StgInt start, + const StgInt end, + float * __restrict__ real, + float * __restrict__ imag, + const complex float * __restrict__ cplx +) +{ + StgInt i; + for (i = start; i < end; ++i) { + const complex float c = cplx[i]; + + real[i] = crealf(c); + imag[i] = cimagf(c); + } +} + +#ifdef __cplusplus +} +#endif + diff --git a/cbits/twine_f64.c b/cbits/twine_f64.c new file mode 100644 --- /dev/null +++ b/cbits/twine_f64.c @@ -0,0 +1,60 @@ +/* + * Module : Twine + * Copyright : [2016] Trevor L. McDonell + * License : BSD3 + * + * Maintainer : Trevor L. McDonell + * Stability : experimental + * Portability : non-portable (GHC extensions) + * + * Convert between Accelerate's Struct-of-Array representation of complex + * numbers and the Array-of-Struct representation used by BLAS. + */ + +#include +#include "HsFFI.h" + +#ifdef __cplusplus +extern "C" { +#endif + +void interleave_f64 +( + const StgInt start, + const StgInt end, + complex double * __restrict__ cplx, + const double * __restrict__ real, + const double * __restrict__ imag +) +{ + StgInt i; + for (i = start; i < end; ++i) { + const double re = real[i]; + const double im = imag[i]; + + cplx[i] = re + im * I; + } +} + +void deinterleave_f64 +( + const StgInt start, + const StgInt end, + double * __restrict__ real, + double * __restrict__ imag, + const complex double * __restrict__ cplx +) +{ + StgInt i; + for (i = start; i < end; ++i) { + const complex double c = cplx[i]; + + real[i] = creal(c); + imag[i] = cimag(c); + } +} + +#ifdef __cplusplus +} +#endif + diff --git a/cubits/twine_f32.ptx b/cubits/twine_f32.ptx new file mode 100644 --- /dev/null +++ b/cubits/twine_f32.ptx @@ -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; +} + + diff --git a/cubits/twine_f64.ptx b/cubits/twine_f64.ptx new file mode 100644 --- /dev/null +++ b/cubits/twine_f64.ptx @@ -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; +} + + diff --git a/test/Backend.hs b/test/Backend.hs new file mode 100644 --- /dev/null +++ b/test/Backend.hs @@ -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) + diff --git a/test/Hedgehog/Gen/Array.hs b/test/Hedgehog/Gen/Array.hs new file mode 100644 --- /dev/null +++ b/test/Hedgehog/Gen/Array.hs @@ -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 + diff --git a/test/Hedgehog/Gen/Shape.hs b/test/Hedgehog/Gen/Shape.hs new file mode 100644 --- /dev/null +++ b/test/Hedgehog/Gen/Shape.hs @@ -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 + diff --git a/test/Level2.hs b/test/Level2.hs new file mode 100644 --- /dev/null +++ b/test/Level2.hs @@ -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 + + + diff --git a/test/Level3.hs b/test/Level3.hs new file mode 100644 --- /dev/null +++ b/test/Level3.hs @@ -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 + diff --git a/test/Main.hs b/test/Main.hs new file mode 100644 --- /dev/null +++ b/test/Main.hs @@ -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 + ] + diff --git a/test/Similar.hs b/test/Similar.hs new file mode 100644 --- /dev/null +++ b/test/Similar.hs @@ -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 +