There are several methods to measure Pi which is a constant that is used to calculate the circumference of a circle from that circle's radius or diameter. One of these methods is the Leibniz formula for Pi:
Follow the lecture / obtain the corresponding code and write the CUDA code by replacing the integration method with the Leibniz formula (Leibniz_pi.c). More information about the Leibniz formula for Pi.
- Your code should follow the lecture through the "Reduction across Blocks" slide.
- Your global_reduce kernel should be able to handle the case where the number of threads is not a power of 2. (This was not covered in the slides, think about it!)
- Your code does not need to have MPI
- CUDA only works on compute nodes, not the login node, so be sure to get a compute node first with srun.
- In your p2.Makefile you may use -Wno-deprecated-gpu-targets to suppress a
warning about deprecated architectures. Do not worry about this warning, everything is fine.
nvcc p2.cu -o p2 -O3 -lm -Wno-deprecated-gpu-targets
Turn in p2.cu (CUDA file), p2.Makefile and p2.README (Explaining the implementation).
We will extend the methods of the last HW into two dimensions.
Download, extract, compile the code lake.tar
This program models the surface of a lake, where some pebbles have been thrown onto the surface. The program works as follows. In the spatial domain, a centralized finite difference is used to inform a zone of how to update itself using the information from its neighbors
The time domain does something similarly, but here using information from the previous two times
The program runs two versions of the algorithm, a CPU version, and a skeleton GPU version. Your
task is to fill in the GPU algorithm to solve the same problem.
Instructions:
V0:
- Reserve one node with one processor. (Each node in the ARC cluster has its own nVidia GPU.)
srun -N1 -n1 --pty /bin/bash
- Run the lake program
./lake {npoints} {npebbles} {end_time} {nthreads}npoints defines the grid size (npoints x npoints), npebbles is the number of pebbles that are generated in the program, end_time is the final time of the simulation, and nthreads will be used withe the GPU implementation.The following runs on a grid of (128 x 128), with 5 pebbles, for 1.0 seconds, using 8 GPU threads (implemented later):
./lake 128 5 1.0 8 Running ./lake with (128 x 128) grid, until 1.000000, with 8 threads CPU took 0.284713 seconds GPU computation: 0.003168 msec GPU took 0.327409 seconds
- View the output in a heatmap with gnuplot:
The output files
lake_i.dat lake_f.dat
can be converted into a .png image using the gnuplot script heatmap.gnu. Rungnuplot heatmap.gnu
This will create the files lake_i.png(the initial configuration), lake_f.png (the final configuration) in the directory.
V1:
- So far you used 5-point stencil for discretizing the wave equation (look at the evolve
method in lake.cu). In this section, you extend that and use 13-point stencil. Under this
configuration, 8 more neighbors are involved in the communication:
- The pseudocode for 5-point is as following:
un[idx] = 2*uc[idx] - uo[idx] + VSQR *(dt * dt) *(( WEST + EAST + NORTH + SOUTH - 4 * uc[idx])/(h * h) + f(pebbles[idx],t)); - The pseudocode for 13-point is as following:
un[idx] = 2*uc[idx] - uo[idx] + VSQR *(dt * dt) * (( WEST + EAST + NORTH + SOUTH + 0.25*(NORTHWEST + NORTHEAST + SOUTHWEST + SOUTHEAST) + 0.125*(WESTWEST + EASTEAST + NORTHNORTH + SOUTHSOUTH) - 6 * uc[idx])/(h * h) + f(pebbles[idx],t)); - You need to calculate the indices for neighbors (NORTHWEST, NORTHEAST, SOUTHWEST, SOUTHEAST, NORTHNORTH, EASTEAST, SOUTHSOUTH, WESTWEST). Note that a lake surface is a grid (2-dimension), and it is represented as a vector (1-dimension) in the code. Make sure the calculated indices do not go outside the bounds of the vector!
- Implement evolve13pt method that performs 13-point stencil under CPU.
- Perform experiments and compare and contrast the two versions (5-point vs 13-point). What are the tradeoffs? Note that the 13-point version evolves faster and you need to run it with shorter end_time when comparing with 5-point version.
- Turn in plots and a description of your observation.
- Use 13-point stencil for the rest of homework.
- Fill in the function run_gpu in the file lakegpu.cu to run the same algorithm
as the cpu version, but using CUDA kernels. The grid will be decomposed on the GPU into 2D blocks.
The program takes as an argument nthreads. This will be the number of threads per block used on the GPU. So, for instance, with nthreads=8, and a domain of grid points (npoints=128 x 128), you will create (npoints/nthreads)x(npoints/nthreads) = (16 x 16) blocks, with (8 x 8) threads on each block.
- You will time your CUDA implementation using cudaEventXXX() API. Be sure to start timing before the first memcpy to the GPU and stop after the last memcpy off of the GPU.
- Compare the CPU/GPU runs for varying grid sizes (16, 32, 64, 128, ..., 1024, etc.)
- Create an MPI version of your program the further decomposes your grid based on processor rank.
- Use 4 nodes (i.e. 4 GPU's) in your implementation.
srun -N4 -n4 --pty /bin/bash
- Each node should communicate boundary information to the appropriate neighbor, then run the CUDA kernel during a time-step (one iteration of evolve).
- Have each node output it's own data file, labeled as
lake_f_0.dat //node 0 lake_f_1.dat //node 1 //etc.
Include in p3.README a discussion of your results. Your discussion should include answering the following questions:
- How well does your algorithm scale on the GPU? Do you find cases (grid size, thread number, etc.) where the GPU implementation does not scale well? Why?
- In the serial code, compare your CPU and GPU runtimes for different grid sizes. When is the GPU better, and when is it worse?
- Integrating CUDA and MPI involves more sophisticated code. What problems did you encounter? How did you fix them?
Hints:
- Here is a sample p3.Makefile for both V2 and V3.
- If you are interested in the particular math behind this algorithm, here is a good introduction. In particular, we are solving the 2D wave equation with sources using finite differencing.
- The code stores the initial configurations of u^0 and u^1 in the variables
double *u_i0; //u^0 double *u_i1; //u^1
These are passed to both the run_cpu and run_gpu routines; both routines should produce the same results. - Create the proper device memory spaces; uc, un, uo should all exist on the device. Between iterations, they should be copied on the device (not sent back to the host). In addition, all data in u_i0, u_i1, as well as pebbles should exist on the device, along with anything the algorithm needs for each iteration.
- In the MPI Implementation, consider ways to speed up each iteration. Do you have to wait on the GPU kernel to finish executing before you exchange boundary information? Consider how you will update the boundary grid zones - do you do it on the CPU after the kernel updates the inner zones or copy over that information to the GPU before the kernel runs? You are encouraged to experiment with different methods to find the most efficient.
Turn in p3.README, lake.cu/lakegpu.cu (version V2), lake_mpi.cu/lakegpu_mpi.cu (version V3), p3.Makefile