Ben Hirsch

Here are some benchmarking notes on CUDA vs the Eigen Library on the two most common operations in my research. Those two operations are the SAXPY operation, which is Y = a * X + Y where X, Y are vectors and a is a scalar. The other main operation is a dot product, which is sum(X * Y) where X and Y are arrays using element-wise array multiplication.

SAXPY

Lots of data here to parse. The main takeaway here is that CUDA equals the performance of Eigen starting at vectors of length 100,000 and above. After that, CUDA roughly outperforms Eigen by an order of magnitude.

I tried a bunch of different parameters, but the main thing was that 128 threads per block with 4 loads in each thread was clearly the winner in terms of performance overall. 128 threads per block obviously creates more blocks; more blocks is essential to maximizing the GPUs performance as they fully occupy all the streaming multiprocessors (SMs). 4 loads per thread means that each thread works on 4 elements, with each element being strided by the width of the block. Doing work in this way allows warps to be swapped out as they request data to be loaded from global memory so that the addition/multiplication instructions can be performed as soon as the data is ready, pipelining everything in a way.

Code sample for 2 loads/thread:

__global__ instruction_level_parallelism(float *x, float *y, float scale, unsigned numThreads) {
   unsigned idx = threadIdx.x + blockIdx.x * blockDim.x;
   if (idx < numThreads)
      x[idx] += y[idx] * scale;

   if (2 * idx < numThreads)
      x[idx + blockDim.x] += y[idx + blockDim.x] * scale;
}

Obviously numThreads should be the number of elements divided by 2 when launching the kernel. But basically there are two loads strided across by the block size.

And finally the difference between 1 and 2 vectors. Normally, you might have vectors stored as:

X = X₁, X₂, X₃, …, X_N
Y = Y₁, Y₂, Y₃, …, Y_N

and SAXPYs those together. But instead, the one vector approach lays them out as:

XY = X₁, Y₁, X₂, Y₂, X₃, Y₃, …, X_N, Y_N

The CUDA cards are able to coalesce memory accesses when consectutive threads access memory locations contiguously. By laying the memory out like this, thread 1 can access memory location 1, thread 2 memory location 2, etc contiguously which maximizes load efficiency.

Dot Products

Here it’s even more interesting. You actually don’t get a similar magnitude of performance until the 100,000 length mark, but even then Eigen is slightly faster, about 1E-5 to 2E-5 faster, which isn’t a lot but (in GenSel) 50,000 markers and a Markov Chain length of 40,000, that’s 2,000,000 dot products which can add up (to like whole seconds). After 1,000,000 CUDA can be 2-5x faster, not an order of magnitude but still faster which is important.

Here it appears the 2 Vector dot product is actually faster. The “Templated Dot Product” is just Mark Harris’ Optimizing Parallel Reduction in CUDA with two vectors. The other dot product is a simple loop which performs a single sum per thread, reducing the vector in half per accumulation iteration. The Product Kernel then Templated Reduce tries to perform an element-wise multiplication followed by the templated reduce. It’s odd that that’s slower in my opinion, but whatever. Also interesting to note is that Eigen is lazy; I had to print the values of the dot products at the end and then remove that last field when analyzing the data just to get the timings right.

Conclusion

My main conclusion was that once you get to vectors at or above lengths of 100,000, CUDA is an excellent choice for matrix math. Below that though you can stick to libraries like Eigen which optimize code very very well. There’s lots of different optimizations you can do, but the main thing is optimizing memory bandwidth and distributing enough work across all of the multiprocessors.