accelerate-cublas 0.0 → 0.1
raw patch · 6 files changed
+283/−39 lines, 6 filesdep +hmatrixdep +pooled-iodep +timeitdep ~acceleratedep ~accelerate-arithmeticdep ~accelerate-cuda
Dependencies added: hmatrix, pooled-io, timeit
Dependency ranges changed: accelerate, accelerate-arithmetic, accelerate-cuda, accelerate-utility, base, cublas, cuda, random, utility-ht, vector
Files
- accelerate-cublas.cabal +33/−13
- benchmark/NewtonInverse.hs +225/−0
- example/Main.hs +5/−7
- src/Data/Array/Accelerate/CUBLAS/Level2/Batched.hs +3/−5
- src/Data/Array/Accelerate/CUBLAS/Level3/Batched.hs +6/−7
- src/Data/Array/Accelerate/CUBLAS/Level3/Batched/Foreign.hs +11/−7
accelerate-cublas.cabal view
@@ -1,10 +1,10 @@ Name: accelerate-cublas-Version: 0.0+Version: 0.1 License: BSD3 License-File: LICENSE Author: Henning Thielemann <haskell@henning-thielemann.de> Maintainer: Henning Thielemann <haskell@henning-thielemann.de>-Homepage: http://code.haskell.org/~thielema/accelerate-cublas/+Homepage: http://hub.darcs.net/thielema/accelerate-cublas/ Category: Math Synopsis: Basic Linear Algebra using native CUBLAS library Description:@@ -19,26 +19,26 @@ default: False Source-Repository this- Tag: 0.0+ Tag: 0.1 Type: darcs- Location: http://code.haskell.org/~thielema/accelerate-cublas/+ Location: http://hub.darcs.net/thielema/accelerate-cublas/ Source-Repository head Type: darcs- Location: http://code.haskell.org/~thielema/accelerate-cublas/+ Location: http://hub.darcs.net/thielema/accelerate-cublas/ Library Build-Depends:- accelerate-arithmetic >=0.0.1 && <0.1,+ accelerate-arithmetic >=0.1 && <0.2, accelerate-utility >=0.1 && <0.2,- accelerate-cuda >=0.15 && <0.16,+ accelerate-cuda >=0.16 && <0.17, accelerate-io >=0.15 && <0.16, accelerate >=0.15 && <0.16, cublas >=0.2.0.2 && <0.3,- cuda >=0.5 && <0.7,- vector >=0.10.11 && <0.11,+ cuda >=0.5 && <0.8,+ vector >=0.10.11 && <0.12, utility-ht >=0.0.8 && <0.1,- base >=4.5 && <4.8+ base >=4.5 && <5 GHC-Options: -Wall -fwarn-missing-import-lists Hs-Source-Dirs: src@@ -60,6 +60,26 @@ accelerate-arithmetic, accelerate-utility, accelerate,- cublas,- random >=1.0 && <1.1,- base >=4.5 && <4.8+ random >=1.0 && <1.2,+ base++Benchmark newton-inverse+ Type: exitcode-stdio-1.0+ Main-Is: NewtonInverse.hs+ Hs-Source-Dirs: benchmark+ Default-Language: Haskell98+ GHC-Options: -Wall -threaded+ GHC-Prof-Options: -fprof-auto+ Build-Depends:+ accelerate-cublas,+ accelerate-arithmetic,+ accelerate-utility,+ accelerate-cuda,+ cuda >=0.5 && <0.8,+ accelerate,+ pooled-io >=0.0 && <0.1,+ timeit >=1.0 && <1.1,+ hmatrix >=0.15.2 && <0.16,+ random >=1.0.1 && <1.2,+ utility-ht,+ base
+ benchmark/NewtonInverse.hs view
@@ -0,0 +1,225 @@+{-# LANGUAGE FlexibleContexts #-}+module Main where++import qualified Data.Array.Accelerate.CUBLAS.Level3.Batched as Batched++import qualified Data.Array.Accelerate.LinearAlgebra as ALinAlg+import qualified Data.Array.Accelerate.Utility.Loop as Loop+import qualified Data.Array.Accelerate.CUDA as CUDA+import qualified Data.Array.Accelerate as A+import Data.Array.Accelerate (All(All), Z(Z), (:.)((:.)))++import qualified Control.Concurrent.PooledIO.Independent as Pooled++import qualified Data.Packed.Matrix as Matrix+import qualified Data.Packed.Vector as Vector+import qualified Numeric.Container as Container+import qualified Numeric.LinearAlgebra.Algorithms as HMLinAlg++import Numeric.Container (Container, (<>))+import Data.Packed.Matrix (Matrix)+import Data.Packed.Vector (Vector)++import qualified System.Random as Rnd+import System.TimeIt (timeIt)++import Text.Printf (printf)++import qualified Data.List.HT as ListHT+import Data.Function.HT (nest)+import Data.Tuple.HT (mapPair)++++newtonInverseStep ::+ (Num a, Container Vector a, Container.Product a) =>+ Matrix a -> Matrix a -> Matrix a+newtonInverseStep a x =+ Container.sub+ (Container.scale 2 x)+ (x <> a <> x)++newtonInverse ::+ (Num a, Container Vector a, Container.Product a) =>+ Int -> Matrix a -> Matrix a -> Matrix a+newtonInverse count start a =+ nest count (newtonInverseStep a) start+++newtonInverseCUBLASStep, newtonInverseCUBLASStepMul ::+ (A.Shape ix, A.Slice ix, Eq ix, Batched.Element a, A.IsNum a, A.Elt a) =>+ Batched.Handle ->+ ALinAlg.Matrix ix a ->+ ALinAlg.Matrix ix a ->+ ALinAlg.Matrix ix a+newtonInverseCUBLASStep h a x =+ Batched.mac h (-1) x (Batched.mul h 1 a x) 2 x++newtonInverseCUBLASStepMul h a x =+ A.zipWith (-) (A.map (2*) x) $+ Batched.mul h 1 x $ Batched.mul h 1 a x++newtonInverseCUBLAS ::+ (A.Shape ix, A.Slice ix, Eq ix, Batched.Element a, A.IsNum a, A.Elt a) =>+ Batched.Handle ->+ A.Exp Int ->+ ALinAlg.Matrix ix a ->+ ALinAlg.Matrix ix a ->+ ALinAlg.Matrix ix a+newtonInverseCUBLAS h n seed a =+ Loop.nest n (newtonInverseCUBLASStep h a) seed+++randomMatrixInv :: Int -> (Matrix Double, Matrix Double)+randomMatrixInv size =+ let x =+ Matrix.fromLists $ take size $ ListHT.sliceVertical size $+ Rnd.randomRs (-1,1::Double) $ Rnd.mkStdGen 42+ in (x, HMLinAlg.inv x)+++parallel :: [a] -> (Int -> a -> IO ()) -> IO ()+parallel xs f = Pooled.run $ zipWith f [0 ..] xs++disturbedMatrices ::+ (Container Vector a) =>+ Matrix a -> [a] -> [Matrix a]+disturbedMatrices x yelems =+ let size = Matrix.rows x+ in map (Container.add x . Matrix.fromLists) $+ ListHT.sliceVertical size $+ ListHT.sliceVertical size $+ yelems++mainHMatrixDirect ::+ (Show a, Container Vector a, HMLinAlg.Field a) =>+ String ->+ Int -> (Matrix a, Matrix a) -> [a] -> IO ()+mainHMatrixDirect typ numberOfMatrices (x, _xinv) yelems = do+ let yinvs = map HMLinAlg.inv $ disturbedMatrices x yelems+ putStrLn $ "hmatrix-direct-" ++ typ+ timeIt $ parallel (take numberOfMatrices yinvs) $ \ n y ->+ writeFile (printf "/tmp/hmatrix-direct-%s%03d.txt" typ n) $ show y++{-+das wird durch mehrere Kerne gar nicht schneller,+vielleicht st�rt das Serialisieren der Zufallszahlen?+-}+mainHMatrix ::+ (Show a, Container Vector a, Container.Product a) =>+ String ->+ Int -> Int -> (Matrix a, Matrix a) -> [a] -> IO ()+mainHMatrix typ numberOfMatrices newtonIts (x, xinv) yelems = do+ let yinvs = map (newtonInverse newtonIts xinv) $ disturbedMatrices x yelems+ putStrLn $ "hmatrix-" ++ typ+ timeIt $ parallel (take numberOfMatrices yinvs) $ \ n y ->+ writeFile (printf "/tmp/hmatrix-%s%03d.txt" typ n) $ show y+++mainCUDA ::+ (A.Elt a, A.IsNum a, Container.Element a) =>+ String ->+ Int -> Int -> (Matrix a, Matrix a) -> [a] -> IO ()+mainCUDA typ numberOfMatrices newtonIts (xm, xinvm) yelems = do+ let size = Matrix.rows xm+ matrixAccFromHM =+ A.fromList (Z :. size :. size) .+ Vector.toList . Matrix.flatten+ xarr = matrixAccFromHM xm+ xinvarr = matrixAccFromHM xinvm++ let ysarr =+ A.fromList (Z :. numberOfMatrices :. size :. size) yelems+ rep = A.replicate (A.lift $ Z :. numberOfMatrices :. All :. All)+ yinvs =+ CUDA.run1+ (\args ->+ case A.unlift args of+ (x, xinv, ys) ->+ ALinAlg.newtonInverse (A.constant newtonIts) (rep xinv) $+ A.zipWith (+) ys (rep x))+ (xarr, xinvarr, ysarr)++ putStrLn $ "cuda-" ++ typ+ timeIt $ writeFile ("/tmp/cuda-"++typ++".txt") $ show yinvs+++mainCUBLASDirect ::+ (Batched.Element a, Container.Element a, A.IsNum a, A.Elt a) =>+ String ->+ Int -> (Matrix a, Matrix a) -> [a] -> IO ()+mainCUBLASDirect typ numberOfMatrices (xm, _xinvm) yelems = do+ let size = Matrix.rows xm+ matrixAccFromHM =+ A.fromList (Z :. size :. size) .+ Vector.toList . Matrix.flatten+ xarr = matrixAccFromHM xm++ handle <- Batched.create+ let ysarr =+ A.fromList (Z :. numberOfMatrices :. size :. size) yelems+ rep = A.replicate (A.lift $ Z :. numberOfMatrices :. All :. All)+ yinvs =+ CUDA.run1+ (\args ->+ case A.unlift args of+ (x, ys) ->+ fst $ Batched.inv handle $ A.zipWith (+) ys (rep x))+ (xarr, ysarr)++ putStrLn $ "cublas-direct-" ++ typ+ timeIt $ writeFile ("/tmp/cublas-direct-"++typ++".txt") $ show yinvs+++mainCUBLAS ::+ (Batched.Element a, Container.Element a, A.IsNum a, A.Elt a) =>+ String ->+ Int -> Int -> (Matrix a, Matrix a) -> [a] -> IO ()+mainCUBLAS typ numberOfMatrices newtonIts (xm, xinvm) yelems = do+ let size = Matrix.rows xm+ matrixAccFromHM =+ A.fromList (Z :. size :. size) .+ Vector.toList . Matrix.flatten+ xarr = matrixAccFromHM xm+ xinvarr = matrixAccFromHM xinvm++ handle <- Batched.create+ let ysarr =+ A.fromList (Z :. numberOfMatrices :. size :. size) yelems+ rep = A.replicate (A.lift $ Z :. numberOfMatrices :. All :. All)+ yinvs =+ CUDA.run1+ (\args ->+ case A.unlift args of+ (x, xinv, ys) ->+ newtonInverseCUBLAS handle (A.constant newtonIts) (rep xinv) $+ A.zipWith (+) ys (rep x))+ (xarr, xinvarr, ysarr)++ putStrLn $ "cublas-" ++ typ+ timeIt $ writeFile ("/tmp/cublas-"++typ++".txt") $ show yinvs+++main :: IO ()+main = do+ let n = 96+ let sz = 50+ let its = 20+ let xmsDouble = randomMatrixInv sz+ ysDouble = Rnd.randomRs (-0.01,0.01) $ Rnd.mkStdGen 23+ let xmsFloat =+ mapPair+ (Container.cmap realToFrac, Container.cmap realToFrac)+ xmsDouble+ ysFloat :: [Float]+ ysFloat = map realToFrac ysDouble+ mainHMatrixDirect "double" n xmsDouble ysDouble+ -- mainHMatrixDirect "float" n xmsFloat ysFloat+ mainHMatrix "double" n its xmsDouble ysDouble+ mainHMatrix "float" n its xmsFloat ysFloat+ mainCUBLASDirect "double" n xmsDouble ysDouble+ mainCUBLASDirect "float" n xmsFloat ysFloat+ mainCUBLAS "double" n its xmsDouble ysDouble+ mainCUBLAS "float" n its xmsFloat ysFloat+ mainCUDA "double" n its xmsDouble ysDouble+ mainCUDA "float" n its xmsFloat ysFloat
example/Main.hs view
@@ -2,7 +2,7 @@ import qualified Data.Array.Accelerate.CUBLAS.Level3.Batched as Batched -import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as ALinAlg+import qualified Data.Array.Accelerate.LinearAlgebra as ALinAlg import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp import Data.Array.Accelerate.Utility.Lift.Exp (expr)@@ -11,8 +11,6 @@ import qualified Data.Array.Accelerate.CUDA as AC import qualified Data.Array.Accelerate as A -import qualified Foreign.CUDA.Cublas as Cublas- import System.Random (randomRs, mkStdGen) @@ -25,7 +23,7 @@ in (A.generate (A.constant $ Z :. numMats :. 3 :. 4) f, A.generate (A.constant $ Z :. numMats :. 4 :. 2) f) -mainMul :: Cublas.Handle -> IO ()+mainMul :: Batched.Handle -> IO () mainMul handle = do print factorMatrices print $ AC.run $@@ -78,7 +76,7 @@ let ix = A.index2 i j in k==*0 ? (x!ix, y!ix) -mainLU :: Cublas.Handle -> IO ()+mainLU :: Batched.Handle -> IO () mainLU handle = do print luMatrices let mat =@@ -90,7 +88,7 @@ Batched.mul handle 1 mat $ A.replicate (A.lift $ Z :. (2::Int) :. A.All :. A.All) rhsMatrix -mainInv :: Cublas.Handle -> IO ()+mainInv :: Batched.Handle -> IO () mainInv handle = do let dim = 4 mat =@@ -103,7 +101,7 @@ main :: IO () main = do- handle <- Cublas.create+ handle <- Batched.create mainMul handle mainLU handle mainInv handle
src/Data/Array/Accelerate/CUBLAS/Level2/Batched.hs view
@@ -1,8 +1,6 @@-{-# LANGUAGE TypeOperators #-}-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE ConstraintKinds #-}-{-# LANGUAGE FlexibleContexts #-} module Data.Array.Accelerate.CUBLAS.Level2.Batched (+ Cublas.Handle,+ Cublas.create, Level3.Element, mul, mac,@@ -11,7 +9,7 @@ import qualified Data.Array.Accelerate.CUBLAS.Level3.Batched as Level3 import Data.Array.Accelerate.CUBLAS.Level3.Batched (Element) -import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as ALinAlg+import qualified Data.Array.Accelerate.LinearAlgebra as ALinAlg import qualified Data.Array.Accelerate as A import Data.Array.Accelerate (Exp)
src/Data/Array/Accelerate/CUBLAS/Level3/Batched.hs view
@@ -1,8 +1,7 @@ {-# LANGUAGE TypeOperators #-}-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE ConstraintKinds #-}-{-# LANGUAGE FlexibleContexts #-} module Data.Array.Accelerate.CUBLAS.Level3.Batched (+ Cublas.Handle,+ Cublas.create, Element, mul, mac,@@ -18,7 +17,7 @@ import qualified Data.Array.Accelerate.CUBLAS.Level3.Batched.Foreign as Foreign import Data.Array.Accelerate.CUBLAS.Level3.Batched.Foreign (Element, Vector) -import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as ALinAlg+import qualified Data.Array.Accelerate.LinearAlgebra as ALinAlg import qualified Data.Array.Accelerate.Utility.Sliced1 as Sliced1 import qualified Data.Array.Accelerate.Utility.Sliced as Sliced@@ -54,7 +53,7 @@ ALinAlg.Matrix ix a mul handle alpha a b = A.foreignAcc- (AF.CUDAForeignAcc "mul" $ uncurry3 $ Foreign.mul handle)+ (AF.CUDAForeignAcc "mul" $ const $ uncurry3 $ Foreign.mul handle) (Acc.modify (expr,acc,acc) $ \(alpha0, a0, b0) -> A.map (alpha0 *) $ ALinAlg.multiplyMatrixMatrix a0 b0)@@ -69,7 +68,7 @@ ALinAlg.Matrix ix a mac handle alpha a b beta c = A.foreignAcc- (AF.CUDAForeignAcc "mac" $+ (AF.CUDAForeignAcc "mac" $ const $ \((alpha0, a0, b0), (beta0, c0)) -> Foreign.mac handle alpha0 a0 b0 beta0 c0) (Acc.modify ((expr,acc,acc),(expr,acc)) $@@ -262,5 +261,5 @@ String -> (acc -> AF.CIO res) -> A.Acc acc -> A.Acc res cudaAcc name f = A.foreignAcc- (AF.CUDAForeignAcc name f)+ (AF.CUDAForeignAcc name $ const f) (error $ name ++ ": requires CUDA backend")
src/Data/Array/Accelerate/CUBLAS/Level3/Batched/Foreign.hs view
@@ -1,6 +1,5 @@ {-# LANGUAGE TypeOperators #-} {-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE ConstraintKinds #-} {-# LANGUAGE FlexibleContexts #-} module Data.Array.Accelerate.CUBLAS.Level3.Batched.Foreign where @@ -191,12 +190,17 @@ return (b, info) -type Element a =- (AF.DevicePtrs (EltRepr a) ~ ((), DevicePtr a),- Fractional (StorableOf a),- Cublas.Cublas (StorableOf a),- Storable (StorableOf a),- Real a)+class+ (AF.DevicePtrs (EltRepr a) ~ ((), DevicePtr a),+ Fractional (StorableOf a),+ Cublas.Cublas (StorableOf a),+ Storable (StorableOf a),+ Real a) =>+ Element a where++instance Element Float where+instance Element Double where+ type family StorableOf float type instance StorableOf Float = CFloat