cudamat/cudamat — pull request #15 by f0k (2013-09-11 22:20:44+00:00)
Title: Choose block count dynamically

Nearly all kernels in `cudamat_kernels.cu` start with the following code:
```C
const unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;
const unsigned int numThreads = blockDim.x * gridDim.x;
for (unsigned int i = idx; i < len; i += numThreads) {
//do something
}
```
All these kernels are run from `cudamat.cu` with grid dimensions `<<>>`, which defaults to 4096 blocks of 512 threads, or 2,097,152 threads in total.

For operations on matrices with fewer than 2 million elements, more threads will be spawned than needed. Although they will never execute the body of the for loop, they still incur a severe performance penalty (see #2).

We should instead only spawn `ceil(len / 512.0) = (len + 511) / 512` blocks of 512 threads, and even reduce the number of threads for matrices with fewer than 512 elements.
(In theory this would allow us to spawn one thread per element, so we could replace the loop with an if statement, just to disable the extra threads of the last block. However, with a one-dimensional grid of thread blocks we can only have 65535 blocks in total and still need the loop to support matrices larger than 65535*512 elements. And actually [it is recommended](https://developer.nvidia.com/content/cuda-pro-tip-write-flexible-kernels-grid-stride-loops) to use a loop of this style for scalability and portability.)

Any comments? Shall I work on a pull request for the change outlined above or do you see a better solution? Setting the number of blocks to a multiple of the number of multiprocessors or CUDA cores might be another way...