accelerate-arithmetic 0.0.1 → 0.1

raw patch · 14 files changed

+533/−919 lines, 14 filesdep −accelerate-cudadep −cublasdep −cudadep ~QuickCheckdep ~acceleratedep ~basePVP ok

version bump matches the API change (PVP)

Dependencies removed: accelerate-cuda, cublas, cuda, hmatrix, pooled-io, random, timeit

Dependency ranges changed: QuickCheck, accelerate, base, utility-ht

API changes (from Hackage documentation)

Files

• accelerate-arithmetic.cabal +15/−33
• benchmark/CUBLASBatched.hs +0/−269
• benchmark/NewtonInverse.hs +0/−224
• private/Data/Array/Accelerate/LinearAlgebra/Private.hs +254/−0
• src/Data/Array/Accelerate/Arithmetic/Example.hs +8/−8
• src/Data/Array/Accelerate/Arithmetic/Interpolation.hs +5/−5
• src/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs +0/−254
• src/Data/Array/Accelerate/Arithmetic/Sparse.hs +0/−115
• src/Data/Array/Accelerate/LinearAlgebra.hs +35/−0
• src/Data/Array/Accelerate/LinearAlgebra/Matrix/Banded.hs +26/−0
• src/Data/Array/Accelerate/LinearAlgebra/Matrix/Sparse.hs +148/−0
• test/Test.hs +1/−0
• test/Test/Data/Array/Accelerate/Arithmetic/LinearAlgebra.hs +2/−2
• test/Test/Data/Array/Accelerate/Arithmetic/Sparse.hs +39/−9

@@ -1,10 +1,10 @@ Name:             accelerate-arithmetic-Version:          0.0.1+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-arithmetic/+Homepage:         http://hub.darcs.net/thielema/accelerate-arithmetic/ Category:         Math Synopsis:         Linear algebra and interpolation using the Accelerate framework Description:@@ -17,66 +17,48 @@ Build-Type:       Simple  Source-Repository this-  Tag:         0.0.1+  Tag:         0.1   Type:        darcs-  Location:    http://code.haskell.org/~thielema/accelerate-arithmetic/+  Location:    http://hub.darcs.net/thielema/accelerate-arithmetic/  Source-Repository head   Type:        darcs-  Location:    http://code.haskell.org/~thielema/accelerate-arithmetic/+  Location:    http://hub.darcs.net/thielema/accelerate-arithmetic/  Library   Build-Depends:     accelerate-utility >=0.1 && <0.2,     accelerate >=0.15 && <0.16,     utility-ht >=0.0.8 && <0.1,-    QuickCheck >=2.4 && <2.8,-    base >=4.5 && <4.8+    QuickCheck >=2.4 && <3,+    base >=4.5 && <4.10    GHC-Options:      -Wall -fwarn-missing-import-lists-  Hs-Source-Dirs:   src+  Hs-Source-Dirs:   src, private   Default-Language: Haskell98   Exposed-Modules:-    Data.Array.Accelerate.Arithmetic.LinearAlgebra-    Data.Array.Accelerate.Arithmetic.Sparse     Data.Array.Accelerate.Arithmetic.Interpolation+    Data.Array.Accelerate.LinearAlgebra+    Data.Array.Accelerate.LinearAlgebra.Matrix.Sparse+    Data.Array.Accelerate.LinearAlgebra.Matrix.Banded   Other-Modules:     Data.Array.Accelerate.Arithmetic.Example+    Data.Array.Accelerate.LinearAlgebra.Private  Test-Suite test   Type: exitcode-stdio-1.0   Main-Is:          Test.hs   GHC-Options:      -Wall -fwarn-missing-import-lists-  Hs-Source-Dirs:   test+  Hs-Source-Dirs:   test, private   Default-Language: Haskell98   Build-Depends:     accelerate-arithmetic,+    accelerate-utility,     accelerate,     QuickCheck,     base   Other-Modules:+    Data.Array.Accelerate.LinearAlgebra.Private     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

@@ -1,269 +0,0 @@-{-# LANGUAGE TypeOperators #-}-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE ConstraintKinds #-}-{-# LANGUAGE FlexibleContexts #-}-module CUBLASBatched where--import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as ALinAlg--import qualified Data.Array.Accelerate.Utility.Lift.Acc as Acc-import Data.Array.Accelerate.Utility.Lift.Acc (acc, expr)--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, A.Slice ix, Eq ix, Element a, A.Elt a, A.IsNum 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)-      (Acc.modify (expr,acc,acc) $ \(alpha0, a0, b0) ->-         A.map (alpha0 *) $-         ALinAlg.multiplyMatrixMatrix a0 b0)-   $-   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 $-      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, A.Slice ix, Eq ix, Element a, A.Elt a, A.IsNum 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" $-       \((alpha0, a0, b0), (beta0, c0)) ->-          macPlain handle alpha0 a0 b0 beta0 c0)-      (Acc.modify ((expr,acc,acc),(expr,acc)) $-       \((alpha0, a0, b0), (beta0, c0)) ->-         A.zipWith (+)-            (A.map (alpha0 *) $-             ALinAlg.multiplyMatrixMatrix a0 b0)-            (A.map (beta0 *) c0))-   $-   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 $-      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

@@ -1,224 +0,0 @@-{-# 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

@@ -0,0 +1,254 @@+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE FlexibleContexts #-}+module Data.Array.Accelerate.LinearAlgebra.Private where++import qualified Data.Array.Accelerate.Utility.Loop as Loop+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+          (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++++columnFromVector ::+   (A.Shape ix, A.Slice ix, A.Elt a) =>+   Vector ix a -> Matrix ix a+columnFromVector a = A.reshape (Exp.indexCons (A.shape a) 1) a++{- |+input must be a matrix with exactly one column+-}+vectorFromColumn ::+   (A.Shape ix, A.Slice ix, A.Elt a) =>+   Matrix ix a -> Vector ix a+vectorFromColumn a = A.reshape (A.indexTail $ A.shape a) a++++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++zipExtrudedVectorWith ::+   (A.Slice ix, A.Shape ix, A.Elt a, A.Elt b, A.Elt c) =>+   (Exp a -> Exp b -> Exp c) ->+   Vector Z a ->+   Vector ix b ->+   Vector ix c+zipExtrudedVectorWith f x y =+   A.zipWith f (extrudeVector (A.indexTail $ A.shape y) x) y++zipExtrudedMatrixWith ::+   (A.Slice ix, A.Shape ix, A.Elt a, A.Elt b, A.Elt c) =>+   (Exp a -> Exp b -> Exp c) ->+   Matrix Z a ->+   Matrix ix b ->+   Matrix ix c+zipExtrudedMatrixWith f x y =+   A.zipWith f (extrudeMatrix (A.indexTail $ A.indexTail $ A.shape y) x) 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)

@@ -1,8 +1,8 @@ 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.LinearAlgebra.Matrix.Sparse as Sparse+import Data.Array.Accelerate.LinearAlgebra (Vector, )  import qualified Data.Array.Accelerate.Interpreter as AI import qualified Data.Array.Accelerate as A@@ -11,9 +11,9 @@  exampleSparseColumnMatrix :: IO () exampleSparseColumnMatrix = do-   let m :: Sparse.ColumnMatrix Z Double+   let m :: Sparse.Columns Z Double        m =-          Sparse.ColumnMatrix (A.lift (3::Int)) $+          Sparse.Columns (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) :@@ -22,13 +22,13 @@        v :: Vector Z Double        v = A.use $ A.fromList (Z :. 5) [1,10,100,1000,10000] -   print $ AI.run $ Sparse.multiplyColumnMatrixVector m v+   print $ AI.run $ Sparse.multiplyColumnsVector m v  exampleSparseRowMatrix :: IO () exampleSparseRowMatrix = do-   let m :: Sparse.RowMatrix Z Double+   let m :: Sparse.Rows Z Double        m =-          Sparse.RowMatrix (A.lift (5::Int)) $+          Sparse.Rows (A.lift (5::Int)) $           A.use $ A.fromList (Z :. 3 :. 2) $           (0,1) : (0,2) :           (3,3) : (1,4) :@@ -38,7 +38,7 @@        v :: Vector Z Double        v = A.use $ A.fromList (Z :. 5) [1,10,100,1000,10000] -   print $ AI.run $ Sparse.multiplyRowMatrixVector m v+   print $ AI.run $ Sparse.multiplyRowsVector m v  exampleLookup :: IO () exampleLookup = do

@@ -4,12 +4,12 @@    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.LinearAlgebra.Matrix.Sparse as Sparse+import qualified Data.Array.Accelerate.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+import Data.Array.Accelerate.LinearAlgebra           (Scalar, Vector, numElems, extrudeVector, )  import qualified Data.Array.Accelerate as A@@ -68,9 +68,9 @@ 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+   Vector Z a -> Vector ix a -> Sparse.Rows ix a sampleBasisFunctions13 interpolate nodes zs =-   Sparse.RowMatrix (numElems nodes) $+   Sparse.Rows (numElems nodes) $    let indices = lookupInterval (extrudeVector (A.shape zs) nodes) zs        minIx = 1        maxIx = numElems nodes - 3

@@ -1,254 +0,0 @@-{-# 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.Lift.Exp as Exp-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----columnFromVector ::-   (A.Shape ix, A.Slice ix, A.Elt a) =>-   Vector ix a -> Matrix ix a-columnFromVector a = A.reshape (Exp.indexCons (A.shape a) 1) a--{- |-input must be a matrix with exactly one column--}-vectorFromColumn ::-   (A.Shape ix, A.Slice ix, A.Elt a) =>-   Matrix ix a -> Vector ix a-vectorFromColumn a = A.reshape (A.indexTail $ A.shape a) a----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--zipExtrudedVectorWith ::-   (A.Slice ix, A.Shape ix, A.Elt a, A.Elt b, A.Elt c) =>-   (Exp a -> Exp b -> Exp c) ->-   Vector Z a ->-   Vector ix b ->-   Vector ix c-zipExtrudedVectorWith f x y =-   A.zipWith f (extrudeVector (A.indexTail $ A.shape y) x) y--zipExtrudedMatrixWith ::-   (A.Slice ix, A.Shape ix, A.Elt a, A.Elt b, A.Elt c) =>-   (Exp a -> Exp b -> Exp c) ->-   Matrix Z a ->-   Matrix ix b ->-   Matrix ix c-zipExtrudedMatrixWith f x y =-   A.zipWith f (extrudeMatrix (A.indexTail $ A.indexTail $ A.shape y) x) 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)

@@ -1,115 +0,0 @@-{-# 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 (expr, )--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 (expr,expr) (expr,expr) $-             \(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

@@ -0,0 +1,35 @@+{-# LANGUAGE TypeOperators #-}+{-# LANGUAGE FlexibleContexts #-}+module Data.Array.Accelerate.LinearAlgebra (+   LinAlg.Scalar,+   LinAlg.Vector,+   LinAlg.Matrix,+   LinAlg.transpose,+   LinAlg.numElems,+   LinAlg.numRows,+   LinAlg.numCols,+   LinAlg.vectorShape,+   LinAlg.matrixShape,+   LinAlg.withVectorIndex,+   LinAlg.withMatrixIndex,+   LinAlg.outer,+   LinAlg.multiplyMatrixVector,+   LinAlg.multiplyMatrixMatrix,+   LinAlg.newtonInverse,+   LinAlg.newtonInverseStep,+   LinAlg.identity,+   LinAlg.scaleRows,+   LinAlg.zipScalarVectorWith,+   LinAlg.zipScalarMatrixWith,+   LinAlg.columnFromVector,+   LinAlg.vectorFromColumn,+   LinAlg.flattenMatrix,+   LinAlg.restoreMatrix,+   LinAlg.extrudeVector,+   LinAlg.extrudeMatrix,+   LinAlg.zipExtrudedVectorWith,+   LinAlg.zipExtrudedMatrixWith,+   LinAlg.gatherFromVector,+   ) where++import qualified Data.Array.Accelerate.LinearAlgebra.Private as LinAlg

@@ -0,0 +1,26 @@+module Data.Array.Accelerate.LinearAlgebra.Matrix.Banded (+   Symmetric(..),+   flattenSymmetric,+   ) where++import Data.Array.Accelerate.LinearAlgebra (Matrix, matrixShape)++import qualified Data.Array.Accelerate.Utility.Lift.Exp as Exp+import qualified Data.Array.Accelerate as A+import Data.Array.Accelerate.Utility.Lift.Exp (expr)+import Data.Array.Accelerate ((:.)((:.)), (>*), (!), (?))+++newtype Symmetric ix a = Symmetric (Matrix ix a)++flattenSymmetric ::+   (A.Slice ix, A.Shape ix, A.Elt a, A.IsNum a) =>+   Symmetric ix a -> Matrix ix a+flattenSymmetric (Symmetric m) =+   case matrixShape m of+      (sh :. rows :. width) ->+         A.generate (A.lift $ sh :. rows :. rows) $+         Exp.modify (expr:.expr:.expr) $ \(ix:.k0:.j0) ->+            let k = min k0 j0+                j = max k0 j0 - k+            in  width >* j ? (m ! A.lift(ix:.k:.j), 0)

@@ -0,0 +1,148 @@+{-# LANGUAGE TypeOperators #-}+module Data.Array.Accelerate.LinearAlgebra.Matrix.Sparse (+   Columns(..),+   multiplyColumnsVector,+   transposeColumns,+   Rows(..),+   multiplyRowsVector,+   transposeRows,+   multiplyColumnsRows,+   realBandedGramian,+   scaleRowRows,+   ) where++import qualified Data.Array.Accelerate.LinearAlgebra.Matrix.Banded as BandMatrix+import qualified Data.Array.Accelerate.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 (expr, )++import Data.Array.Accelerate.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 Columns ix a =+        Columns {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)++multiplyColumnsVector ::+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+   Columns ix a ->+   Vector ix a ->+   Vector ix a+multiplyColumnsVector (Columns 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)++transposeColumns ::+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+   Columns ix a ->+   Rows ix a+transposeColumns (Columns n x) =+   Rows n $ LinAlg.transpose x+++{- |+Sparse matrix with a definite number of non-zero entries per row.+-}+data Rows ix a =+        Rows {numCols :: Exp Int, rowMatrix :: Matrix ix (Int, a)}++multiplyRowsVector ::+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+   Rows ix a ->+   Vector ix a ->+   Vector ix a+multiplyRowsVector (Rows _cols m) v =+   A.fold1 (+) $+   A.zipWith (*) (A.map A.snd m) $+   Arrange.gather (realIndex m) v++transposeRows ::+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+   Rows ix a ->+   Columns ix a+transposeRows (Rows n x) =+   (Columns n $ LinAlg.transpose x)++multiplyColumnsRows ::+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+   Columns ix a ->+   Rows ix a ->+   Matrix ix a+multiplyColumnsRows (Columns rows x) (Rows cols y) =+   let (ixs,prods) = A.unzip $ matchMatrices x y+       global = A.indexTail . A.indexTail . A.indexTail+   in  Arrange.scatter (+)+          (Arrange.mapWithIndex+             (Exp.modify2 expr (expr,expr) $ \mix (k,j) ->+                global mix :. k :. j) $+           ixs)+          (A.fill (A.lift $ global (A.shape prods) :. rows :. cols) 0)+          prods++{- |+Compute x^T*x, given that it has a band structure.+You must pass the band-width as parameter+and you must make sure that the Gramian stays within this band.+Otherwise you cause out-of-bounds array accesses.+So far, only correct for real matrices.+-}+realBandedGramian ::+   (A.Shape ix, A.Slice ix, A.IsNum a, A.Elt a) =>+   Exp Int ->+   Rows ix a ->+   BandMatrix.Symmetric ix a+realBandedGramian width (Rows cols y) =+   let (ixs,prods) = A.unzip $ matchMatrices (LinAlg.transpose y) y+       global = A.indexTail . A.indexTail . A.indexTail+   in  BandMatrix.Symmetric $+       Arrange.scatter (+)+          (Arrange.mapWithIndex+             (Exp.modify2 expr (expr,expr) $ \mix (k,j) ->+                k>*j ? (A.ignore, A.lift $ global mix :. k :. j-k)) $+           ixs)+          (A.fill (A.lift $ global (A.shape prods) :. cols :. width) 0)+          prods++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 (expr,expr) (expr,expr) $+             \(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 -> Rows ix a -> Rows ix a+scaleRowRows s (Rows n x) =+   Rows n $+   LinAlg.zipScalarVectorWith (\si xi -> Exp.mapSnd (si*) xi) s x

@@ -10,6 +10,7 @@ test :: IO () test = mapM_ (\(msg,act) -> putStr (msg++": ") >> act) $    ("sparseMatrix", quickCheck (\(Mod.Blind x) -> Sparse.multiplication x)) :+   ("bandedGramian", quickCheck (\(Mod.Blind x) -> Sparse.bandedGramian 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)) :

@@ -1,9 +1,9 @@ module Test.Data.Array.Accelerate.Arithmetic.LinearAlgebra where -import qualified Data.Array.Accelerate.Arithmetic.LinearAlgebra as LinAlg+import qualified Data.Array.Accelerate.LinearAlgebra.Private as LinAlg import qualified Data.Array.Accelerate as A -import Data.Array.Accelerate.Arithmetic.LinearAlgebra (Matrix, numCols, )+import Data.Array.Accelerate.LinearAlgebra.Private (Matrix, numCols, ) import Data.Array.Accelerate (Z(Z), (:.)((:.)),)  import Test.Data.Array.Accelerate.Arithmetic.Utility (arbitraryArray, (=!=), )

@@ -4,11 +4,12 @@  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.LinearAlgebra.Matrix.Banded as BandMatrix+import qualified Data.Array.Accelerate.LinearAlgebra.Matrix.Sparse as Sparse+import qualified Data.Array.Accelerate.LinearAlgebra as LinAlg  import qualified Data.Array.Accelerate as A-import Data.Array.Accelerate (Z(Z), (:.)((:.)))+import Data.Array.Accelerate (Exp, Z(Z), (:.)((:.)))  import qualified Test.QuickCheck as QC @@ -20,8 +21,8 @@ data    CRVTriple a =       CRVTriple-         (Sparse.ColumnMatrix Z a)-         (Sparse.RowMatrix Z a)+         (Sparse.Columns Z a)+         (Sparse.Rows Z a)          (LinAlg.Vector Z a)  instance (QC.Arbitrary a, A.Elt a) => QC.Arbitrary (CRVTriple a) where@@ -40,13 +41,42 @@       v <- arbitraryArray (Z :. nr) QC.arbitrary       return $          CRVTriple-            (Sparse.ColumnMatrix (A.lift nc) (A.use mc))-            (Sparse.RowMatrix (A.lift nr) (A.use mr))+            (Sparse.Columns (A.lift nc) (A.use mc))+            (Sparse.Rows (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+   LinAlg.multiplyMatrixVector (Sparse.multiplyColumnsRows mc mr) v    =!=-   Sparse.multiplyColumnMatrixVector mc (Sparse.multiplyRowMatrixVector mr v)+   Sparse.multiplyColumnsVector mc (Sparse.multiplyRowsVector mr v)++++data BandGramian a = BandGramian (Exp Int) (Sparse.Rows Z a)++instance (QC.Arbitrary a, A.Elt a) => QC.Arbitrary (BandGramian a) where+   arbitrary = do+      width <- QC.choose (1,10)+      rows <- QC.choose (1,100)+      cols <- QC.choose (width,100)++      m <-+         fmap (A.fromList (Z :. rows :. width) . concat) $+         QC.vectorOf rows $+         liftM2+            (\start row -> zip [start..] row)+            (QC.choose (0,cols-width))+            (QC.vectorOf width QC.arbitrary)++      return $+         BandGramian (A.lift width)+            (Sparse.Rows (A.lift cols) (A.use m))+++bandedGramian :: BandGramian Word32 -> Bool+bandedGramian (BandGramian width m) =+   Sparse.multiplyColumnsRows (Sparse.transposeRows m) m+   =!=+   BandMatrix.flattenSymmetric (Sparse.realBandedGramian width m)