Efficient Successive Over Relaxation using Cuda?

Question

I recently implemented Successive Over Relaxation using Cuda as a part of my course project and was curious to know how I can make the code more efficient.

I'm using Red/Black SOR scheme which is popular version of parallel SOR. To update the value at a grid point, I need to access 4 adjacent values which are generally not contiguous in the memory and hence the GPU is not able to use memory coalescing feature of cuda which makes the algorithm memory bound. To overcome this I have reordered the storage of odd and even grid points in a separate vectors as mentioned in the paper here: http://link.springer.com/chapter/10.1007%2F978-3-642-31464-3_60 which takes advantage of memory coalescing.

Implementing the technique mentioned in the paper I see a speedup of around ~15 (compared to ~10 with traditional Red/Black SOR) when compared with 1 CPU. How much performance speed up I should be expecting? I feel speedup of ~15 is quite low?

I'm calling the kernel functions to update odd values and then even values sequentially from host CPU in following way :

for(size_t it = 0; it < itmax; it++){
      cuda_odd_update <<<dimGrid_odd,dimBlock>>> (odd,even);
      cuda_even_update <<<dimGrid_even,dimBlock>>> (odd,even);
}

where odd and even are the pointers to reordered vectors in device memory. I'm using this approach as I'm not aware of any way to synchronize all the threads in the device. Is there any better way to implement the execution loop?

The code below implements the technique mentioned in the paper and not a traditional Red/Black SOR:

__global__ void cuda_odd_update(double* odd,double* even){

  size_t tx = blockIdx.x*blockDim.x + threadIdx.x;
  size_t ty = blockIdx.y*blockDim.y + threadIdx.y;
  size_t odd_index = tx*height_odd+ty;
  size_t even_index = tx*height_even+ty;

  if (( (ty == 0 && tx%2 == 0) || (ty == height_odd-1 && tx%2 == 1) || (ty > 0 && ty < height_odd-1) ) && (tx > 0 && tx < width-1 && ty < height_odd) ){

      odd[odd_index] = (1.0 - omega)*odd[odd_index] + omega/(2*(1+beta))
                     * ( (tx%2)*even[even_index]+(1-tx%2)*even[even_index+1]
                       + (tx%2)*even[even_index-1]+(1-tx%2)*even[even_index]
                       + beta * ( even[even_index-height_even]
                                + even[even_index+height_even]   )) ;

  }

}


__global__ void cuda_even_update(double* odd,double* even){

  size_t tx = blockIdx.x*blockDim.x + threadIdx.x;
  size_t ty = blockIdx.y*blockDim.y + threadIdx.y;
  size_t odd_index = tx*height_odd+ty;
  size_t even_index = tx*height_even+ty;

  if(( (ty == 0 & tx%2 == 1) || (ty == height_even-1 && tx%2 == 0) || (ty > 0 && ty < height_even-1) ) && ( tx > 0 && tx < width-1 && ty < height_even)){


      even[even_index] = (1.0 - omega)*even[even_index] + omega/(2*(1+beta))
                       * ( (1-tx%2)*odd[odd_index] + (tx%2)*odd[odd_index+1]
                         + (1-tx%2)*odd[odd_index-1] + (tx%2)*odd[odd_index]
                         + beta * ( odd[odd_index-height_odd]
                                  + odd[odd_index+height_odd]  )) ;


  }

}

I'm currently using Intel Xeon CPU E5-2660-2.20GHz (8 cores) and NVIDIA TESLA M2090. The nvprof gives the following results (for 1000 iterations):

CPU time : 62.69
==27665== NVPROF is profiling process 27665, command: ./sor
GPU time : 3.73
End!
==27665== Profiling application: ./sor
==27665== Profiling result:
Time(%)      Time     Calls       Avg       Min       Max  Name
 49.81%  1.83752s      1000  1.8375ms  1.8333ms  1.8422ms  cuda_even_update(double*, double*)
 49.67%  1.83226s      1000  1.8323ms  1.8290ms  1.8364ms  cuda_odd_update(double*, double*)
  0.26%  9.6281ms         2  4.8141ms  4.6306ms  4.9975ms  [CUDA memcpy DtoH]
  0.26%  9.4569ms         2  4.7284ms  4.7065ms  4.7504ms  [CUDA memcpy HtoD]

Answer 1

Synchronization. You can't synchronize threads across different blocks in the same grid. On the other hand, if you execute kernels the way you do, they are guaranteed to execute in sequence, so all the threads in even_update will see all updates made by odd_update in the same iteration. (See NVIDIA's own documentation.)

There is no good way that I know of to synchronize threads across a whole kernel, e.g. see this discussion), so having a main loop execute kernels in sequence seems fine. You can measure

Speedup. Speedup relative to the CPU is not such a useful performance measure. I know it is always reported, and so on, but it can be a little bit misleading. For example, someone using the same GPU and the same memory, but a faster CPU will get a lower speedup, but that says nothing about the performance of their code. NVIDIA's best practices guide has some things to say about how to measure performance of your code. Since your algorithm is memory-bandwidth-limited (it performs a relatively small number of operations per memory access), the quantity I'd say is more important than the speedup is memory throughput, the number of bytes per second your program actually reads/writes (effective bandwidth), especially as a proportion of theoretical peak bandwidth of your GPU. Also helpful is to know how much more memory your kernel reads/writes due to coalescing requirements than it strictly requires.

Have you tried measuring these things? Fraction of theoretical peak bandwidth is a more useful guidance for optimization than speedup relative to CPU.

Performance. As @GodricSeer points out, executing lots of small kernels might be a problem.

One immediate thing to try is to assign more work to each thread. So instead of each thread doing work for one index $j$, make each thread do the work for $m$ indices $j_{0:m-1}$, suitably spread out so that at every step contiguous threads to the work of contiguous indices. This reduces the amount of work the scheduler has to do by a factor of $m$.

Another issue might be with things like $\mathit{even}[\mathit{even\_index}\pm1]$, where the current threads might be reading memory that was already read by a thread in the same block. I haven't run your code, but one thing to try is to load parts of $\mathit{even}$ into shared memory so that threads with $\mathrm{threadIdx}.y$ differing by $1$ can use each other's global memory reads. There might be an issue with how much shared memory this uses is the size of a block is too large, so this needs careful explicit testing.

When you run NVIDIA's profiler on your code, does it highlight any interesting issues?