Cuda kernels are written inefficiently

[nikitinvv]

Cuda kernels in usfft.cu, convolution.cu are not optimal.

They are heavy and use a huge number of registers which may significantly slowdown the code. They should be split into smaller ones.

1. If there exist 'if,else' statement like for fwd/adj operator then it is better to split it into 2, like it was at the beginning.

2. Loops with many iterations inside the kernels should be avoided by adding new threads, this way allowing the scheduler to switch threads more efficiently, like it was at the beginning

3. Non-sequential and uncoallesced memory access in the loops make the code slower, the following loop structure is unacceptable from cuda optimization point of view. It also looks unreadable.

   // for each image
   for (int ti = blockIdx.z; ti < nimage; ti += gridDim.z) {
   // for each scan position
   for (int ts = blockIdx.y; ts < nscan; ts += gridDim.y) {
     ....
     // for x,y coords in patch
     for (int py = blockIdx.x; py < patch_shape; py += gridDim.x) {
       for (int px = threadIdx.x; px < patch_shape; px += blockDim.x) {

   Use 3d thread blocks and grids associated with each array dimension, like it was at the beginning. This will give you natural coalesced memory access.

4. Block size should be a power of 2, optimal combinations (1024,1,1), (256,4,4), (32,32,1), (16,16,4).. What we have in the code:

   block = (min(self.scatter_kernel.max_threads_per_block, (2 * m)**3),)

   Say m=3, then block = 216 - not optimal.

Would It be better to return to my initial kernels implementations?

[carterbox]

Related to https://github.com/tomography/tike/pull/69#pullrequestreview-421461953

[carterbox]

for (int py = blockIdx.x; py < patch_shape; py += gridDim.x) {
        for (int px = threadIdx.x; px < patch_shape; px += blockDim.x) {

3.

• This memory access pattern already is both sequential and coalesced for the block and grid sizes that I have chosen. Each block processes one row of one patch (a row is sequential in memory). Because the maximum number of threads-per-block is 1024, you cannot process a patch larger than 32 x 32 with a single block. Instead I assign each block a single row and then the maximum patch size becomes 1024x1024 for coalesced memory access (in both the patch and the full image). The following code chunk is equivalent to the grid-stride-loops above when gridDim.x >= patch_shape and blockDim.x >= patch_shape. The difference is that grid-stride-loops will execute correctly for any block and grid shape; you can still tune the kernel launches by adjusting the block and grid shapes.

  int py = blockIdx.x;
  if (py >= patch_shape) return;
  int px = threadIdx.x;
  if (px >= patch_shape) return;

• The outer loops (below) control the block indexes, but since the maximum length of the z,y dimensions of a grid are 65535 (256 x 256), we should move the scan dimension to the x dimension of the grid (max size 2**31). Reiterating that the grid-stride-loops don't behave differently than an if () return; statement if the problem size is smaller than the grid.

  // for each image
  for (int ti = blockIdx.z; ti < nimage; ti += gridDim.z) {
  // for each scan position
  for (int ts = blockIdx.y; ts < nscan; ts += gridDim.y) {

[carterbox]

4. We can make the block-size always a power of 2. We just need to write a function that finds the next power of two after or equal to a number.