In this exercise, we will work with a C implementation of the Jacobi method for solving Laplace's equation in 2D. The Jacobi method is an iterative scheme that repeatedly updates each interior grid point to the average of its four immediate neighbors, while the boundary values of the grid stay fixed. (Though this method is very simple to implement, it is not actually recommended for this problem, as it converges very slowly.) The code below is a correct implementation of the algorithm paired with a hasty attempt at OpenMP offload on the two large loops. We will diagnose and address some of the issues with this implementation. Can you see where to improve it?...

Note that in the first loop, the average of the four neighboring grid values of A[j][i] is computed and stored in Anew[j][i], and in the second loop, the values of A are updated from Anew. Offload works in each case because the default data mapping for arrays is tofrom, which means that A and Anew are transferred to/from the GPU before/after each loop. The scalar maxdiff, which is used to test convergence, is mapped as tofrom as well because it is involved in a reduction clause. (If maxdiff were not involved in a reduction clause, it would use the normal default mapping for scalars, firstprivate, and could not to return a result to the host.) Overall, the code as it appears above is correct, but its performance is extremely poor.

In the following three problems, we will incrementally improve this code to make full use of the OpenMP offload techniques we have encountered. In doing so, we will not change what the program computes; we will only add or modify OpenMP directives. If you want to try out your modifications, the commands to compile and run your codes on Vista would look like this:

Problem 1: Distribute the Workload

As noted in the previous topic, simply using #pragma omp target does not automatically distribute the workload across the GPU. Fix the code by adding the appropriate worksharing directives and constructs. You should attempt this problem before looking at the solution below.

Explanation: In Teams and Parallel and Distribute and For/Do, we saw that an OpenMP offload region created without the teams directive causes the entire region to be executed on a single GPU SM. We can do much better than that. The teams and parallel directives, together, create a league of teams of threads (equivalent to CUDA blocks of threads), and the distribute and for directives then distribute the loop iterations across the teams and the threads, respectively. In addition, notice that each offload region consists of a nested loop of depth two, so we can add the collapse(2) clause to merge the two loops into one larger loop that is more efficiently distributed across the GPU threads.

Problem 2: Add Explicit Mapping

In the two OpenMP directives, we are implicitly copying A and Anew to the GPU and back on every iteration. If you read the code carefully, you will notice this is not necessary. The next challenge is to reduce the amount of data being transferred by adding directional data movers to the OpenMP directives. You should attempt this problem before looking at the solution below.

Explanation: In the first loop, the values of A are only read and the values of Anew are only written, so we can use the to and from mapping types, respectively. In the second loop, the roles are reversed, and the directions of the transfers are as well. Although the boundary values in Anew are never used, we transfer the entire array every time for the sake of simplicity. We could alternatively restrict the mapping of Anew to just its interior values, by specifying map(from: Anew[1:N-2][1:N-2]) and map(to: Anew[1:N-2][1:N-2]). (Remember, when mapping in C/C++, the second number inside the brackets is the element count, not the final index.) However, the second mapping for A needs to be be map(tofrom: A[0:N][0:N]) rather than map(from: A[1:N-2][1:N-2]) in order preserve the boundary values of A.

If you compile and run this code, you will find it runs slower, even though some unnecessary data transfers have been eliminated. This is because specifying more precise data movement may come at some cost in efficiency, as it is generally faster to transfer data in large, contiguous chunks.

Problem 3: Create a Data Region

Now that we've explicitly mapped data in the offload directive, you might wonder, is it even necessary to map data in every offload region? Can we just map A and Anew once and be done? For the final step, use what we know about persistent data mapping, and implement the solution. You should attempt this problem before looking at the solution below.

Explanation: The solution is to create a data region that maps A and Anew before any iterations are performed. We can then remove the explicit mappings from the loops, since A and Anew are already present on the device. Note that A takes tofrom only because we need A as the output at the end, while Anew can take the alloc mapping type because its initial values are not important; it only stores intermediate values of A.

Footnote: Each of the codes presented above prints out the final value of the Jacobi iterations at the exact center of the grid, A[N/2][N/2]. Analytically, it can be shown that this value should be 0.25, yet the value reported is 0.000005, which is not even close. The reason for the discrepancy is that the Jacobi method is extremely slow to converge, and an accurate answer would require millions of iterations. A better method for this problem is Red-Black Successive Over-Relaxation (SOR), which achieves a far closer approximation without doing much more work. If you are interested in trying out an implementation of this method that has been parallelized using OpenMP offload to the GPU, then download, compile, and run rbsor_off.c.

Credit: This exercise was inspired by a presentation by NVIDIA at SC17.