title: Parallel dimension reduction (3D to 2D with sum) using Cuda
body: <p>In a CUDA application, I have an <code>N x N x D</code> matrix that I want to reduce to <code>N x D</code> by summing over the entire first (or second) axis. How do I do this most efficiently?</p>
  
  <p>Typically, N is greater than 10000 and D is 2 or 3.</p>
  
  <p>A quick and naive solution using atomicAdd would be the following:</p>
  
  <pre><code>namespace kernel {
      __global__ void sumNND(float* devPtrIn, float* devPtrOut, const int N, const int D) {
          int index = blockIdx.x * blockDim.x + threadIdx.x;
          int stride = blockDim.x * gridDim.x;
  
          for (int id = index; id &lt; N * N * D; id += stride) {
              const unsigned int d = id % D;
              const unsigned int i = (id - d) / D;
              const unsigned int n = i / N;
              const unsigned int m = i % N;
  
              atomicAdd(&amp;devPtrOut[d + D * n], devPtrIn[d + D * n + N * m]);
          }
      }
  }
  
  void sumNND(const int numBlocks, const int blockSize, float* devPtrIn, float* devPtrOut, const int N, const int D) {
      HANDLE_ERROR(cudaMemset(devPtrOut, 0, N * D * sizeof(float)));
      kernel::sumNND&lt;&lt;&lt;numBlocks, blockSize&gt;&gt;&gt;(devPtrIn, devPtrOut, N, D);
      HANDLE_ERROR(cudaDeviceSynchronize());
  }
  </code></pre>
  
  <p>where <code>sumNND</code> is being called with</p>
  
  <p><code>loopSize = N * N * D</code>,
  <code>blockSize = 768</code> and
  <code>numBlocks = (loopSize + blockSize - 1) / blockSize</code>.</p>
  
  <p>This is (not surprisingly) a bottleneck in my timeline, but I can't figure out how to effectively parallelize the dimension reduction. Any pointers?</p>