diff --git a/README.md b/README.md index a82ea0f..c826381 100644 --- a/README.md +++ b/README.md @@ -3,211 +3,70 @@ CUDA Stream Compaction **University of Pennsylvania, CIS 565: GPU Programming and Architecture, Project 2** -* (TODO) YOUR NAME HERE -* Tested on: (TODO) Windows 22, i7-2222 @ 2.22GHz 22GB, GTX 222 222MB (Moore 2222 Lab) +* Megan Moore +* Tested on: Windows 7, i7-4770 @ 3.40GHz 16GB (Moore 100 Lab C) + +``` + +**************** +** SCAN TESTS ** +**************** + [ 3 29 33 19 0 16 10 40 39 50 44 30 9 ... 4 -858993460 ] +==== cpu scan, power-of-two ==== + [ 0 3 32 65 84 84 100 110 150 189 239 283 313 ... 6684 6688 ] +==== cpu scan, non-power-of-two ==== + [ 0 3 32 65 84 84 100 110 150 189 239 283 313 ... 6613 6626 ] + passed +==== naive scan, power-of-two ==== + passed +==== naive scan, non-power-of-two ==== + passed +==== work-efficient scan, power-of-two ==== + passed +==== work-efficient scan, non-power-of-two ==== + passed +==== thrust scan, power-of-two ==== + passed +==== thrust scan, non-power-of-two ==== + passed + +***************************** +** STREAM COMPACTION TESTS ** +***************************** + [ 4 3 0 3 4 2 3 2 3 1 1 1 4 ... 3 -858993460 ] +==== cpu compact without scan, power-of-two ==== + [ 4 3 3 4 2 3 2 3 1 1 1 4 3 ... 3 -858993460 ] + passed +==== cpu compact without scan, non-power-of-two ==== + [ 4 3 3 4 2 3 2 3 1 1 1 4 3 ... 4 4 ] + passed +==== cpu compact with scan ==== + [ 4 3 3 4 2 3 2 3 1 1 1 4 3 ... 3 -858993460 ] + passed +==== work-efficient compact, power-of-two ==== + passed +==== work-efficient compact, non-power-of-two ==== + passed +Press any key to continue . . . + +``` -### (TODO: Your README) - -Include analysis, etc. (Remember, this is public, so don't put -anything here that you don't want to share with the world.) - -Instructions (delete me) -======================== - -This is due Sunday, September 13 at midnight. - -**Summary:** In this project, you'll implement GPU stream compaction in CUDA, -from scratch. This algorithm is widely used, and will be important for -accelerating your path tracer project. - -Your stream compaction implementations in this project will simply remove `0`s -from an array of `int`s. In the path tracer, you will remove terminated paths -from an array of rays. - -In addition to being useful for your path tracer, this project is meant to -reorient your algorithmic thinking to the way of the GPU. On GPUs, many -algorithms can benefit from massive parallelism and, in particular, data -parallelism: executing the same code many times simultaneously with different -data. - -You'll implement a few different versions of the *Scan* (*Prefix Sum*) -algorithm. First, you'll implement a CPU version of the algorithm to reinforce -your understanding. Then, you'll write a few GPU implementations: "naive" and -"work-efficient." Finally, you'll use some of these to implement GPU stream -compaction. - -**Algorithm overview & details:** There are two primary references for details -on the implementation of scan and stream compaction. - -* The [slides on Parallel Algorithms](https://github.com/CIS565-Fall-2015/cis565-fall-2015.github.io/raw/master/lectures/2-Parallel-Algorithms.pptx) - for Scan, Stream Compaction, and Work-Efficient Parallel Scan. -* GPU Gems 3, Chapter 39 - [Parallel Prefix Sum (Scan) with CUDA](http://http.developer.nvidia.com/GPUGems3/gpugems3_ch39.html). - -Your GPU stream compaction implementation will live inside of the -`stream_compaction` subproject. This way, you will be able to easily copy it -over for use in your GPU path tracer. - - -## Part 0: The Usual - -This project (and all other CUDA projects in this course) requires an NVIDIA -graphics card with CUDA capability. Any card with Compute Capability 2.0 -(`sm_20`) or greater will work. Check your GPU on this -[compatibility table](https://developer.nvidia.com/cuda-gpus). -If you do not have a personal machine with these specs, you may use those -computers in the Moore 100B/C which have supported GPUs. - -**HOWEVER**: If you need to use the lab computer for your development, you will -not presently be able to do GPU performance profiling. This will be very -important for debugging performance bottlenecks in your program. - -### Useful existing code - -* `stream_compaction/common.h` - * `checkCUDAError` macro: checks for CUDA errors and exits if there were any. - * `ilog2ceil(x)`: computes the ceiling of log2(x), as an integer. -* `main.cpp` - * Some testing code for your implementations. - - -## Part 1: CPU Scan & Stream Compaction - -This stream compaction method will remove `0`s from an array of `int`s. - -In `stream_compaction/cpu.cu`, implement: - -* `StreamCompaction::CPU::scan`: compute an exclusive prefix sum. -* `StreamCompaction::CPU::compactWithoutScan`: stream compaction without using - the `scan` function. -* `StreamCompaction::CPU::compactWithScan`: stream compaction using the `scan` - function. Map the input array to an array of 0s and 1s, scan it, and use - scatter to produce the output. You will need a **CPU** scatter implementation - for this (see slides or GPU Gems chapter for an explanation). - -These implementations should only be a few lines long. - - -## Part 2: Naive GPU Scan Algorithm - -In `stream_compaction/naive.cu`, implement `StreamCompaction::Naive::scan` - -This uses the "Naive" algorithm from GPU Gems 3, Section 39.2.1. We haven't yet -taught shared memory, and you **shouldn't use it yet**. Example 39-1 uses -shared memory, but is limited to operating on very small arrays! Instead, write -this using global memory only. As a result of this, you will have to do -`ilog2ceil(n)` separate kernel invocations. - -Beware of errors in Example 39-1 in the book; both the pseudocode and the CUDA -code in the online version of Chapter 39 are known to have a few small errors -(in superscripting, missing braces, bad indentation, etc.) - -Since the parallel scan algorithm operates on a binary tree structure, it works -best with arrays with power-of-two length. Make sure your implementation works -on non-power-of-two sized arrays (see `ilog2ceil`). This requires extra memory -- your intermediate array sizes will need to be rounded to the next power of -two. - - -## Part 3: Work-Efficient GPU Scan & Stream Compaction - -### 3.1. Scan - -In `stream_compaction/efficient.cu`, implement -`StreamCompaction::Efficient::scan` - -All of the text in Part 2 applies. - -* This uses the "Work-Efficient" algorithm from GPU Gems 3, Section 39.2.2. -* Beware of errors in Example 39-2. -* Test non-power-of-two sized arrays. - -### 3.2. Stream Compaction - -This stream compaction method will remove `0`s from an array of `int`s. - -In `stream_compaction/efficient.cu`, implement -`StreamCompaction::Efficient::compact` - -For compaction, you will also need to implement the scatter algorithm presented -in the slides and the GPU Gems chapter. - -In `stream_compaction/common.cu`, implement these for use in `compact`: - -* `StreamCompaction::Common::kernMapToBoolean` -* `StreamCompaction::Common::kernScatter` - - -## Part 4: Using Thrust's Implementation - -In `stream_compaction/thrust.cu`, implement: - -* `StreamCompaction::Thrust::scan` - -This should be a very short function which wraps a call to the Thrust library -function `thrust::exclusive_scan(first, last, result)`. - -To measure timing, be sure to exclude memory operations by passing -`exclusive_scan` a `thrust::device_vector` (which is already allocated on the -GPU). You can create a `thrust::device_vector` by creating a -`thrust::host_vector` from the given pointer, then casting it. - - -## Part 5: Radix Sort (Extra Credit) (+10) - -Add an additional module to the `stream_compaction` subproject. Implement radix -sort using one of your scan implementations. Add tests to check its correctness. - - -## Write-up - -1. Update all of the TODOs at the top of this README. -2. Add a description of this project including a list of its features. -3. Add your performance analysis (see below). - -All extra credit features must be documented in your README, explaining its -value (with performance comparison, if applicable!) and showing an example how -it works. For radix sort, show how it is called and an example of its output. - -Always profile with Release mode builds and run without debugging. ### Questions * Roughly optimize the block sizes of each of your implementations for minimal run time on your GPU. - * (You shouldn't compare unoptimized implementations to each other!) + + * Four different block sizes (128, 256, 512, 1024) were tested against four different array sizes (128, 256, 512, 1024). Based on the cudaEvent timing, none of the different combinations led to a notable difference in times. When the times did differ, it was only by a few tenths of a millisecond. Also, the speed ups that occured with a blocksize for one of the scan functions, caused other scan functions to slow down. Therefore, I used a consistent blocksize of 128 for all scan functions. * Compare all of these GPU Scan implementations (Naive, Work-Efficient, and Thrust) to the serial CPU version of Scan. Plot a graph of the comparison (with array size on the independent axis). - * You should use CUDA events for timing. Be sure **not** to include any - explicit memory operations in your performance measurements, for - comparability. - * To guess at what might be happening inside the Thrust implementation, take - a look at the Nsight timeline for its execution. + ![](images/Graph.png "Array size analysis") + * I tried to use the chrono implimentation to time the CPU scan function. However, I continuously got 0 nanoseconds for each test. + * The thrust application is taking a much longer time than all the other implimentations. Seeing as the cudaEvent_t's, start and stop, occur right before and after the function call to thrust::exclusive_scan(), it could be possible that these events are picking up memory transfer timing. Whereas, with all the other function, I was able to record the time without the memory transfer time being included. One other interesting thing about the thrust application, is that the time for the non-power of 2 array is significantly lower than the time for the power of 2 array. This, in theory, makes sense because a smaller array should take less time. However, with all the other applications, they take approximately the same amount of time because the kernel is 2^(ilog2ceil(n)) times (for an array of size n). This can be seen in the graph, as the naive/work-efficient power of 2 and non-power of 2 lines are almost identicle. In the thrust application, they must be allocating their memory better, as they clearly do not have to call the kernel 2^(ilog2ceil(n)) times. * Write a brief explanation of the phenomena you see here. - * Can you find the performance bottlenecks? Is it memory I/O? Computation? Is - it different for each implementation? - -* Paste the output of the test program into a triple-backtick block in your - README. - * If you add your own tests (e.g. for radix sort or to test additional corner - cases), be sure to mention it explicitly. - -These questions should help guide you in performance analysis on future -assignments, as well. - -## Submit - -If you have modified any of the `CMakeLists.txt` files at all (aside from the -list of `SOURCE_FILES`), you must test that your project can build in Moore -100B/C. Beware of any build issues discussed on the Google Group. + * It is interesting that the naive implementation is running faster than the work-efficient implementation. This seems to go against what we think should be happening. An idea of what could possibly be causing this is that the work-efficient implementation requires two function calls to kernels. If the array is not so large, that it requires multiple blocks, then it could finish the naive scan implementation in the same time it finishes just one of the functions (up/down sweep) in the work-efficient scan implementation. + * I would guess that the bottlenecks are occuring at the memory transfers. Especially in the work-efficient implementation, because it has so many different arrays (idata, odata, bools, and indices) it requires more memory allocation and transfers between the host and device. I would not think the computation is a large part of the time, considering it is only addition, and the amount of addition that is necessary decreases throughout the function (in all of the implementations). -1. Open a GitHub pull request so that we can see that you have finished. - The title should be "Submission: YOUR NAME". -2. Send an email to the TA (gmail: kainino1+cis565@) with: - * **Subject**: in the form of `[CIS565] Project 2: PENNKEY` - * Direct link to your pull request on GitHub - * In the form of a grade (0-100+) with comments, evaluate your own - performance on the project. - * Feedback on the project itself, if any. diff --git a/images/Graph.png b/images/Graph.png new file mode 100644 index 0000000..d618525 Binary files /dev/null and b/images/Graph.png differ diff --git a/src/main.cpp b/src/main.cpp index efc8c06..c29b215 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -16,7 +16,7 @@ #include "testing_helpers.hpp" int main(int argc, char* argv[]) { - const int SIZE = 1 << 8; + const int SIZE = 1 << 8; const int NPOT = SIZE - 3; int a[SIZE], b[SIZE], c[SIZE]; diff --git a/src/testing_helpers.hpp b/src/testing_helpers.hpp index 69ac2e9..20bb30b 100644 --- a/src/testing_helpers.hpp +++ b/src/testing_helpers.hpp @@ -17,7 +17,7 @@ void printDesc(const char *desc) { template void printCmpResult(int n, T *a, T *b) { - printf(" %s \n", + printf(" %s \n", cmpArrays(n, a, b) ? "FAIL VALUE" : "passed"); } diff --git a/stream_compaction/common.cu b/stream_compaction/common.cu index fe872d4..e03178e 100644 --- a/stream_compaction/common.cu +++ b/stream_compaction/common.cu @@ -23,7 +23,13 @@ namespace Common { * which map to 0 will be removed, and elements which map to 1 will be kept. */ __global__ void kernMapToBoolean(int n, int *bools, const int *idata) { - // TODO + int thrId = threadIdx.x + (blockIdx.x * blockDim.x); + if (idata[thrId] == 0) { + bools[thrId] = 0; + } + else { + bools[thrId] = 1; + } } /** @@ -32,7 +38,12 @@ __global__ void kernMapToBoolean(int n, int *bools, const int *idata) { */ __global__ void kernScatter(int n, int *odata, const int *idata, const int *bools, const int *indices) { - // TODO + int thrId = threadIdx.x + (blockIdx.x * blockDim.x); + if (thrId < n) { + if (bools[thrId] == 1) { + odata[indices[thrId]] = idata[thrId]; + } + } } } diff --git a/stream_compaction/cpu.cu b/stream_compaction/cpu.cu index e600c29..00be952 100644 --- a/stream_compaction/cpu.cu +++ b/stream_compaction/cpu.cu @@ -1,4 +1,6 @@ #include +#include +#include #include "cpu.h" namespace StreamCompaction { @@ -9,7 +11,11 @@ namespace CPU { */ void scan(int n, int *odata, const int *idata) { // TODO - printf("TODO\n"); + odata[0] = 0; + for (int i = 1; i < n; i++) { + odata[i] = idata[i-1] + odata[i-1]; + } + } /** @@ -18,8 +24,19 @@ void scan(int n, int *odata, const int *idata) { * @returns the number of elements remaining after compaction. */ int compactWithoutScan(int n, int *odata, const int *idata) { - // TODO - return -1; + using namespace std; + int sum = 0; + auto begin = std::chrono::high_resolution_clock::now(); + for (int i = 0; i < n; i++) { + if (idata[i] != 0) { + odata[sum] = idata[i]; + sum++; + } + + } + auto end = std::chrono::high_resolution_clock::now(); + //std::cout << std::chrono::duration_cast(end-begin).count() << "ns" << std::endl; + return sum; } /** @@ -28,8 +45,23 @@ int compactWithoutScan(int n, int *odata, const int *idata) { * @returns the number of elements remaining after compaction. */ int compactWithScan(int n, int *odata, const int *idata) { - // TODO - return -1; + int* bools = new int[n]; + for (int i = 0; i < n; i++) { + if (idata[i] == 0) { + bools[i] = 0; + } + else { + bools[i] = 1; + } + } + int* scanArray = new int[n]; + scan(n, scanArray, bools); + for (int i = 0; i < n; i++) { + if (bools[i] == 1) { + odata[scanArray[i]] = idata[i]; + } + } + return scanArray[n - 1] + bools[n - 1]; } } diff --git a/stream_compaction/efficient.cu b/stream_compaction/efficient.cu index b2f739b..298afa1 100644 --- a/stream_compaction/efficient.cu +++ b/stream_compaction/efficient.cu @@ -6,14 +6,91 @@ namespace StreamCompaction { namespace Efficient { -// TODO: __global__ +int* g_odata; +int* g_idata; +int* dev_bools; +int* dev_indices; +__global__ void generate_zeros(int *data) { + int i = threadIdx.x; + data[i] = 0; +} + +__global__ void set_zero(int size, int n, int *data) { + int i = threadIdx.x; + if (i >= n - 1) { + data[i] = 0; + } +} +__global__ void kern_up_sweep(int n, int *odata, const int *idata, int layer) { + int thrId = threadIdx.x + (blockIdx.x * blockDim.x); + if ((thrId < n) && (thrId%layer == 0)) { + odata[thrId + layer - 1] += idata[thrId + (layer / 2) - 1]; + } +} + +__global__ void kern_down_sweep(int n, int *odata, const int *idata, int layer) { + int thrId = threadIdx.x + (blockIdx.x * blockDim.x); + if ((thrId < n) && (thrId%layer == 0)) { + int temp = idata[thrId + (layer / 2) - 1]; + odata[thrId + (layer / 2) - 1] = idata[thrId + layer - 1]; + odata[thrId + layer - 1] += temp; + } + + +} /** * Performs prefix-sum (aka scan) on idata, storing the result into odata. */ void scan(int n, int *odata, const int *idata) { - // TODO - printf("TODO\n"); + cudaEvent_t start, stop; + cudaEventCreate(&start); + cudaEventCreate(&stop); + + int blockSize = 128; + int numBlocks = ceil((float)n / (float)blockSize); + int powTwo = pow(2, ilog2ceil(n)); + dim3 fullBlocksPerGrid((powTwo + blockSize - 1) / blockSize); + cudaMalloc((void**)&g_odata, powTwo * sizeof(int)); + cudaMalloc((void**)&g_idata, powTwo * sizeof(int)); + + generate_zeros<<<1, powTwo>>>(g_odata); + generate_zeros<<<1, powTwo>>>(g_idata); + + int* scanArray = new int[n]; + //scanArray[0] = 0; + for (int i = 0; i < n; i++) { + scanArray[i] = idata[i]; + } + + cudaMemcpy(g_odata, odata, n*sizeof(int), cudaMemcpyHostToDevice); + cudaMemcpy(g_idata, scanArray, n*sizeof(int), cudaMemcpyHostToDevice); + + cudaEventRecord(start); + for (int d = 0; d <= ilog2ceil(n) - 1; d++) { + int layer = pow(2, d + 1); + g_odata = g_idata; + + kern_up_sweep<<>>(powTwo, g_odata, g_idata, layer); + g_idata = g_odata; + } + + + set_zero<<<1, powTwo>>>(powTwo, n, g_idata); + + for (int d = ilog2ceil(n) - 1; d >= 0; d--) { + int layer = pow(2, d + 1); + g_odata = g_idata; + kern_down_sweep<<>>(powTwo, g_odata, g_idata, layer); + g_idata = g_odata; + } + cudaEventRecord(stop); + + cudaEventSynchronize(stop); + float milliseconds = 0; + cudaEventElapsedTime(&milliseconds, start, stop); + //printf("%f - ", milliseconds); + cudaMemcpy(odata, g_odata, n*sizeof(int), cudaMemcpyDeviceToHost); } /** @@ -25,9 +102,35 @@ void scan(int n, int *odata, const int *idata) { * @param idata The array of elements to compact. * @returns The number of elements remaining after compaction. */ -int compact(int n, int *odata, const int *idata) { - // TODO - return -1; +int compact(int n, int *odata, const int *idata) { + int blockSize = 128; + int numBlocks = ceil((float)n / (float)blockSize); + int powTwo = pow(2, ilog2ceil(n)); + dim3 fullBlocksPerGrid((powTwo + blockSize - 1) / blockSize); + + cudaMalloc((void**)&g_odata, powTwo * sizeof(int)); + cudaMalloc((void**)&g_idata, powTwo * sizeof(int)); + cudaMalloc((void**)&dev_bools, powTwo * sizeof(int)); + cudaMalloc((void**)&dev_indices, powTwo * sizeof(int)); + + int* indices = new int[n]; + int* bools = new int[n]; + + cudaMemcpy(g_idata, idata, n*sizeof(int), cudaMemcpyHostToDevice); + + Common::kernMapToBoolean<<>>(powTwo, dev_bools, g_idata); + cudaMemcpy(bools, dev_bools, n*sizeof(int), cudaMemcpyDeviceToHost); + + scan(n, indices, bools); + + cudaMemcpy(dev_indices, indices, n*sizeof(int), cudaMemcpyHostToDevice); + + cudaMemcpy(g_idata, idata, n*sizeof(int), cudaMemcpyHostToDevice); + Common::kernScatter<<>>(powTwo, g_odata, g_idata, dev_bools, dev_indices); + + cudaMemcpy(odata, g_odata, n*sizeof(int), cudaMemcpyDeviceToHost); + + return indices[n-1] + bools[n-1]; } } diff --git a/stream_compaction/naive.cu b/stream_compaction/naive.cu index 3d86b60..e778738 100644 --- a/stream_compaction/naive.cu +++ b/stream_compaction/naive.cu @@ -6,14 +6,56 @@ namespace StreamCompaction { namespace Naive { -// TODO: __global__ +int *g_odata; +int *g_idata; +__global__ void kern_scan(int n, int *odata, const int *idata, int layer) { + int thrId = threadIdx.x + (blockIdx.x * blockDim.x); + if (thrId >= layer) { + + odata[thrId] = idata[thrId - layer] + idata[thrId]; + + } + else { + odata[thrId] = idata[thrId]; + } + +} /** * Performs prefix-sum (aka scan) on idata, storing the result into odata. */ void scan(int n, int *odata, const int *idata) { - // TODO - printf("TODO\n"); + cudaMalloc((void**)&g_odata, n * sizeof(int)); + cudaMalloc((void**)&g_idata, n * sizeof(int)); + cudaEvent_t start, stop; + cudaEventCreate(&start); + cudaEventCreate(&stop); + int blockSize = 128; + dim3 numBlocks = (int)ceil((float)n / (float)blockSize); + int powTwo = pow(2, ilog2ceil(n)); + dim3 fullBlocksPerGrid((powTwo + blockSize - 1) / blockSize); + + int* scanArray = new int[n]; + scanArray[0] = 0; + for (int i = 1; i < n; i++) { + scanArray[i] = idata[i - 1]; + } + + cudaMemcpy(g_odata, odata, n*sizeof(int), cudaMemcpyHostToDevice); + cudaMemcpy(g_idata, scanArray, n*sizeof(int), cudaMemcpyHostToDevice); + cudaEventRecord(start); + for (int d = 1; d <= ilog2ceil(n); d++) { + int layer = pow(2, d - 1); + kern_scan<<>>(n, g_odata, g_idata, layer); + g_idata = g_odata; + } + cudaEventRecord(stop); + + cudaEventSynchronize(stop); + float milliseconds = 0; + cudaEventElapsedTime(&milliseconds, start, stop); + //printf("%f - ", milliseconds); + cudaMemcpy(odata, g_odata, n*sizeof(int), cudaMemcpyDeviceToHost); } } diff --git a/stream_compaction/thrust.cu b/stream_compaction/thrust.cu index d8dbb32..d182d1b 100644 --- a/stream_compaction/thrust.cu +++ b/stream_compaction/thrust.cu @@ -16,6 +16,24 @@ void scan(int n, int *odata, const int *idata) { // TODO use `thrust::exclusive_scan` // example: for device_vectors dv_in and dv_out: // thrust::exclusive_scan(dv_in.begin(), dv_in.end(), dv_out.begin()); + thrust::device_vector dv_in(idata, idata + n); + thrust::device_vector dv_out(odata, odata + n); + + cudaEvent_t start, stop; + cudaEventCreate(&start); + cudaEventCreate(&stop); + + cudaEventRecord(start); + thrust::exclusive_scan(dv_in.begin(), dv_in.end(), dv_out.begin()); + cudaEventRecord(stop); + + cudaEventSynchronize(stop); + float milliseconds = 0; + cudaEventElapsedTime(&milliseconds, start, stop); + //printf("%f - ", milliseconds); + + thrust::copy(dv_out.begin(), dv_out.end(), odata); + } }