GPGPU Programming in Haskell with Accelerate

Video of the presentation: http://youtu.be/ARqE4yT2Z0o
Workshop presentation on using Accelerate effectively: https://speakerdeck.com/tmcdonell/gpgpu-programming-in-haskell-with-accelerate-workshop

Current graphics cards are massively parallel multicore processors optimised for workloads with a large degree of SIMD parallelism. Peak performance of these devices is far greater than that of traditional CPUs, however this is difficult to realise because good performance requires highly idiomatic programs, whose development is work intensive and requires expert knowledge. To raise the level of abstraction we are developing a domain-specific high-level language in Haskell for programming these devices. Computations are expressed in the form of parameterised collective operations —such as maps, reductions, and permutations— over multi-dimensional arrays. These computations are online compiled and executed on the graphics processor.

In this talk, I will introduce the Accelerate project; the language and its embedding in Haskell, as well as the online code generator and runtime system targeting CUDA GPUs.

Transcript

1. GPGPU Programming in Haskell
   with Accelerate
   Trevor L. McDonell
   
   University of New South Wales
   
   !
   @tlmcdonell
   
   tmcdonell@cse.unsw.edu.au
   
   !
   https://github.com/AccelerateHS
   

   View Slide

2. What is GPGPU Programming?
   • General Purpose Programming on Graphics Processing Units (GPUs)
   
   • Using your graphics card for something other than playing games
   
   • GPUs have many more cores than a CPU
   
   - GeForce GTX Titan
   
   - 2688 cores @ 837 MHz
   
   - 6 GB memory @ 288 GB/s
   

   View Slide

3. What is GPGPU Programming?
   • Main differences:
   
   - Single program multiple data (SPMD / SIMD), or just data-parallelism
   
   - All the cores run the same program, but on different data
   
   • We can’t program these in the same way as a CPU
   
   - Different instruction sets: can’t run a Haskell program directly
   
   - More restrictive hardware designs, limited control structures
   
   • GPUs have their own memory
   
   - Data has to be explicitly moved back and forth
   

   View Slide

4. Dot product four ways
   • Dot-product: pair-wise multiply two arrays and sum the result
   
   • C (sequential):
   float  dotp(float  *xs,  float  *ys,  int  size)  
   {  
          int      i;  
          float  sum  =  0;  
   !
          for  (i  =  0;  i  <  size;  ++i)  {  
                  sum  =  sum  +  xs[i]  *  ys[i];  
          }  
   !
          return  sum;  
   }
   

   View Slide

5. Dot product four ways
   • Haskell (sequential):
   dotp  ::  [Float]  -­‐>  [Float]  -­‐>  Float  
   dotp  xs  ys  =  foldl  (+)  0  
                        (  zipWith  (*)  xs  ys  )
   

   View Slide

6. Dot product four ways
   • Haskell (sequential):
   - [Float] is a list of floating point numbers
   dotp  ::  [Float]  -­‐>  [Float]  -­‐>  Float  
   dotp  xs  ys  =  foldl  (+)  0  
                        (  zipWith  (*)  xs  ys  )
   

   View Slide

7. Dot product four ways
   • Haskell (sequential):
   - [Float] is a list of floating point numbers
   - zipWith applies the function (*) element-wise to the two input lists
   dotp  ::  [Float]  -­‐>  [Float]  -­‐>  Float  
   dotp  xs  ys  =  foldl  (+)  0  
                        (  zipWith  (*)  xs  ys  )
   

   View Slide

8. Dot product four ways
   • Haskell (sequential):
   - [Float] is a list of floating point numbers
   - zipWith applies the function (*) element-wise to the two input lists
   - foldl sums the result of zipWith
   dotp  ::  [Float]  -­‐>  [Float]  -­‐>  Float  
   dotp  xs  ys  =  foldl  (+)  0  
                        (  zipWith  (*)  xs  ys  )
   

   View Slide

9. Dot product four ways
   • CUDA (parallel):
   - Step 1: element-wise multiplication
   __global__  
   void  zipWithMult(float  *xs,  float  *ys,  float  *zs,  int  size)  
   {  
          int  i  =  blockDim.x  *  blockIdx.x  +  threadIdx.x;  
   !
          if  (i  <  size)  {  
                  zs[i]  =  xs[i]  *  ys[i];  
          }  
   }
   

   View Slide

10. Dot product four ways
    • CUDA (parallel):
    - Step 1: element-wise multiplication
    - __global__ indicates this is a function that runs on the GPU
    __global__  
    void  zipWithMult(float  *xs,  float  *ys,  float  *zs,  int  size)  
    {  
           int  i  =  blockDim.x  *  blockIdx.x  +  threadIdx.x;  
    !
           if  (i  <  size)  {  
                   zs[i]  =  xs[i]  *  ys[i];  
           }  
    }
    

    View Slide

11. Dot product four ways
    • CUDA (parallel):
    - Step 1: element-wise multiplication
    - __global__ indicates this is a function that runs on the GPU
    - spawn one thread for each element in the vector: unique for each thread
    __global__  
    void  zipWithMult(float  *xs,  float  *ys,  float  *zs,  int  size)  
    {  
           int  i  =  blockDim.x  *  blockIdx.x  +  threadIdx.x;  
    !
           if  (i  <  size)  {  
                   zs[i]  =  xs[i]  *  ys[i];  
           }  
    }
    

    View Slide

12. Dot product four ways
    • CUDA (parallel):
    
    - Step 2: vector reduction … is somewhat complex
    struct  SharedMemory  
    {  
           __device__  inline  operator  float  *()  
           {  
                   extern  __shared__  int  __smem[];  
                   return  (float  *)__smem;  
           }  
    !
           __device__  inline  operator  const  float  *()  const  
           {  
                   extern  __shared__  int  __smem[];  
                   return  (float  *)__smem;  
           }  
    };  
    !
    template    
    __global__  void  
    reduce_kernel(float  *g_idata,  float  *g_odata,  unsigned  int  n)  
    {  
           float  *sdata  =  SharedMemory();  
    !
           unsigned  int  tid            =  threadIdx.x;  
           unsigned  int  i                =  blockIdx.x*blockSize*2  +  threadIdx.x;  
           unsigned  int  gridSize  =  blockSize*2*gridDim.x;  
    !
           float  sum  =  0;  
    !
           while  (i  <  n)  {  
                   sum  +=  g_idata[i];  
    !
                   if  (nIsPow2  ||  i  +  blockSize  <  n)  
                           sum  +=  g_idata[i+blockSize];  
    !
                   i  +=  gridSize;  
           }  
    !
           sdata[tid]  =  sum;  
           __syncthreads();  
    !
           if  (blockSize  >=  512)  {  
                   if  (tid  <  256)  {  
                           sdata[tid]  =  sum  =  sum  +  sdata[tid  +  256];  
                   }  
    !
                   __syncthreads();  
           }  
    !
           if  (blockSize  >=  256)  {  
                   if  (tid  <  128)  {  
                           sdata[tid]  =  sum  =  sum  +  sdata[tid  +  128];  
                   }  
    !
                   __syncthreads();  
           }  
    !
           if  (blockSize  >=  128)  {  
                   if  (tid  <    64)  {  
                           sdata[tid]  =  sum  =  sum  +  sdata[tid  +    64];  
                   }  
    !
                   __syncthreads();  
           }  
    !
           if  (tid  <  32)  
           {  
                   volatile  float  *smem  =  sdata;  
    !
                   if  (blockSize  >=    64)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +  32];  
                   }  
    !
                   if  (blockSize  >=    32)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +  16];  
                   }  
    !
                   if  (blockSize  >=    16)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +    8];  
                   }  
    !
                   if  (blockSize  >=      8)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +    4];  
                   }  
    !
                   if  (blockSize  >=      4)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +    2];  
                   }  
    !
                   if  (blockSize  >=      2)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +    1];  
                   }  
           }  
    !
           if  (tid  ==  0)  
                   g_odata[blockIdx.x]  =  sdata[0];  
    }  
    !
    void  getNumBlocksAndThreads(int  n,  int  maxBlocks,  int  maxThreads,  int  &blocks,  int  &threads)  
    {  
           cudaDeviceProp  prop;  
           int  device;  
           checkCudaErrors(cudaGetDevice(&device));  
           checkCudaErrors(cudaGetDeviceProperties(&prop,  device));  
    !
           threads  =  (n  <  maxThreads*2)  ?  nextPow2((n  +  1)/  2)  :  maxThreads;  
           blocks    =  (n  +  (threads  *  2  -­‐  1))  /  (threads  *  2);  
    !
           if  (blocks  >  prop.maxGridSize[0])  
           {  
                   blocks    /=  2;  
                   threads  *=  2;  
           }  
    !
           blocks  =  min(maxBlocks,  blocks);  
    }  
    !
    float  
    reduce(int  n,  float  *d_idata,  float  *d_odata)  
    {  
           int  threads        =  0;  
           int  blocks          =  0;  
           int  maxThreads  =  256;  
           int  maxBlocks    =  64;  
           int  size              =  n  
    !
           while  (size  >  1)  {  
                   getNumBlocksAndThreads(size,  maxBlocks,  maxThreads,  blocks,  threads);  
    !
                   int  smemSize  =  (threads  <=  32)  ?  2  *  threads  *  sizeof(float)  :  threads  *  sizeof(float);  
                   dim3  dimBlock(threads,  1,  1);  
                   dim3  dimGrid(blocks,  1,  1);  
    !
                   if  (isPow2(size))  
                   {  
                           switch  (threads)  
                           {  
                                   case  512:  
                                           reduce_kernel<512,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  256:  
                                           reduce_kernel<256,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  128:  
                                           reduce_kernel<128,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  64:  
                                           reduce_kernel<  64,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  32:  
                                           reduce_kernel<  32,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  16:  
                                           reduce_kernel<  16,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    8:  
                                           reduce_kernel<    8,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    4:  
                                           reduce_kernel<    4,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    2:  
                                           reduce_kernel<    2,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    1:  
                                           reduce_kernel<    1,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                           }  
                   }  
                   else  
                   {  
                           switch  (threads)  
                           {  
                                   case  512:  
                                           reduce_kernel<512,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  256:  
                                           reduce_kernel<256,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  128:  
                                           reduce_kernel<128,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  64:  
                                           reduce_kernel<  64,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  32:  
                                           reduce_kernel<  32,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  16:  
                                           reduce_kernel<  16,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    8:  
                                           reduce_kernel<    8,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    4:  
                                           reduce_kernel<    4,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    2:  
                                           reduce_kernel<    2,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    1:  
                                           reduce_kernel<    1,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                           }  
                   }  
    !
                   size  =  (size  +  (threads*2-­‐1))  /  (threads*2);  
           }  
    !
           float  sum;  
           checkCudaErrors(cudaMemcpy(&sum,  d_odata,  sizeof(float),  cudaMemcpyDeviceToHost));  
    !
    }  
    !
    

    View Slide

13. Dot product four ways
    • CUDA (parallel):
    
    - Step 2: vector reduction … is somewhat complex
    struct  SharedMemory  
    {  
           __device__  inline  operator  float  *()  
           {  
                   extern  __shared__  int  __smem[];  
                   return  (float  *)__smem;  
           }  
    !
           __device__  inline  operator  const  float  *()  const  
           {  
                   extern  __shared__  int  __smem[];  
                   return  (float  *)__smem;  
           }  
    };  
    !
    template    
    __global__  void  
    reduce_kernel(float  *g_idata,  float  *g_odata,  unsigned  int  n)  
    {  
           float  *sdata  =  SharedMemory();  
    !
           unsigned  int  tid            =  threadIdx.x;  
           unsigned  int  i                =  blockIdx.x*blockSize*2  +  threadIdx.x;  
           unsigned  int  gridSize  =  blockSize*2*gridDim.x;  
    !
           float  sum  =  0;  
    !
           while  (i  <  n)  {  
                   sum  +=  g_idata[i];  
    !
                   if  (nIsPow2  ||  i  +  blockSize  <  n)  
                           sum  +=  g_idata[i+blockSize];  
    !
                   i  +=  gridSize;  
           }  
    !
           sdata[tid]  =  sum;  
           __syncthreads();  
    !
           if  (blockSize  >=  512)  {  
                   if  (tid  <  256)  {  
                           sdata[tid]  =  sum  =  sum  +  sdata[tid  +  256];  
                   }  
    !
                   __syncthreads();  
           }  
    !
           if  (blockSize  >=  256)  {  
                   if  (tid  <  128)  {  
                           sdata[tid]  =  sum  =  sum  +  sdata[tid  +  128];  
                   }  
    !
                   __syncthreads();  
           }  
    !
           if  (blockSize  >=  128)  {  
                   if  (tid  <    64)  {  
                           sdata[tid]  =  sum  =  sum  +  sdata[tid  +    64];  
                   }  
    !
                   __syncthreads();  
           }  
    !
           if  (tid  <  32)  
           {  
                   volatile  float  *smem  =  sdata;  
    !
                   if  (blockSize  >=    64)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +  32];  
                   }  
    !
                   if  (blockSize  >=    32)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +  16];  
                   }  
    !
                   if  (blockSize  >=    16)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +    8];  
                   }  
    !
                   if  (blockSize  >=      8)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +    4];  
                   }  
    !
                   if  (blockSize  >=      4)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +    2];  
                   }  
    !
                   if  (blockSize  >=      2)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +    1];  
                   }  
           }  
    !
           if  (tid  ==  0)  
                   g_odata[blockIdx.x]  =  sdata[0];  
    }  
    !
    void  getNumBlocksAndThreads(int  n,  int  maxBlocks,  int  maxThreads,  int  &blocks,  int  &threads)  
    {  
           cudaDeviceProp  prop;  
           int  device;  
           checkCudaErrors(cudaGetDevice(&device));  
           checkCudaErrors(cudaGetDeviceProperties(&prop,  device));  
    !
           threads  =  (n  <  maxThreads*2)  ?  nextPow2((n  +  1)/  2)  :  maxThreads;  
           blocks    =  (n  +  (threads  *  2  -­‐  1))  /  (threads  *  2);  
    !
           if  (blocks  >  prop.maxGridSize[0])  
           {  
                   blocks    /=  2;  
                   threads  *=  2;  
           }  
    !
           blocks  =  min(maxBlocks,  blocks);  
    }  
    !
    float  
    reduce(int  n,  float  *d_idata,  float  *d_odata)  
    {  
           int  threads        =  0;  
           int  blocks          =  0;  
           int  maxThreads  =  256;  
           int  maxBlocks    =  64;  
           int  size              =  n  
    !
           while  (size  >  1)  {  
                   getNumBlocksAndThreads(size,  maxBlocks,  maxThreads,  blocks,  threads);  
    !
                   int  smemSize  =  (threads  <=  32)  ?  2  *  threads  *  sizeof(float)  :  threads  *  sizeof(float);  
                   dim3  dimBlock(threads,  1,  1);  
                   dim3  dimGrid(blocks,  1,  1);  
    !
                   if  (isPow2(size))  
                   {  
                           switch  (threads)  
                           {  
                                   case  512:  
                                           reduce_kernel<512,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  256:  
                                           reduce_kernel<256,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  128:  
                                           reduce_kernel<128,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  64:  
                                           reduce_kernel<  64,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  32:  
                                           reduce_kernel<  32,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  16:  
                                           reduce_kernel<  16,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    8:  
                                           reduce_kernel<    8,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    4:  
                                           reduce_kernel<    4,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    2:  
                                           reduce_kernel<    2,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    1:  
                                           reduce_kernel<    1,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                           }  
                   }  
                   else  
                   {  
                           switch  (threads)  
                           {  
                                   case  512:  
                                           reduce_kernel<512,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  256:  
                                           reduce_kernel<256,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  128:  
                                           reduce_kernel<128,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  64:  
                                           reduce_kernel<  64,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  32:  
                                           reduce_kernel<  32,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  16:  
                                           reduce_kernel<  16,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    8:  
                                           reduce_kernel<    8,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    4:  
                                           reduce_kernel<    4,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    2:  
                                           reduce_kernel<    2,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    1:  
                                           reduce_kernel<    1,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                           }  
                   }  
    !
                   size  =  (size  +  (threads*2-­‐1))  /  (threads*2);  
           }  
    !
           float  sum;  
           checkCudaErrors(cudaMemcpy(&sum,  d_odata,  sizeof(float),  cudaMemcpyDeviceToHost));  
    !
    }  
    !
    

    View Slide

14. Dot product four ways
    • CUDA (parallel):
    
    - Step 2: vector reduction … is somewhat complex
    struct  SharedMemory  
    {  
           __device__  inline  operator  float  *()  
           {  
                   extern  __shared__  int  __smem[];  
                   return  (float  *)__smem;  
           }  
    !
           __device__  inline  operator  const  float  *()  const  
           {  
                   extern  __shared__  int  __smem[];  
                   return  (float  *)__smem;  
           }  
    };  
    !
    template    
    __global__  void  
    reduce_kernel(float  *g_idata,  float  *g_odata,  unsigned  int  n)  
    {  
           float  *sdata  =  SharedMemory();  
    !
           unsigned  int  tid            =  threadIdx.x;  
           unsigned  int  i                =  blockIdx.x*blockSize*2  +  threadIdx.x;  
           unsigned  int  gridSize  =  blockSize*2*gridDim.x;  
    !
           float  sum  =  0;  
    !
           while  (i  <  n)  {  
                   sum  +=  g_idata[i];  
    !
                   if  (nIsPow2  ||  i  +  blockSize  <  n)  
                           sum  +=  g_idata[i+blockSize];  
    !
                   i  +=  gridSize;  
           }  
    !
           sdata[tid]  =  sum;  
           __syncthreads();  
    !
           if  (blockSize  >=  512)  {  
                   if  (tid  <  256)  {  
                           sdata[tid]  =  sum  =  sum  +  sdata[tid  +  256];  
                   }  
    !
                   __syncthreads();  
           }  
    !
           if  (blockSize  >=  256)  {  
                   if  (tid  <  128)  {  
                           sdata[tid]  =  sum  =  sum  +  sdata[tid  +  128];  
                   }  
    !
                   __syncthreads();  
           }  
    !
           if  (blockSize  >=  128)  {  
                   if  (tid  <    64)  {  
                           sdata[tid]  =  sum  =  sum  +  sdata[tid  +    64];  
                   }  
    !
                   __syncthreads();  
           }  
    !
           if  (tid  <  32)  
           {  
                   volatile  float  *smem  =  sdata;  
    !
                   if  (blockSize  >=    64)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +  32];  
                   }  
    !
                   if  (blockSize  >=    32)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +  16];  
                   }  
    !
                   if  (blockSize  >=    16)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +    8];  
                   }  
    !
                   if  (blockSize  >=      8)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +    4];  
                   }  
    !
                   if  (blockSize  >=      4)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +    2];  
                   }  
    !
                   if  (blockSize  >=      2)  {  
                           smem[tid]  =  sum  =  sum  +  smem[tid  +    1];  
                   }  
           }  
    !
           if  (tid  ==  0)  
                   g_odata[blockIdx.x]  =  sdata[0];  
    }  
    !
    void  getNumBlocksAndThreads(int  n,  int  maxBlocks,  int  maxThreads,  int  &blocks,  int  &threads)  
    {  
           cudaDeviceProp  prop;  
           int  device;  
           checkCudaErrors(cudaGetDevice(&device));  
           checkCudaErrors(cudaGetDeviceProperties(&prop,  device));  
    !
           threads  =  (n  <  maxThreads*2)  ?  nextPow2((n  +  1)/  2)  :  maxThreads;  
           blocks    =  (n  +  (threads  *  2  -­‐  1))  /  (threads  *  2);  
    !
           if  (blocks  >  prop.maxGridSize[0])  
           {  
                   blocks    /=  2;  
                   threads  *=  2;  
           }  
    !
           blocks  =  min(maxBlocks,  blocks);  
    }  
    !
    float  
    reduce(int  n,  float  *d_idata,  float  *d_odata)  
    {  
           int  threads        =  0;  
           int  blocks          =  0;  
           int  maxThreads  =  256;  
           int  maxBlocks    =  64;  
           int  size              =  n  
    !
           while  (size  >  1)  {  
                   getNumBlocksAndThreads(size,  maxBlocks,  maxThreads,  blocks,  threads);  
    !
                   int  smemSize  =  (threads  <=  32)  ?  2  *  threads  *  sizeof(float)  :  threads  *  sizeof(float);  
                   dim3  dimBlock(threads,  1,  1);  
                   dim3  dimGrid(blocks,  1,  1);  
    !
                   if  (isPow2(size))  
                   {  
                           switch  (threads)  
                           {  
                                   case  512:  
                                           reduce_kernel<512,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  256:  
                                           reduce_kernel<256,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  128:  
                                           reduce_kernel<128,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  64:  
                                           reduce_kernel<  64,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  32:  
                                           reduce_kernel<  32,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  16:  
                                           reduce_kernel<  16,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    8:  
                                           reduce_kernel<    8,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    4:  
                                           reduce_kernel<    4,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    2:  
                                           reduce_kernel<    2,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    1:  
                                           reduce_kernel<    1,  true><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                           }  
                   }  
                   else  
                   {  
                           switch  (threads)  
                           {  
                                   case  512:  
                                           reduce_kernel<512,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  256:  
                                           reduce_kernel<256,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  128:  
                                           reduce_kernel<128,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  64:  
                                           reduce_kernel<  64,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  32:  
                                           reduce_kernel<  32,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case  16:  
                                           reduce_kernel<  16,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    8:  
                                           reduce_kernel<    8,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    4:  
                                           reduce_kernel<    4,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    2:  
                                           reduce_kernel<    2,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                                   case    1:  
                                           reduce_kernel<    1,  false><<<  dimGrid,  dimBlock,  smemSize  >>>(d_idata,  d_odata,  size);  
                                           break;  
                           }  
                   }  
    !
                   size  =  (size  +  (threads*2-­‐1))  /  (threads*2);  
           }  
    !
           float  sum;  
           checkCudaErrors(cudaMemcpy(&sum,  d_odata,  sizeof(float),  cudaMemcpyDeviceToHost));  
    !
    }  
    !
    o_O
    

    View Slide

15. Dot product four ways
    • Accelerate (parallel):
    dotp  ::  Acc  (Vector  Float)  -­‐>  Acc  (Vector  Float)  -­‐>  Acc  (Scalar  Float)  
    dotp  xs  ys  =  fold  (+)  0  
                         (  zipWith  (*)  xs  ys  )
    

    View Slide

16. Dot product four ways
    • Accelerate (parallel):
    • Recall the sequential Haskell version:
    dotp  ::  Acc  (Vector  Float)  -­‐>  Acc  (Vector  Float)  -­‐>  Acc  (Scalar  Float)  
    dotp  xs  ys  =  fold  (+)  0  
                         (  zipWith  (*)  xs  ys  )
    dotp  ::  [Float]  -­‐>  [Float]  -­‐>  Float  
    dotp  xs  ys  =  foldl  (+)  0  
                         (  zipWith  (*)  xs  ys  )
    

    View Slide

17. Dot product four ways
    • Accelerate (parallel):
    • Recall the sequential Haskell version:
    dotp  ::  Acc  (Vector  Float)  -­‐>  Acc  (Vector  Float)  -­‐>  Acc  (Scalar  Float)  
    dotp  xs  ys  =  fold  (+)  0  
                         (  zipWith  (*)  xs  ys  )
    dotp  ::  [Float]  -­‐>  [Float]  -­‐>  Float  
    dotp  xs  ys  =  foldl  (+)  0  
                         (  zipWith  (*)  xs  ys  )
    

    View Slide

18. Dot product four ways
    • Accelerate (parallel):
    • Recall the sequential Haskell version:
    dotp  ::  Acc  (Vector  Float)  -­‐>  Acc  (Vector  Float)  -­‐>  Acc  (Scalar  Float)  
    dotp  xs  ys  =  fold  (+)  0  
                         (  zipWith  (*)  xs  ys  )
    dotp  ::  [Float]  -­‐>  [Float]  -­‐>  Float  
    dotp  xs  ys  =  foldl  (+)  0  
                         (  zipWith  (*)  xs  ys  )
    

    View Slide

19. Dot product four ways
    • Accelerate (parallel):
    • Recall the sequential Haskell version:
    dotp  ::  Acc  (Vector  Float)  -­‐>  Acc  (Vector  Float)  -­‐>  Acc  (Scalar  Float)  
    dotp  xs  ys  =  fold  (+)  0  
                         (  zipWith  (*)  xs  ys  )
    dotp  ::  [Float]  -­‐>  [Float]  -­‐>  Float  
    dotp  xs  ys  =  foldl  (+)  0  
                         (  zipWith  (*)  xs  ys  )
    left-to-right
    traversal
    

    View Slide

20. Dot product four ways
    • Accelerate (parallel):
    • Recall the sequential Haskell version:
    dotp  ::  Acc  (Vector  Float)  -­‐>  Acc  (Vector  Float)  -­‐>  Acc  (Scalar  Float)  
    dotp  xs  ys  =  fold  (+)  0  
                         (  zipWith  (*)  xs  ys  )
    dotp  ::  [Float]  -­‐>  [Float]  -­‐>  Float  
    dotp  xs  ys  =  foldl  (+)  0  
                         (  zipWith  (*)  xs  ys  )
    left-to-right
    traversal
    neither left nor right:
    happens in parallel (tree-like)
    

    View Slide

21. Dot product four ways
    • Accelerate (parallel):
    • Recall the sequential Haskell version:
    dotp  ::  Acc  (Vector  Float)  -­‐>  Acc  (Vector  Float)  -­‐>  Acc  (Scalar  Float)  
    dotp  xs  ys  =  fold  (+)  0  
                         (  zipWith  (*)  xs  ys  )
    dotp  ::  [Float]  -­‐>  [Float]  -­‐>  Float  
    dotp  xs  ys  =  foldl  (+)  0  
                         (  zipWith  (*)  xs  ys  )
    left-to-right
    traversal
    neither left nor right:
    happens in parallel (tree-like)
    But… how does it perform?
    

    View Slide

22. Dot product four ways
    0.1
    1
    10
    100
    2 4 6 8 10 12 14 16 18 20
    Run Time (ms)
    Elements (millions)
    Dot product
    sequential
    Accelerate
    CUBLAS
    Tesla T10 (240 cores @ 1.3 GHz) vs. Xenon E5405 (2GHz)
    

    View Slide

23. Dot product four ways
    0.1
    1
    10
    100
    2 4 6 8 10 12 14 16 18 20
    Run Time (ms)
    Elements (millions)
    Dot product
    sequential
    Accelerate
    CUBLAS
    1.2x
    Tesla T10 (240 cores @ 1.3 GHz) vs. Xenon E5405 (2GHz)
    

    View Slide

24. Dot product four ways
    0.1
    1
    10
    100
    2 4 6 8 10 12 14 16 18 20
    Run Time (ms)
    Elements (millions)
    Dot product
    sequential
    Accelerate
    CUBLAS
    30x
    1.2x
    Tesla T10 (240 cores @ 1.3 GHz) vs. Xenon E5405 (2GHz)
    

    View Slide

25. View Slide

26. Mandelbrot fractal
    

    View Slide

27. Mandelbrot fractal
    n-body gravitational simulation
    

    View Slide

28. Mandelbrot fractal
    n-body gravitational simulation
    Canny edge detection
    

    View Slide

29. Mandelbrot fractal
    n-body gravitational simulation
    Canny edge detection
    SmoothLife cellular automata
    

    View Slide

30. Mandelbrot fractal
    n-body gravitational simulation
    Canny edge detection
    SmoothLife cellular automata
    stable fluid flow
    

    View Slide

31. Mandelbrot fractal
    n-body gravitational simulation
    Canny edge detection
    SmoothLife cellular automata
    stable fluid flow
    ...  
    d6b821d937a4170b3c4f8ad93495575d:  saitek1  
    d0e52829bf7962ee0aa90550ffdcccaa:  laura1230  
    494a8204b800c41b2da763f9bbbcc462:  lina03  
    d8ff07c52a95b30800809758f84ce28c:  Jenny10  
    e81bed02faa9892f8360c705241191ae:  carmen89  
    46f7d75718029de99dd81fd907034bc9:  mellon22  
    0dd3c176cf34486ec00b526b6920b782:  helena04  
    9351c4bc8c8ba17b58d5a6a1f839f356:  85548554  
    9c36c5599f40d08f874559ac824d091a:  585123456  
    4b4dce6c91b429e8360aa65f97342e90:  5678go  
    3aa561d4c17d9d58443fc15d10cc86ae:  momo55  
    !
    Recovered  150/1000  (15.00  %)  digests  in  59.45  s,  185.03  MHash/sec
    Password “recovery” (MD5 dictionary attack)
    

    View Slide

32. Accelerate
    • Accelerate is a Domain-Specific Language for GPU programming
    Haskell/Accelerate
    program
    CUDA code
    Compile with NVIDIA’s
    compiler & load onto the GPU
    Copy result back to Haskell
    Transform Accelerate
    program into CUDA program
    

    View Slide

33. Accelerate
    • Accelerate is a Domain-Specific Language for GPU programming
    
    - This process may happen several times during program execution
    
    - Code and data fragments get cached and reused
    
    • An Accelerate program is a Haskell program that generates a CUDA program
    
    - However, in many respects this still looks like a Haskell program
    
    - Shares various concepts with Repa, a Haskell array library for CPUs
    

    View Slide

34. Accelerate
    • To execute an Accelerate computation (on the GPU):
    - run comes from whichever backend we have chosen (CUDA)
    run  ::  Arrays  a  =>  Acc  a  -­‐>  a
    

    View Slide

35. Accelerate
    • To execute an Accelerate computation (on the GPU):
    - run comes from whichever backend we have chosen (CUDA)
    - Arrays constrains the result to be an Array, or tuple thereof
    run  ::  Arrays  a  =>  Acc  a  -­‐>  a
    

    View Slide

36. Accelerate
    • To execute an Accelerate computation (on the GPU):
    - run comes from whichever backend we have chosen (CUDA)
    - Arrays constrains the result to be an Array, or tuple thereof
    • What is Acc?
    run  ::  Arrays  a  =>  Acc  a  -­‐>  a
    

    View Slide

37. Accelerate
    • To execute an Accelerate computation (on the GPU):
    - run comes from whichever backend we have chosen (CUDA)
    - Arrays constrains the result to be an Array, or tuple thereof
    • What is Acc?
    - This is our DSL type
    - A data structure (Abstract Syntax Tree) representing a computation that
    once executed will yield a result of type ‘a’
    run  ::  Arrays  a  =>  Acc  a  -­‐>  a
    

    View Slide

38. Accelerate
    • To execute an Accelerate computation (on the GPU):
    - run comes from whichever backend we have chosen (CUDA)
    - Arrays constrains the result to be an Array, or tuple thereof
    • What is Acc?
    - This is our DSL type
    - A data structure (Abstract Syntax Tree) representing a computation that
    once executed will yield a result of type ‘a’
    run  ::  Arrays  a  =>  Acc  a  -­‐>  a
    Accelerate is a library of
    collective operations over
    arrays of type Acc  a
    

    View Slide

39. Accelerate
    • Accelerate computations take place on arrays
    - Parallelism is introduced in the form of collective operations over arrays
    Accelerate
    computation
    Arrays in Arrays out
    run  ::  Arrays  a  =>  Acc  a  -­‐>  a
    

    View Slide

40. Accelerate
    • Accelerate computations take place on arrays
    - Parallelism is introduced in the form of collective operations over arrays
    • Arrays have two type parameters
    data  Array  sh  e
    Accelerate
    computation
    Arrays in Arrays out
    run  ::  Arrays  a  =>  Acc  a  -­‐>  a
    

    View Slide

41. Accelerate
    • Accelerate computations take place on arrays
    - Parallelism is introduced in the form of collective operations over arrays
    • Arrays have two type parameters
    - The shape of the array, or dimensionality
    data  Array  sh  e
    Accelerate
    computation
    Arrays in Arrays out
    run  ::  Arrays  a  =>  Acc  a  -­‐>  a
    

    View Slide

42. Accelerate
    • Accelerate computations take place on arrays
    - Parallelism is introduced in the form of collective operations over arrays
    • Arrays have two type parameters
    - The shape of the array, or dimensionality
    - The element type of the array: Int, Float, etc.
    data  Array  sh  e
    Accelerate
    computation
    Arrays in Arrays out
    run  ::  Arrays  a  =>  Acc  a  -­‐>  a
    

    View Slide

43. Arrays
    • Supported array element types are members of the Elt class:
    
    - ()
    
    - Int, Int32, Int64, Word, Word32, Word64...
    
    - Float, Double
    
    - Char
    
    - Bool
    
    - Tuples up to 9-tuples of these, including nested tuples
    
    • Note that Array itself is not an allowable element type. There are no nested
    arrays in Accelerate, regular arrays only!
    data  Array  sh  e
    

    View Slide

44. Accelerate by example
    Mandelbrot fractal
    

    View Slide

45. Mandelbrot set generator
    • Basics
    
    - Pick a window onto the complex plane & divide into pixels
    
    - A point is in the set if the value of does not diverge to infinity
    
    !
    - Each pixel has a value given by its coordinates in the complex plane
    
    - Colour depends on number of iterations before divergence
    
    • Each pixel is independent: lots of data parallelism
    c
    n
    |z|
    zn+1 = c + z2
    n
    

    View Slide

46. Mandelbrot set generator
    • First, some types:
    - A pair of floating point numbers for the real and imaginary parts
    data  Complex     =  (Float,  Float)  
    data  Array  sh  e
    

    View Slide

47. Mandelbrot set generator
    • First, some types:
    - A pair of floating point numbers for the real and imaginary parts
    - DIM2 is a type synonym for a two dimensional Shape
    data  Complex     =  (Float,  Float)  
    data  ComplexPlane  
    =  Array  DIM2  Complex
    data  Array  sh  e
    

    View Slide

48. Shapes
    • Shapes determines the dimensions of the array and the type of the index
    
    - Z represents a rank-zero array (singleton array with one element)
    
    - (:.) increases the rank by adding a new dimension on the right
    data  Z                        =  Z  
    data  tail  :.  head  =  tail  :.  head
    

    View Slide

49. Shapes
    • Shapes determines the dimensions of the array and the type of the index
    
    - Z represents a rank-zero array (singleton array with one element)
    
    - (:.) increases the rank by adding a new dimension on the right
    data  Z                        =  Z  
    data  tail  :.  head  =  tail  :.  head
    

    View Slide

50. Shapes
    • Shapes determines the dimensions of the array and the type of the index
    
    - Z represents a rank-zero array (singleton array with one element)
    
    - (:.) increases the rank by adding a new dimension on the right
    data  Z                        =  Z  
    data  tail  :.  head  =  tail  :.  head
    

    View Slide

51. Shapes
    • Shapes determines the dimensions of the array and the type of the index
    
    - Z represents a rank-zero array (singleton array with one element)
    
    - (:.) increases the rank by adding a new dimension on the right
    • Examples:
    
    - One-dimensional array (Vector) indexed by Int: (Z  :.  Int)
    
    - Two-dimensional array, indexed by Int: (Z  :.  Int  :.  Int)
    data  Z                        =  Z  
    data  tail  :.  head  =  tail  :.  head
    

    View Slide

52. Shapes
    • Shapes determines the dimensions of the array and the type of the index
    
    - Z represents a rank-zero array (singleton array with one element)
    
    - (:.) increases the rank by adding a new dimension on the right
    • Examples:
    
    - One-dimensional array (Vector) indexed by Int: (Z  :.  Int)
    
    - Two-dimensional array, indexed by Int: (Z  :.  Int  :.  Int)
    • This style is used at both the type and value level:
    data  Z                        =  Z  
    data  tail  :.  head  =  tail  :.  head
    sh  ::  Z  :.  Int  
    sh  =    Z  :.  10
    

    View Slide

53. Shapes
    • Shapes determines the dimensions of the array and the type of the index
    
    - Z represents a rank-zero array (singleton array with one element)
    
    - (:.) increases the rank by adding a new dimension on the right
    • Examples:
    
    - One-dimensional array (Vector) indexed by Int: (Z  :.  Int)
    
    - Two-dimensional array, indexed by Int: (Z  :.  Int  :.  Int)
    • This style is used at both the type and value level:
    data  Z                        =  Z  
    data  tail  :.  head  =  tail  :.  head
    sh  ::  Z  :.  Int  
    sh  =    Z  :.  10
    

    View Slide

54. Shapes
    • Shapes determines the dimensions of the array and the type of the index
    
    - Z represents a rank-zero array (singleton array with one element)
    
    - (:.) increases the rank by adding a new dimension on the right
    • Examples:
    
    - One-dimensional array (Vector) indexed by Int: (Z  :.  Int)
    
    - Two-dimensional array, indexed by Int: (Z  :.  Int  :.  Int)
    • This style is used at both the type and value level:
    data  Z                        =  Z  
    data  tail  :.  head  =  tail  :.  head
    sh  ::  Z  :.  Int  
    sh  =    Z  :.  10
    

    View Slide

55. Mandelbrot set generator
    • The initial complex plane:
    • generate is a collective operation that yields a value for an index in the array
    generate  ::  (Shape  sh,  Elt  a)  
                     =>  Exp  sh  
                     -­‐>  (Exp  sh  -­‐>  Exp  a)  
                     -­‐>  Acc  (Array  sh  a)
    z0 = c
    

    View Slide

56. Mandelbrot set generator
    • The initial complex plane:
    • generate is a collective operation that yields a value for an index in the array
    - Supported shape and element types: we will use DIM2 and Complex
    generate  ::  (Shape  sh,  Elt  a)  
                     =>  Exp  sh  
                     -­‐>  (Exp  sh  -­‐>  Exp  a)  
                     -­‐>  Acc  (Array  sh  a)
    z0 = c
    

    View Slide

57. Mandelbrot set generator
    • The initial complex plane:
    • generate is a collective operation that yields a value for an index in the array
    - Supported shape and element types: we will use DIM2 and Complex
    - Size of the result array: number of pixels in the image
    generate  ::  (Shape  sh,  Elt  a)  
                     =>  Exp  sh  
                     -­‐>  (Exp  sh  -­‐>  Exp  a)  
                     -­‐>  Acc  (Array  sh  a)
    z0 = c
    

    View Slide

58. Mandelbrot set generator
    • The initial complex plane:
    • generate is a collective operation that yields a value for an index in the array
    - Supported shape and element types: we will use DIM2 and Complex
    - Size of the result array: number of pixels in the image
    - A function to apply at every index: generate the values of at each pixel
    generate  ::  (Shape  sh,  Elt  a)  
                     =>  Exp  sh  
                     -­‐>  (Exp  sh  -­‐>  Exp  a)  
                     -­‐>  Acc  (Array  sh  a)
    c
    z0 = c
    

    View Slide

59. Mandelbrot set generator
    • The initial complex plane:
    • generate is a collective operation that yields a value for an index in the array
    - Supported shape and element types: we will use DIM2 and Complex
    - Size of the result array: number of pixels in the image
    - A function to apply at every index: generate the values of at each pixel
    generate  ::  (Shape  sh,  Elt  a)  
                     =>  Exp  sh  
                     -­‐>  (Exp  sh  -­‐>  Exp  a)  
                     -­‐>  Acc  (Array  sh  a)
    c
    z0 = c
    

    View Slide

60. Mandelbrot set generator
    • The initial complex plane:
    • generate is a collective operation that yields a value for an index in the array
    - Supported shape and element types: we will use DIM2 and Complex
    - Size of the result array: number of pixels in the image
    - A function to apply at every index: generate the values of at each pixel
    generate  ::  (Shape  sh,  Elt  a)  
                     =>  Exp  sh  
                     -­‐>  (Exp  sh  -­‐>  Exp  a)  
                     -­‐>  Acc  (Array  sh  a)
    c
    z0 = c
    If Acc is our DSL type,
    what is Exp?
    

    View Slide

61. A Stratified Language
    • Accelerate is split into two worlds: Acc and Exp
    
    - Acc represents collective operations over instances of Arrays
    
    - Exp is a scalar computation on things of type Elt
    

    View Slide

62. A Stratified Language
    • Accelerate is split into two worlds: Acc and Exp
    
    - Acc represents collective operations over instances of Arrays
    
    - Exp is a scalar computation on things of type Elt
    • Collective operations in Acc comprise many scalar operations in Exp,
    executed in parallel over Arrays
    
    - Scalar operations can not contain collective operations
    

    View Slide

63. A Stratified Language
    • Accelerate is split into two worlds: Acc and Exp
    
    - Acc represents collective operations over instances of Arrays
    
    - Exp is a scalar computation on things of type Elt
    • Collective operations in Acc comprise many scalar operations in Exp,
    executed in parallel over Arrays
    
    - Scalar operations can not contain collective operations
    • This stratification excludes nested data parallelism
    

    View Slide

64. Collective Operations
    • Collective operations comprise many scalar operations applied in parallel
    example  ::  Acc  (Vector  Int)  
    example  =  generate  (constant  (Z:.10))  
                                         (\ix  -­‐>  f  ix)
    

    View Slide

65. Collective Operations
    • Collective operations comprise many scalar operations applied in parallel
    - constant lifts a plain value into Exp land of scalar expressions
    example  ::  Acc  (Vector  Int)  
    example  =  generate  (constant  (Z:.10))  
                                         (\ix  -­‐>  f  ix)
    constant  ::  Elt  e  =>  e  -­‐>  Exp  e
    

    View Slide

66. Collective Operations
    • Collective operations comprise many scalar operations applied in parallel
    - constant lifts a plain value into Exp land of scalar expressions
    example  ::  Acc  (Vector  Int)  
    example  =  generate  (constant  (Z:.10))  
                                         (\ix  -­‐>  f  ix)
    generate  (constant  (Z:.10))  f
    constant  ::  Elt  e  =>  e  -­‐>  Exp  e
    

    View Slide

67. Collective Operations
    • Collective operations comprise many scalar operations applied in parallel
    - constant lifts a plain value into Exp land of scalar expressions
    example  ::  Acc  (Vector  Int)  
    example  =  generate  (constant  (Z:.10))  
                                         (\ix  -­‐>  f  ix)
    f  0 f  1 f  2 f  3 f  4 f  5 f  6 f  7 f  8 f  9
    generate  (constant  (Z:.10))  f
    f  ::  Exp  DIM1  -­‐>  Exp  Int
    constant  ::  Elt  e  =>  e  -­‐>  Exp  e
    

    View Slide

68. Collective Operations
    • Collective operations comprise many scalar operations applied in parallel
    - constant lifts a plain value into Exp land of scalar expressions
    example  ::  Acc  (Vector  Int)  
    example  =  generate  (constant  (Z:.10))  
                                         (\ix  -­‐>  f  ix)
    f  0 f  1 f  2 f  3 f  4 f  5 f  6 f  7 f  8 f  9
    Acc  (Vector  Int)
    constant  ::  Elt  e  =>  e  -­‐>  Exp  e
    

    View Slide

69. Mandelbrot set generator
    genPlane  ::  Exp  (Float,Float,Float,Float)  -­‐>  Acc  ComplexPlane  
    genPlane  view  =  
       generate  (constant  (Z:.600:.800))  
                         (\ix  -­‐>  let  Z  :.  y  :.  x  =  unlift  ix  
                                         in  lift  (  xmin  +  (fromIntegral  x  *  viewx)  /  800  
                                                         ,  ymin  +  (fromIntegral  y  *  viewy)  /  600  ))  
       where  
           (xmin,ymin,xmax,ymax)  =  unlift  view  
           viewx                                  =  xmax  -­‐  xmin  
           viewy                                  =  ymax  -­‐  ymin
    

    View Slide

70. Mandelbrot set generator
    genPlane  ::  Exp  (Float,Float,Float,Float)  -­‐>  Acc  ComplexPlane  
    genPlane  view  =  
       generate  (constant  (Z:.600:.800))  
                         (\ix  -­‐>  let  Z  :.  y  :.  x  =  unlift  ix  
                                         in  lift  (  xmin  +  (fromIntegral  x  *  viewx)  /  800  
                                                         ,  ymin  +  (fromIntegral  y  *  viewy)  /  600  ))  
       where  
           (xmin,ymin,xmax,ymax)  =  unlift  view  
           viewx                                  =  xmax  -­‐  xmin  
           viewy                                  =  ymax  -­‐  ymin
    

    View Slide

71. - unlift has a funky type, but just unpacks “stuff”. lift does the reverse.
    Mandelbrot set generator
    genPlane  ::  Exp  (Float,Float,Float,Float)  -­‐>  Acc  ComplexPlane  
    genPlane  view  =  
       generate  (constant  (Z:.600:.800))  
                         (\ix  -­‐>  let  Z  :.  y  :.  x  =  unlift  ix  
                                         in  lift  (  xmin  +  (fromIntegral  x  *  viewx)  /  800  
                                                         ,  ymin  +  (fromIntegral  y  *  viewy)  /  600  ))  
       where  
           (xmin,ymin,xmax,ymax)  =  unlift  view  
           viewx                                  =  xmax  -­‐  xmin  
           viewy                                  =  ymax  -­‐  ymin
    unlift  ::  Exp  (Z  :.  Int  :.  Int)  
                 -­‐>  Z  :.  Exp  Int  :.  Exp  Int
    

    View Slide

72. - unlift has a funky type, but just unpacks “stuff”. lift does the reverse.
    Mandelbrot set generator
    genPlane  ::  Exp  (Float,Float,Float,Float)  -­‐>  Acc  ComplexPlane  
    genPlane  view  =  
       generate  (constant  (Z:.600:.800))  
                         (\ix  -­‐>  let  Z  :.  y  :.  x  =  unlift  ix  
                                         in  lift  (  xmin  +  (fromIntegral  x  *  viewx)  /  800  
                                                         ,  ymin  +  (fromIntegral  y  *  viewy)  /  600  ))  
       where  
           (xmin,ymin,xmax,ymax)  =  unlift  view  
           viewx                                  =  xmax  -­‐  xmin  
           viewy                                  =  ymax  -­‐  ymin
    unlift  ::  Exp  (F,F,F,F)  
                 -­‐>  (Exp  F,  Exp  F,  Exp  F,  Exp  F)
    

    View Slide

73. - unlift has a funky type, but just unpacks “stuff”. lift does the reverse.
    Mandelbrot set generator
    genPlane  ::  Exp  (Float,Float,Float,Float)  -­‐>  Acc  ComplexPlane  
    genPlane  view  =  
       generate  (constant  (Z:.600:.800))  
                         (\ix  -­‐>  let  Z  :.  y  :.  x  =  unlift  ix  
                                         in  lift  (  xmin  +  (fromIntegral  x  *  viewx)  /  800  
                                                         ,  ymin  +  (fromIntegral  y  *  viewy)  /  600  ))  
       where  
           (xmin,ymin,xmax,ymax)  =  unlift  view  
           viewx                                  =  xmax  -­‐  xmin  
           viewy                                  =  ymax  -­‐  ymin
    

    View Slide

74. - unlift has a funky type, but just unpacks “stuff”. lift does the reverse.
    - Even though operations are in Exp, can still use standard Haskell
    operations like (*) and (-­‐)
    Mandelbrot set generator
    genPlane  ::  Exp  (Float,Float,Float,Float)  -­‐>  Acc  ComplexPlane  
    genPlane  view  =  
       generate  (constant  (Z:.600:.800))  
                         (\ix  -­‐>  let  Z  :.  y  :.  x  =  unlift  ix  
                                         in  lift  (  xmin  +  (fromIntegral  x  *  viewx)  /  800  
                                                         ,  ymin  +  (fromIntegral  y  *  viewy)  /  600  ))  
       where  
           (xmin,ymin,xmax,ymax)  =  unlift  view  
           viewx                                  =  xmax  -­‐  xmin  
           viewy                                  =  ymax  -­‐  ymin
    

    View Slide

75. Mandelbrot set generator
    • Let’s define the function we will iterate
    zn+1 = c + z2
    n
    next  ::  Exp  Complex  -­‐>  Exp  Complex  -­‐>  Exp  Complex  
    next  c  z  =  plus  c  (times  z  z)
    

    View Slide

76. Mandelbrot set generator
    • Let’s define the function we will iterate
    - Use lift and unlift as before to unpack the tuples
    zn+1 = c + z2
    n
    next  ::  Exp  Complex  -­‐>  Exp  Complex  -­‐>  Exp  Complex  
    next  c  z  =  plus  c  (times  z  z)
    plus  ::  Exp  Complex  -­‐>  Exp  Complex  -­‐>  Exp  Complex  
    plus  a  b  =  lift  (ax+bx,  ay+by)  
       where  
           (ax,ay)  =  unlift  a    ::  (Exp  Float,  Exp  Float)  
           (bx,by)  =  unlift  b    ::  (Exp  Float,  Exp  Float)
    

    View Slide

77. Mandelbrot set generator
    • Let’s define the function we will iterate
    - Use lift and unlift as before to unpack the tuples
    zn+1 = c + z2
    n
    next  ::  Exp  Complex  -­‐>  Exp  Complex  -­‐>  Exp  Complex  
    next  c  z  =  plus  c  (times  z  z)
    plus  ::  Exp  Complex  -­‐>  Exp  Complex  -­‐>  Exp  Complex  
    plus  a  b  =  lift  (ax+bx,  ay+by)  
       where  
           (ax,ay)  =  unlift  a    ::  (Exp  Float,  Exp  Float)  
           (bx,by)  =  unlift  b    ::  (Exp  Float,  Exp  Float)
    lift  ::  (Exp  Float,  Exp  Float)  
             -­‐>  Exp  (Float,  Float)
    

    View Slide

78. Mandelbrot set generator
    • Let’s define the function we will iterate
    - Use lift and unlift as before to unpack the tuples
    - Note that we had to add some type signatures to unlift
    zn+1 = c + z2
    n
    next  ::  Exp  Complex  -­‐>  Exp  Complex  -­‐>  Exp  Complex  
    next  c  z  =  plus  c  (times  z  z)
    plus  ::  Exp  Complex  -­‐>  Exp  Complex  -­‐>  Exp  Complex  
    plus  a  b  =  lift  (ax+bx,  ay+by)  
       where  
           (ax,ay)  =  unlift  a    ::  (Exp  Float,  Exp  Float)  
           (bx,by)  =  unlift  b    ::  (Exp  Float,  Exp  Float)
    lift  ::  (Exp  Float,  Exp  Float)  
             -­‐>  Exp  (Float,  Float)
    

    View Slide

79. Mandelbrot set generator
    • Complication: GPUs must do the same thing to lots of different data
    - Keep a pair (z,i) for each pixel, where i is the point iteration diverged
    zn+1 = c + z2
    n
    

    View Slide

80. Mandelbrot set generator
    • Complication: GPUs must do the same thing to lots of different data
    - Keep a pair (z,i) for each pixel, where i is the point iteration diverged
    zn+1 = c + z2
    n
    step  ::  Exp  Complex  -­‐>  Exp  (Complex,  Int)  -­‐>  Exp  (Complex,  Int)  
    step  c  zi  =  f  (fst  zi)  (snd  zi)  
       where  
           f  z  i  =  let  z'  =  next  c  z  
                           in    dot  z'  >*  4  ?  (  zi,  lift  (z',  i+1)  )
    

    View Slide

81. Mandelbrot set generator
    • Complication: GPUs must do the same thing to lots of different data
    - Keep a pair (z,i) for each pixel, where i is the point iteration diverged
    - fst and snd to extract individual tuple components
    zn+1 = c + z2
    n
    step  ::  Exp  Complex  -­‐>  Exp  (Complex,  Int)  -­‐>  Exp  (Complex,  Int)  
    step  c  zi  =  f  (fst  zi)  (snd  zi)  
       where  
           f  z  i  =  let  z'  =  next  c  z  
                           in    dot  z'  >*  4  ?  (  zi,  lift  (z',  i+1)  )
    

    View Slide

82. Mandelbrot set generator
    • Complication: GPUs must do the same thing to lots of different data
    - Keep a pair (z,i) for each pixel, where i is the point iteration diverged
    - fst and snd to extract individual tuple components
    - dot is the magnitude of a complex number squared
    zn+1 = c + z2
    n
    step  ::  Exp  Complex  -­‐>  Exp  (Complex,  Int)  -­‐>  Exp  (Complex,  Int)  
    step  c  zi  =  f  (fst  zi)  (snd  zi)  
       where  
           f  z  i  =  let  z'  =  next  c  z  
                           in    dot  z'  >*  4  ?  (  zi,  lift  (z',  i+1)  )
    

    View Slide

83. Mandelbrot set generator
    • Complication: GPUs must do the same thing to lots of different data
    - Keep a pair (z,i) for each pixel, where i is the point iteration diverged
    - fst and snd to extract individual tuple components
    - dot is the magnitude of a complex number squared
    - Conditionals can lead to SIMD divergence, so use sparingly
    zn+1 = c + z2
    n
    (?)  ::  Elt  t  =>  Exp  Bool  -­‐>  (Exp  t,  Exp  t)  -­‐>  Exp  t
    step  ::  Exp  Complex  -­‐>  Exp  (Complex,  Int)  -­‐>  Exp  (Complex,  Int)  
    step  c  zi  =  f  (fst  zi)  (snd  zi)  
       where  
           f  z  i  =  let  z'  =  next  c  z  
                           in    dot  z'  >*  4  ?  (  zi,  lift  (z',  i+1)  )
    

    View Slide

84. Mandelbrot set generator
    • Complication: GPUs must do the same thing to lots of different data
    - Keep a pair (z,i) for each pixel, where i is the point iteration diverged
    - fst and snd to extract individual tuple components
    - dot is the magnitude of a complex number squared
    - Conditionals can lead to SIMD divergence, so use sparingly
    zn+1 = c + z2
    n
    (?)  ::  Elt  t  =>  Exp  Bool  -­‐>  (Exp  t,  Exp  t)  -­‐>  Exp  t
    step  ::  Exp  Complex  -­‐>  Exp  (Complex,  Int)  -­‐>  Exp  (Complex,  Int)  
    step  c  zi  =  f  (fst  zi)  (snd  zi)  
       where  
           f  z  i  =  let  z'  =  next  c  z  
                           in    dot  z'  >*  4  ?  (  zi,  lift  (z',  i+1)  )
    

    View Slide

85. Mandelbrot set generator
    • Complication: GPUs must do the same thing to lots of different data
    - Keep a pair (z,i) for each pixel, where i is the point iteration diverged
    - fst and snd to extract individual tuple components
    - dot is the magnitude of a complex number squared
    - Conditionals can lead to SIMD divergence, so use sparingly
    zn+1 = c + z2
    n
    (?)  ::  Elt  t  =>  Exp  Bool  -­‐>  (Exp  t,  Exp  t)  -­‐>  Exp  t
    step  ::  Exp  Complex  -­‐>  Exp  (Complex,  Int)  -­‐>  Exp  (Complex,  Int)  
    step  c  zi  =  f  (fst  zi)  (snd  zi)  
       where  
           f  z  i  =  let  z'  =  next  c  z  
                           in    dot  z'  >*  4  ?  (  zi,  lift  (z',  i+1)  )
    

    View Slide

86. Mandelbrot set generator
    • Complication: GPUs must do the same thing to lots of different data
    - Keep a pair (z,i) for each pixel, where i is the point iteration diverged
    - fst and snd to extract individual tuple components
    - dot is the magnitude of a complex number squared
    - Conditionals can lead to SIMD divergence, so use sparingly
    zn+1 = c + z2
    n
    (?)  ::  Elt  t  =>  Exp  Bool  -­‐>  (Exp  t,  Exp  t)  -­‐>  Exp  t
    step  ::  Exp  Complex  -­‐>  Exp  (Complex,  Int)  -­‐>  Exp  (Complex,  Int)  
    step  c  zi  =  f  (fst  zi)  (snd  zi)  
       where  
           f  z  i  =  let  z'  =  next  c  z  
                           in    dot  z'  >*  4  ?  (  zi,  lift  (z',  i+1)  )
    

    View Slide

87. Mandelbrot set generator
    • Accelerate is a meta programming language
    - So, just use regular Haskell to unroll the loop a fixed number of times
    zn+1 = c + z2
    n
    mandel  ::  Exp  (Float,Float,Float,Float)  -­‐>  Acc  (Array  DIM2  (Complex,  Int))  
    mandel  view  =  
       Prelude.foldl1  (.)  (Prelude.replicate  255  go)  zs0  
       where  
           cs          =  genPlane  view  
           zs0        =  zip  cs  (fill  (shape  cs)  0)  
           go  zs    =  zipWith  step  cs  zs
    

    View Slide

88. Mandelbrot set generator
    • Accelerate is a meta programming language
    - So, just use regular Haskell to unroll the loop a fixed number of times
    - zipWith applies its function step pairwise to elements of the two arrays
    zn+1 = c + z2
    n
    mandel  ::  Exp  (Float,Float,Float,Float)  -­‐>  Acc  (Array  DIM2  (Complex,  Int))  
    mandel  view  =  
       Prelude.foldl1  (.)  (Prelude.replicate  255  go)  zs0  
       where  
           cs          =  genPlane  view  
           zs0        =  zip  cs  (fill  (shape  cs)  0)  
           go  zs    =  zipWith  step  cs  zs
    

    View Slide

89. Mandelbrot set generator
    • Accelerate is a meta programming language
    - So, just use regular Haskell to unroll the loop a fixed number of times
    - zipWith applies its function step pairwise to elements of the two arrays
    - Replicate the transition step from to
    zn+1 = c + z2
    n
    mandel  ::  Exp  (Float,Float,Float,Float)  -­‐>  Acc  (Array  DIM2  (Complex,  Int))  
    mandel  view  =  
       Prelude.foldl1  (.)  (Prelude.replicate  255  go)  zs0  
       where  
           cs          =  genPlane  view  
           zs0        =  zip  cs  (fill  (shape  cs)  0)  
           go  zs    =  zipWith  step  cs  zs
    zn
    zn+1
    

    View Slide

90. Mandelbrot set generator
    • Accelerate is a meta programming language
    - So, just use regular Haskell to unroll the loop a fixed number of times
    - zipWith applies its function step pairwise to elements of the two arrays
    - Replicate the transition step from to
    - Applies the steps in sequence, beginning with the initial data zs0
    zn+1 = c + z2
    n
    mandel  ::  Exp  (Float,Float,Float,Float)  -­‐>  Acc  (Array  DIM2  (Complex,  Int))  
    mandel  view  =  
       Prelude.foldl1  (.)  (Prelude.replicate  255  go)  zs0  
       where  
           cs          =  genPlane  view  
           zs0        =  zip  cs  (fill  (shape  cs)  0)  
           go  zs    =  zipWith  step  cs  zs
    zn
    zn+1
    

    View Slide

91. Mandelbrot set generator
    • Accelerate is a meta programming language
    - So, just use regular Haskell to unroll the loop a fixed number of times
    - zipWith applies its function step pairwise to elements of the two arrays
    - Replicate the transition step from to
    - Applies the steps in sequence, beginning with the initial data zs0
    zn+1 = c + z2
    n
    mandel  ::  Exp  (Float,Float,Float,Float)  -­‐>  Acc  (Array  DIM2  (Complex,  Int))  
    mandel  view  =  
       Prelude.foldl1  (.)  (Prelude.replicate  255  go)  zs0  
       where  
           cs          =  genPlane  view  
           zs0        =  zip  cs  (fill  (shape  cs)  0)  
           go  zs    =  zipWith  step  cs  zs
    zn
    zn+1
    f (g x) ≣ (f . g) x
    

    View Slide

92. In the workshop…
    • More example code, computational kernels
    
    - Common pitfalls
    
    - Tips for good performance
    
    - Figuring out what went wrong (or knowing who to blame)
    

    View Slide

93. In the workshop…
    • More example code, computational kernels
    
    - Common pitfalls
    
    - Tips for good performance
    
    - Figuring out what went wrong (or knowing who to blame)
    me
    

    View Slide

94. Questions?
    https://github.com/AccelerateHS/
    http://xkcd.com/365/
    

    View Slide

95. Extra Slides…
    

    View Slide

96. Seriously?
    

    View Slide

97. Arrays
    • Create an array from a list:
    
    !
    - Generates a multidimensional array by consuming elements from the list
    and adding them to the array in row-major order
    
    • Example:
    data  Array  sh  e
    fromList  ::  (Shape  sh,  Elt  e)  =>  sh  -­‐>  [e]  -­‐>  Array  sh  e
    ghci>  fromList  (Z:.10)  [1..10]
    

    View Slide

98. Arrays
    • Create an array from a list:
    data  Array  sh  e
    >  fromList  (Z:.10)  [1..10]  ::  Vector  Float  
    Array  (Z  :.  10)  [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0]
    

    View Slide

99. Arrays
    • Create an array from a list:
    • Multidimensional arrays are similar:
    
    - Elements are filled along the right-most dimension first
    data  Array  sh  e
    >  fromList  (Z:.10)  [1..10]  ::  Vector  Float  
    Array  (Z  :.  10)  [1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0]
    >  fromList  (Z:.3:.5)  [1..]  ::  Array  DIM2  Int  
    Array  (Z  :.  3  :.  5)  [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]
    1 2 3 4 5
    6 7 8 9 10
    11 12 13 14 15
    

    View Slide

100. Arrays
     • Array indices start counting from zero
     data  Array  sh  e
     >  let  mat  =  fromList  (Z:.3:.5)  [1..]  ::  Array  DIM2  Int  
     >  indexArray  mat  (Z:.2:.1)  
     12
     1 2 3 4 5
     6 7 8 9 10
     11 12 13 14 15
     

     View Slide

101. Arrays
     • Similarly, an array of (nested) tuples:
     
     - This is just a trick: internally converted into a tuple of arrays
     data  Array  sh  e
     >  fromList  (Z:.2:.3)  $  P.zip  [1..]  ['a'..]  ::  Array  DIM2  (Int,Char)  
     Array  (Z  :.  2  :.  3)  [(1,'a'),(2,'b'),(3,'c'),(4,'d'),(5,'e'),(6,'f')]
     1 2 3
     4 5 6
     a b c
     d e f
     ( )
     ,
     

     View Slide

102. Data.Array.Accelerate
     • Need to import both the base library as well as a backend
     
     - There is also an interpreter available for testing
     
     - Runs without using the GPU (much more slowly of course)
     import  Prelude                                        as  P  
     import  Data.Array.Accelerate            as  A  
     import  Data.Array.Accelerate.CUDA  as  CUDA
     

     View Slide

103. Data.Array.Accelerate
     

     View Slide

104. Data.Array.Accelerate
     • To get arrays into the Acc world:
     
     !
     - This may involve copying data to the GPU
     use  ::  Arrays  a  =>  a  -­‐>  Acc  a
     

     View Slide

105. Data.Array.Accelerate
     • To get arrays into the Acc world:
     
     !
     - This may involve copying data to the GPU
     • use injects arrays into our DSL
     • run executes the computation to get arrays out
     • Using Accelerate focuses on everything in between
     use  ::  Arrays  a  =>  a  -­‐>  Acc  a
     

     View Slide

106. Collective Operations
     • Example: add one to each element of an array
     >  let  arr  =  fromList  (Z:.3:.5)  [1..]  ::  Array  DIM2  Int  
     >  run  $  A.map  (+1)  (use  arr)  
     Array  (Z  :.  3  :.  5)  [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
     

     View Slide

107. Collective Operations
     • Example: add one to each element of an array
     • What is the type of map?
     >  let  arr  =  fromList  (Z:.3:.5)  [1..]  ::  Array  DIM2  Int  
     >  run  $  A.map  (+1)  (use  arr)  
     Array  (Z  :.  3  :.  5)  [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
     map  ::  (Shape  sh,  Elt  a,  Elt  b)  
            =>  (Exp  a  -­‐>  Exp  b)  
            -­‐>  Acc  (Array  sh  a)  
            -­‐>  Acc  (Array  sh  b)
     

     View Slide

108. Collective Operations
     • Example: add one to each element of an array
     • What is the type of map?
     >  let  arr  =  fromList  (Z:.3:.5)  [1..]  ::  Array  DIM2  Int  
     >  run  $  A.map  (+1)  (use  arr)  
     Array  (Z  :.  3  :.  5)  [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
     map  ::  (Shape  sh,  Elt  a,  Elt  b)  
            =>  (Exp  a  -­‐>  Exp  b)  
            -­‐>  Acc  (Array  sh  a)  
            -­‐>  Acc  (Array  sh  b)
     Supported shape &
     element types
     

     View Slide

109. Collective Operations
     • Example: add one to each element of an array
     • What is the type of map?
     >  let  arr  =  fromList  (Z:.3:.5)  [1..]  ::  Array  DIM2  Int  
     >  run  $  A.map  (+1)  (use  arr)  
     Array  (Z  :.  3  :.  5)  [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
     map  ::  (Shape  sh,  Elt  a,  Elt  b)  
            =>  (Exp  a  -­‐>  Exp  b)  
            -­‐>  Acc  (Array  sh  a)  
            -­‐>  Acc  (Array  sh  b)
     DSL array
     Supported shape &
     element types
     

     View Slide

110. Collective Operations
     • Example: add one to each element of an array
     • What is the type of map?
     >  let  arr  =  fromList  (Z:.3:.5)  [1..]  ::  Array  DIM2  Int  
     >  run  $  A.map  (+1)  (use  arr)  
     Array  (Z  :.  3  :.  5)  [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
     map  ::  (Shape  sh,  Elt  a,  Elt  b)  
            =>  (Exp  a  -­‐>  Exp  b)  
            -­‐>  Acc  (Array  sh  a)  
            -­‐>  Acc  (Array  sh  b)
     DSL array
     Function to apply at every
     element. But what is Exp?
     Supported shape &
     element types
     

     View Slide

111. • The type class overloading trick is used for standard Haskell classes
     
     !
     • Standard boolean operations are available with slightly different names
     
     - The standard names can not be overloaded
     
     !
     • Conditionals
     
     - Use sparingly: leads to SIMD divergence
     (+1)  ::  (Elt  a,  IsNum  a)  =>  Exp  a  -­‐>  Exp  a
     Scalar Expressions
     (==*)  ::  (Elt  t,  IsScalar  t)  =>  Exp  t  -­‐>  Exp  t  -­‐>  Exp  Bool  
     (/=*),  (<*),  (>*),  min,  max,  (||*),  (&&*)      -­‐-­‐  and  so  on...
     (?)  ::  Elt  t  =>  Exp  Bool  -­‐>  (Exp  t,  Exp  t)  -­‐>  Exp  t
     

     View Slide

112. Scalar Expressions
     • Bring a Haskell value into Exp land
     
     !
     • Lift an expression into a singleton array
     
     !
     • Extract the element from a singleton array
     constant  ::  Elt  e  =>  e  -­‐>  Exp  e
     unit  ::  Exp  e  -­‐>  Acc  (Scalar  e)
     the  ::  Acc  (Scalar  e)  -­‐>  Exp  e
     

     View Slide

113. Reductions
     • Folding (+) over a vector produces a sum
     
     - The result is a one-element array (scalar). Why?
     
     !
     >  let  xs  =  fromList  (Z:.10)  [1..]  ::  Vector  Int  
     >  run  $  A.fold  (+)  0  (use  xs)  
     Array  (Z)  [55]
     

     View Slide

114. Reductions
     • Folding (+) over a vector produces a sum
     
     - The result is a one-element array (scalar). Why?
     
     !
     • Fold has an interesting type:
     >  let  xs  =  fromList  (Z:.10)  [1..]  ::  Vector  Int  
     >  run  $  A.fold  (+)  0  (use  xs)  
     Array  (Z)  [55]
     fold  ::  (Shape  sh,  Elt  a)  
              =>  (Exp  a  -­‐>  Exp  a  -­‐>  Exp  a)  
              -­‐>  Exp  a  
              -­‐>  Acc  (Array  (sh:.Int)  a)  
              -­‐>  Acc  (Array  sh                a)
     

     View Slide

115. Reductions
     • Folding (+) over a vector produces a sum
     
     - The result is a one-element array (scalar). Why?
     
     !
     • Fold has an interesting type:
     >  let  xs  =  fromList  (Z:.10)  [1..]  ::  Vector  Int  
     >  run  $  A.fold  (+)  0  (use  xs)  
     Array  (Z)  [55]
     fold  ::  (Shape  sh,  Elt  a)  
              =>  (Exp  a  -­‐>  Exp  a  -­‐>  Exp  a)  
              -­‐>  Exp  a  
              -­‐>  Acc  (Array  (sh:.Int)  a)  
              -­‐>  Acc  (Array  sh                a)
     input array
     

     View Slide

116. Reductions
     • Folding (+) over a vector produces a sum
     
     - The result is a one-element array (scalar). Why?
     
     !
     • Fold has an interesting type:
     >  let  xs  =  fromList  (Z:.10)  [1..]  ::  Vector  Int  
     >  run  $  A.fold  (+)  0  (use  xs)  
     Array  (Z)  [55]
     fold  ::  (Shape  sh,  Elt  a)  
              =>  (Exp  a  -­‐>  Exp  a  -­‐>  Exp  a)  
              -­‐>  Exp  a  
              -­‐>  Acc  (Array  (sh:.Int)  a)  
              -­‐>  Acc  (Array  sh                a)
     outer dimension removed
     input array
     

     View Slide

117. Reductions
     • Fold occurs over the outer dimension of the array
     >  let  mat  =  fromList  (Z:.3:.5)  [1..]  ::  Array  DIM2  Int  
     >  run  $  A.fold  (+)  0  (use  mat)  
     Array  (Z  :.  3)  [15,40,65]
     1 2 3 4 5
     6 7 8 9 10
     11 12 13 14 15
     15
     40
     65
     

     View Slide

118. Reductions
     • Is this a left-fold or a right-fold?
     
     - Neither! The fold happens in parallel, tree-like
     
     - Therefore the function must be associative: (Exp  a  -­‐>  Exp  a  -­‐>  Exp  a)
     
     - (We pretend that floating point operations are associative, though strictly
     speaking they are not)
     

     View Slide

119. Stencils
     • A stencil is a map with access to the neighbourhood around each element
     - Useful in many scientific & image processing algorithms
     laplace  ::  Stencil3x3  Int  -­‐>  Exp  Int  
     laplace  ((_,t,_)  
                    ,(l,c,r)  
                    ,(_,b,_))  =  t  +  b  +  l  +  r  -­‐  4*c
     t
     l c r
     b
     

     View Slide

120. Stencils
     • A stencil is a map with access to the neighbourhood around each element
     - Useful in many scientific & image processing algorithms
     - Boundary conditions specify how to handle out-of-bounds neighbours
     laplace  ::  Stencil3x3  Int  -­‐>  Exp  Int  
     laplace  ((_,t,_)  
                    ,(l,c,r)  
                    ,(_,b,_))  =  t  +  b  +  l  +  r  -­‐  4*c
     >  let  mat  =  fromList  (Z:.3:.5)  [1..]  ::  Array  DIM2  Int  
     >  run  $  stencil  laplace  (Constant  0)  (use  mat)  
     Array  (Z  :.  3  :.  5)  [4,3,2,1,-­‐6,-­‐5,0,0,0,-­‐11,-­‐26,-­‐17,-­‐18,-­‐19,-­‐36]
     t
     l c r
     b
     

     View Slide

121. Stencils
     • A stencil is a map with access to the neighbourhood around each element
     
     - Useful in many scientific & image processing algorithms
     

     View Slide

122. Index Transforms
     • Index transforms change the order of elements, not their values
     
     - We usually want to push such operations into the consumer
     
     • backpermute specifies which element to read from a source array
     backpermute  ::  (Shape  ix,  Shape  ix',  Elt  a)  
                            =>  Exp  ix'  
                            -­‐>  (Exp  ix'  -­‐>  Exp  ix)  
                            -­‐>  Acc  (Array  ix    a)  
                            -­‐>  Acc  (Array  ix'  a)
     

     View Slide

123. Index Transforms
     • Index transforms change the order of elements, not their values
     
     - We usually want to push such operations into the consumer
     
     • backpermute specifies which element to read from a source array
     backpermute  ::  (Shape  ix,  Shape  ix',  Elt  a)  
                            =>  Exp  ix'  
                            -­‐>  (Exp  ix'  -­‐>  Exp  ix)  
                            -­‐>  Acc  (Array  ix    a)  
                            -­‐>  Acc  (Array  ix'  a)
     shape of result
     

     View Slide

124. Index Transforms
     • Index transforms change the order of elements, not their values
     
     - We usually want to push such operations into the consumer
     
     • backpermute specifies which element to read from a source array
     backpermute  ::  (Shape  ix,  Shape  ix',  Elt  a)  
                            =>  Exp  ix'  
                            -­‐>  (Exp  ix'  -­‐>  Exp  ix)  
                            -­‐>  Acc  (Array  ix    a)  
                            -­‐>  Acc  (Array  ix'  a)
     shape of result
     source data
     

     View Slide

125. Index Transforms
     • Index transforms change the order of elements, not their values
     
     - We usually want to push such operations into the consumer
     
     • backpermute specifies which element to read from a source array
     backpermute  ::  (Shape  ix,  Shape  ix',  Elt  a)  
                            =>  Exp  ix'  
                            -­‐>  (Exp  ix'  -­‐>  Exp  ix)  
                            -­‐>  Acc  (Array  ix    a)  
                            -­‐>  Acc  (Array  ix'  a)
     shape of result
     source data
     index mapping from
     destination array to
     source
     

     View Slide

126. Index Transforms
     • Index transforms change the order of elements, not their values
     transpose  mat  =  
        let  swap  =  lift1  $  \(Z:.j:.i)  -­‐>  Z:.i:.j  ::  Z  :.  Exp  Int  :.  Exp  Int  
        in    backpermute  (swap  $  shape  mat)  swap  mat
     1 2 3 4 5
     6 7 8 9 10
     11 12 13 14 15
     1 6 11
     2 7 12
     3 8 13
     4 9 14
     5 10 15
     

     View Slide