diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Setup.lhs b/Setup.lhs new file mode 100644 --- /dev/null +++ b/Setup.lhs @@ -0,0 +1,3 @@ +#! /usr/bin/env runhaskell +> import Distribution.Simple +> main = defaultMain diff --git a/accelerate-arithmetic.cabal b/accelerate-arithmetic.cabal new file mode 100644 --- /dev/null +++ b/accelerate-arithmetic.cabal @@ -0,0 +1,82 @@ +Name: accelerate-arithmetic +Version: 0.0 +License: BSD3 +License-File: LICENSE +Author: Henning Thielemann +Maintainer: Henning Thielemann +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 diff --git a/benchmark/CUBLASBatched.hs b/benchmark/CUBLASBatched.hs new file mode 100644 --- /dev/null +++ b/benchmark/CUBLASBatched.hs @@ -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 diff --git a/benchmark/NewtonInverse.hs b/benchmark/NewtonInverse.hs new file mode 100644 --- /dev/null +++ b/benchmark/NewtonInverse.hs @@ -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 diff --git a/src/Data/Array/Accelerate/Arithmetic/Example.hs b/src/Data/Array/Accelerate/Arithmetic/Example.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Arithmetic/Example.hs @@ -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 diff --git a/src/Data/Array/Accelerate/Arithmetic/Interpolation.hs b/src/Data/Array/Accelerate/Arithmetic/Interpolation.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Arithmetic/Interpolation.hs @@ -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))]) diff --git a/src/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs b/src/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs @@ -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) diff --git a/src/Data/Array/Accelerate/Arithmetic/Sparse.hs b/src/Data/Array/Accelerate/Arithmetic/Sparse.hs new file mode 100644 --- /dev/null +++ b/src/Data/Array/Accelerate/Arithmetic/Sparse.hs @@ -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 diff --git a/test/Test.hs b/test/Test.hs new file mode 100644 --- /dev/null +++ b/test/Test.hs @@ -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 diff --git a/test/Test/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs b/test/Test/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs new file mode 100644 --- /dev/null +++ b/test/Test/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs @@ -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) diff --git a/test/Test/Data/Array/Accelerate/Arithmetic/Sparse.hs b/test/Test/Data/Array/Accelerate/Arithmetic/Sparse.hs new file mode 100644 --- /dev/null +++ b/test/Test/Data/Array/Accelerate/Arithmetic/Sparse.hs @@ -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) diff --git a/test/Test/Data/Array/Accelerate/Arithmetic/Utility.hs b/test/Test/Data/Array/Accelerate/Arithmetic/Utility.hs new file mode 100644 --- /dev/null +++ b/test/Test/Data/Array/Accelerate/Arithmetic/Utility.hs @@ -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