# Efficient Successive Over Relaxation using Cuda?

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]

• For iterative solvers, about an order of magnitude speedup sounds reasonable using a GPU (see www-users.cs.umn.edu/~saad/PDF/umsi-2010-112.pdf for benchmarks). You may also be suffering from kernel startup time (depending on your problem size) since you are doing so twice per iteration. – Godric Seer Dec 10 '14 at 18:28
• is there any way to include iteration loop within the kernel functions? – Pranav Dec 10 '14 at 18:37
• It has been a long time since I have programmed in cuda, so I won't make a recommendation. I just know that for small problems kernel startup can cause errors in benchmarking since it tends to be long (many microseconds). If each kernel execution takes more than a millisecond or so it likely isn't a problem. – Godric Seer Dec 10 '14 at 18:39
• Can you give more details on what performance you are getting? Speedup of 15 relative to your CPU isn't much to go on, to be honest. – Kirill Dec 10 '14 at 20:26
• Is your CPU program using all the cores fully? It's not really a fair comparison otherwise. – Kirill Dec 10 '14 at 22:02

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

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?

• you are right, I should probably do bandwidth related studies to measure performance more accurately. Also I have pasted the nvprof results in the post. – Pranav Dec 10 '14 at 21:57
• @Pranav Can you make nvprof output memory/instruction throughput? I definitely remember the visual profiler can do that; it can also give you warnings on performance problems. – Kirill Dec 10 '14 at 22:00
• well, I don't have access to visual profiler as I am running them remotely..and the other output flags prints info for every iterations.. – Pranav Dec 10 '14 at 22:21
• @Pranav I think it is impossible for me to say anything more concrete without much more detailed information. – Kirill Dec 10 '14 at 22:32