diff --git a/Data/Array/Accelerate/AST.hs b/Data/Array/Accelerate/AST.hs --- a/Data/Array/Accelerate/AST.hs +++ b/Data/Array/Accelerate/AST.hs @@ -1,5 +1,5 @@ {-# LANGUAGE GADTs, EmptyDataDecls, FlexibleContexts, TypeFamilies, TypeOperators #-} -{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, TypeSynonymInstances #-} +{-# LANGUAGE CPP, FlexibleInstances, MultiParamTypeClasses, TypeSynonymInstances #-} -- Module : Data.Array.Accelerate.AST -- Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee @@ -84,7 +84,9 @@ import Data.Array.Accelerate.Array.Sugar import Data.Array.Accelerate.Tuple +#include "accelerate.h" + -- Typed de Bruijn indices -- ----------------------- @@ -111,7 +113,7 @@ prj ZeroIdx (Push _ v) = v prj (SuccIdx idx) (Push val _) = prj idx val prj _ _ = - error "Data.Array.Accelerate: prj: inconsistent valuation" + INTERNAL_ERROR(error) "prj" "inconsistent valuation" -- Array expressions diff --git a/Data/Array/Accelerate/Analysis/Shape.hs b/Data/Array/Accelerate/Analysis/Shape.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/Analysis/Shape.hs @@ -0,0 +1,53 @@ +{-# LANGUAGE ScopedTypeVariables, GADTs #-} +-- | +-- Module : Data.Array.Accelerate.Analysis.Shape +-- Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Manuel M T Chakravarty +-- Stability : experimental +-- Portability : non-portable (GHC extensions) +-- + +module Data.Array.Accelerate.Analysis.Shape (accDim, accDim2) + where + +import Data.Array.Accelerate.AST +import Data.Array.Accelerate.Type +import Data.Array.Accelerate.Array.Sugar + + +-- |Reify the dimensionality of the result type of an array computation +-- +accDim :: forall aenv dim e. OpenAcc aenv (Array dim e) -> Int +accDim (Let _ acc) = accDim acc +accDim (Let2 _ acc) = accDim acc +accDim (Avar _) = ndim (elemType (undefined::dim)) +accDim (Use (Array _ _)) = ndim (elemType (undefined::dim)) +accDim (Unit _) = 0 +accDim (Reshape _ _) = ndim (elemType (undefined::dim)) +accDim (Replicate _ _ _) = ndim (elemType (undefined::dim)) +accDim (Index _ _ _) = ndim (elemType (undefined::dim)) +accDim (Map _ acc) = accDim acc +accDim (ZipWith _ _ acc) = accDim acc +accDim (Fold _ _ _) = 0 +accDim (FoldSeg _ _ _ _) = 1 +accDim (Permute _ acc _ _) = accDim acc +accDim (Backpermute _ _ _) = ndim (elemType (undefined::dim)) +accDim (Stencil _ _ acc) = accDim acc +accDim (Stencil2 _ _ acc _ _) = accDim acc + +-- |Reify the dimensionality of the results of a computation that yields two +-- arrays +-- +accDim2 :: OpenAcc aenv (Array dim1 e1, Array dim2 e2) -> (Int,Int) +accDim2 (Scanl _ _ acc) = (accDim acc, 0) +accDim2 (Scanr _ _ acc) = (accDim acc, 0) + +-- Count the number of components to a tuple type +-- +ndim :: TupleType a -> Int +ndim UnitTuple = 0 +ndim (SingleTuple _) = 1 +ndim (PairTuple a b) = ndim a + ndim b + diff --git a/Data/Array/Accelerate/Analysis/Type.hs b/Data/Array/Accelerate/Analysis/Type.hs --- a/Data/Array/Accelerate/Analysis/Type.hs +++ b/Data/Array/Accelerate/Analysis/Type.hs @@ -17,7 +17,7 @@ module Data.Array.Accelerate.Analysis.Type ( -- * Query AST types - accType, accType2, accShapeType, accShapeType2, expType, sizeOf + accType, accType2, expType, sizeOf ) where @@ -62,35 +62,6 @@ -> (TupleType (ElemRepr e1), TupleType (ElemRepr e2)) accType2 (Scanl _ e acc) = (accType acc, expType e) accType2 (Scanr _ e acc) = (accType acc, expType e) - --- |Reify the shape type of the result of an array computation --- -accShapeType :: forall aenv dim e. - OpenAcc aenv (Array dim e) -> TupleType (ElemRepr dim) -accShapeType (Let _ acc) = accShapeType acc -accShapeType (Let2 _ acc) = accShapeType acc -accShapeType (Avar _) = elemType (undefined::dim) -accShapeType (Use (Array _ _)) = elemType (undefined::dim) -accShapeType (Unit _) = elemType (undefined::DIM0) -accShapeType (Reshape _ _) = elemType (undefined::dim) -accShapeType (Replicate _ _ _) = elemType (undefined::dim) -accShapeType (Index _ _ _) = elemType (undefined::dim) -accShapeType (Map _ acc) = accShapeType acc -accShapeType (ZipWith _ _ acc) = accShapeType acc -accShapeType (Fold _ _ _) = elemType (undefined::DIM0) -accShapeType (FoldSeg _ _ _ _) = elemType (undefined::DIM1) -accShapeType (Permute _ acc _ _) = accShapeType acc -accShapeType (Backpermute _ _ _) = elemType (undefined::dim) -accShapeType (Stencil _ _ acc) = accShapeType acc -accShapeType (Stencil2 _ _ acc _ _) = accShapeType acc - --- |Reify the shape type of the results of a computation that yields two arrays --- -accShapeType2 :: OpenAcc aenv (Array dim1 e1, Array dim2 e2) - -> (TupleType (ElemRepr dim1), TupleType (ElemRepr dim2)) -accShapeType2 (Scanl _ _ acc) = (accShapeType acc, elemType (undefined::())) -accShapeType2 (Scanr _ _ acc) = (accShapeType acc, elemType (undefined::())) - -- |Reify the result type of a scalar expression. -- diff --git a/Data/Array/Accelerate/Array/Delayed.hs b/Data/Array/Accelerate/Array/Delayed.hs --- a/Data/Array/Accelerate/Array/Delayed.hs +++ b/Data/Array/Accelerate/Array/Delayed.hs @@ -40,7 +40,7 @@ data Delayed (Array dim e) = (Ix dim, Elem e) => DelayedArray { shapeDA :: ElemRepr dim - , repfDA :: (ElemRepr dim -> ElemRepr e) + , repfDA :: ElemRepr dim -> ElemRepr e } delay arr@(Array sh _) = DelayedArray sh (fromElem . (arr!) . toElem) force (DelayedArray sh f) = newArray (toElem sh) (toElem . f . fromElem) diff --git a/Data/Array/Accelerate/Array/Representation.hs b/Data/Array/Accelerate/Array/Representation.hs --- a/Data/Array/Accelerate/Array/Representation.hs +++ b/Data/Array/Accelerate/Array/Representation.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE GADTs, TypeFamilies, FlexibleContexts, FlexibleInstances #-} +{-# LANGUAGE CPP, GADTs, TypeFamilies, FlexibleContexts, FlexibleInstances #-} -- |Embedded array processing language: array representation -- @@ -19,7 +19,9 @@ -- friends import Data.Array.Accelerate.Type +#include "accelerate.h" + -- |Index representation -- @@ -75,10 +77,8 @@ size (sh, sz) = size sh * sz (sh1, sz1) `intersect` (sh2, sz2) = (sh1 `intersect` sh2, sz1 `min` sz2) ignore = (ignore, -1) - index (sh, sz) (ix, i) - | i >= 0 && i < sz = index sh ix + size sh * i - | otherwise - = error "Data.Array.Accelerate.Array: index out of bounds" + index (sh, sz) (ix, i) = BOUNDS_CHECK(checkIndex) "index" i sz + $ index sh ix + size sh * i bound (sh, sz) (ix, i) bndy | i < 0 = case bndy of Clamp -> bound sh ix bndy `addDim` 0 diff --git a/Data/Array/Accelerate/CUDA.hs b/Data/Array/Accelerate/CUDA.hs --- a/Data/Array/Accelerate/CUDA.hs +++ b/Data/Array/Accelerate/CUDA.hs @@ -20,6 +20,7 @@ import Prelude hiding (catch) import Control.Exception +import System.IO.Unsafe import Foreign.CUDA.Driver.Error @@ -36,13 +37,15 @@ -- Accelerate: CUDA -- ~~~~~~~~~~~~~~~~ --- | Compiles and run a complete embedded array program using the CUDA backend +-- | Compile and run a complete embedded array program using the CUDA backend -- -run :: Arrays a => Sugar.Acc a -> IO a -run acc = - evalCUDA (execute (Sugar.convertAcc acc) >>= collect) - `catch` - \e -> INTERNAL_ERROR(error) "unhandled" (show (e :: CUDAException)) +{-# NOINLINE run #-} +run :: Arrays a => Sugar.Acc a -> a +run acc + = unsafePerformIO + $ evalCUDA (execute (Sugar.convertAcc acc) >>= collect) + `catch` + \e -> INTERNAL_ERROR(error) "unhandled" (show (e :: CUDAException)) execute :: Arrays a => Acc a -> CIO a execute acc = compileAcc acc >> executeAcc acc diff --git a/Data/Array/Accelerate/CUDA/Analysis/Device.hs b/Data/Array/Accelerate/CUDA/Analysis/Device.hs new file mode 100644 --- /dev/null +++ b/Data/Array/Accelerate/CUDA/Analysis/Device.hs @@ -0,0 +1,41 @@ +-- | +-- Module : Data.Array.Accelerate.CUDA.Analysis.Device +-- Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell +-- License : BSD3 +-- +-- Maintainer : Manuel M T Chakravarty +-- Stability : experimental +-- Portability : non-partable (GHC extensions) +-- + +module Data.Array.Accelerate.CUDA.Analysis.Device + where + +import Data.Ord +import Data.List +import Data.Function +import Foreign.CUDA.Driver.Device +import qualified Foreign.CUDA.Driver as CUDA + + +-- Select the best of the available CUDA capable devices. This prefers devices +-- with higher compute capability, followed by maximum throughput. This does not +-- take into account any other factors, such as whether the device is currently +-- in use by another process. +-- +-- Ignore the possibility of emulation-mode devices, as this has been deprecated +-- as of CUDA v3.0 (compute-capability == 9999.9999) +-- +selectBestDevice :: IO (Device, DeviceProperties) +selectBestDevice = do + dev <- mapM CUDA.device . enumFromTo 0 . subtract 1 . fromIntegral =<< CUDA.count + prop <- mapM CUDA.props dev + return . head . sortBy (cmp `on` snd) $ zip dev prop + where + compute = computeCapability + flops d = multiProcessorCount d * clockRate d * cores d + cores d | compute d < 2 = 8 + | otherwise = 32 + cmp x y | compute x == compute y = comparing flops x y + | otherwise = comparing compute x y + diff --git a/Data/Array/Accelerate/CUDA/Analysis/Hash.hs b/Data/Array/Accelerate/CUDA/Analysis/Hash.hs --- a/Data/Array/Accelerate/CUDA/Analysis/Hash.hs +++ b/Data/Array/Accelerate/CUDA/Analysis/Hash.hs @@ -22,6 +22,8 @@ import Data.Array.Accelerate.Pretty () import Data.Array.Accelerate.Analysis.Type import Data.Array.Accelerate.CUDA.CodeGen +import Data.Array.Accelerate.Array.Representation +import qualified Data.Array.Accelerate.Array.Sugar as Sugar #include "accelerate.h" @@ -30,8 +32,8 @@ -- Generate a unique key for each kernel computation (not extensively tested...) -- accToKey :: OpenAcc aenv a -> String -accToKey (Replicate _ e a) = chr 17 : showTy (accType a) ++ showExp e -accToKey (Index _ a e) = chr 29 : showTy (accType a) ++ showExp e +accToKey r@(Replicate s e a) = chr 17 : showTy (accType a) ++ showExp e ++ showSI s e a r +accToKey r@(Index s a e) = chr 29 : showTy (accType a) ++ showExp e ++ showSI s e r a accToKey (Map f a) = chr 41 : showTy (accType a) ++ showFun f accToKey (ZipWith f x y) = chr 53 : showTy (accType x) ++ showTy (accType y) ++ showFun f accToKey (Fold f e a) = chr 67 : showTy (accType a) ++ showFun f ++ showExp e @@ -54,6 +56,19 @@ showExp :: OpenExp env aenv a -> String showExp e = render . hcat . map pretty $ evalState (codeGenExp e) [] + +showSI :: SliceIndex (Sugar.ElemRepr slix) (Sugar.ElemRepr sl) co (Sugar.ElemRepr dim) + -> Exp aenv slix {- dummy -} + -> OpenAcc aenv (Sugar.Array sl e) {- dummy -} + -> OpenAcc aenv (Sugar.Array dim e) {- dummy -} + -> String +showSI sl _ _ _ = slice sl 0 + where + slice :: SliceIndex slix sl co dim -> Int -> String + slice (SliceNil) _ = [] + slice (SliceAll sliceIdx) n = '_' : slice sliceIdx n + slice (SliceFixed sliceIdx) n = show n ++ slice sliceIdx (n+1) + {- -- hash function from the dragon book pp437; assumes 7 bit characters and needs diff --git a/Data/Array/Accelerate/CUDA/Analysis/Launch.hs b/Data/Array/Accelerate/CUDA/Analysis/Launch.hs --- a/Data/Array/Accelerate/CUDA/Analysis/Launch.hs +++ b/Data/Array/Accelerate/CUDA/Analysis/Launch.hs @@ -36,19 +36,25 @@ -- TLM: this could probably be stored in the KernelEntry -- launchConfig :: OpenAcc aenv a -> Int -> CUDA.Fun -> CIO (Int, Int, Integer) -launchConfig acc@(Fold _ _ _) n _ = - return (256, gridSize undefined acc n 256, toInteger (sharedMem undefined acc 256)) - launchConfig acc n fn = do regs <- liftIO $ CUDA.requires fn CUDA.NumRegs stat <- liftIO $ CUDA.requires fn CUDA.SharedSizeBytes -- static memory only prop <- getM deviceProps let dyn = sharedMem prop acc - (cta, occ) = CUDA.optimalBlockSize prop (const regs) ((stat+) . dyn) + (cta, occ) = blockSize prop acc regs ((stat+) . dyn) mbk = CUDA.multiProcessorCount prop * CUDA.activeThreadBlocks occ return (cta, mbk `min` gridSize prop acc n cta, toInteger (dyn cta)) + + +-- | +-- Determine the optimal thread block size for a given array computation. Fold +-- requires blocks with a power-of-two number of threads. +-- +blockSize :: CUDA.DeviceProperties -> OpenAcc aenv a -> Int -> (Int -> Int) -> (Int, CUDA.Occupancy) +blockSize p (Fold _ _ _) r s = CUDA.optimalBlockSizeBy p CUDA.incPow2 (const r) s +blockSize p _ r s = CUDA.optimalBlockSizeBy p CUDA.incWarp (const r) s -- | diff --git a/Data/Array/Accelerate/CUDA/CodeGen.hs b/Data/Array/Accelerate/CUDA/CodeGen.hs --- a/Data/Array/Accelerate/CUDA/CodeGen.hs +++ b/Data/Array/Accelerate/CUDA/CodeGen.hs @@ -28,6 +28,7 @@ import Data.Array.Accelerate.Tuple import Data.Array.Accelerate.Pretty () import Data.Array.Accelerate.Analysis.Type +import Data.Array.Accelerate.Analysis.Shape import Data.Array.Accelerate.Array.Representation import qualified Data.Array.Accelerate.AST as AST import qualified Data.Array.Accelerate.Array.Sugar as Sugar @@ -74,18 +75,18 @@ codeGenAcc' (AST.Scanr f1 e1 _) = mkScanr (codeGenExpType e1) <$> codeGenExp e1 <*> codeGenFun f1 codeGenAcc' op@(AST.Map f1 a1) = mkMap (codeGenAccType op) (codeGenAccType a1) <$> codeGenFun f1 codeGenAcc' op@(AST.ZipWith f1 a1 a0) - = mkZipWith (codeGenAccType op) (codeGenShapeType op) - (codeGenAccType a1) (codeGenShapeType a1) - (codeGenAccType a0) (codeGenShapeType a0) + = mkZipWith (codeGenAccType op) (accDim op) + (codeGenAccType a1) (accDim a1) + (codeGenAccType a0) (accDim a0) <$> codeGenFun f1 codeGenAcc' op@(AST.Permute f1 _ f2 a1) - = mkPermute (codeGenAccType a1) (codeGenShapeType op) (codeGenShapeType a1) + = mkPermute (codeGenAccType a1) (accDim op) (accDim a1) <$> codeGenFun f1 <*> codeGenFun f2 codeGenAcc' op@(AST.Backpermute _ f1 a1) - = mkBackpermute (codeGenAccType a1) (codeGenShapeType op) (codeGenShapeType a1) + = mkBackpermute (codeGenAccType a1) (accDim op) (accDim a1) <$> codeGenFun f1 codeGenAcc' x = @@ -109,8 +110,14 @@ -- computations must be hoisted out of scalar expressions before code generation -- or execution: kernel functions can not invoke other array computations. -- --- TLM 2010-06-24: Shape for free array variables?? +-- TLM 2010-06-24: +-- Shape for free array variables?? -- +-- TLM 2010-09-03: +-- We push a conditional expression into each component of a tuple, in order +-- to support tuples as the result of a branch. However, I doubt nvcc is +-- clever enough to do CSE on this expression. +-- codeGenExp :: forall env aenv t. AST.OpenExp env aenv t -> CG [CExpr] codeGenExp (AST.Shape _) = return . unit $ CVar (internalIdent "shape") internalNode codeGenExp (AST.PrimConst c) = return . unit $ codeGenPrimConst c @@ -133,9 +140,7 @@ codeGenExp (AST.Cond p e1 e2) = do [a] <- codeGenExp p - [b] <- codeGenExp e1 - [c] <- codeGenExp e2 - return [CCond a (Just b) c internalNode] + zipWith (\b c -> CCond a (Just b) c internalNode) <$> codeGenExp e1 <*> codeGenExp e2 codeGenExp (AST.IndexScalar a1 e1) = do n <- length <$> get @@ -187,12 +192,12 @@ -> AST.Exp aenv slix -> CG CUTranslSkel codeGenIndex sl acc acc' slix = - return . mkIndex ty dimSl dimCo dimIn0 $ restrict sl (length dimCo-1,length dimSl-1) + return . mkIndex ty dimSl dimCo dimIn0 $ restrict sl (dimCo-1,dimSl-1) where ty = codeGenAccType acc - dimCo = codeGenExpType slix - dimSl = codeGenShapeType acc' - dimIn0 = codeGenShapeType acc + dimCo = length (codeGenExpType slix) + dimSl = accDim acc' + dimIn0 = accDim acc restrict :: SliceIndex slix sl co dim -> (Int,Int) -> [CExpr] restrict (SliceNil) _ = [] @@ -207,30 +212,27 @@ (Sugar.ElemRepr sl) co (Sugar.ElemRepr dim) - -> AST.Exp aenv slix + -> AST.Exp aenv slix -- dummy to fix type variables -> AST.OpenAcc aenv (Sugar.Array sl e) -> AST.OpenAcc aenv (Sugar.Array dim e) -> CG CUTranslSkel codeGenReplicate sl _slix acc acc' = - return . mkReplicate ty dimSl dimOut . post $ extend sl (length dimOut-1) + return . mkReplicate ty dimSl dimOut $ extend sl (dimOut-1) where ty = codeGenAccType acc - dimSl = codeGenShapeType acc - dimOut = codeGenShapeType acc' + dimSl = accDim acc + dimOut = accDim acc' extend :: SliceIndex slix sl co dim -> Int -> [CExpr] extend (SliceNil) _ = [] extend (SliceAll sliceIdx) n = mkPrj dimOut "dim" n : extend sliceIdx (n-1) extend (SliceFixed sliceIdx) n = extend sliceIdx (n-1) - post [] = [CConst (CIntConst (cInteger 0) internalNode)] - post xs = xs - -mkPrj :: [CType] -> String -> Int -> CExpr -mkPrj ty var c - | length ty <= 1 = CVar (internalIdent var) internalNode - | otherwise = CMember (CVar (internalIdent var) internalNode) (internalIdent ('a':show c)) False internalNode +mkPrj :: Int -> String -> Int -> CExpr +mkPrj ndim var c + | ndim <= 1 = CVar (internalIdent var) internalNode + | otherwise = CMember (CVar (internalIdent var) internalNode) (internalIdent ('a':show c)) False internalNode -- Types @@ -243,9 +245,6 @@ codeGenExpType :: AST.OpenExp aenv env t -> [CType] codeGenExpType = codeGenTupleType . expType - -codeGenShapeType :: AST.OpenAcc aenv (Sugar.Array dim e) -> [CType] -codeGenShapeType = codeGenTupleType . accShapeType -- Implementation diff --git a/Data/Array/Accelerate/CUDA/CodeGen/Skeleton.hs b/Data/Array/Accelerate/CUDA/CodeGen/Skeleton.hs --- a/Data/Array/Accelerate/CUDA/CodeGen/Skeleton.hs +++ b/Data/Array/Accelerate/CUDA/CodeGen/Skeleton.hs @@ -33,7 +33,8 @@ skel = "fold.inl" code = CTranslUnit ( mkTupleTypeAsc 2 ty ++ - [ mkIdentity identity + [ mkTuplePartition ty + , mkIdentity identity , mkApply 2 apply ]) (mkNodeInfo (initPos skel) (Name 0)) @@ -43,7 +44,8 @@ skel = "fold_segmented.inl" code = CTranslUnit ( mkTupleTypeAsc 2 ty ++ - [ mkTypedef "Int" False (head int) + [ mkTuplePartition ty + , mkTypedef "Int" False (head int) , mkIdentity identity , mkApply 2 apply ]) (mkNodeInfo (initPos skel) (Name 0)) @@ -64,8 +66,8 @@ (mkNodeInfo (initPos skel) (Name 0)) -mkZipWith :: [CType] -> [CType] -> [CType] -> [CType] -> [CType] -> [CType] -> [CExpr] -> CUTranslSkel -mkZipWith tyOut shOut tyIn1 shIn1 tyIn0 shIn0 apply = CUTranslSkel code skel +mkZipWith :: [CType] -> Int -> [CType] -> Int -> [CType] -> Int -> [CExpr] -> CUTranslSkel +mkZipWith tyOut dimOut tyIn1 dimIn1 tyIn0 dimIn0 apply = CUTranslSkel code skel where skel = "zipWith.inl" code = CTranslUnit @@ -73,9 +75,9 @@ mkTupleType (Just 1) tyIn1 ++ mkTupleType (Just 0) tyIn0 ++ [ mkApply 2 apply - , mkDim "DimOut" shOut - , mkDim "DimIn1" shIn1 - , mkDim "DimIn0" shIn0 ]) + , mkDim "DimOut" dimOut + , mkDim "DimIn1" dimIn1 + , mkDim "DimIn0" dimIn0 ]) (mkNodeInfo (initPos skel) (Name 0)) @@ -119,7 +121,7 @@ -- Permutation -------------------------------------------------------------------------------- -mkPermute :: [CType] -> [CType] -> [CType] -> [CExpr] -> [CExpr] -> CUTranslSkel +mkPermute :: [CType] -> Int -> Int -> [CExpr] -> [CExpr] -> CUTranslSkel mkPermute ty dimOut dimIn0 combinefn indexfn = CUTranslSkel code skel where skel = "permute.inl" @@ -127,19 +129,19 @@ ( mkTupleTypeAsc 2 ty ++ [ mkDim "DimOut" dimOut , mkDim "DimIn0" dimIn0 - , mkProject indexfn + , mkProject Forward indexfn , mkApply 2 combinefn ]) (mkNodeInfo (initPos skel) (Name 0)) -mkBackpermute :: [CType] -> [CType] -> [CType] -> [CExpr] -> CUTranslSkel -mkBackpermute ty dimOut dimIn0 index = CUTranslSkel code skel +mkBackpermute :: [CType] -> Int -> Int -> [CExpr] -> CUTranslSkel +mkBackpermute ty dimOut dimIn0 indexFn = CUTranslSkel code skel where skel = "backpermute.inl" code = CTranslUnit ( mkTupleTypeAsc 1 ty ++ [ mkDim "DimOut" dimOut , mkDim "DimIn0" dimIn0 - , mkProject index ]) + , mkProject Backward indexFn ]) (mkNodeInfo (initPos skel) (Name 0)) @@ -147,7 +149,7 @@ -- Multidimensional Index and Replicate -------------------------------------------------------------------------------- -mkIndex :: [CType] -> [CType] -> [CType] -> [CType] -> [CExpr] -> CUTranslSkel +mkIndex :: [CType] -> Int -> Int -> Int -> [CExpr] -> CUTranslSkel mkIndex ty dimSl dimCo dimIn0 slix = CUTranslSkel code skel where skel = "slice.inl" @@ -160,7 +162,7 @@ (mkNodeInfo (initPos skel) (Name 0)) -mkReplicate :: [CType] -> [CType] -> [CType] -> [CExpr] -> CUTranslSkel +mkReplicate :: [CType] -> Int -> Int -> [CExpr] -> CUTranslSkel mkReplicate ty dimSl dimOut slix = CUTranslSkel code skel where skel = "replicate.inl" diff --git a/Data/Array/Accelerate/CUDA/CodeGen/Tuple.hs b/Data/Array/Accelerate/CUDA/CodeGen/Tuple.hs --- a/Data/Array/Accelerate/CUDA/CodeGen/Tuple.hs +++ b/Data/Array/Accelerate/CUDA/CodeGen/Tuple.hs @@ -8,7 +8,7 @@ -- Portability : non-partable (GHC extensions) -- -module Data.Array.Accelerate.CUDA.CodeGen.Tuple (mkTupleType, mkTupleTypeAsc) +module Data.Array.Accelerate.CUDA.CodeGen.Tuple (mkTupleType, mkTupleTypeAsc, mkTuplePartition) where import Language.C @@ -84,4 +84,23 @@ | n <= 1 = [CBlockStmt (CExpr (Just (CAssign CAssignOp (CIndex (CVar (internalIdent "d_out") internalNode) (CVar (internalIdent "idx") internalNode) internalNode) (CVar (internalIdent "val") internalNode) internalNode)) internalNode)] | otherwise = reverse . take n . flip map (enumFrom 0 :: [Int]) $ \v -> CBlockStmt (CExpr (Just (CAssign CAssignOp (CIndex (CMember (CVar (internalIdent "d_out") internalNode) (internalIdent ('a':show v)) False internalNode) (CVar (internalIdent "idx") internalNode) internalNode) (CMember (CVar (internalIdent "val") internalNode) (internalIdent ('a':show v)) False internalNode) internalNode)) internalNode) + + +mkTuplePartition :: [CType] -> CExtDecl +mkTuplePartition ty = + CFDefExt + (CFunDef + [CStorageSpec (CStatic internalNode),CTypeQual (CInlineQual internalNode),CTypeQual (CAttrQual (CAttr (internalIdent "device") [] internalNode)),CTypeSpec (CTypeDef (internalIdent "ArrOut") internalNode)] + (CDeclr (Just (internalIdent "partition")) [CFunDeclr (Right ([CDecl [CTypeQual (CConstQual internalNode),CTypeSpec (CVoidType internalNode)] [(Just (CDeclr (Just (internalIdent "s_data")) [CPtrDeclr [] internalNode] Nothing [] internalNode),Nothing,Nothing)] internalNode,CDecl [CTypeQual (CConstQual internalNode),CTypeSpec (CIntType internalNode)] [(Just (CDeclr (Just (internalIdent "n")) [] Nothing [] internalNode),Nothing,Nothing)] internalNode],False)) [] internalNode] Nothing [] internalNode) + [] + (CCompound [] (stmts ++ [CBlockDecl (CDecl [CTypeSpec (CTypeDef (internalIdent "ArrOut") internalNode)] [(Just (CDeclr (Just (internalIdent "r")) [] Nothing [] internalNode),Just initp,Nothing)] internalNode) ,CBlockStmt (CReturn (Just (CVar (internalIdent "r") internalNode)) internalNode)]) internalNode) + internalNode) + where + n = length ty + var s = CVar (internalIdent s) internalNode + names = map (('p':) . show) [n-1,n-2..0] + initp = mkInitList (map var names) + stmts = zipWith (\l r -> CBlockDecl (CDecl (map CTypeSpec l) r internalNode)) ty + . zipWith3 (\p t s -> [(Just (CDeclr (Just (internalIdent p)) [CPtrDeclr [] internalNode] Nothing [] internalNode),Just (CInitExpr (CCast (CDecl (map CTypeSpec t) [(Just (CDeclr Nothing [CPtrDeclr [] internalNode] Nothing [] internalNode),Nothing,Nothing)] internalNode) s internalNode) internalNode),Nothing)]) names ty + $ var "s_data" : map (\v -> (CUnary CAdrOp (CIndex (var v) (CVar (internalIdent "n") internalNode) internalNode) internalNode)) names diff --git a/Data/Array/Accelerate/CUDA/CodeGen/Util.hs b/Data/Array/Accelerate/CUDA/CodeGen/Util.hs --- a/Data/Array/Accelerate/CUDA/CodeGen/Util.hs +++ b/Data/Array/Accelerate/CUDA/CodeGen/Util.hs @@ -14,6 +14,7 @@ import Language.C import Data.Array.Accelerate.CUDA.CodeGen.Data +data Direction = Forward | Backward -------------------------------------------------------------------------------- -- Common device functions @@ -27,8 +28,9 @@ = mkDeviceFun "apply" (typename "TyOut") $ map (\n -> (typename ("TyIn"++show n), 'x':show n)) [argc-1,argc-2..0] -mkProject :: [CExpr] -> CExtDecl -mkProject = mkDeviceFun "project" (typename "DimOut") [(typename "DimIn0","x0")] +mkProject :: Direction -> [CExpr] -> CExtDecl +mkProject Forward = mkDeviceFun "project" (typename "DimOut") [(typename "DimIn0","x0")] +mkProject Backward = mkDeviceFun "project" (typename "DimIn0") [(typename "DimOut","x0")] mkSliceIndex :: [CExpr] -> CExtDecl mkSliceIndex = @@ -50,9 +52,9 @@ fromBool True = CConst $ CIntConst (cInteger 1) internalNode fromBool False = CConst $ CIntConst (cInteger 0) internalNode -mkDim :: String -> [CType] -> CExtDecl -mkDim name ty = - mkTypedef name False [CTypeDef (internalIdent ("DIM" ++ show (length ty))) internalNode] +mkDim :: String -> Int -> CExtDecl +mkDim name n = + mkTypedef name False [CTypeDef (internalIdent ("DIM" ++ show n)) internalNode] mkTypedef :: String -> Bool -> CType -> CExtDecl mkTypedef var ptr ty = @@ -62,6 +64,7 @@ internalNode mkInitList :: [CExpr] -> CInit +mkInitList [] = CInitExpr (CConst (CIntConst (cInteger 0) internalNode)) internalNode mkInitList [x] = CInitExpr x internalNode mkInitList xs = CInitList (map (\e -> ([],CInitExpr e internalNode)) xs) internalNode diff --git a/Data/Array/Accelerate/CUDA/Compile.hs b/Data/Array/Accelerate/CUDA/Compile.hs --- a/Data/Array/Accelerate/CUDA/Compile.hs +++ b/Data/Array/Accelerate/CUDA/Compile.hs @@ -149,7 +149,7 @@ compileFlags cufile = do arch <- computeCapability <$> getM deviceProps ddir <- liftIO getDataDir - return [ "-I.", "-I", ddir + return [ "-I.", "-I", ddir "cubits" , "-O2", "--compiler-options", "-fno-strict-aliasing" , "-arch=sm_" ++ show (round (arch * 10) :: Int) , "-DUNIX" @@ -174,9 +174,11 @@ outputName :: OpenAcc aenv a -> FilePath -> CIO FilePath outputName acc cufile = do n <- freshVar - x <- liftIO $ doesFileExist (base ++ show n <.> suffix) + x <- liftIO $ doesFileExist (filename n) if x then outputName acc cufile - else return (base ++ show n <.> suffix) + else return (filename n) where (base,suffix) = splitExtension cufile + filename n = base ++ pad (show n) <.> suffix + pad s = replicate (4-length s) '0' ++ s diff --git a/Data/Array/Accelerate/CUDA/Execute.hs b/Data/Array/Accelerate/CUDA/Execute.hs --- a/Data/Array/Accelerate/CUDA/Execute.hs +++ b/Data/Array/Accelerate/CUDA/Execute.hs @@ -109,13 +109,18 @@ (ax1,ax2) <- executeOpenAcc x env executeOpenAcc y (env `Push` ax1 `Push` ax2) -executeOpenAcc (Unit e) env = do +executeOpenAcc (Reshape e a) env = do + ix <- executeExp e env + (Array sh ad) <- executeOpenAcc a env + BOUNDS_CHECK(check) "reshape" "shape mismatch" (Sugar.size ix == size sh) + $ return (Array (Sugar.fromElem ix) ad) + +executeOpenAcc (Unit e) env = do v <- executeExp e env let ad = fst . runArrayData $ (,undefined) <$> do - arr <- newArrayData 1 - writeArrayData arr 0 (Sugar.fromElem v) + arr <- newArrayData 1024 -- FIXME: small arrays moved by the GC + writeArrayData arr 0 (Sugar.fromElem v) return arr - mallocArray ad 1 pokeArrayAsync ad 1 Nothing return (Array (Sugar.fromElem ()) ad) @@ -268,10 +273,59 @@ release f_arr return res +-- Unified dispatch handler for left/right scan. +-- +-- TLM 2010-07-02: +-- This is a little awkward. At its core we have an inclusive scan routine, +-- whereas accelerate actually returns an exclusive scan result. Multi-block +-- arrays require the inclusive scan when calculating the partial block sums, +-- which are then added to every element of an interval. Single-block arrays +-- on the other hand need to be "converted" to an exclusive result in the +-- second pass. +-- +-- This optimised could be implemented by enabling some some extra code to the +-- skeleton, and dispatching accordingly. +-- +-- TLM 2010-08-19: +-- On the other hand, the inclusive scan core provides a way in which we could +-- add support for non-identic seed elements. +-- +dispatch (Scanr f x ad) env mdl = dispatch (Scanl f x ad) env mdl +dispatch acc@(Scanl f x ad) env mdl = do + fscan <- liftIO $ CUDA.getFun mdl "inclusive_scan" + fadd <- liftIO $ CUDA.getFun mdl "exclusive_update" + (Array sh in0) <- executeOpenAcc ad env + (cta,grid,smem) <- launchConfig acc (size sh) fscan + let a_out@(Array _ out) = newArray (Sugar.toElem sh) + a_sum@(Array _ sum) = newArray () + a_bks@(Array _ bks) = newArray grid + n = size sh + interval = (n + grid - 1) `div` grid -dispatch acc@(Scanl _ _ _) env mdl = dispatchScan acc env mdl -dispatch acc@(Scanr _ _ _) env mdl = dispatchScan acc env mdl + unify :: Array dim e -> Array dim e -> CIO () + unify _ _ = return () + unify a_out a_bks -- TLM: *cough* + + mallocArray out n + mallocArray sum 1 + mallocArray bks grid + d_out <- devicePtrs out + d_in0 <- devicePtrs in0 + d_bks <- devicePtrs bks + d_sum <- devicePtrs sum + f_arr <- (++) <$> liftExp x env <*> liftFun f env + t_var <- bind mdl f_arr + + launch' (cta,grid,smem) fscan (d_out ++ d_in0 ++ d_bks ++ t_var ++ map CUDA.IArg [n,interval]) + launch' (cta,1,smem) fscan (d_bks ++ d_bks ++ d_sum ++ map CUDA.IArg [grid,interval]) + launch' (cta,grid,smem) fadd (d_out ++ d_bks ++ map CUDA.IArg [n,interval]) + + freeArray in0 + freeArray bks + release f_arr + return (a_out, a_sum) + dispatch acc@(Permute f1 df f2 ad) env mdl = do fn <- liftIO $ CUDA.getFun mdl "permute" (Array sh def) <- executeOpenAcc df env @@ -366,60 +420,6 @@ (unlines ["unsupported array primitive", render . nest 2 $ text (show x)]) --- Unified dispatch handler for left/right scan. --- --- TLM 2010-07-02: --- This is a little awkward. At its core we have an inclusive scan routine, --- whereas accelerate actually returns an exclusive scan result. Multi-block --- arrays require the inclusive scan when calculating the partial block sums, --- which are then added to every element of an interval. Single-block arrays --- on the other hand need to be "converted" to an exclusive result in the --- second pass. --- --- This optimised could be implemented by enabling some some extra code to the --- skeleton, and dispatching accordingly. --- -dispatchScan :: OpenAcc aenv a -> Val aenv -> CUDA.Module -> CIO a -dispatchScan (Scanr f x ad) env mdl = dispatchScan (Scanl f x ad) env mdl -dispatchScan acc@(Scanl f x ad) env mdl = do - fscan <- liftIO $ CUDA.getFun mdl "inclusive_scan" - fadd <- liftIO $ CUDA.getFun mdl "exclusive_update" - (Array sh in0) <- executeOpenAcc ad env - (cta,grid,smem) <- launchConfig acc (size sh) fscan - let a_out@(Array _ out) = newArray (Sugar.toElem sh) - a_sum@(Array _ sum) = newArray () - a_bks@(Array _ bks) = newArray grid - n = size sh - interval = (n + grid - 1) `div` grid - - unify :: Array dim e -> Array dim e -> CIO () - unify _ _ = return () - - unify a_out a_bks -- TLM: *cough* - - mallocArray out n - mallocArray sum 1 - mallocArray bks grid - d_out <- devicePtrs out - d_in0 <- devicePtrs in0 - d_bks <- devicePtrs bks - d_sum <- devicePtrs sum - f_arr <- (++) <$> liftExp x env <*> liftFun f env - t_var <- bind mdl f_arr - - launch' (cta,grid,smem) fscan (d_out ++ d_in0 ++ d_bks ++ t_var ++ map CUDA.IArg [n,interval]) - launch' (cta,1,smem) fscan (d_bks ++ d_bks ++ d_sum ++ map CUDA.IArg [grid,interval]) - launch' (cta,grid,smem) fadd (d_out ++ d_bks ++ map CUDA.IArg [n,interval]) - - freeArray in0 - freeArray bks - release f_arr - return (a_out, a_sum) - -dispatchScan _ _ _ = - error "we can never get here" - - -- Initiate the device computation. The first version selects launch parameters -- automatically, the second requires them explicitly. This tuple contains -- threads per block, grid size, and dynamic shared memory, respectively. @@ -447,14 +447,18 @@ -- FIXME: small arrays are relocated by the GC ad = fst . runArrayData $ (,undefined) <$> newArrayData (1024 `max` Sugar.size sh) --- Extract shape dimensions as a list of function parameters. Not that this will --- convert to the base integer width of the device, namely, 32-bits. Singleton --- dimensions are considered to be of unit size. +-- Extract shape dimensions as a list of function parameters. Note that this +-- will convert to the base integer width of the device, namely, 32-bits. +-- Singleton dimensions are considered to be of unit size. -- +-- Internally, Accelerate uses snoc-based tuple projection, while the data +-- itself is stored in reading order. Ensure we match the behaviour of regular +-- tuples and code generation thereof. +-- convertIx :: Ix dim => dim -> [CUDA.FunParam] convertIx = post . map CUDA.IArg . shapeToList where post [] = [CUDA.IArg 1] - post xs = xs + post xs = reverse xs -- Convert a slice specification into storable index projection components. -- Note: implicit conversion Int -> Int32 diff --git a/Data/Array/Accelerate/CUDA/State.hs b/Data/Array/Accelerate/CUDA/State.hs --- a/Data/Array/Accelerate/CUDA/State.hs +++ b/Data/Array/Accelerate/CUDA/State.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE TemplateHaskell, TupleSections, TypeOperators #-} +{-# LANGUAGE CPP, TemplateHaskell, TupleSections, TypeOperators #-} -- | -- Module : Data.Array.Accelerate.CUDA.State -- Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell @@ -15,7 +15,7 @@ module Data.Array.Accelerate.CUDA.State ( - evalCUDA, runCUDA, CIO, unique, outputDir, deviceProps, memoryTable, kernelTable, + evalCUDA, runCUDA, CIO, unique, outputDir, deviceProps, deviceContext, memoryTable, kernelTable, AccTable, KernelEntry(KernelEntry), kernelName, kernelStatus, MemTable, MemoryEntry(MemoryEntry), refcount, memsize, arena, @@ -28,25 +28,33 @@ import Control.Category import Data.Int -import Data.Maybe -import Data.Binary +import Data.IORef import Data.Record.Label -import Control.Arrow import Control.Applicative -import Control.Exception -import Control.Monad.State (StateT(..), liftM) +import Control.Monad +import Control.Monad.State (StateT(..)) import Data.HashTable (HashTable) +import Foreign.Ptr import qualified Data.HashTable as HT +import qualified Foreign.CUDA.Driver as CUDA import System.Directory import System.FilePath -import System.Posix.Process import System.Posix.Types (ProcessID) +import System.Mem.Weak +import System.IO.Unsafe -import Foreign.Ptr -import qualified Foreign.CUDA.Driver as CUDA +import Data.Array.Accelerate.CUDA.Analysis.Device +#ifdef ACCELERATE_CUDA_PERSISTENT_CACHE +import Data.Binary (encodeFile, decodeFile) +import Control.Arrow (second) +import Paths_accelerate (getDataDir) +#else +import System.Posix.Process (getProcessID) +#endif + -- Types -- ~~~~~ @@ -58,11 +66,12 @@ type CIO = StateT CUDAState IO data CUDAState = CUDAState { - _unique :: Int, - _outputDir :: FilePath, - _deviceProps :: CUDA.DeviceProperties, - _memoryTable :: MemTable, - _kernelTable :: AccTable + _unique :: Int, + _outputDir :: FilePath, + _deviceProps :: CUDA.DeviceProperties, + _deviceContext :: CUDA.Context, + _memoryTable :: MemTable, + _kernelTable :: AccTable } -- | @@ -90,78 +99,110 @@ _arena :: WordPtr } +$(mkLabels [''CUDAState, ''MemoryEntry, ''KernelEntry]) --- The CUDA State Monad --- ~~~~~~~~~~~~~~~~~~~~ +-- Execution State +-- ~~~~~~~~~~~~~~~ -- Return the output directory for compilation by-products, creating if it does --- not exist. This currently maps to a temporary directory, but could be mode to --- point towards a persistent cache (eg: getAppUserDataDirectory) +-- not exist. -- getOutputDir :: IO FilePath getOutputDir = do +#ifdef ACCELERATE_CUDA_PERSISTENT_CACHE + tmp <- getDataDir + dir <- canonicalizePath $ tmp "cache" +#else tmp <- getTemporaryDirectory pid <- getProcessID dir <- canonicalizePath $ tmp "ac" ++ show pid +#endif createDirectoryIfMissing True dir return dir --- Store the kernel module map to file to the given directory +-- Store the kernel module map to file -- -save :: FilePath -> AccTable -> IO () -save f m = encodeFile f . map (second _kernelName) . filter compiled =<< HT.toList m - where - compiled (_,KernelEntry _ (Right _)) = True - compiled _ = False +saveIndexFile :: CUDAState -> IO () +#ifdef ACCELERATE_CUDA_PERSISTENT_CACHE +saveIndexFile s = encodeFile (_outputDir s "_index") . map (second _kernelName) =<< HT.toList (_kernelTable s) +#else +saveIndexFile _ = return () +#endif -- Read the kernel index map file (if it exists), loading modules into the -- current context -- -load :: FilePath -> IO (AccTable, Int) -load f = do +loadIndexFile :: FilePath -> IO (AccTable, Int) +#ifdef ACCELERATE_CUDA_PERSISTENT_CACHE +loadIndexFile f = do x <- doesFileExist f e <- if x then mapM reload =<< decodeFile f else return [] - (,length e) <$> HT.fromList HT.hashString e where - reload (k,n) = - (k,) . KernelEntry n . Right <$> CUDA.loadFile (n `replaceExtension` ".cubin") + reload (k,n) = (k,) . KernelEntry n . Right <$> CUDA.loadFile (n `replaceExtension` ".cubin") +#else +loadIndexFile _ = (,0) <$> HT.new (==) HT.hashString +#endif +-- Select and initialise the CUDA device, and create a new execution context. +-- This will be done only once per program execution, as initialising the CUDA +-- context is relatively expensive. +-- +-- Would like to put the finaliser on the state token, since finalising the +-- context affects the various hash tables. However, this places the finaliser +-- on the CUDAState "box", and the box is removed by optimisations causing the +-- finaliser to fire prematurely. +-- +initialise :: IO CUDAState +initialise = do + CUDA.initialise [] + (d,prp) <- selectBestDevice + ctx <- CUDA.create d [CUDA.SchedAuto] + dir <- getOutputDir + mem <- HT.new (==) fromIntegral + (knl,n) <- loadIndexFile (dir "_index") + addFinalizer ctx (CUDA.destroy ctx) + return $ CUDAState n dir prp ctx mem knl + + -- | -- Evaluate a CUDA array computation under a newly initialised environment, -- discarding the final state. -- evalCUDA :: CIO a -> IO a -evalCUDA = liftM fst . runCUDA +evalCUDA = liftM fst . runCUDA runCUDA :: CIO a -> IO (a, CUDAState) runCUDA acc = - bracket (initialise Nothing) finalise $ \(dev,_ctx) -> do - props <- CUDA.props dev - dir <- getOutputDir - tab <- HT.new (==) fromIntegral - (k,n) <- load (dir "_index") - (a,s) <- runStateT acc (CUDAState n dir props tab k) - -- - -- TLM 2010-06-05: assert all memory has been released ?? - -- does CUDA.destroy release memory ?? - save (dir "_index") (_kernelTable s) + let + {-# NOINLINE ref #-} -- hic sunt dracones: truly unsafe use of unsafePerformIO + ref = unsafePerformIO (initialise >>= newIORef) + in do + state <- readIORef ref + clearMemTable state -- ugly kludge + (a,s) <- runStateT acc state + saveIndexFile s + writeIORef ref s return (a,s) - where - finalise = CUDA.destroy . snd - initialise n = do - CUDA.initialise [] - dev <- CUDA.device (fromMaybe 0 n) -- TLM: select the "best" device ?? - ctx <- CUDA.create dev [CUDA.SchedAuto] - return (dev, ctx) +-- In case of memory leaks, which we should fix, manually release any lingering +-- device arrays. These would otherwise remain until the program exits. +-- +clearMemTable :: CUDAState -> IO () +clearMemTable st = do + CUDA.sync -- TLM: not blocking?? + entries <- HT.toList (_memoryTable st) + forM_ entries $ \(k,v) -> do + HT.delete (_memoryTable st) k + CUDA.free (CUDA.wordPtrToDevPtr (_arena v)) -$(mkLabels [''CUDAState, ''MemoryEntry, ''KernelEntry]) +-- Utility +-- ~~~~~~~ -- | A unique name supply -- diff --git a/Data/Array/Accelerate/Internal/Check.hs b/Data/Array/Accelerate/Internal/Check.hs --- a/Data/Array/Accelerate/Internal/Check.hs +++ b/Data/Array/Accelerate/Internal/Check.hs @@ -57,7 +57,8 @@ error file line kind loc msg = P.error . unlines $ (if kind == Internal - then (["*** Internal error in package accelerate ***" + then (["" + ,"*** Internal error in package accelerate ***" ,"*** Please submit a bug report at http://trac.haskell.org/accelerate"]++) else id) [ file ++ ":" ++ show line ++ " (" ++ loc ++ "): " ++ msg ] diff --git a/Data/Array/Accelerate/Interpreter.hs b/Data/Array/Accelerate/Interpreter.hs --- a/Data/Array/Accelerate/Interpreter.hs +++ b/Data/Array/Accelerate/Interpreter.hs @@ -1,5 +1,5 @@ {-# OPTIONS_HADDOCK prune #-} -{-# LANGUAGE GADTs, BangPatterns, PatternGuards #-} +{-# LANGUAGE CPP, GADTs, BangPatterns, PatternGuards #-} {-# LANGUAGE TypeFamilies, ScopedTypeVariables, FlexibleContexts #-} -- | -- Module : Data.Array.Accelerate.Interpreter @@ -50,7 +50,9 @@ import qualified Data.Array.Accelerate.Smart as Sugar import qualified Data.Array.Accelerate.Array.Sugar as Sugar +#include "accelerate.h" + -- Program execution -- ----------------- @@ -144,12 +146,10 @@ reshapeOp :: Sugar.Ix dim => dim -> Delayed (Array dim' e) -> Delayed (Array dim e) reshapeOp newShape darr@(DelayedArray {shapeDA = oldShape}) - | Sugar.size newShape == size oldShape = let Array _ adata = force darr in - delay $ Array (Sugar.fromElem newShape) adata - | otherwise - = error "Data.Array.Accelerate.Interpreter.reshape: shape mismatch" + BOUNDS_CHECK(check) "reshape" "shape mismatch" (Sugar.size newShape == size oldShape) + $ delay $ Array (Sugar.fromElem newShape) adata replicateOp :: (Sugar.Ix dim, Sugar.Elem slix) => SliceIndex (Sugar.ElemRepr slix) @@ -201,11 +201,9 @@ in ((sl', sz), \(ix, i) -> (pf' ix, i)) restrict (SliceFixed sliceIndex) (slix, i) (sh, sz) - | i < sz = let (sl', pf') = restrict sliceIndex slix sh in - (sl', \ix -> (pf' ix, i)) - | otherwise = error "Index out of bounds" + BOUNDS_CHECK(checkIndex) "index" i sz $ (sl', \ix -> (pf' ix, i)) mapOp :: Sugar.Elem e' => (e -> e') diff --git a/Data/Array/Accelerate/Language.hs b/Data/Array/Accelerate/Language.hs --- a/Data/Array/Accelerate/Language.hs +++ b/Data/Array/Accelerate/Language.hs @@ -338,7 +338,35 @@ tuple = tup5 untuple = untup5 +instance (Elem a, Elem b, Elem c, Elem d, Elem e, Elem f) + => Tuple (Exp a, Exp b, Exp c, Exp d, Exp e, Exp f) where + type TupleT (Exp a, Exp b, Exp c, Exp d, Exp e, Exp f) + = Exp (a, b, c, d, e, f) + tuple = tup6 + untuple = untup6 +instance (Elem a, Elem b, Elem c, Elem d, Elem e, Elem f, Elem g) + => Tuple (Exp a, Exp b, Exp c, Exp d, Exp e, Exp f, Exp g) where + type TupleT (Exp a, Exp b, Exp c, Exp d, Exp e, Exp f, Exp g) + = Exp (a, b, c, d, e, f, g) + tuple = tup7 + untuple = untup7 + +instance (Elem a, Elem b, Elem c, Elem d, Elem e, Elem f, Elem g, Elem h) + => Tuple (Exp a, Exp b, Exp c, Exp d, Exp e, Exp f, Exp g, Exp h) where + type TupleT (Exp a, Exp b, Exp c, Exp d, Exp e, Exp f, Exp g, Exp h) + = Exp (a, b, c, d, e, f, g, h) + tuple = tup8 + untuple = untup8 + +instance (Elem a, Elem b, Elem c, Elem d, Elem e, Elem f, Elem g, Elem h, Elem i) + => Tuple (Exp a, Exp b, Exp c, Exp d, Exp e, Exp f, Exp g, Exp h, Exp i) where + type TupleT (Exp a, Exp b, Exp c, Exp d, Exp e, Exp f, Exp g, Exp h, Exp i) + = Exp (a, b, c, d, e, f, g, h, i) + tuple = tup9 + untuple = untup9 + + -- |Extract the first component of a pair -- fst :: forall a b. (Elem a, Elem b) => Exp (a, b) -> Exp a @@ -454,20 +482,20 @@ -- toInteger = -- makes no sense instance (Elem t, IsFloating t) => Floating (Exp t) where - pi = mkPi - sin = mkSin - cos = mkCos - tan = mkTan - asin = mkAsin - acos = mkAcos - atan = mkAtan - asinh = mkAsinh - acosh = mkAcosh - atanh = mkAtanh - exp = mkExpFloating - sqrt = mkSqrt - log = mkLog - (**) = mkFPow + pi = mkPi + sin = mkSin + cos = mkCos + tan = mkTan + asin = mkAsin + acos = mkAcos + atan = mkAtan + asinh = mkAsinh + acosh = mkAcosh + atanh = mkAtanh + exp = mkExpFloating + sqrt = mkSqrt + log = mkLog + (**) = mkFPow logBase = mkLogBase -- FIXME: add other ops diff --git a/Data/Array/Accelerate/Pretty.hs b/Data/Array/Accelerate/Pretty.hs --- a/Data/Array/Accelerate/Pretty.hs +++ b/Data/Array/Accelerate/Pretty.hs @@ -104,7 +104,7 @@ prettyBoundary _ Clamp = text "Clamp" prettyBoundary _ Mirror = text "Mirror" prettyBoundary _ Wrap = text "Wrap" -prettyBoundary _ (Constant e) = parens $ text "Wrap" <+> text (show (toElem e :: e)) +prettyBoundary _ (Constant e) = parens $ text "Constant" <+> text (show (toElem e :: e)) prettyArrOp :: String -> [Doc] -> Doc prettyArrOp name docs = hang (text name) 2 $ sep docs diff --git a/Data/Array/Accelerate/Smart.hs b/Data/Array/Accelerate/Smart.hs --- a/Data/Array/Accelerate/Smart.hs +++ b/Data/Array/Accelerate/Smart.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE GADTs, TypeFamilies, ScopedTypeVariables, FlexibleContexts #-} +{-# LANGUAGE CPP, GADTs, TypeFamilies, ScopedTypeVariables, FlexibleContexts #-} {-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, TypeSynonymInstances #-} -- Module : Data.Array.Accelerate.Smart @@ -29,7 +29,8 @@ constant, -- * Smart constructors and destructors for tuples - tup2, tup3, tup4, tup5, untup2, untup3, untup4, untup5, + tup2, tup3, tup4, tup5, tup6, tup7, tup8, tup9, + untup2, untup3, untup4, untup5, untup6, untup7, untup8, untup9, -- * Smart constructors for constants mkMinBound, mkMaxBound, mkPi, @@ -62,7 +63,9 @@ import qualified Data.Array.Accelerate.AST as AST import Data.Array.Accelerate.Pretty () +#include "accelerate.h" + -- Monadic array computations -- -------------------------- @@ -139,7 +142,7 @@ -> Boundary a -> Acc (Array dim a) -> Acc (Array dim b) - Stencil2 :: (Ix dim, Elem a, Elem b, Elem c, + Stencil2 :: (Ix dim, Elem a, Elem b, Elem c, Stencil dim a stencil1, Stencil dim b stencil2) => (stencil1 -> stencil2 -> Exp c) -> Boundary a @@ -273,8 +276,8 @@ prjIdx 0 (PushLayout _ ix) = fromJust (gcast ix) -- can't go wrong unless the library is wrong! prjIdx n (PushLayout l _) = prjIdx (n - 1) l -prjIdx _ EmptyLayout - = error "Data.Array.Accelerate.Smart.prjIdx: internal error" +prjIdx _ EmptyLayout = + INTERNAL_ERROR(error) "prjIdx" "inconsistent valuation" -- |Convert an open expression with given environment layouts. -- @@ -643,6 +646,36 @@ = Tuple $ NilTup `SnocTup` x1 `SnocTup` x2 `SnocTup` x3 `SnocTup` x4 `SnocTup` x5 +tup6 :: (Elem a, Elem b, Elem c, Elem d, Elem e, Elem f) + => (Exp a, Exp b, Exp c, Exp d, Exp e, Exp f) -> Exp (a, b, c, d, e, f) +tup6 (x1, x2, x3, x4, x5, x6) + = Tuple $ + NilTup `SnocTup` x1 `SnocTup` x2 `SnocTup` x3 `SnocTup` x4 `SnocTup` x5 `SnocTup` x6 + +tup7 :: (Elem a, Elem b, Elem c, Elem d, Elem e, Elem f, Elem g) + => (Exp a, Exp b, Exp c, Exp d, Exp e, Exp f, Exp g) + -> Exp (a, b, c, d, e, f, g) +tup7 (x1, x2, x3, x4, x5, x6, x7) + = Tuple $ + NilTup `SnocTup` x1 `SnocTup` x2 `SnocTup` x3 + `SnocTup` x4 `SnocTup` x5 `SnocTup` x6 `SnocTup` x7 + +tup8 :: (Elem a, Elem b, Elem c, Elem d, Elem e, Elem f, Elem g, Elem h) + => (Exp a, Exp b, Exp c, Exp d, Exp e, Exp f, Exp g, Exp h) + -> Exp (a, b, c, d, e, f, g, h) +tup8 (x1, x2, x3, x4, x5, x6, x7, x8) + = Tuple $ + NilTup `SnocTup` x1 `SnocTup` x2 `SnocTup` x3 `SnocTup` x4 + `SnocTup` x5 `SnocTup` x6 `SnocTup` x7 `SnocTup` x8 + +tup9 :: (Elem a, Elem b, Elem c, Elem d, Elem e, Elem f, Elem g, Elem h, Elem i) + => (Exp a, Exp b, Exp c, Exp d, Exp e, Exp f, Exp g, Exp h, Exp i) + -> Exp (a, b, c, d, e, f, g, h, i) +tup9 (x1, x2, x3, x4, x5, x6, x7, x8, x9) + = Tuple $ + NilTup `SnocTup` x1 `SnocTup` x2 `SnocTup` x3 `SnocTup` x4 + `SnocTup` x5 `SnocTup` x6 `SnocTup` x7 `SnocTup` x8 `SnocTup` x9 + untup2 :: (Elem a, Elem b) => Exp (a, b) -> (Exp a, Exp b) untup2 e = ((SuccTupIdx ZeroTupIdx) `Prj` e, ZeroTupIdx `Prj` e) @@ -665,6 +698,48 @@ SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx)) `Prj` e, SuccTupIdx (SuccTupIdx ZeroTupIdx) `Prj` e, SuccTupIdx ZeroTupIdx `Prj` e, + ZeroTupIdx `Prj` e) + +untup6 :: (Elem a, Elem b, Elem c, Elem d, Elem e, Elem f) + => Exp (a, b, c, d, e, f) -> (Exp a, Exp b, Exp c, Exp d, Exp e, Exp f) +untup6 e = (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx)))) `Prj` e, + SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx))) `Prj` e, + SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx)) `Prj` e, + SuccTupIdx (SuccTupIdx ZeroTupIdx) `Prj` e, + SuccTupIdx ZeroTupIdx `Prj` e, + ZeroTupIdx `Prj` e) + +untup7 :: (Elem a, Elem b, Elem c, Elem d, Elem e, Elem f, Elem g) + => Exp (a, b, c, d, e, f, g) -> (Exp a, Exp b, Exp c, Exp d, Exp e, Exp f, Exp g) +untup7 e = (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx))))) `Prj` e, + SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx)))) `Prj` e, + SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx))) `Prj` e, + SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx)) `Prj` e, + SuccTupIdx (SuccTupIdx ZeroTupIdx) `Prj` e, + SuccTupIdx ZeroTupIdx `Prj` e, + ZeroTupIdx `Prj` e) + +untup8 :: (Elem a, Elem b, Elem c, Elem d, Elem e, Elem f, Elem g, Elem h) + => Exp (a, b, c, d, e, f, g, h) -> (Exp a, Exp b, Exp c, Exp d, Exp e, Exp f, Exp g, Exp h) +untup8 e = (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx)))))) `Prj` e, + SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx))))) `Prj` e, + SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx)))) `Prj` e, + SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx))) `Prj` e, + SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx)) `Prj` e, + SuccTupIdx (SuccTupIdx ZeroTupIdx) `Prj` e, + SuccTupIdx ZeroTupIdx `Prj` e, + ZeroTupIdx `Prj` e) + +untup9 :: (Elem a, Elem b, Elem c, Elem d, Elem e, Elem f, Elem g, Elem h, Elem i) + => Exp (a, b, c, d, e, f, g, h, i) -> (Exp a, Exp b, Exp c, Exp d, Exp e, Exp f, Exp g, Exp h, Exp i) +untup9 e = (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx))))))) `Prj` e, + SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx)))))) `Prj` e, + SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx))))) `Prj` e, + SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx)))) `Prj` e, + SuccTupIdx (SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx))) `Prj` e, + SuccTupIdx (SuccTupIdx (SuccTupIdx ZeroTupIdx)) `Prj` e, + SuccTupIdx (SuccTupIdx ZeroTupIdx) `Prj` e, + SuccTupIdx ZeroTupIdx `Prj` e, ZeroTupIdx `Prj` e) -- Smart constructor for constants diff --git a/accelerate.cabal b/accelerate.cabal --- a/accelerate.cabal +++ b/accelerate.cabal @@ -1,5 +1,5 @@ Name: accelerate -Version: 0.8.0.0 +Version: 0.8.1.0 Cabal-version: >= 1.6 Tested-with: GHC >= 6.12.3 Build-type: Simple @@ -19,8 +19,10 @@ and 'Bool' arrays. . Known bugs in this version: - + . + * New in 0.8.1.0: bug fixes and some performance tweaks + . * New in 0.8.0.0: 'replicate', 'slice' and 'foldSeg' supported in the CUDA backend; frontend and interpreter support for 'stencil'; bug fixes . @@ -40,18 +42,20 @@ -- Should be in the Library stanza, and only enabled for the CUDA backend, -- but Cabal does not support that. -Data-dir: cubits -Data-files: accelerate_cuda_extras.h - accelerate_cuda_utils.h - backpermute.inl - fold.inl - fold_segmented.inl - map.inl - permute.inl - replicate.inl - slice.inl - zipWith.inl - thrust/scan_safe.inl +Data-files: cubits/accelerate_cuda_extras.h + cubits/accelerate_cuda_function.h + cubits/accelerate_cuda_shape.h + cubits/accelerate_cuda_texture.h + cubits/accelerate_cuda_util.h + cubits/backpermute.inl + cubits/fold.inl + cubits/fold_segmented.inl + cubits/map.inl + cubits/permute.inl + cubits/replicate.inl + cubits/slice.inl + cubits/zipWith.inl + cubits/thrust/scan_safe.inl Extra-source-files: INSTALL include/accelerate.h @@ -80,6 +84,10 @@ Description: Enable the CUDA parallel backend for NVIDIA GPUs Default: True +Flag pcache + Description: Enable the persistent caching of the compiled CUDA modules (experimental) + Default: False + Flag test-suite Description: Export extra test modules Default: False @@ -112,7 +120,7 @@ Build-depends: binary, bytestring, containers, - cuda >= 0.2 && < 0.3, + cuda >= 0.2.2 && < 0.3, directory, fclabels >= 0.9 && < 1.0, filepath, @@ -134,7 +142,7 @@ If flag(test-suite) Exposed-modules: Data.Array.Accelerate.Test - Data.Array.Accelerate.Test.QuickCheck + Other-modules: Data.Array.Accelerate.Test.QuickCheck Data.Array.Accelerate.Test.QuickCheck.Arbitrary Other-modules: Data.Array.Accelerate.Internal.Check @@ -143,6 +151,7 @@ Data.Array.Accelerate.Array.Representation Data.Array.Accelerate.Array.Sugar Data.Array.Accelerate.Analysis.Type + Data.Array.Accelerate.Analysis.Shape Data.Array.Accelerate.AST Data.Array.Accelerate.Debug Data.Array.Accelerate.Language @@ -156,7 +165,8 @@ If flag(cuda) CPP-options: -DACCELERATE_CUDA_BACKEND - Other-modules: Data.Array.Accelerate.CUDA.Analysis.Hash + Other-modules: Data.Array.Accelerate.CUDA.Analysis.Device + Data.Array.Accelerate.CUDA.Analysis.Hash Data.Array.Accelerate.CUDA.Analysis.Launch Data.Array.Accelerate.CUDA.Array.Data Data.Array.Accelerate.CUDA.Array.Device @@ -169,6 +179,9 @@ Data.Array.Accelerate.CUDA.Execute Data.Array.Accelerate.CUDA.Smart Data.Array.Accelerate.CUDA.State + + if flag(pcache) + CPP-options: -DACCELERATE_CUDA_PERSISTENT_CACHE if flag(bounds-checks) cpp-options: -DACCELERATE_BOUNDS_CHECKS diff --git a/cubits/accelerate_cuda_extras.h b/cubits/accelerate_cuda_extras.h --- a/cubits/accelerate_cuda_extras.h +++ b/cubits/accelerate_cuda_extras.h @@ -1,344 +1,21 @@ /* ----------------------------------------------------------------------------- * - * Module : Extras - * Copyright : (c) [2009..2010] Trevor L. McDonell - * License : BSD + * Module : Extras + * Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell + * License : BSD3 * + * Maintainer : Trevor L. McDonell + * Stability : experimental + * * ---------------------------------------------------------------------------*/ #ifndef __ACCELERATE_CUDA_EXTRAS_H__ #define __ACCELERATE_CUDA_EXTRAS_H__ -#include -#include - -#include - -/* ----------------------------------------------------------------------------- - * Textures - * ----------------------------------------------------------------------------- - * - * CUDA texture definitions and access functions are defined in terms of - * templates, and hence only available through the C++ interface. Expose some - * dummy wrappers to enable parsing with language-c. - * - * We don't have a definition for `Int' or `Word', since the bitwidth of the - * Haskell and C types may be different. - * - * NOTE ON 64-BIT TYPES - * The CUDA device uses little-endian arithmetic. We haven't accounted for the - * fact that this may be different on the host, for both initial data transfer - * and the unpacking below. - */ -#if defined(__cplusplus) && defined(__CUDACC__) - -typedef texture TexWord64; -typedef texture TexWord32; -typedef texture TexWord16; -typedef texture TexWord8; -typedef texture TexInt64; -typedef texture TexInt32; -typedef texture TexInt16; -typedef texture TexInt8; -typedef texture TexFloat; -typedef texture TexCChar; - -static __inline__ __device__ uint8_t indexArray(TexWord8 t, const int x) { return tex1Dfetch(t,x); } -static __inline__ __device__ uint16_t indexArray(TexWord16 t, const int x) { return tex1Dfetch(t,x); } -static __inline__ __device__ uint32_t indexArray(TexWord32 t, const int x) { return tex1Dfetch(t,x); } -static __inline__ __device__ uint64_t indexArray(TexWord64 t, const int x) -{ - union { uint2 x; uint64_t y; } v; - v.x = tex1Dfetch(t,x); - return v.y; -} - -static __inline__ __device__ int8_t indexArray(TexInt8 t, const int x) { return tex1Dfetch(t,x); } -static __inline__ __device__ int16_t indexArray(TexInt16 t, const int x) { return tex1Dfetch(t,x); } -static __inline__ __device__ int32_t indexArray(TexInt32 t, const int x) { return tex1Dfetch(t,x); } -static __inline__ __device__ int64_t indexArray(TexInt64 t, const int x) -{ - union { int2 x; int64_t y; } v; - v.x = tex1Dfetch(t,x); - return v.y; -} - -static __inline__ __device__ float indexArray(TexFloat t, const int x) { return tex1Dfetch(t,x); } -static __inline__ __device__ char indexArray(TexCChar t, const int x) { return tex1Dfetch(t,x); } - -#if defined(__LP64__) -typedef TexInt64 TexCLong; -typedef TexWord64 TexCULong; -#else -typedef TexInt32 TexCLong; -typedef TexWord32 TexCULong; -#endif - -/* - * Synonyms for C-types. NVCC will force Ints to be 32-bits. - */ -typedef TexInt8 TexCSChar; -typedef TexWord8 TexCUChar; -typedef TexInt16 TexCShort; -typedef TexWord16 TexCUShort; -typedef TexInt32 TexCInt; -typedef TexWord32 TexCUInt; -typedef TexInt64 TexCLLong; -typedef TexWord64 TexCULLong; -typedef TexFloat TexCFloat; - -/* - * Doubles, only available when compiled for Compute 1.3 and greater - */ -#if !defined(CUDA_NO_SM_13_DOUBLE_INTRINSICS) -typedef texture TexDouble; -typedef TexDouble TexCDouble; - -static __inline__ __device__ double indexDArray(TexDouble t, const int x) -{ - int2 v = tex1Dfetch(t,x); - return __hiloint2double(v.y,v.x); -} -#endif - -#else - -typedef void* TexWord64; -typedef void* TexWord32; -typedef void* TexWord16; -typedef void* TexWord8; -typedef void* TexInt64; -typedef void* TexInt32; -typedef void* TexInt16; -typedef void* TexInt8; -typedef void* TexCShort; -typedef void* TexCUShort; -typedef void* TexCInt; -typedef void* TexCUInt; -typedef void* TexCLong; -typedef void* TexCULong; -typedef void* TexCLLong; -typedef void* TexCULLong; -typedef void* TexFloat; -typedef void* TexDouble; -typedef void* TexCFloat; -typedef void* TexCDouble; -typedef void* TexCChar; -typedef void* TexCSChar; -typedef void* TexCUChar; - -void* indexArray(const void*, const int); - -#endif - -/* ----------------------------------------------------------------------------- - * Shapes - * -------------------------------------------------------------------------- */ - -typedef int32_t Ix; -typedef Ix DIM0; -typedef Ix DIM1; -typedef struct { Ix a1,a0; } DIM2; -typedef struct { Ix a2,a1,a0; } DIM3; -typedef struct { Ix a3,a2,a1,a0; } DIM4; -typedef struct { Ix a4,a3,a2,a1,a0; } DIM5; - -#ifdef __cplusplus - -/* - * Number of dimensions of a shape - */ -static __inline__ __device__ int dim(DIM1 sh) { return 1; } -static __inline__ __device__ int dim(DIM2 sh) { return 2; } -static __inline__ __device__ int dim(DIM3 sh) { return 3; } -static __inline__ __device__ int dim(DIM4 sh) { return 4; } -static __inline__ __device__ int dim(DIM5 sh) { return 5; } - -/* - * Yield the total number of elements in a shape - */ -static __inline__ __device__ int size(DIM1 sh) { return sh; } -static __inline__ __device__ int size(DIM2 sh) { return sh.a0 * sh.a1; } -static __inline__ __device__ int size(DIM3 sh) { return sh.a0 * sh.a1 * sh.a2; } -static __inline__ __device__ int size(DIM4 sh) { return sh.a0 * sh.a1 * sh.a2 * sh.a3; } -static __inline__ __device__ int size(DIM5 sh) { return sh.a0 * sh.a1 * sh.a2 * sh.a3 * sh.a4; } - -/* - * Convert the individual dimensions of a linear array into a shape - */ -static __inline__ __device__ DIM1 shape(Ix a) -{ - return a; -} - -static __inline__ __device__ DIM2 shape(Ix b, Ix a) -{ - DIM2 sh = { b, a }; - return sh; -} - -static __inline__ __device__ DIM3 shape(Ix c, Ix b, Ix a) -{ - DIM3 sh = { c, b, a }; - return sh; -} - -static __inline__ __device__ DIM4 shape(Ix d, Ix c, Ix b, Ix a) -{ - DIM4 sh = { d, c, b, a }; - return sh; -} - -static __inline__ __device__ DIM5 shape(Ix e, Ix d, Ix c, Ix b, Ix a) -{ - DIM5 sh = { e, d, c, b, a }; - return sh; -} - -/* - * Yield the index position in a linear, row-major representation of the array. - * First argument is the shape of the array, the second the index - */ -static __inline__ __device__ Ix toIndex(DIM1 sh, DIM1 ix) -{ - return ix; -} - -static __inline__ __device__ Ix toIndex(DIM2 sh, DIM2 ix) -{ - DIM1 sh_ = sh.a1; - DIM1 ix_ = ix.a1; - return toIndex(sh_, ix_) * sh.a0 + ix.a0; -} - -static __inline__ __device__ Ix toIndex(DIM3 sh, DIM3 ix) -{ - DIM2 sh_ = { sh.a2, sh.a1 }; - DIM2 ix_ = { ix.a2, ix.a1 }; - return toIndex(sh_, ix_) * sh.a0 + ix.a0; -} - -static __inline__ __device__ Ix toIndex(DIM4 sh, DIM4 ix) -{ - DIM3 sh_ = { sh.a3, sh.a2, sh.a1 }; - DIM3 ix_ = { ix.a3, ix.a2, ix.a1 }; - return toIndex(sh_, ix_) * sh.a0 + ix.a0; -} - -static __inline__ __device__ Ix toIndex(DIM5 sh, DIM5 ix) -{ - DIM4 sh_ = { sh.a4, sh.a3, sh.a2, sh.a1 }; - DIM4 ix_ = { ix.a4, ix.a3, ix.a2, ix.a1 }; - return toIndex(sh_, ix_) * sh.a0 + ix.a0; -} - -/* - * Inverse of 'toIndex' - */ -static __inline__ __device__ DIM1 fromIndex(DIM1 sh, Ix ix) -{ - return ix; -} - -static __inline__ __device__ DIM2 fromIndex(DIM2 sh, Ix ix) -{ - DIM1 sh_ = shape(sh.a1); - DIM1 ix_ = fromIndex(sh_, ix / sh.a0); - return shape(ix_, ix % sh.a0); -} - -static __inline__ __device__ DIM3 fromIndex(DIM3 sh, Ix ix) -{ - DIM2 sh_ = shape(sh.a2, sh.a1); - DIM2 ix_ = fromIndex(sh_, ix / sh.a0); - return shape(ix_.a1, ix_.a0, ix % sh.a0); -} - -static __inline__ __device__ DIM4 fromIndex(DIM4 sh, Ix ix) -{ - DIM3 sh_ = shape(sh.a3, sh.a2, sh.a1); - DIM3 ix_ = fromIndex(sh_, ix / sh.a0); - return shape(ix_.a2, ix_.a1, ix_.a0, ix % sh.a0); -} - -static __inline__ __device__ DIM5 fromIndex(DIM5 sh, Ix ix) -{ - DIM4 sh_ = shape(sh.a4, sh.a3, sh.a2, sh.a1); - DIM4 ix_ = fromIndex(sh_, ix / sh.a0); - return shape(ix_.a3, ix_.a2, ix_.a1, ix_.a0, ix % sh.a0); -} - - -/* - * Test for the magic index `ignore` - */ -static __inline__ __device__ int ignore(DIM1 ix) -{ - return ix == -1; -} - -static __inline__ __device__ int ignore(DIM2 ix) -{ - return ix.a0 == -1 && ix.a1 == -1; -} - -static __inline__ __device__ int ignore(DIM3 ix) -{ - return ix.a0 == -1 && ix.a1 == -1 && ix.a2 == -1; -} - -static __inline__ __device__ int ignore(DIM4 ix) -{ - return ix.a0 == -1 && ix.a1 == -1 && ix.a2 == -1 && ix.a3 == -1; -} - -static __inline__ __device__ int ignore(DIM5 ix) -{ - return ix.a0 == -1 && ix.a1 == -1 && ix.a2 == -1 && ix.a3 == -1 && ix.a4 == -1; -} - - -#else - -int dim(Ix sh); -int size(Ix sh); -int shape(Ix sh); -int toIndex(Ix sh, Ix ix); -int fromIndex(Ix sh, Ix ix); +#include "accelerate_cuda_function.h" +#include "accelerate_cuda_shape.h" +#include "accelerate_cuda_texture.h" +#include "accelerate_cuda_util.h" #endif -/* ----------------------------------------------------------------------------- - * Functions - * -------------------------------------------------------------------------- */ - -#ifdef __cplusplus -template -static __inline__ __device__ T rotateL(const T x, int32_t i) -{ - return (i &= 8 * sizeof(x) - 1) == 0 ? x : x << i | x >> 8 * sizeof(x) - i; -} - -template -static __inline__ __device__ T rotateR(const T x, int32_t i) -{ - return (i &= 8 * sizeof(x) - 1) == 0 ? x : x >> i | x << 8 * sizeof(x) - i; -} - -template -static __inline__ __device__ T idiv(const T x, const T y) -{ - return x > 0 && y < 0 ? (x - y - 1) / y : (x < 0 && y > 0 ? (x - y + 1) / y : x / y); -} - -template -static __inline__ __device__ T mod(const T x, const T y) -{ - const T r = x % y; - return x > 0 && y < 0 || x < 0 && y > 0 ? (r != 0 ? r + y : 0) : r; -} -#endif - -#endif // __ACCELERATE_CUDA_EXTRAS_H__ - -// vim: filetype=cuda.c diff --git a/cubits/accelerate_cuda_function.h b/cubits/accelerate_cuda_function.h new file mode 100644 --- /dev/null +++ b/cubits/accelerate_cuda_function.h @@ -0,0 +1,129 @@ +/* ----------------------------------------------------------------------------- + * + * Module : Function + * Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell + * License : BSD3 + * + * Maintainer : Trevor L. McDonell + * Stability : experimental + * + * ---------------------------------------------------------------------------*/ + +#ifndef __ACCELERATE_CUDA_FUNCTION_H__ +#define __ACCELERATE_CUDA_FUNCTION_H__ + +#include +#include + +#ifdef __cplusplus + +/* ----------------------------------------------------------------------------- + * Device functions required to support generated code + * -------------------------------------------------------------------------- */ + +/* + * Left/Right bitwise rotation + */ +template +static __inline__ __device__ T rotateL(const T x, const int32_t i) +{ + const int32_t i8 = i & 8 * sizeof(x) - 1; + return i8 == 0 ? x : x << i8 | x >> 8 * sizeof(x) - i8; +} + +template +static __inline__ __device__ T rotateR(const T x, const int32_t i) +{ + const int32_t i8 = i & 8 * sizeof(x) - 1; + return i8 == 0 ? x : x >> i8 | x << 8 * sizeof(x) - i8; +} + +/* + * Integer division, truncated towards negative infinity + */ +template +static __inline__ __device__ T idiv(const T x, const T y) +{ + return x > 0 && y < 0 ? (x - y - 1) / y : (x < 0 && y > 0 ? (x - y + 1) / y : x / y); +} + +/* + * Integer modulus, Haskell style + */ +template +static __inline__ __device__ T mod(const T x, const T y) +{ + const T r = x % y; + return x > 0 && y < 0 || x < 0 && y > 0 ? (r != 0 ? r + y : 0) : r; +} + + +/* ----------------------------------------------------------------------------- + * Additional helper functions + * -------------------------------------------------------------------------- */ + +/* + * Determine if the input is a power of two + */ +template +static __inline__ __host__ __device__ T isPow2(const T x) +{ + return ((x&(x-1)) == 0); +} + +/* + * Compute the next highest power of two + */ +template +static __inline__ __host__ T ceilPow2(const T x) +{ +#if 0 + --x; + x |= x >> 1; + x |= x >> 2; + x |= x >> 4; + x |= x >> 8; + x |= x >> 16; + return ++x; +#endif + + return (isPow2(x)) ? x : 1u << (int) ceil(log2((double)x)); +} + +/* + * Compute the next lowest power of two + */ +template +static __inline__ __host__ T floorPow2(const T x) +{ +#if 0 + float nf = (float) n; + return 1 << (((*(int*)&nf) >> 23) - 127); +#endif + + int exp; + frexp(x, &exp); + return 1 << (exp - 1); +} + +/* + * computes next highest multiple of f from x + */ +template +static __inline__ __host__ T multiple(const T x, const T f) +{ + return ((x + (f-1)) / f); +} + +/* + * MS Excel-style CEIL() function. Rounds x up to nearest multiple of f + */ +template +static __inline__ __host__ T ceiling(const T x, const T f) +{ + return multiple(x, f) * f; +} + +#endif // __cplusplus +#endif // __ACCELERATE_CUDA_FUNCTION_H__ + diff --git a/cubits/accelerate_cuda_shape.h b/cubits/accelerate_cuda_shape.h new file mode 100644 --- /dev/null +++ b/cubits/accelerate_cuda_shape.h @@ -0,0 +1,304 @@ +/* ----------------------------------------------------------------------------- + * + * Module : Shape + * Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell + * License : BSD3 + * + * Maintainer : Trevor L. McDonell + * Stability : experimental + * + * ---------------------------------------------------------------------------*/ + +#ifndef __ACCELERATE_CUDA_SHAPE_H__ +#define __ACCELERATE_CUDA_SHAPE_H__ + +#include +#include + +typedef int32_t Ix; +typedef Ix DIM0; +typedef Ix DIM1; +typedef struct { Ix a1,a0; } DIM2; +typedef struct { Ix a2,a1,a0; } DIM3; +typedef struct { Ix a3,a2,a1,a0; } DIM4; +typedef struct { Ix a4,a3,a2,a1,a0; } DIM5; +typedef struct { Ix a5,a4,a3,a2,a1,a0; } DIM6; +typedef struct { Ix a6,a5,a4,a3,a2,a1,a0; } DIM7; +typedef struct { Ix a7,a6,a5,a4,a3,a2,a1,a0; } DIM8; +typedef struct { Ix a8,a7,a6,a5,a4,a3,a2,a1,a0; } DIM9; + +#ifdef __cplusplus + +/* + * Number of dimensions of a shape + */ +static __inline__ __device__ int dim(DIM1 sh) { return 1; } +static __inline__ __device__ int dim(DIM2 sh) { return 2; } +static __inline__ __device__ int dim(DIM3 sh) { return 3; } +static __inline__ __device__ int dim(DIM4 sh) { return 4; } +static __inline__ __device__ int dim(DIM5 sh) { return 5; } +static __inline__ __device__ int dim(DIM6 sh) { return 6; } +static __inline__ __device__ int dim(DIM7 sh) { return 7; } +static __inline__ __device__ int dim(DIM8 sh) { return 8; } +static __inline__ __device__ int dim(DIM9 sh) { return 9; } + +/* + * Yield the total number of elements in a shape + */ +static __inline__ __device__ int size(DIM1 sh) { return sh; } +static __inline__ __device__ int size(DIM2 sh) { return sh.a0 * sh.a1; } +static __inline__ __device__ int size(DIM3 sh) { return sh.a0 * sh.a1 * sh.a2; } +static __inline__ __device__ int size(DIM4 sh) { return sh.a0 * sh.a1 * sh.a2 * sh.a3; } +static __inline__ __device__ int size(DIM5 sh) { return sh.a0 * sh.a1 * sh.a2 * sh.a3 * sh.a4; } +static __inline__ __device__ int size(DIM6 sh) { return sh.a0 * sh.a1 * sh.a2 * sh.a3 * sh.a4 * sh.a5; } +static __inline__ __device__ int size(DIM7 sh) { return sh.a0 * sh.a1 * sh.a2 * sh.a3 * sh.a4 * sh.a5 * sh.a6; } +static __inline__ __device__ int size(DIM8 sh) { return sh.a0 * sh.a1 * sh.a2 * sh.a3 * sh.a4 * sh.a5 * sh.a6 * sh.a7; } +static __inline__ __device__ int size(DIM9 sh) { return sh.a0 * sh.a1 * sh.a2 * sh.a3 * sh.a4 * sh.a5 * sh.a6 * sh.a7 * sh.a8; } + +/* + * Convert the individual dimensions of a linear array into a shape + */ +static __inline__ __device__ DIM1 shape(Ix a) +{ + return a; +} + +static __inline__ __device__ DIM2 shape(Ix b, Ix a) +{ + DIM2 sh = { b, a }; + return sh; +} + +static __inline__ __device__ DIM3 shape(Ix c, Ix b, Ix a) +{ + DIM3 sh = { c, b, a }; + return sh; +} + +static __inline__ __device__ DIM4 shape(Ix d, Ix c, Ix b, Ix a) +{ + DIM4 sh = { d, c, b, a }; + return sh; +} + +static __inline__ __device__ DIM5 shape(Ix e, Ix d, Ix c, Ix b, Ix a) +{ + DIM5 sh = { e, d, c, b, a }; + return sh; +} + +static __inline__ __device__ DIM6 shape(Ix f, Ix e, Ix d, Ix c, Ix b, Ix a) +{ + DIM6 sh = { f, e, d, c, b, a }; + return sh; +} + +static __inline__ __device__ DIM7 shape(Ix g, Ix f, Ix e, Ix d, Ix c, Ix b, Ix a) +{ + DIM7 sh = { g, f, e, d, c, b, a }; + return sh; +} + +static __inline__ __device__ DIM8 shape(Ix h, Ix g, Ix f, Ix e, Ix d, Ix c, Ix b, Ix a) +{ + DIM8 sh = { h, g, f, e, d, c, b, a }; + return sh; +} + +static __inline__ __device__ DIM9 shape(Ix i, Ix h, Ix g, Ix f, Ix e, Ix d, Ix c, Ix b, Ix a) +{ + DIM9 sh = { i, h, g, f, e, d, c, b, a }; + return sh; +} + + +/* + * Yield the index position in a linear, row-major representation of the array. + * First argument is the shape of the array, the second the index + */ +static __inline__ __device__ Ix toIndex(DIM1 sh, DIM1 ix) +{ + return ix; +} + +static __inline__ __device__ Ix toIndex(DIM2 sh, DIM2 ix) +{ + DIM1 sh_ = sh.a0; + DIM1 ix_ = ix.a0; + return toIndex(sh_, ix_) * sh.a1 + ix.a1; +} + +static __inline__ __device__ Ix toIndex(DIM3 sh, DIM3 ix) +{ + DIM2 sh_ = { sh.a1, sh.a0 }; + DIM2 ix_ = { ix.a1, ix.a0 }; + return toIndex(sh_, ix_) * sh.a2 + ix.a2; +} + +static __inline__ __device__ Ix toIndex(DIM4 sh, DIM4 ix) +{ + DIM3 sh_ = { sh.a2, sh.a1, sh.a0 }; + DIM3 ix_ = { ix.a2, ix.a1, ix.a0 }; + return toIndex(sh_, ix_) * sh.a3 + ix.a3; +} + +static __inline__ __device__ Ix toIndex(DIM5 sh, DIM5 ix) +{ + DIM4 sh_ = { sh.a3, sh.a2, sh.a1, sh.a0 }; + DIM4 ix_ = { ix.a3, ix.a2, ix.a1, ix.a0 }; + return toIndex(sh_, ix_) * sh.a4 + ix.a4; +} + +static __inline__ __device__ Ix toIndex(DIM6 sh, DIM6 ix) +{ + DIM5 sh_ = { sh.a4, sh.a3, sh.a2, sh.a1, sh.a0 }; + DIM5 ix_ = { ix.a4, ix.a3, ix.a2, ix.a1, ix.a0 }; + return toIndex(sh_, ix_) * sh.a5 + ix.a5; +} + +static __inline__ __device__ Ix toIndex(DIM7 sh, DIM7 ix) +{ + DIM6 sh_ = { sh.a5, sh.a4, sh.a3, sh.a2, sh.a1, sh.a0 }; + DIM6 ix_ = { ix.a5, ix.a4, ix.a3, ix.a2, ix.a1, ix.a0 }; + return toIndex(sh_, ix_) * sh.a6 + ix.a6; +} + +static __inline__ __device__ Ix toIndex(DIM8 sh, DIM8 ix) +{ + DIM7 sh_ = { sh.a6, sh.a5, sh.a4, sh.a3, sh.a2, sh.a1, sh.a0 }; + DIM7 ix_ = { ix.a6, ix.a5, ix.a4, ix.a3, ix.a2, ix.a1, ix.a0 }; + return toIndex(sh_, ix_) * sh.a7 + ix.a7; +} + +static __inline__ __device__ Ix toIndex(DIM9 sh, DIM9 ix) +{ + DIM8 sh_ = { sh.a7, sh.a6, sh.a5, sh.a4, sh.a3, sh.a2, sh.a1, sh.a0 }; + DIM8 ix_ = { ix.a7, ix.a6, ix.a5, ix.a4, ix.a3, ix.a2, ix.a1, ix.a0 }; + return toIndex(sh_, ix_) * sh.a8 + ix.a8; +} + + +/* + * Inverse of 'toIndex' + */ +static __inline__ __device__ DIM1 fromIndex(DIM1 sh, Ix ix) +{ + return ix; +} + +static __inline__ __device__ DIM2 fromIndex(DIM2 sh, Ix ix) +{ + DIM1 sh_ = shape(sh.a0); + DIM1 ix_ = fromIndex(sh_, ix / sh.a1); + return shape(ix % sh.a1, ix_); +} + +static __inline__ __device__ DIM3 fromIndex(DIM3 sh, Ix ix) +{ + DIM2 sh_ = shape(sh.a1, sh.a0); + DIM2 ix_ = fromIndex(sh_, ix / sh.a2); + return shape(ix % sh.a2, ix_.a1, ix_.a0); +} + +static __inline__ __device__ DIM4 fromIndex(DIM4 sh, Ix ix) +{ + DIM3 sh_ = shape(sh.a2, sh.a1, sh.a0); + DIM3 ix_ = fromIndex(sh_, ix / sh.a3); + return shape(ix % sh.a3, ix_.a2, ix_.a1, ix_.a0); +} + +static __inline__ __device__ DIM5 fromIndex(DIM5 sh, Ix ix) +{ + DIM4 sh_ = shape(sh.a3, sh.a2, sh.a1, sh.a0); + DIM4 ix_ = fromIndex(sh_, ix / sh.a4); + return shape(ix % sh.a4, ix_.a3, ix_.a2, ix_.a1, ix_.a0); +} + +static __inline__ __device__ DIM6 fromIndex(DIM6 sh, Ix ix) +{ + DIM5 sh_ = shape(sh.a4, sh.a3, sh.a2, sh.a1, sh.a0); + DIM5 ix_ = fromIndex(sh_, ix / sh.a5); + return shape(ix % sh.a5, ix_.a4, ix_.a3, ix_.a2, ix_.a1, ix_.a0); +} + +static __inline__ __device__ DIM7 fromIndex(DIM7 sh, Ix ix) +{ + DIM6 sh_ = shape(sh.a5, sh.a4, sh.a3, sh.a2, sh.a1, sh.a0); + DIM6 ix_ = fromIndex(sh_, ix / sh.a6); + return shape(ix % sh.a6, ix_.a5, ix_.a4, ix_.a3, ix_.a2, ix_.a1, ix_.a0); +} + +static __inline__ __device__ DIM8 fromIndex(DIM8 sh, Ix ix) +{ + DIM7 sh_ = shape(sh.a6, sh.a5, sh.a4, sh.a3, sh.a2, sh.a1, sh.a0); + DIM7 ix_ = fromIndex(sh_, ix / sh.a7); + return shape(ix % sh.a7, ix_.a6, ix_.a5, ix_.a4, ix_.a3, ix_.a2, ix_.a1, ix_.a0); +} + +static __inline__ __device__ DIM9 fromIndex(DIM9 sh, Ix ix) +{ + DIM8 sh_ = shape(sh.a7, sh.a6, sh.a5, sh.a4, sh.a3, sh.a2, sh.a1, sh.a0); + DIM8 ix_ = fromIndex(sh_, ix / sh.a8); + return shape(ix % sh.a8, ix_.a7, ix_.a6, ix_.a5, ix_.a4, ix_.a3, ix_.a2, ix_.a1, ix_.a0); +} + + +/* + * Test for the magic index `ignore` + */ +static __inline__ __device__ int ignore(DIM1 ix) +{ + return ix == -1; +} + +static __inline__ __device__ int ignore(DIM2 ix) +{ + return ix.a0 == -1 && ix.a1 == -1; +} + +static __inline__ __device__ int ignore(DIM3 ix) +{ + return ix.a0 == -1 && ix.a1 == -1 && ix.a2 == -1; +} + +static __inline__ __device__ int ignore(DIM4 ix) +{ + return ix.a0 == -1 && ix.a1 == -1 && ix.a2 == -1 && ix.a3 == -1; +} + +static __inline__ __device__ int ignore(DIM5 ix) +{ + return ix.a0 == -1 && ix.a1 == -1 && ix.a2 == -1 && ix.a3 == -1 && ix.a4 == -1; +} + +static __inline__ __device__ int ignore(DIM6 ix) +{ + return ix.a0 == -1 && ix.a1 == -1 && ix.a2 == -1 && ix.a3 == -1 && ix.a4 == -1 && ix.a5 == -1; +} + +static __inline__ __device__ int ignore(DIM7 ix) +{ + return ix.a0 == -1 && ix.a1 == -1 && ix.a2 == -1 && ix.a3 == -1 && ix.a4 == -1 && ix.a5 == -1 && ix.a6 == -1; +} + +static __inline__ __device__ int ignore(DIM8 ix) +{ + return ix.a0 == -1 && ix.a1 == -1 && ix.a2 == -1 && ix.a3 == -1 && ix.a4 == -1 && ix.a5 == -1 && ix.a6 == -1 && ix.a7 == -1; +} + +static __inline__ __device__ int ignore(DIM9 ix) +{ + return ix.a0 == -1 && ix.a1 == -1 && ix.a2 == -1 && ix.a3 == -1 && ix.a4 == -1 && ix.a5 == -1 && ix.a6 == -1 && ix.a7 == -1 && ix.a8 == -1; +} + +#else + +int dim(Ix sh); +int size(Ix sh); +int shape(Ix sh); +int toIndex(Ix sh, Ix ix); +int fromIndex(Ix sh, Ix ix); + +#endif // __cplusplus +#endif // __ACCELERATE_CUDA_SHAPE_H__ + diff --git a/cubits/accelerate_cuda_texture.h b/cubits/accelerate_cuda_texture.h new file mode 100644 --- /dev/null +++ b/cubits/accelerate_cuda_texture.h @@ -0,0 +1,132 @@ +/* ----------------------------------------------------------------------------- + * + * Module : Texture + * Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell + * License : BSD3 + * + * Maintainer : Trevor L. McDonell + * Stability : experimental + * + * CUDA texture definitions and access functions are defined in terms of + * templates, and hence only available through the C++ interface. Expose some + * dummy wrappers to enable parsing with language-c. + * + * We don't have a definition for `Int' or `Word', since the bitwidth of the + * Haskell and C types may be different. + * + * NOTE ON 64-BIT TYPES + * The CUDA device uses little-endian arithmetic. We haven't accounted for the + * fact that this may be different on the host, for both initial data transfer + * and the unpacking below. + * + * ---------------------------------------------------------------------------*/ + +#ifndef __ACCELERATE_CUDA_TEXTURE_H__ +#define __ACCELERATE_CUDA_TEXTURE_H__ + +#include +#include + +#if defined(__cplusplus) && defined(__CUDACC__) + +typedef texture TexWord64; +typedef texture TexWord32; +typedef texture TexWord16; +typedef texture TexWord8; +typedef texture TexInt64; +typedef texture TexInt32; +typedef texture TexInt16; +typedef texture TexInt8; +typedef texture TexFloat; +typedef texture TexCChar; + +static __inline__ __device__ uint8_t indexArray(TexWord8 t, const int x) { return tex1Dfetch(t,x); } +static __inline__ __device__ uint16_t indexArray(TexWord16 t, const int x) { return tex1Dfetch(t,x); } +static __inline__ __device__ uint32_t indexArray(TexWord32 t, const int x) { return tex1Dfetch(t,x); } +static __inline__ __device__ uint64_t indexArray(TexWord64 t, const int x) +{ + union { uint2 x; uint64_t y; } v; + v.x = tex1Dfetch(t,x); + return v.y; +} + +static __inline__ __device__ int8_t indexArray(TexInt8 t, const int x) { return tex1Dfetch(t,x); } +static __inline__ __device__ int16_t indexArray(TexInt16 t, const int x) { return tex1Dfetch(t,x); } +static __inline__ __device__ int32_t indexArray(TexInt32 t, const int x) { return tex1Dfetch(t,x); } +static __inline__ __device__ int64_t indexArray(TexInt64 t, const int x) +{ + union { int2 x; int64_t y; } v; + v.x = tex1Dfetch(t,x); + return v.y; +} + +static __inline__ __device__ float indexArray(TexFloat t, const int x) { return tex1Dfetch(t,x); } +static __inline__ __device__ char indexArray(TexCChar t, const int x) { return tex1Dfetch(t,x); } + +#if defined(__LP64__) +typedef TexInt64 TexCLong; +typedef TexWord64 TexCULong; +#else +typedef TexInt32 TexCLong; +typedef TexWord32 TexCULong; +#endif + +/* + * Synonyms for C-types. NVCC will force Ints to be 32-bits. + */ +typedef TexInt8 TexCSChar; +typedef TexWord8 TexCUChar; +typedef TexInt16 TexCShort; +typedef TexWord16 TexCUShort; +typedef TexInt32 TexCInt; +typedef TexWord32 TexCUInt; +typedef TexInt64 TexCLLong; +typedef TexWord64 TexCULLong; +typedef TexFloat TexCFloat; + +/* + * Doubles, only available when compiled for Compute 1.3 and greater + */ +typedef texture TexDouble; +typedef TexDouble TexCDouble; + +#if !defined(CUDA_NO_SM_13_DOUBLE_INTRINSICS) +static __inline__ __device__ double indexDArray(TexDouble t, const int x) +{ + int2 v = tex1Dfetch(t,x); + return __hiloint2double(v.y,v.x); +} +#endif + +#else + +typedef void* TexWord64; +typedef void* TexWord32; +typedef void* TexWord16; +typedef void* TexWord8; +typedef void* TexInt64; +typedef void* TexInt32; +typedef void* TexInt16; +typedef void* TexInt8; +typedef void* TexCShort; +typedef void* TexCUShort; +typedef void* TexCInt; +typedef void* TexCUInt; +typedef void* TexCLong; +typedef void* TexCULong; +typedef void* TexCLLong; +typedef void* TexCULLong; +typedef void* TexFloat; +typedef void* TexDouble; +typedef void* TexCFloat; +typedef void* TexCDouble; +typedef void* TexCChar; +typedef void* TexCSChar; +typedef void* TexCUChar; + +void* indexArray(const void*, const int); +void* indexDArray(const void*, const int); + +#endif // defined(__cplusplus) && defined(__CUDACC__) +#endif // __ACCELERATE_CUDA_TEXTURE_H__ + diff --git a/cubits/accelerate_cuda_util.h b/cubits/accelerate_cuda_util.h new file mode 100644 --- /dev/null +++ b/cubits/accelerate_cuda_util.h @@ -0,0 +1,69 @@ +/* ----------------------------------------------------------------------------- + * + * Module : Util + * Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell + * License : BSD3 + * + * Maintainer : Trevor L. McDonell + * Stability : experimental + * + * ---------------------------------------------------------------------------*/ + +#ifndef __ACCELERATE_CUDA_UTIL_H__ +#define __ACCELERATE_CUDA_UTIL_H__ + +#include +#include + +/* + * Core assert function. Don't let this escape... + */ +#if defined(__CUDACC__) || !defined(__DEVICE_EMULATION__) +#define __assert(e, file, line) ((void)0) +#else +#define __assert(e, file, line) \ + ((void) fprintf (stderr, "%s:%u: failed assertion `%s'\n", file, line, e), abort()) +#endif + +/* + * Test the given expression, and abort the program if it evaluates to false. + * Only available in debug mode. + */ +#ifndef _DEBUG +#define assert(e) ((void)0) +#else +#define assert(e) \ + ((void) ((e) ? (void(0)) : __assert (#e, __FILE__, __LINE__))) +#endif + +/* + * Macro to insert __syncthreads() in device emulation mode + */ +#ifdef __DEVICE_EMULATION__ +#define __EMUSYNC __syncthreads() +#else +#define __EMUSYNC +#endif + +/* + * Check the return status of CUDA API calls, and abort with an appropriate + * error string on failure. + */ +#define CUDA_SAFE_CALL_NO_SYNC(call) \ + do { \ + cudaError err = call; \ + if(cudaSuccess != err) { \ + const char *str = cudaGetErrorString(err); \ + __assert(str, __FILE__, __LINE__); \ + } \ + } while (0) + +#define CUDA_SAFE_CALL(call) \ + do { \ + CUDA_SAFE_CALL_NO_SYNC(call); \ + CUDA_SAFE_CALL_NO_SYNC(cudaThreadSynchronize()); \ + } while (0) + +#undef __assert +#endif // __ACCELERATE_CUDA_UTIL_H__ + diff --git a/cubits/accelerate_cuda_utils.h b/cubits/accelerate_cuda_utils.h deleted file mode 100644 --- a/cubits/accelerate_cuda_utils.h +++ /dev/null @@ -1,136 +0,0 @@ -/* ----------------------------------------------------------------------------- - * - * Module : Utils - * Copyright : (c) [2009..2010] Trevor L. McDonell - * License : BSD - * - * ---------------------------------------------------------------------------*/ - -#ifndef __UTILS_H__ -#define __UTILS_H__ - -#include -#include - -/* - * Core assert function. Don't let this escape... - */ -#if defined(__CUDACC__) || !defined(__DEVICE_EMULATION__) -#define __assert(e, file, line) ((void)0) -#else -#define __assert(e, file, line) \ - ((void) fprintf (stderr, "%s:%u: failed assertion `%s'\n", file, line, e), abort()) -#endif - -/* - * Test the given expression, and abort the program if it evaluates to false. - * Only available in debug mode. - */ -#ifndef _DEBUG -#define assert(e) ((void)0) -#else -#define assert(e) \ - ((void) ((e) ? (void(0)) : __assert (#e, __FILE__, __LINE__))) -#endif - -/* - * Macro to insert __syncthreads() in device emulation mode - */ -#ifdef __DEVICE_EMULATION__ -#define __EMUSYNC __syncthreads() -#else -#define __EMUSYNC -#endif - -/* - * Check the return status of CUDA API calls, and abort with an appropriate - * error string on failure. - */ -#define CUDA_SAFE_CALL_NO_SYNC(call) \ - do { \ - cudaError err = call; \ - if(cudaSuccess != err) { \ - const char *str = cudaGetErrorString(err); \ - __assert(str, __FILE__, __LINE__); \ - } \ - } while (0) - -#define CUDA_SAFE_CALL(call) \ - do { \ - CUDA_SAFE_CALL_NO_SYNC(call); \ - CUDA_SAFE_CALL_NO_SYNC(cudaThreadSynchronize()); \ - } while (0) - - -#ifdef __cplusplus -extern "C" { -#endif - -/* - * Determine if the input is a power of two - */ -inline int -isPow2(unsigned int x) -{ - return ((x&(x-1)) == 0); -} - -/* - * Compute the next highest power of two - */ -inline unsigned int -ceilPow2(unsigned int x) -{ -#if 0 - --x; - x |= x >> 1; - x |= x >> 2; - x |= x >> 4; - x |= x >> 8; - x |= x >> 16; - return ++x; -#endif - - return (isPow2(x)) ? x : 1u << (int) ceil(log2((double)x)); -} - -/* - * Compute the next lowest power of two - */ -inline unsigned int -floorPow2(unsigned int x) -{ -#if 0 - float nf = (float) n; - return 1 << (((*(int*)&nf) >> 23) - 127); -#endif - - int exp; - frexp(x, &exp); - return 1 << (exp - 1); -} - -/* - * computes next highest multiple of f from x - */ -inline unsigned int -multiple(unsigned int x, unsigned int f) -{ - return ((x + (f-1)) / f); -} - -/* - * MS Excel-style CEIL() function. Rounds x up to nearest multiple of f - */ -inline unsigned int -ceiling(unsigned int x, unsigned int f) -{ - return multiple(x, f) * f; -} - -#undef __assert -#ifdef __cplusplus -} -#endif -#endif - diff --git a/cubits/backpermute.inl b/cubits/backpermute.inl --- a/cubits/backpermute.inl +++ b/cubits/backpermute.inl @@ -1,9 +1,12 @@ /* ----------------------------------------------------------------------------- * - * Module : Backpermute - * Copyright : (c) [2009..2010] Trevor L. McDonell - * License : BSD + * Kernel : Backpermute + * Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell + * License : BSD3 * + * Maintainer : Trevor L. McDonell + * Stability : experimental + * * Backwards permutation (gather) an array according to the permutation * function. The input `shape' is that of the output array. * @@ -22,8 +25,8 @@ const DimIn0 shIn0 ) { - Ix shapeSize = size(shOut); - const Ix gridSize = __umul24(blockDim.x, gridDim.x); + const Ix shapeSize = size(shOut); + const Ix gridSize = __umul24(blockDim.x, gridDim.x); for (Ix ix = __umul24(blockDim.x, blockIdx.x) + threadIdx.x; ix < shapeSize; ix += gridSize) { diff --git a/cubits/fold.inl b/cubits/fold.inl --- a/cubits/fold.inl +++ b/cubits/fold.inl @@ -1,21 +1,16 @@ /* ----------------------------------------------------------------------------- * - * Module : Fold - * Copyright : (c) [2009..2010] Trevor L. McDonell - * License : BSD + * Kernel : Fold + * Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell + * License : BSD3 * + * Maintainer : Trevor L. McDonell + * Stability : experimental + * * Reduce an array to a single value with a binary associative function * * ---------------------------------------------------------------------------*/ -/* - * We require block sizes are a power of two. This value gives maximum occupancy - * for both 1.x and 2.x class devices. - */ -#ifndef BLOCK_SIZE -#define BLOCK_SIZE 256 -#endif - #ifndef LENGTH_IS_POW_2 #define LENGTH_IS_POW_2 0 #endif @@ -35,7 +30,8 @@ const Ix shape ) { - extern volatile __shared__ TyOut s_data[]; + extern volatile __shared__ void* s_ptr[]; + ArrOut s_data = partition(s_ptr, blockDim.x); /* * Calculate first level of reduction reading into shared memory @@ -43,7 +39,7 @@ Ix i; TyOut sum = identity(); const Ix tid = threadIdx.x; - const Ix gridSize = BLOCK_SIZE * 2 * gridDim.x; + const Ix gridSize = blockDim.x * 2 * gridDim.x; /* * Reduce multiple elements per thread. The number is determined by the @@ -52,7 +48,7 @@ * * The loop stride of `gridSize' is used to maintain coalescing. */ - for (i = blockIdx.x * BLOCK_SIZE * 2 + tid; i < shape; i += gridSize) + for (i = blockIdx.x * blockDim.x * 2 + tid; i < shape; i += gridSize) { sum = apply(sum, get0(d_in0, i)); @@ -60,20 +56,20 @@ * Ensure we don't read out of bounds. This is optimised away if the * input length is a power of two */ - if (LENGTH_IS_POW_2 || i + BLOCK_SIZE < shape) - sum = apply(sum, get0(d_in0, i+BLOCK_SIZE)); + if (LENGTH_IS_POW_2 || i + blockDim.x < shape) + sum = apply(sum, get0(d_in0, i+blockDim.x)); } /* * Each thread puts its local sum into shared memory, then threads * cooperatively reduce the shared array to a single value. */ - s_data[tid] = sum; + set(s_data, tid, sum); __syncthreads(); - if (BLOCK_SIZE >= 512) { if (tid < 256) { s_data[tid] = sum = apply(sum, s_data[tid+256]); } __syncthreads(); } - if (BLOCK_SIZE >= 256) { if (tid < 128) { s_data[tid] = sum = apply(sum, s_data[tid+128]); } __syncthreads(); } - if (BLOCK_SIZE >= 128) { if (tid < 64) { s_data[tid] = sum = apply(sum, s_data[tid+ 64]); } __syncthreads(); } + if (blockDim.x >= 512) { if (tid < 256) { sum = apply(sum, get0(s_data, tid+256)); set(s_data, tid, sum); } __syncthreads(); } + if (blockDim.x >= 256) { if (tid < 128) { sum = apply(sum, get0(s_data, tid+128)); set(s_data, tid, sum); } __syncthreads(); } + if (blockDim.x >= 128) { if (tid < 64) { sum = apply(sum, get0(s_data, tid+ 64)); set(s_data, tid, sum); } __syncthreads(); } if (tid < 32) { @@ -81,12 +77,12 @@ * Use an extra warps worth of elements of shared memory, to let threads * index beyond the input data without using any branch instructions. */ - s_data[tid] = sum = apply(sum, s_data[tid + 32]); - s_data[tid] = sum = apply(sum, s_data[tid + 16]); - s_data[tid] = sum = apply(sum, s_data[tid + 8]); - s_data[tid] = sum = apply(sum, s_data[tid + 4]); - s_data[tid] = sum = apply(sum, s_data[tid + 2]); - s_data[tid] = sum = apply(sum, s_data[tid + 1]); + sum = apply(sum, get0(s_data, tid+32)); set(s_data, tid, sum); + sum = apply(sum, get0(s_data, tid+16)); set(s_data, tid, sum); + sum = apply(sum, get0(s_data, tid+ 8)); set(s_data, tid, sum); + sum = apply(sum, get0(s_data, tid+ 4)); set(s_data, tid, sum); + sum = apply(sum, get0(s_data, tid+ 2)); set(s_data, tid, sum); + sum = apply(sum, get0(s_data, tid+ 1)); } /* diff --git a/cubits/fold_segmented.inl b/cubits/fold_segmented.inl --- a/cubits/fold_segmented.inl +++ b/cubits/fold_segmented.inl @@ -1,32 +1,32 @@ /* ----------------------------------------------------------------------------- * - * Module : FoldSeg - * Copyright : (c) [2009..2010] Trevor L. McDonell, Rami G. Mukhtar - * License : BSD + * Kernel : FoldSeg + * Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell + * License : BSD3 * + * Maintainer : Trevor L. McDonell + * Stability : experimental + * * Reduction an array to a single value for each segment, using a binary * associative function * * ---------------------------------------------------------------------------*/ -#define WARP_SIZE 32 - /* * Cooperatively reduce an array to a single value. The computation requires an * extra half-warp worth of elements of shared memory per block, to let threads * index beyond the input data without using any branch instructions. */ -template static __inline__ __device__ -T reduce_warp(volatile T* s_data, T sum) +TyOut reduce_warp(ArrOut s_data, TyOut sum) { - s_data[threadIdx.x] = sum; - s_data[threadIdx.x] = sum = apply(sum, s_data[threadIdx.x + 16]); - s_data[threadIdx.x] = sum = apply(sum, s_data[threadIdx.x + 8]); - s_data[threadIdx.x] = sum = apply(sum, s_data[threadIdx.x + 4]); - s_data[threadIdx.x] = sum = apply(sum, s_data[threadIdx.x + 2]); - s_data[threadIdx.x] = sum = apply(sum, s_data[threadIdx.x + 1]); + set(s_data, threadIdx.x, sum); + sum = apply(sum, get0(s_data, threadIdx.x + 16)); set(s_data, threadIdx.x, sum); + sum = apply(sum, get0(s_data, threadIdx.x + 8)); set(s_data, threadIdx.x, sum); + sum = apply(sum, get0(s_data, threadIdx.x + 4)); set(s_data, threadIdx.x, sum); + sum = apply(sum, get0(s_data, threadIdx.x + 2)); set(s_data, threadIdx.x, sum); + sum = apply(sum, get0(s_data, threadIdx.x + 1)); return sum; } @@ -63,18 +63,18 @@ const Ix length // of input array d_in0 ) { - const Ix vectors_per_block = blockDim.x / WARP_SIZE; + const Ix vectors_per_block = blockDim.x / warpSize; const Ix num_vectors = vectors_per_block * gridDim.x; const Ix thread_id = blockDim.x * blockIdx.x + threadIdx.x; - const Ix vector_id = thread_id / WARP_SIZE; - const Ix thread_lane = threadIdx.x & (WARP_SIZE - 1); - const Ix vector_lane = threadIdx.x / WARP_SIZE; + const Ix vector_id = thread_id / warpSize; + const Ix thread_lane = threadIdx.x & (warpSize - 1); + const Ix vector_lane = threadIdx.x / warpSize; /* * Manually partition (dynamically-allocated) shared memory */ - extern __shared__ Ix s_ptrs[][2]; - TyOut* s_data = (TyOut*) &s_ptrs[vectors_per_block][2]; + extern volatile __shared__ Ix s_ptrs[][2]; + ArrOut s_data = partition((void*) &s_ptrs[vectors_per_block][2], blockDim.x); for (Ix seg = vector_id; seg < num_segments; seg += num_vectors) { @@ -105,7 +105,7 @@ * local sum. This is then reduced cooperatively in shared memory. */ TyOut sum = identity(); - for (Ix i = start + thread_lane; i < end; i += WARP_SIZE) + for (Ix i = start + thread_lane; i < end; i += warpSize) sum = apply(sum, get0(d_in0, i)); sum = reduce_warp(s_data, sum); diff --git a/cubits/map.inl b/cubits/map.inl --- a/cubits/map.inl +++ b/cubits/map.inl @@ -1,8 +1,11 @@ /* ----------------------------------------------------------------------------- * - * Module : Map - * Copyright : (c) 2010 Trevor L. McDonell - * License : BSD + * Kernel : Map + * Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell + * License : BSD3 + * + * Maintainer : Trevor L. McDonell + * Stability : experimental * * Apply the function to each element of the array. Each thread processes * multiple elements, striding the array by the grid size. diff --git a/cubits/permute.inl b/cubits/permute.inl --- a/cubits/permute.inl +++ b/cubits/permute.inl @@ -1,9 +1,12 @@ /* ----------------------------------------------------------------------------- * - * Module : Permute - * Copyright : (c) [2009..2010] Trevor L. McDonell - * License : BSD + * Kernel : Permute + * Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell + * License : BSD3 * + * Maintainer : Trevor L. McDonell + * Stability : experimental + * * Forward permutation, characterised by a function that determines for each * element in the source array where it should go in the target. The output * array should be initialised with a default value, as the permutation may be @@ -24,8 +27,8 @@ const DimIn0 shIn0 ) { - Ix shapeSize = size(shIn0); - const Ix gridSize = __umul24(blockDim.x, gridDim.x); + const Ix shapeSize = size(shIn0); + const Ix gridSize = __umul24(blockDim.x, gridDim.x); for (Ix ix = __umul24(blockDim.x, blockIdx.x) + threadIdx.x; ix < shapeSize; ix += gridSize) { diff --git a/cubits/replicate.inl b/cubits/replicate.inl --- a/cubits/replicate.inl +++ b/cubits/replicate.inl @@ -1,9 +1,12 @@ /* ----------------------------------------------------------------------------- * - * Module : Replicate - * Copyright : (c) [2009..2010] Trevor L. McDonell - * License : BSD + * Kernel : Replicate + * Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell + * License : BSD3 * + * Maintainer : Trevor L. McDonell + * Stability : experimental + * * ---------------------------------------------------------------------------*/ extern "C" @@ -16,8 +19,8 @@ const SliceDim sliceDim ) { - Ix shapeSize = size(sliceDim); - const Ix gridSize = __umul24(blockDim.x, gridDim.x); + const Ix shapeSize = size(sliceDim); + const Ix gridSize = __umul24(blockDim.x, gridDim.x); for (Ix ix = __umul24(blockDim.x, blockIdx.x) + threadIdx.x; ix < shapeSize; ix += gridSize) { diff --git a/cubits/slice.inl b/cubits/slice.inl --- a/cubits/slice.inl +++ b/cubits/slice.inl @@ -1,9 +1,12 @@ /* ----------------------------------------------------------------------------- * - * Module : Slice - * Copyright : (c) [2009..2010] Trevor L. McDonell - * License : BSD + * Kernel : Slice + * Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell + * License : BSD3 * + * Maintainer : Trevor L. McDonell + * Stability : experimental + * * ---------------------------------------------------------------------------*/ extern "C" @@ -17,8 +20,8 @@ const SliceDim sliceDim ) { - Ix shapeSize = size(slice); - const Ix gridSize = __umul24(blockDim.x, gridDim.x); + const Ix shapeSize = size(slice); + const Ix gridSize = __umul24(blockDim.x, gridDim.x); for (Ix ix = __umul24(blockDim.x, blockIdx.x) + threadIdx.x; ix < shapeSize; ix += gridSize) { diff --git a/cubits/zipWith.inl b/cubits/zipWith.inl --- a/cubits/zipWith.inl +++ b/cubits/zipWith.inl @@ -1,9 +1,12 @@ /* ----------------------------------------------------------------------------- * - * Module : ZipWith - * Copyright : (c) 2010 Trevor L. McDonell - * License : BSD + * Kernel : ZipWith + * Copyright : [2008..2010] Manuel M T Chakravarty, Gabriele Keller, Sean Lee, Trevor L. McDonell + * License : BSD3 * + * Maintainer : Trevor L. McDonell + * Stability : experimental + * * Combine two arrays using the given binary operator. Each thread processes * multiple elements, striding the array by the grid size. * @@ -21,8 +24,8 @@ const DimIn0 shIn0 ) { - Ix shapeSize = size(shOut); - const Ix gridSize = __umul24(blockDim.x, gridDim.x); + const Ix shapeSize = size(shOut); + const Ix gridSize = __umul24(blockDim.x, gridDim.x); for (Ix ix = __umul24(blockDim.x, blockIdx.x) + threadIdx.x; ix < shapeSize; ix += gridSize) { diff --git a/examples/simple/src/Main.hs b/examples/simple/src/Main.hs --- a/examples/simple/src/Main.hs +++ b/examples/simple/src/Main.hs @@ -23,7 +23,7 @@ instance (Ix dim, IArray UArray e) => NFData (UArray dim e) where - rnf a = a ! (head (indices a)) `seq` () + rnf a = a ! head (indices a) `seq` () -- Generate a benchmark test iff the reference and accelerate tests succeed. @@ -37,13 +37,13 @@ -> IO Benchmark benchmark name sim ref acc = do - putStr "Interpreter : " ; v1 <- validate sim (ref ()) (Acc.toIArray $ Interp.run (acc ())) - putStr "CUDA : " ; v2 <- validate sim (ref ()) . Acc.toIArray =<< (CUDA.run (acc ())) + putStr "Interpreter : " ; v1 <- validate sim (ref ()) (Acc.toIArray $ Interp.run (acc ())) + putStr "CUDA : " ; v2 <- validate sim (ref ()) (Acc.toIArray $ CUDA.run (acc ())) if not (v1 && v2) then return $ bgroup "" [] else return $ bgroup name [ bench "ref" $ nf ref () - , bench "cuda" $ (CUDA.run . acc) () + , bench "cuda" $ whnf (CUDA.run . acc) () ] @@ -58,6 +58,7 @@ ys' <- convertVector ys benchmark "dotp" similar (run_ref xs ys) (run_acc xs' ys') where + {-# NOINLINE run_ref #-} run_ref x y () = dotp_ref x y run_acc x y () = dotp x y @@ -72,6 +73,7 @@ alpha <- uniform gen benchmark "saxpy" similar (run_ref alpha xs ys) (run_acc alpha xs' ys') where + {-# NOINLINE run_ref #-} run_ref alpha x y () = saxpy_ref alpha x y run_acc alpha x y () = saxpy alpha x y @@ -83,6 +85,7 @@ xs' <- convertVector xs benchmark "filter" similar (run_ref xs) (run_acc xs') where + {-# NOINLINE run_ref #-} run_ref x () = filter_ref (< 0.5) x run_acc x () = filter (Acc.<* 0.5) x @@ -102,6 +105,7 @@ in evaluate (v `Acc.indexArray` 0) >> return v benchmark "smvm" similar (run_ref segd inds vals vec) (run_acc segd' mat' vec') where + {-# NOINLINE run_ref #-} run_ref d i x v () = smvm_ref (d, (i,x)) v run_acc d x v () = smvm (d,x) v @@ -117,7 +121,7 @@ defaultMain =<< sequence [ test_dotp gen 100000 , test_saxpy gen 100000 - , test_filter gen 1800 + , test_filter gen 10000 , test_smvm gen (0,42) (2400,400) ]