accelerate-arithmetic (empty) → 0.0
raw patch · 13 files changed
+1235/−0 lines, 13 filesdep +QuickCheckdep +acceleratedep +accelerate-arithmeticsetup-changed
Dependencies added: QuickCheck, accelerate, accelerate-arithmetic, accelerate-cuda, accelerate-utility, base, cublas, cuda, hmatrix, pooled-io, random, timeit, utility-ht
Files
- LICENSE +27/−0
- Setup.lhs +3/−0
- accelerate-arithmetic.cabal +82/−0
- benchmark/CUBLASBatched.hs +259/−0
- benchmark/NewtonInverse.hs +224/−0
- src/Data/Array/Accelerate/Arithmetic/Example.hs +49/−0
- src/Data/Array/Accelerate/Arithmetic/Interpolation.hs +97/−0
- src/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs +218/−0
- src/Data/Array/Accelerate/Arithmetic/Sparse.hs +115/−0
- test/Test.hs +20/−0
- test/Test/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs +41/−0
- test/Test/Data/Array/Accelerate/Arithmetic/Sparse.hs +52/−0
- test/Test/Data/Array/Accelerate/Arithmetic/Utility.hs +48/−0
+ LICENSE view
@@ -0,0 +1,27 @@+Copyright (c) Henning Thielemann 2014++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions+are met:+1. Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.+2. Redistributions in binary form must reproduce the above copyright+ notice, this list of conditions and the following disclaimer in the+ documentation and/or other materials provided with the distribution.+3. Neither the name of the author nor the names of his contributors+ may be used to endorse or promote products derived from this software+ without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE+ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS+OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)+HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT+LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY+OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF+SUCH DAMAGE.
+ Setup.lhs view
@@ -0,0 +1,3 @@+#! /usr/bin/env runhaskell+> import Distribution.Simple+> main = defaultMain
+ accelerate-arithmetic.cabal view
@@ -0,0 +1,82 @@+Name: accelerate-arithmetic+Version: 0.0+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-arithmetic/+Category: Math+Synopsis: Linear algebra and interpolation using the Accelerate framework+Description:+ Linear algebra and interpolation via the @accelerate@ package.+ This can be used for computations on GPUs+ but it does not contain processor optimizations+ or optimizations for CUDA.+Tested-With: GHC==7.8.2+Cabal-Version: >=1.14+Build-Type: Simple++Source-Repository this+ Tag: 0.0+ Type: darcs+ Location: http://code.haskell.org/~thielema/accelerate-arithmetic/++Source-Repository head+ Type: darcs+ Location: http://code.haskell.org/~thielema/accelerate-arithmetic/++Library+ Build-Depends:+ accelerate-utility >=0.0 && <0.1,+ accelerate >=0.15 && <0.16,+ utility-ht >=0.0.8 && <0.1,+ QuickCheck >=2.4 && <2.8,+ base >=4.5 && <4.8++ GHC-Options: -Wall -fwarn-missing-import-lists+ Hs-Source-Dirs: src+ Default-Language: Haskell98+ Exposed-Modules:+ Data.Array.Accelerate.Arithmetic.LinearAlgebra+ Data.Array.Accelerate.Arithmetic.Sparse+ Data.Array.Accelerate.Arithmetic.Interpolation+ Other-Modules:+ Data.Array.Accelerate.Arithmetic.Example++Test-Suite test+ Type: exitcode-stdio-1.0+ Main-Is: Test.hs+ GHC-Options: -Wall -fwarn-missing-import-lists+ Hs-Source-Dirs: test+ Default-Language: Haskell98+ Build-Depends:+ accelerate-arithmetic,+ accelerate,+ QuickCheck,+ base+ Other-Modules:+ Test.Data.Array.Accelerate.Arithmetic.LinearAlgebra+ Test.Data.Array.Accelerate.Arithmetic.Sparse+ Test.Data.Array.Accelerate.Arithmetic.Utility++Benchmark newton-inverse+ Type: exitcode-stdio-1.0+ Main-Is: NewtonInverse.hs+ Hs-Source-Dirs: benchmark+ Other-Modules: CUBLASBatched+ Default-Language: Haskell98+ GHC-Options: -Wall -threaded+ GHC-Prof-Options: -fprof-auto+ Build-Depends:+ accelerate-arithmetic,+ accelerate-utility,+ accelerate-cuda >=0.15 && <0.16,+ cublas >=0.2.0.2 && <0.3,+ cuda >=0.5 && <0.7,+ accelerate,+ pooled-io >=0.0 && <0.1,+ timeit >=1.0 && <1.1,+ hmatrix >=0.15.2 && <0.16,+ random >=1.0.1 && <1.1,+ utility-ht,+ base
+ benchmark/CUBLASBatched.hs view
@@ -0,0 +1,259 @@+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE FlexibleContexts #-}+module CUBLASBatched where++import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as ALinAlg++import Data.Array.Accelerate.Array.Sugar (EltRepr)+import Data.Array.Accelerate (Array, DIM3, Acc, Z (..), (:.) (..), Exp)+import qualified Data.Array.Accelerate.CUDA.Foreign as AF+import qualified Data.Array.Accelerate.CUDA as AC+import qualified Data.Array.Accelerate as A++import qualified Foreign.CUDA.Cublas as Cublas+import Foreign.CUDA.Ptr (DevicePtr, castDevPtr, advanceDevPtr)++import Foreign.C.Types (CFloat, CDouble)+import Foreign.Storable (Storable)++import Data.Tuple.HT (uncurry3)+++type Matrix ix = Array (ix :. Int :. Int)+type Vector ix = Array (ix :. Int)+type Scalar ix = Array ix++mul ::+ (A.Shape ix, Eq ix, Element a, A.Elt a) =>+ Cublas.Handle ->+ Exp a ->+ ALinAlg.Matrix ix a -> ALinAlg.Matrix ix a ->+ ALinAlg.Matrix ix a+mul handle alpha a b =+ A.foreignAcc+ (AF.CUDAForeignAcc "mul" $ uncurry3 $ mulPlain handle)+ (error "Requires CUDA backend")+ $+ A.lift (A.unit alpha, a, b)++mulPlain ::+ (A.Shape ix, Eq ix, Element a, A.Elt a) =>+ Cublas.Handle ->+ A.Scalar a -> Matrix ix a -> Matrix ix a ->+ AF.CIO (Matrix ix a)+mulPlain handle alpha a b = do+ let (aNumMatrices :. n :. k) = A.arrayShape a+ let (bNumMatrices :. _k :. m) = A.arrayShape b+ let numMatrices =+ if aNumMatrices == bNumMatrices+ then aNumMatrices+ else error "mul: mismatching shapes of matrix arrays"+ c <- AF.allocateArray (numMatrices :. n :. m)+ (pas, lda) <- arrayPtrs a+ (pbs, ldb) <- arrayPtrs b+ (pcs, ldc) <- arrayPtrs c+ AF.liftIO $ do+ Cublas.gemmBatched handle Cublas.N Cublas.N m n k+ (storableFromScalar alpha)+ pbs ldb+ pas lda+ 0+ pcs ldc+ (A.arraySize numMatrices)+ return c++mac ::+ (A.Shape ix, Eq ix, Element a, A.Elt a) =>+ Cublas.Handle ->+ Exp a -> ALinAlg.Matrix ix a -> ALinAlg.Matrix ix a ->+ Exp a -> ALinAlg.Matrix ix a ->+ ALinAlg.Matrix ix a+mac handle alpha a b beta c =+ A.foreignAcc+ (AF.CUDAForeignAcc "mac" $+ \(aalpha, aa, ab, abeta, ac) ->+ macPlain handle aalpha aa ab abeta ac)+ (error "Requires CUDA backend")+ $+ A.lift (A.unit alpha, a, b, A.unit beta, c)++macPlain ::+ (A.Shape ix, Eq ix, Element a, A.Elt a) =>+ Cublas.Handle ->+ A.Scalar a -> Matrix ix a -> Matrix ix a ->+ A.Scalar a -> Matrix ix a ->+ AF.CIO (Matrix ix a)+macPlain handle alpha a b beta c = do+ let (aNumMatrices :. n :. k ) = A.arrayShape a+ let (bNumMatrices :. _k :. m ) = A.arrayShape b+ let (cNumMatrices :. n' :. m') = A.arrayShape c+ let numMatrices =+ if aNumMatrices == bNumMatrices+ &&+ aNumMatrices == cNumMatrices+ then aNumMatrices+ else error "mac: mismatching shapes of matrix arrays"+ d <- AF.allocateArray (numMatrices :. n' :. m')+ AF.copyArray c d+ (pas, lda) <- arrayPtrs a+ (pbs, ldb) <- arrayPtrs b+ (pds, ldd) <- arrayPtrs d+ AF.liftIO $ do+ Cublas.gemmBatched handle Cublas.N Cublas.N m n k+ (storableFromScalar alpha)+ pbs ldb+ pas lda+ (storableFromScalar beta)+ pds ldd+ (A.arraySize numMatrices)+ return d++lu ::+ (A.Shape ix, Eq ix, Element a, A.Elt a) =>+ Cublas.Handle ->+ ALinAlg.Matrix ix a ->+ (ALinAlg.Matrix ix a, ALinAlg.Vector ix Int, ALinAlg.Scalar ix Int)+lu handle =+ A.unlift+ .+ A.foreignAcc+ (AF.CUDAForeignAcc "lu" $ luPlain handle)+ (error "Requires CUDA backend")++luPlain ::+ (A.Shape ix, Eq ix, Element a, A.Elt a) =>+ Cublas.Handle ->+ Matrix ix a ->+ AF.CIO (Matrix ix a, Vector ix Int, Scalar ix Int)+luPlain handle a = do+ let sh@(numMatrices :. n :. k) = A.arrayShape a+ let size =+ if n == k+ then n+ else error "lu: matrices must have square shape"+ b <- AF.allocateArray sh+ AF.copyArray a b+ (pbs, ldb) <- arrayPtrs b++ pivot <- AF.allocateArray (numMatrices :. size)+ pivotPtr <- fmap (castDevPtr . snd) $ AF.devicePtrsOfArray pivot++ info <- AF.allocateArray numMatrices+ infoPtr <- fmap (castDevPtr . snd) $ AF.devicePtrsOfArray info++ AF.liftIO $+ Cublas.getrfBatched handle size+ pbs ldb+ pivotPtr infoPtr+ (A.arraySize numMatrices)+ return (b, pivot, info)+++luInv ::+ (A.Shape ix, Eq ix, Element a, A.Elt a) =>+ Cublas.Handle ->+ (ALinAlg.Matrix ix a, ALinAlg.Vector ix Int, ALinAlg.Scalar ix Int) ->+ ALinAlg.Matrix ix a+luInv handle =+ A.foreignAcc+ (AF.CUDAForeignAcc "luInv" $ luInvPlain handle)+ (error "Requires CUDA backend")+ .+ A.lift++luInvPlain ::+ (A.Shape ix, Eq ix, Element a, A.Elt a) =>+ Cublas.Handle ->+ (Matrix ix a, Vector ix Int, Scalar ix Int) ->+ AF.CIO (Matrix ix a)+luInvPlain handle (a, pivot, info) = do+ let sh@(numMatrices :. n :. k) = A.arrayShape a+ let size =+ if n == k+ then n+ else error "luInv: matrices must have square shape"+ c <- AF.allocateArray sh+ AF.copyArray a c+ (pas, lda) <- arrayPtrs a+ (pcs, ldc) <- arrayPtrs c++ pivotPtr <- fmap (castDevPtr . snd) $ AF.devicePtrsOfArray pivot+ infoPtr <- fmap (castDevPtr . snd) $ AF.devicePtrsOfArray info++ AF.liftIO $+ Cublas.getriBatched handle size+ pas lda+ pivotPtr+ pcs ldc+ infoPtr+ (A.arraySize numMatrices)+ return c+++inv ::+ (A.Shape ix, Eq ix, Element a, A.Elt a) =>+ Cublas.Handle ->+ ALinAlg.Matrix ix a ->+ (ALinAlg.Matrix ix a, ALinAlg.Scalar ix Int)+inv handle a =+ let sol@(_,_,info) = lu handle a+ in (luInv handle sol, info)+++type Element a =+ (AF.DevicePtrs (EltRepr a) ~ ((), DevicePtr a),+ Fractional (StorableOf a),+ Cublas.Cublas (StorableOf a),+ Storable (StorableOf a),+ Real a)++type family StorableOf float+type instance StorableOf Float = CFloat+type instance StorableOf Double = CDouble++storableFromScalar ::+ (Real a, StorableOf a ~ b, Fractional b) => A.Scalar a -> b+storableFromScalar x = realToFrac $ A.indexArray x Z++arrayPtrs ::+ (Storable a, StorableOf e ~ a,+ A.Shape ix,+ AF.DevicePtrs (EltRepr e) ~ ((), DevicePtr e)) =>+ Array (ix :. Int :. Int) e -> AF.CIO ([DevicePtr a], Int)+arrayPtrs arr = do+ let (numMatrices :. n :. k) = A.arrayShape arr+ pa <- fmap (castDevPtr . snd) $ AF.devicePtrsOfArray arr+ return (genPointers (n*k) pa (A.arraySize numMatrices), k)++genPointers ::+ (Storable a) =>+ Int -> DevicePtr a -> Int -> [DevicePtr a]+genPointers size p n =+ take n $ iterate (flip advanceDevPtr size) p+++genMatrices :: (Acc (Array DIM3 Double), Acc (Array DIM3 Double))+genMatrices = (a,b)+ where+ a = A.generate (A.constant sha) $ \ix ->+ let (Z :. i :. j :. k) = unlift ix+ in A.fromIntegral (i+j+k)+ b = A.generate (A.constant shb) $ \ix ->+ let (Z :. i :. j :. k) = unlift ix+ in A.fromIntegral (i+j+k)+ numMats = 100 :: Int+ sha = Z :. numMats :. (3 :: Int) :. (4 :: Int)+ shb = Z :. numMats :. (4 :: Int) :. (2 :: Int)+ unlift :: Exp (Z :. Int :. Int :. Int)+ -> Z :. Exp Int :. Exp Int :. Exp Int+ unlift = A.unlift++test :: IO ()+test = do+ handle <- Cublas.create+ print genMatrices+ print $ AC.run $+ case genMatrices of+ (a,b) -> mul handle 1 a b
+ benchmark/NewtonInverse.hs view
@@ -0,0 +1,224 @@+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE FlexibleContexts #-}+module Main where++import qualified CUBLASBatched as Batched+import qualified Foreign.CUDA.Cublas as Cublas++import qualified Data.Array.Accelerate.Arithmetic.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) =>+ Cublas.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) =>+ Cublas.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++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 <- Cublas.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 <- Cublas.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::Double) $ 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
+ src/Data/Array/Accelerate/Arithmetic/Example.hs view
@@ -0,0 +1,49 @@+module Data.Array.Accelerate.Arithmetic.Example where++import qualified Data.Array.Accelerate.Arithmetic.Interpolation as Ip+import qualified Data.Array.Accelerate.Arithmetic.Sparse as Sparse+import Data.Array.Accelerate.Arithmetic.LinearAlgebra (Vector, )++import qualified Data.Array.Accelerate.Interpreter as AI+import qualified Data.Array.Accelerate as A+import Data.Array.Accelerate (Array, Z(Z), (:.)((:.)), )+++exampleSparseColumnMatrix :: IO ()+exampleSparseColumnMatrix = do+ let m :: Sparse.ColumnMatrix Z Double+ m =+ Sparse.ColumnMatrix (A.lift (3::Int)) $+ A.use $ A.fromList (Z :. 2 :. 5) $+ (0,1) : (2,2) : (1,3) : (0,4) : (2,5) :+ (1,6) : (2,7) : (0,8) : (2,9) : (1,10) :+ []++ v :: Vector Z Double+ v = A.use $ A.fromList (Z :. 5) [1,10,100,1000,10000]++ print $ AI.run $ Sparse.multiplyColumnMatrixVector m v++exampleSparseRowMatrix :: IO ()+exampleSparseRowMatrix = do+ let m :: Sparse.RowMatrix Z Double+ m =+ Sparse.RowMatrix (A.lift (5::Int)) $+ A.use $ A.fromList (Z :. 3 :. 2) $+ (0,1) : (0,2) :+ (3,3) : (1,4) :+ (3,5) : (4,6) :+ []++ v :: Vector Z Double+ v = A.use $ A.fromList (Z :. 5) [1,10,100,1000,10000]++ print $ AI.run $ Sparse.multiplyRowMatrixVector m v++exampleLookup :: IO ()+exampleLookup = do+ let nodes :: Array A.DIM2 Double+ nodes = A.fromList (Z :. 3 :. 5) [0 ..]+ x :: Array A.DIM1 Double+ x = A.fromList (Z :. 3) [0.2, 6.7, 13.1]+ print $ AI.run1 (Ip.lookupInterval (A.use nodes)) x
+ src/Data/Array/Accelerate/Arithmetic/Interpolation.hs view
@@ -0,0 +1,97 @@+module Data.Array.Accelerate.Arithmetic.Interpolation (+ bisect,+ lookupInterval,+ Interpolator13, sampleBasisFunctions13,+ ) where++import qualified Data.Array.Accelerate.Arithmetic.Sparse as Sparse+import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg+import qualified Data.Array.Accelerate.Utility.Arrange as Arrange+import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp+import qualified Data.Array.Accelerate.Utility.Loop as Loop+import Data.Array.Accelerate.Arithmetic.LinearAlgebra+ (Scalar, Vector, numElems, extrudeVector, )++import qualified Data.Array.Accelerate as A+import Data.Array.Accelerate (Exp, Any(Any), Z(Z), (:.)((:.)), )++import Data.Ord.HT (limit, )+++bisect ::+ (A.Slice ix, A.Shape ix, A.IsScalar a, A.Elt a) =>+ Vector ix a ->+ Scalar ix a ->+ Scalar ix (Int, Int) ->+ Scalar ix (Int, Int)+bisect nodes xs bounds =+ let centers =+ A.map+ (A.uncurry $ \lower upper -> div (lower+upper) 2)+ bounds+ in A.zipWith3+ (\center interval leftBranch ->+ A.cond leftBranch+ (Exp.mapSnd (const center) interval)+ (Exp.mapFst (const center) interval))+ centers bounds $+ A.zipWith (A.<*) xs $+ Arrange.gather (Arrange.mapWithIndex Exp.indexCons centers) nodes++lookupInterval ::+ (A.Slice ix, A.Shape ix, A.IsScalar a, A.Elt a) =>+ Vector ix a ->+ Scalar ix a ->+ Scalar ix Int+lookupInterval nodes x =+ A.map A.fst $+ Loop.nestLog2 (numElems nodes) (bisect nodes x) $+ A.fill (A.shape x) $+ A.lift (0 :: Exp Int, numElems nodes)+++outerVector ::+ (A.Shape ix, A.Slice ix, A.Elt a, A.Elt b, A.Elt c) =>+ (Exp a -> Exp b -> Exp c) ->+ Scalar ix a -> Vector Z b -> Vector ix c+outerVector f x y =+ A.zipWith f+ (A.replicate (A.lift $ Any :. numElems y) x)+ (extrudeVector (A.shape x) y)+++{- |+One node before index 0 and three nodes starting from index 0.+-}+type Interpolator13 a = (a,a) -> (a,a) -> (a,a) -> (a,a) -> a -> a++sampleBasisFunctions13 ::+ (A.Slice ix, A.Shape ix, A.Elt a, A.IsFloating a, Num a) =>+ Interpolator13 (Exp a) ->+ Vector Z a -> Vector ix a -> Sparse.RowMatrix ix a+sampleBasisFunctions13 interpolate nodes zs =+ Sparse.RowMatrix (numElems nodes) $+ let indices = lookupInterval (extrudeVector (A.shape zs) nodes) zs+ minIx = 1+ maxIx = numElems nodes - 3+ limitIndices = A.map (limit (minIx, maxIx)) indices+ gatherFromNodes d =+ LinAlg.gatherFromVector (A.map (d+) limitIndices) nodes+ in outerVector+ (A.lift2 $+ \(n, ln, z, x) (k, y) ->+ case (Exp.unliftQuadruple x, Exp.unliftQuadruple y) of+ ((xm1,x0,x1,x2), (ym1,y0,y1,y2)) ->+ (ln+k :: Exp Int,+ A.cond (n A.<* minIx) y0 $+ A.cond (n A.>* maxIx) y1 $+ interpolate (xm1,ym1) (x0,y0) (x1,y1) (x2,y2) z))+ (A.zip4 indices limitIndices zs+ (A.zip4+ (gatherFromNodes (-1))+ (gatherFromNodes 0)+ (gatherFromNodes 1)+ (gatherFromNodes 2)))+ (A.use $+ A.fromList (Z:.4)+ [(-1, (1,0,0,0)), (0, (0,1,0,0)), (1, (0,0,1,0)), (2, (0,0,0,1))])
+ src/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs view
@@ -0,0 +1,218 @@+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE FlexibleContexts #-}+module Data.Array.Accelerate.Arithmetic.LinearAlgebra where++import qualified Data.Array.Accelerate.Utility.Loop as Loop+import qualified Data.Array.Accelerate.Utility.Arrange as Arrange+import qualified Data.Array.Accelerate as A+import Data.Array.Accelerate+ (Acc, Array, Exp, Any(Any), All(All), Z(Z), (:.)((:.)))++++type Scalar ix a = Acc (Array ix a)+type Vector ix a = Acc (Array (ix :. Int) a)+type Matrix ix a = Acc (Array (ix :. Int :. Int) a)++transpose ::+ (A.Shape ix, A.Slice ix, A.Elt a) =>+ Matrix ix a -> Matrix ix a+transpose m =+ A.backpermute+ (A.lift $ swapIndex $ matrixShape m)+ (A.lift . swapIndex . A.unlift)+ m++swapIndex ::+ Exp ix :. Exp Int :. Exp Int ->+ Exp ix :. Exp Int :. Exp Int+swapIndex (ix :. r :. c) = (ix :. c :. r)+++numElems :: (A.Shape ix, A.Slice ix, A.Elt a) => Vector ix a -> Exp Int+numElems m = case vectorShape m of _ix :. n -> n++numRows :: (A.Shape ix, A.Slice ix, A.Elt a) => Matrix ix a -> Exp Int+numRows m = case matrixShape m of _ix :. rows :. _cols -> rows++numCols :: (A.Shape ix, A.Slice ix, A.Elt a) => Matrix ix a -> Exp Int+numCols m = case matrixShape m of _ix :. _rows :. cols -> cols++vectorShape ::+ (A.Shape ix, A.Slice ix, A.Elt a) =>+ Vector ix a -> Exp ix :. Exp Int+vectorShape m = A.unlift $ A.shape m++matrixShape ::+ (A.Shape ix, A.Slice ix, A.Elt a) =>+ Matrix ix a -> Exp ix :. Exp Int :. Exp Int+matrixShape m = A.unlift $ A.shape m++withVectorIndex ::+ (A.Shape ix, A.Slice ix, A.Lift Exp a) =>+ (Exp ix :. Exp Int -> a) ->+ (Exp (ix :. Int) -> Exp (A.Plain a))+withVectorIndex f = A.lift . f . A.unlift++withMatrixIndex ::+ (A.Shape ix, A.Slice ix, A.Lift Exp a) =>+ (Exp ix :. Exp Int :. Exp Int -> a) ->+ (Exp (ix :. Int :. Int) -> Exp (A.Plain a))+withMatrixIndex f = A.lift . f . A.unlift+++outer ::+ (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+ Vector ix a -> Vector ix a -> Matrix ix a+outer x y =+ A.zipWith (*)+ (A.replicate (A.lift $ Any :. All :. numElems y) x)+ (A.replicate (A.lift $ Any :. numElems x :. All) y)++multiplyMatrixVector ::+ (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+ Matrix ix a ->+ Vector ix a ->+ Vector ix a+multiplyMatrixVector m v =+ case matrixShape m of+ (_ix :. rows :. _cols) ->+ A.fold1 (+) $+ A.zipWith (*) m+ (A.replicate (A.lift $ Any :. rows :. All) v)++multiplyMatrixMatrix ::+ (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+ Matrix ix a ->+ Matrix ix a ->+ Matrix ix a+multiplyMatrixMatrix x y =+ case (matrixShape x, matrixShape y) of+ (_ :. rows :. _cols, _ :. _rows :. cols) ->+ A.fold1 (+) $ transpose $+ A.zipWith (*)+ (A.replicate (A.lift $ Any :. All :. All :. cols) x)+ (A.replicate (A.lift $ Any :. rows :. All :. All) y)++newtonInverseStep ::+ (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+ Matrix ix a ->+ Matrix ix a ->+ Matrix ix a+newtonInverseStep a x =+ A.zipWith (-) (A.map (2*) x) $+ multiplyMatrixMatrix x $ multiplyMatrixMatrix a x++identity ::+ (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+ Exp (ix :. Int :. Int) -> Matrix ix a+identity sh =+ A.generate sh+ (withMatrixIndex $+ \(_ :. r :. c) -> A.fromIntegral $ A.boolToInt (r A.==* c))++newtonInverse ::+ (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+ Exp Int ->+ Matrix ix a ->+ Matrix ix a ->+ Matrix ix a+newtonInverse n seed a =+ Loop.nest n (newtonInverseStep a) seed++++scaleRows ::+ (A.Slice ix, A.Shape ix, A.Elt a, A.IsNum a) =>+ Vector ix a -> Matrix ix a -> Matrix ix a+scaleRows s x =+ zipScalarVectorWith (*) s x++++zipScalarVectorWith ::+ (A.Slice ix, A.Shape ix, A.Elt a, A.Elt b, A.Elt c) =>+ (Exp a -> Exp b -> Exp c) ->+ Scalar ix a -> Vector ix b -> Vector ix c+zipScalarVectorWith f x ys =+ case vectorShape ys of+ _ix :. dim ->+ A.zipWith f (A.replicate (A.lift (Any :. dim)) x) ys++zipScalarMatrixWith ::+ (A.Slice ix, A.Shape ix, A.Elt a, A.Elt b, A.Elt c) =>+ (Exp a -> Exp b -> Exp c) ->+ Scalar ix a -> Matrix ix b -> Matrix ix c+zipScalarMatrixWith f x ys =+ case matrixShape ys of+ _ix :. rows :. cols ->+ A.zipWith f+ (A.replicate (A.lift (Any :. rows :. cols)) x) ys++flattenMatrix, flattenMatrixReshape, flattenMatrixBackPermute ::+ (A.Slice ix, A.Shape ix, A.Elt a) =>+ Matrix ix a -> Vector ix a+flattenMatrix = flattenMatrixBackPermute++flattenMatrixReshape m =+ case matrixShape m of+ ix :. rows :. cols ->+ A.reshape (A.lift $ ix :. rows*cols) m++accDivMod :: Integral a => a -> a -> (a, a)+accDivMod x y = (div x y, mod x y)++flattenMatrixBackPermute m =+ case matrixShape m of+ ix :. rows :. cols ->+ A.backpermute+ (A.lift $ ix :. rows*cols)+ (withVectorIndex $+ \(vix :. n) -> case accDivMod n cols of (r,c) -> vix :. r :. c)+ m+++restoreMatrix, restoreMatrixReshape, restoreMatrixBackPermute ::+ (A.Slice ix, A.Shape ix, A.Elt a) =>+ Exp Int -> Vector ix a -> Matrix ix a+restoreMatrix = restoreMatrixBackPermute++restoreMatrixReshape cols v =+ case vectorShape v of+ ix :. n ->+ A.reshape (A.lift $ ix :. div n cols :. cols) v++restoreMatrixBackPermute cols v =+ case vectorShape v of+ ix :. n ->+ A.backpermute+ (A.lift $ ix :. div n cols :. cols)+ (withMatrixIndex $ \(vix :. k :. j) -> vix :. k*cols+j)+ v++++extrudeVector ::+ (A.Shape ix, A.Slice ix, A.Elt a) =>+ Exp ix -> Vector Z a -> Vector ix a+extrudeVector shape y =+ -- A.replicate (A.lift $ shape :. All) y+ A.backpermute+ (A.lift $ shape :. numElems y)+ (A.index1 . A.indexHead)+ y++extrudeMatrix ::+ (A.Shape ix, A.Slice ix, A.Elt a) =>+ Exp ix -> Matrix Z a -> Matrix ix a+extrudeMatrix shape y =+ A.backpermute+ (A.lift $ shape :. numRows y :. numCols y)+ (withMatrixIndex $ \(_:.r:.c) -> Z:.r:.c)+ y++gatherFromVector ::+ (A.Shape ix, A.Elt a) =>+ Scalar ix Int -> Vector Z a -> Scalar ix a+gatherFromVector indices =+ Arrange.gather (A.map A.index1 indices)
+ src/Data/Array/Accelerate/Arithmetic/Sparse.hs view
@@ -0,0 +1,115 @@+{-# LANGUAGE TypeOperators #-}+module Data.Array.Accelerate.Arithmetic.Sparse where++import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg+import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp+import qualified Data.Array.Accelerate.Utility.Arrange as Arrange+import qualified Data.Array.Accelerate as A+import Data.Array.Accelerate.Utility.Lift.Exp (atom, )++import Data.Array.Accelerate.Arithmetic.LinearAlgebra+ (Matrix, Vector, matrixShape, )+import Data.Array.Accelerate+ (Exp, Any(Any), All(All), (:.)((:.)), )+++{- |+Sparse matrix with a definite number of non-zero entries per column.+-}+data ColumnMatrix ix a =+ ColumnMatrix {numRows :: Exp Int, columnMatrix :: Matrix ix (Int, a)}++realIndex ::+ (A.Shape ix, A.Slice ix, A.Elt a) =>+ Matrix ix (Int, a) ->+ Matrix ix (ix :. Int)+realIndex m =+ A.zipWith Exp.indexCons+ (A.generate (A.shape m) (A.indexTail . A.indexTail))+ (A.map A.fst m)++multiplyColumnMatrixVector ::+ (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+ ColumnMatrix ix a ->+ Vector ix a ->+ Vector ix a+multiplyColumnMatrixVector (ColumnMatrix rows m) v =+ Arrange.scatter (+)+ (realIndex m)+ (case matrixShape m of+ sh :. _rows :. _cols -> A.fill (A.lift $ sh :. rows) 0) $+ A.zipWith (*)+ (A.map A.snd m)+ (A.replicate (A.lift $ Any :. LinAlg.numRows m :. All) v)++transposeColumnMatrix ::+ (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+ ColumnMatrix ix a ->+ RowMatrix ix a+transposeColumnMatrix (ColumnMatrix n x) =+ RowMatrix n $ LinAlg.transpose x+++{- |+Sparse matrix with a definite number of non-zero entries per row.+-}+data RowMatrix ix a =+ RowMatrix {numCols :: Exp Int, rowMatrix :: Matrix ix (Int, a)}++multiplyRowMatrixVector ::+ (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+ RowMatrix ix a ->+ Vector ix a ->+ Vector ix a+multiplyRowMatrixVector (RowMatrix _cols m) v =+ A.fold1 (+) $+ A.zipWith (*) (A.map A.snd m) $+ Arrange.gather (realIndex m) v++transposeRowMatrix ::+ (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+ RowMatrix ix a ->+ ColumnMatrix ix a+transposeRowMatrix (RowMatrix n x) =+ (ColumnMatrix n $ LinAlg.transpose x)++multiplyMatrixMatrix ::+ (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+ ColumnMatrix ix a ->+ RowMatrix ix a ->+ Matrix ix a+multiplyMatrixMatrix+ (ColumnMatrix rows x) (RowMatrix cols y) =+ case matchMatrices x y of+ m ->+ let global = A.indexTail . A.indexTail . A.indexTail+ in Arrange.scatter (+)+ (Arrange.mapWithIndex+ (\mix tix ->+ A.lift $ global mix :. A.fst tix :. A.snd tix) $+ A.map A.fst m)+ (A.fill (A.lift $ global (A.shape m) :. rows :. cols) 0)+ (A.map A.snd m)++matchMatrices ::+ (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+ Matrix ix (Int, a) ->+ Matrix ix (Int, a) ->+ Matrix (ix :. Int) ((Int, Int), a)+matchMatrices x y =+ case (matrixShape x, matrixShape y) of+ (_ :. xRows :. _xCols, _ :. _yRows :. yCols) ->+ -- it must be xCols == yRows+ A.zipWith+ (Exp.modify2 (atom,atom) (atom,atom) $+ \(n,xi) (m,yi) -> ((n, m), xi*yi))+ (A.replicate (A.lift $ Any :. All :. All :. yCols) x)+ (A.replicate (A.lift $ Any :. xRows :. All :. All) y)+++scaleRowRows ::+ (A.Slice ix, A.Shape ix, A.Elt a, A.IsNum a) =>+ Vector ix a -> RowMatrix ix a -> RowMatrix ix a+scaleRowRows s (RowMatrix n x) =+ RowMatrix n $+ LinAlg.zipScalarVectorWith (\si xi -> Exp.mapSnd (si*) xi) s x
+ test/Test.hs view
@@ -0,0 +1,20 @@+module Main where++import qualified Test.Data.Array.Accelerate.Arithmetic.Sparse as Sparse+import qualified Test.Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg++import qualified Test.QuickCheck.Modifiers as Mod+import Test.QuickCheck (quickCheck)+++test :: IO ()+test = mapM_ (\(msg,act) -> putStr (msg++": ") >> act) $+ ("sparseMatrix", quickCheck (\(Mod.Blind x) -> Sparse.multiplication x)) :+ ("flattenMatrix", quickCheck (\(Mod.Blind x) -> LinAlg.flattenMatrix x)) :+ ("restoreMatrix", quickCheck (\(Mod.Blind x) -> LinAlg.restoreMatrix x)) :+ ("flattenRestoreMatrix", quickCheck (\(Mod.Blind x) -> LinAlg.flattenRestoreMatrix x)) :+ []+++main :: IO ()+main = test
+ test/Test/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs view
@@ -0,0 +1,41 @@+module Test.Data.Array.Accelerate.Arithmetic.LinearAlgebra where++import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg+import qualified Data.Array.Accelerate as A++import Data.Array.Accelerate.Arithmetic.LinearAlgebra (Matrix, numCols, )+import Data.Array.Accelerate (Z(Z), (:.)((:.)),)++import Test.Data.Array.Accelerate.Arithmetic.Utility (arbitraryArray, (=!=), )++import qualified Test.QuickCheck as QC++import Data.Word (Word32, )+++newtype ArbMatrix a = ArbMatrix (Matrix Z a)++instance (QC.Arbitrary a, A.Elt a) => QC.Arbitrary (ArbMatrix a) where+ arbitrary = do+ nc <- QC.choose (1,100)+ nr <- QC.choose (1,100)+ fmap (ArbMatrix . A.use) $+ arbitraryArray (Z :. nr :. nc) QC.arbitrary+++flattenMatrix :: ArbMatrix Word32 -> Bool+flattenMatrix (ArbMatrix m) =+ LinAlg.flattenMatrixReshape m+ =!=+ LinAlg.flattenMatrixBackPermute m++restoreMatrix :: ArbMatrix Word32 -> Bool+restoreMatrix (ArbMatrix m) =+ let v = LinAlg.flattenMatrix m+ in LinAlg.restoreMatrixReshape (numCols m) v+ =!=+ LinAlg.restoreMatrixBackPermute (numCols m) v++flattenRestoreMatrix :: ArbMatrix Word32 -> Bool+flattenRestoreMatrix (ArbMatrix m) =+ m =!= LinAlg.restoreMatrix (numCols m) (LinAlg.flattenMatrix m)
+ test/Test/Data/Array/Accelerate/Arithmetic/Sparse.hs view
@@ -0,0 +1,52 @@+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE TypeOperators #-}+module Test.Data.Array.Accelerate.Arithmetic.Sparse where++import Test.Data.Array.Accelerate.Arithmetic.Utility (arbitraryArray, (=!=), )++import qualified Data.Array.Accelerate.Arithmetic.Sparse as Sparse+import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg++import qualified Data.Array.Accelerate as A+import Data.Array.Accelerate (Z(Z), (:.)((:.)))++import qualified Test.QuickCheck as QC++import Control.Monad (liftM2, )++import Data.Word (Word32, )+++data+ CRVTriple a =+ CRVTriple+ (Sparse.ColumnMatrix Z a)+ (Sparse.RowMatrix Z a)+ (LinAlg.Vector Z a)++instance (QC.Arbitrary a, A.Elt a) => QC.Arbitrary (CRVTriple a) where+ arbitrary = do+ k <- QC.choose (1,200)+ nc <- QC.choose (1,100)+ nr <- QC.choose (1,100)+ cc <- QC.choose (1,10)+ cr <- QC.choose (1,10)+ mc <-+ arbitraryArray (Z :. cc :. k) $+ liftM2 (,) (QC.choose (0,nc-1)) QC.arbitrary+ mr <-+ arbitraryArray (Z :. k :. cr) $+ liftM2 (,) (QC.choose (0,nr-1)) QC.arbitrary+ v <- arbitraryArray (Z :. nr) QC.arbitrary+ return $+ CRVTriple+ (Sparse.ColumnMatrix (A.lift nc) (A.use mc))+ (Sparse.RowMatrix (A.lift nr) (A.use mr))+ (A.use v)+++multiplication :: CRVTriple Word32 -> Bool+multiplication (CRVTriple mc mr v) =+ LinAlg.multiplyMatrixVector (Sparse.multiplyMatrixMatrix mc mr) v+ =!=+ Sparse.multiplyColumnMatrixVector mc (Sparse.multiplyRowMatrixVector mr v)
+ test/Test/Data/Array/Accelerate/Arithmetic/Utility.hs view
@@ -0,0 +1,48 @@+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE FlexibleInstances #-}+module Test.Data.Array.Accelerate.Arithmetic.Utility where++import qualified Data.Array.Accelerate.Interpreter as AI+import qualified Data.Array.Accelerate as A+import Data.Array.Accelerate (Acc, Array, Z, (:.)((:.)), )++import qualified Test.QuickCheck as QC+++infix 4 =!=++(=!=) ::+ (Eq sh, Eq e, A.Shape sh, A.Elt e) =>+ Acc (Array sh e) -> Acc (Array sh e) -> Bool+x =!= y =+ let xi = AI.run x+ yi = AI.run y+ in A.arrayShape xi == A.arrayShape yi+ &&+ A.toList xi == A.toList yi+++class A.Shape sh => Shape sh where+ switchShape ::+ f Z ->+ (forall sh1. Shape sh1 => f (sh1 :. Int)) ->+ f sh++instance Shape Z where switchShape f _ = f+instance Shape sh => Shape (sh :. Int) where switchShape _ f = f++newtype ShapeSize sh = ShapeSize {getShapeSize :: sh -> Int}++shapeSize :: Shape sh => sh -> Int+shapeSize =+ getShapeSize $+ switchShape+ (ShapeSize $ const 1)+ (ShapeSize $ \(ix:.n) -> shapeSize ix * n)++arbitraryArray ::+ (Shape sh, A.Elt a) => sh -> QC.Gen a -> QC.Gen (Array sh a)+arbitraryArray sh gen =+ fmap (A.fromList sh) $+ QC.vectorOf (shapeSize sh) gen