I don't mean to be blunt, but have you guys taken a course in linear algebra?
On Mon, Dec 3, 2012 at 9:21 PM, Trevor L. McDonell < tmcdon...@cse.unsw.edu.au> wrote: > As far as I am aware, the only description is in the Repa paper. I you are > right, it really should be explained properly somewhere… > > At a simpler example, here is the outer product of two vectors [1]. > > vvProd :: (IsNum e, Elt e) => Acc (Vector e) -> Acc (Vector e) -> Acc > (Matrix e) > vvProd xs ys = A.zipWith (*) xsRepl ysRepl > where > n = A.size xs > m = A.size ys > > xsRepl = A.replicate (lift (Z :. All :. m )) xs > ysRepl = A.replicate (lift (Z :. n :. All)) ys > > If we then `A.fold (+) 0` the matrix, it would reduce along each row > producing a vector. So the first element of that vector is going to be > calculated as (xs[0] * ys[0] + xs[0] * ys[1] + … xs[0] * ys[m-1]). That's > the idea we want for our matrix multiplication … but I agree, it is > difficult for me to visualise as well. > > I do the same sort of trick with the n-body demo to get all n^2 particle > interactions. > > -Trev > > > [1]: http://en.wikipedia.org/wiki/Outer_product#Vector_multiplication > > > > On 04/12/2012, at 3:41 AM, Clark Gaebel <cgae...@uwaterloo.ca> wrote: > > Ah. I see now. Silly Haskell making inefficient algorithms hard to write > and efficient ones easy. It's actually kind of annoying when learning, but > probably for the best. > > Is there a good write-up of the algorithm you're using somewhere? The Repa > paper was very brief in its explaination, and I'm having trouble > visualizing the mapping of the 2D matricies into 3 dimensions. > > - Clark > > > On Mon, Dec 3, 2012 at 2:06 AM, Trevor L. McDonell < > tmcdon...@cse.unsw.edu.au> wrote: > >> Hi Clark, >> >> The trick is that most accelerate operations work over multidimensional >> arrays, so you can still get around the fact that we are limited to flat >> data-parallelism only. >> >> Here is matrix multiplication in Accelerate, lifted from the first Repa >> paper [1]. >> >> >> import Data.Array.Accelerate as A >> >> type Matrix a = Array DIM2 a >> >> matMul :: (IsNum e, Elt e) => Acc (Matrix e) -> Acc (Matrix e) -> Acc >> (Matrix e) >> matMul arr brr >> = A.fold (+) 0 >> $ A.zipWith (*) arrRepl brrRepl >> where >> Z :. rowsA :. _ = unlift (shape arr) :: Z :. Exp Int :. Exp Int >> Z :. _ :. colsB = unlift (shape brr) :: Z :. Exp Int :. Exp Int >> >> arrRepl = A.replicate (lift $ Z :. All :. colsB :. All) >> arr >> brrRepl = A.replicate (lift $ Z :. rowsA :. All :. All) >> (A.transpose brr) >> >> >> If you use github sources rather than the hackage package, those >> intermediate replicates will get fused away. >> >> >> Cheers, >> -Trev >> >> [1] http://www.cse.unsw.edu.au/~chak/papers/KCLPL10.html >> >> >> >> >> On 03/12/2012, at 5:07 PM, Clark Gaebel <cgae...@uwaterloo.ca> wrote: >> >> Hello cafe, >> >> I've recently started learning about cuda and hetrogenous programming, >> and have been using accelerate [1] to help me out. Right now, I'm running >> into trouble in that I can't call parallel code from sequential code. Turns >> out GPUs aren't exactly like Repa =P. >> >> Here's what I have so far: >> >> import qualified Data.Array.Accelerate as A >> import Data.Array.Accelerate ( (:.)(..) >> , Acc >> , Vector >> , Scalar >> , Elt >> , fold >> , slice >> , constant >> , Array >> , Z(..), DIM1, DIM2 >> , fromList >> , All(..) >> , generate >> , lift, unlift >> , shape >> ) >> import Data.Array.Accelerate.Interpreter ( run ) >> >> dotP :: (Num a, Elt a) => Acc (Vector a) -> Acc (Vector a) -> Acc (Scalar >> a) >> dotP xs ys = fold (+) 0 $ A.zipWith (*) xs ys >> >> type Matrix a = Array DIM2 a >> >> getRow :: Elt a => Int -> Acc (Matrix a) -> Acc (Vector a) >> getRow n mat = slice mat . constant $ Z :. n :. All >> >> -- Naive matrix multiplication: >> -- >> -- index (i, j) is equal to the ith row of 'a' `dot` the jth row of 'b' >> matMul :: A.Acc (Matrix Double) -> A.Acc (Matrix Double) -> A.Acc (Matrix >> Double) >> matMul a b' = A.generate (constant $ Z :. nrows :. ncols) $ >> \ix -> >> let (Z :. i :. j) = unlift ix >> in getRow i a `dotP` getRow j b >> where >> b = A.transpose b' -- I assume row indexing is faster than column >> indexing... >> (Z :. nrows :. _ ) = unlift $ shape a >> (Z :. _ :. ncols) = unlift $ shape b >> >> >> This, of course, gives me errors right now because I'm calling getRow and >> dotP from within the generation function, which expects Exp[ression]s, not >> Acc[elerated computation]s. >> >> So maybe I need to replace that line with an inner for loop? Is there an >> easy way to do that with Accelerate? >> >> Thanks for your help, >> - Clark >> >> [1] http://hackage.haskell.org/package/accelerate >> _______________________________________________ >> Haskell-Cafe mailing list >> Haskell-Cafe@haskell.org >> http://www.haskell.org/mailman/listinfo/haskell-cafe >> >> >> >> _______________________________________________ >> Haskell-Cafe mailing list >> Haskell-Cafe@haskell.org >> http://www.haskell.org/mailman/listinfo/haskell-cafe >> >> > > > _______________________________________________ > Haskell-Cafe mailing list > Haskell-Cafe@haskell.org > http://www.haskell.org/mailman/listinfo/haskell-cafe > >
_______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
