diff --git a/README.md b/README.md index 0e38ddb1..13450701 100644 --- a/README.md +++ b/README.md @@ -1,14 +1,85 @@ -CUDA Stream Compaction -====================== +# University of Pennsylvania, CIS 5650: GPU Programming and Architecture +## Project 2 - Stream Compaction -**University of Pennsylvania, CIS 565: GPU Programming and Architecture, Project 2** +* Zwe Tun + * LinkedIn: https://www.linkedin.com/in/zwe-tun-6b7191256/ +* Tested on: Intel(R) i7-14700HX, 2100 Mhz, RTX 5060 Laptop -* (TODO) YOUR NAME HERE - * (TODO) [LinkedIn](), [personal website](), [twitter](), etc. -* Tested on: (TODO) Windows 22, i7-2222 @ 2.22GHz 22GB, GTX 222 222MB (Moore 2222 Lab) +## Overview +Stream compaction is a widely used algorithm with applications in areas such as image compression, data filtering, and parallel computing. This project explores four different implementations of stream compaction: -### (TODO: Your README) +- CPU compaction without scan -Include analysis, etc. (Remember, this is public, so don't put -anything here that you don't want to share with the world.) +- CPU compaction with scan +- GPU-based compaction + +Each method demonstrates a different approach to optimizing performance and leveraging hardware capabilities. The goal is to analyze the performance trade-offs and efficiency of each implementation. + +## Implementation + +### CPU Compaction Without Scan +This implementation uses an iterative approach that incrementally places non-zero elements from the input array into an output array. An index counter tracks the next available position in the output. + +### CPU Compaction With Scan +In this approach, an exclusive scan is introduced to calculate the output indices of the non-zero elements. The scan is performed using an iterative loop based on the following formula: + +x[i] = x[i - 1] + input[i - 1] + +Once the scan array is generated, we pass over the input. For each element i, if input[i] is non-zero, we place it into the output array at the position defined by scan[i]. + +### GPU-based Compaction +The GPU version builds on the scan-based CPU approach. Its performance depends on two scan implementations: + +**Naive Scan:** +Uses GPU threads to compute the scan in multiple layers. At each layer, threads read from specific indices and write to new indices, all in place, gradually building up the scan. + +**Work-Efficient Scan:** +Uses two phases: + +Up-sweep: +Performs a parallel reduction to build a balanced binary tree of partial sums. + +Down-sweep: +Traverses back down the tree to compute the exclusive scan in-place. At each pass, a node passes its value to its left child, and sets the right child to the sum of the previous left child’s value and its value. + +Once the scan is complete, the result is copied back to the host (CPU), where the compaction is performed + +## Performance Analysis +To evaluate and compare each implementation, a benchmark is used based on the Thrust API, specifically leveraging thrust::exclusive_scan. Performance is measured for both CPU and GPU executions by wrapping each in a performance timer class. The CPU execution time is recorded using std::chrono for high-precision timing, while GPU timing is measured using CUDA events. + +Note, memory allocation and transfer operations (cudaMalloc, cudaMemcpy, etc.) are excluded from the timing results. + + +  + + +![Stream Compaction](img/block1.png) + +*Performance across various GPU block sizes for 256 array size, 128 block size used for futher tests.* + + +  + + +![Stream Compaction](img/Scan.png) +![Stream Compaction](img/compact.png) + +*Performance results for 8,388,608 (2^23) sized array, 128 block size.* + +  + + +![Stream Compaction](img/size.png) + +*Performance results for 8,388,608 (2^23) non-power-two vs power-of-two sized array, 128 block size.* + +  + +## Questions +### To guess at what might be happening inside the Thrust implementation (e.g. allocation, memory copy), take a look at the Nsight timeline for its execution. Your analysis here doesn't have to be detailed, since you aren't even looking at the code for the implementation. Write a brief explanation of the phenomena you see here. + +Looking at the timeline, we can observe numerous memory operations and multiple kernel launches. These kernel launches appear to be synchronized and follow a pattern similar to the scan approaches used in this project's algorithm. + +### Can you find the performance bottlenecks? Is it memory I/O? Computation? Is it different for each implementation? +For the CPU implementation, the bottleneck is predominantly computation based on the size of the array. Since the CPU executes sequentially, the main cost is due to math operations at each interation. For the GPU, memory I/O will make up the majority of the cost. Computation is parallelized through threads so global memory accesses and copies will be limiting. diff --git a/img/Scan.png b/img/Scan.png new file mode 100644 index 00000000..5918c31f Binary files /dev/null and b/img/Scan.png differ diff --git a/img/block.png b/img/block.png new file mode 100644 index 00000000..c317a752 Binary files /dev/null and b/img/block.png differ diff --git a/img/block1.png b/img/block1.png new file mode 100644 index 00000000..b284724b Binary files /dev/null and b/img/block1.png differ diff --git a/img/compact.png b/img/compact.png new file mode 100644 index 00000000..73a34182 Binary files /dev/null and b/img/compact.png differ diff --git a/img/size.png b/img/size.png new file mode 100644 index 00000000..d61d8c17 Binary files /dev/null and b/img/size.png differ diff --git a/src/main.cpp b/src/main.cpp index 3d5c8820..74a574bd 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -15,9 +15,9 @@ const int SIZE = 1 << 8; // feel free to change the size of array const int NPOT = SIZE - 3; // Non-Power-Of-Two -int *a = new int[SIZE]; -int *b = new int[SIZE]; -int *c = new int[SIZE]; +int* a = new int[SIZE]; +int* b = new int[SIZE]; +int* c = new int[SIZE]; int main(int argc, char* argv[]) { // Scan tests @@ -29,7 +29,7 @@ int main(int argc, char* argv[]) { genArray(SIZE - 1, a, 50); // Leave a 0 at the end to test that edge case a[SIZE - 1] = 0; - printArray(SIZE, a, true); + // printArray(SIZE, a, true); // initialize b using StreamCompaction::CPU::scan you implement // We use b for further comparison. Make sure your StreamCompaction::CPU::scan is correct. @@ -44,7 +44,7 @@ int main(int argc, char* argv[]) { printDesc("cpu scan, non-power-of-two"); StreamCompaction::CPU::scan(NPOT, c, a); printElapsedTime(StreamCompaction::CPU::timer().getCpuElapsedTimeForPreviousOperation(), "(std::chrono Measured)"); - printArray(NPOT, c, true); + //printArray(NPOT, c, true); printCmpResult(NPOT, b, c); zeroArray(SIZE, c); @@ -81,18 +81,27 @@ int main(int argc, char* argv[]) { //printArray(NPOT, c, true); printCmpResult(NPOT, b, c); + /* zeroArray(SIZE, c); + int* myTest = new int[SIZE] {0, 1, 2, 3, 4, 5, 6, 7}; + int* myResult = new int[SIZE] {0, 1, 2, 6, 4, 9, 6, 28}; + printDesc("work-efficient scan, upSweep test"); + StreamCompaction::Efficient::scan(SIZE, c, myTest); + printElapsedTime(StreamCompaction::Efficient::timer().getGpuElapsedTimeForPreviousOperation(), "(CUDA Measured)"); + printArray(SIZE, c, true); + printCmpResult(SIZE, myResult, c);*/ + zeroArray(SIZE, c); printDesc("thrust scan, power-of-two"); StreamCompaction::Thrust::scan(SIZE, c, a); printElapsedTime(StreamCompaction::Thrust::timer().getGpuElapsedTimeForPreviousOperation(), "(CUDA Measured)"); - //printArray(SIZE, c, true); + // printArray(SIZE, c, true); printCmpResult(SIZE, b, c); zeroArray(SIZE, c); printDesc("thrust scan, non-power-of-two"); StreamCompaction::Thrust::scan(NPOT, c, a); printElapsedTime(StreamCompaction::Thrust::timer().getGpuElapsedTimeForPreviousOperation(), "(CUDA Measured)"); - //printArray(NPOT, c, true); + // printArray(NPOT, c, true); printCmpResult(NPOT, b, c); printf("\n"); @@ -104,7 +113,7 @@ int main(int argc, char* argv[]) { genArray(SIZE - 1, a, 4); // Leave a 0 at the end to test that edge case a[SIZE - 1] = 0; - printArray(SIZE, a, true); + // printArray(SIZE, a, true); int count, expectedCount, expectedNPOT; @@ -115,7 +124,7 @@ int main(int argc, char* argv[]) { count = StreamCompaction::CPU::compactWithoutScan(SIZE, b, a); printElapsedTime(StreamCompaction::CPU::timer().getCpuElapsedTimeForPreviousOperation(), "(std::chrono Measured)"); expectedCount = count; - printArray(count, b, true); + // printArray(count, b, true); printCmpLenResult(count, expectedCount, b, b); zeroArray(SIZE, c); @@ -123,21 +132,21 @@ int main(int argc, char* argv[]) { count = StreamCompaction::CPU::compactWithoutScan(NPOT, c, a); printElapsedTime(StreamCompaction::CPU::timer().getCpuElapsedTimeForPreviousOperation(), "(std::chrono Measured)"); expectedNPOT = count; - printArray(count, c, true); + // printArray(count, c, true); printCmpLenResult(count, expectedNPOT, b, c); zeroArray(SIZE, c); printDesc("cpu compact with scan"); count = StreamCompaction::CPU::compactWithScan(SIZE, c, a); printElapsedTime(StreamCompaction::CPU::timer().getCpuElapsedTimeForPreviousOperation(), "(std::chrono Measured)"); - printArray(count, c, true); + // printArray(count, c, true); printCmpLenResult(count, expectedCount, b, c); zeroArray(SIZE, c); printDesc("work-efficient compact, power-of-two"); count = StreamCompaction::Efficient::compact(SIZE, c, a); printElapsedTime(StreamCompaction::Efficient::timer().getGpuElapsedTimeForPreviousOperation(), "(CUDA Measured)"); - //printArray(count, c, true); + // printArray(count, c, true); printCmpLenResult(count, expectedCount, b, c); zeroArray(SIZE, c); @@ -151,4 +160,4 @@ int main(int argc, char* argv[]) { delete[] a; delete[] b; delete[] c; -} +} \ No newline at end of file diff --git a/stream_compaction/common.cu b/stream_compaction/common.cu index 2ed6d630..205faa81 100644 --- a/stream_compaction/common.cu +++ b/stream_compaction/common.cu @@ -1,4 +1,5 @@ #include "common.h" +#include void checkCUDAErrorFn(const char *msg, const char *file, int line) { cudaError_t err = cudaGetLastError(); @@ -23,7 +24,12 @@ namespace StreamCompaction { * 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 index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index < n) { + bools[index] = (idata[index] != 0) ? 1 : 0; + } + } /** @@ -32,7 +38,12 @@ namespace StreamCompaction { */ __global__ void kernScatter(int n, int *odata, const int *idata, const int *bools, const int *indices) { - // TODO + + int index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index < n && bools[index]) { + odata[indices[index]] = idata[index]; + } + } } diff --git a/stream_compaction/cpu.cu b/stream_compaction/cpu.cu index 719fa115..fb47d248 100644 --- a/stream_compaction/cpu.cu +++ b/stream_compaction/cpu.cu @@ -19,7 +19,14 @@ namespace StreamCompaction { */ void scan(int n, int *odata, const int *idata) { timer().startCpuTimer(); - // TODO + + odata[0] = 0; + for (int i = 1; i < n; i++) { + + odata[i] = odata[i - 1] + idata[i - 1]; + + } + timer().endCpuTimer(); } @@ -30,9 +37,17 @@ namespace StreamCompaction { */ int compactWithoutScan(int n, int *odata, const int *idata) { timer().startCpuTimer(); - // TODO + + int count = 0; + for (int i = 0; i < n; i++) { + if (idata[i] > 0) { + odata[count] = idata[i]; + count++; + } + } + timer().endCpuTimer(); - return -1; + return count; } /** @@ -42,9 +57,49 @@ namespace StreamCompaction { */ int compactWithScan(int n, int *odata, const int *idata) { timer().startCpuTimer(); - // TODO + + + + int count = 0; + int* temp = new int[n]; + int* scanned = new int[n]; + + for (int i = 0; i < n; i++) { + if (idata[i] != 0) { + odata[i] = 1; + } + else { + odata[i] = 0; + } + + + } + + + scanned[0] = 0; + for (int i = 1; i < n; i++) { + + scanned[i] = scanned[i - 1] + odata[i - 1]; + } + + + + for (int i = 0; i < n; i++) { + if (odata[i] == 1) { + + temp[scanned[i]] = idata[i]; + count++; + } + } + + + + for (int i = 0; i < count; i++) { + odata[i] = temp[i]; + } + timer().endCpuTimer(); - return -1; + return count; } } } diff --git a/stream_compaction/efficient.cu b/stream_compaction/efficient.cu index 2db346ee..dff4b152 100644 --- a/stream_compaction/efficient.cu +++ b/stream_compaction/efficient.cu @@ -2,6 +2,7 @@ #include #include "common.h" #include "efficient.h" +#include namespace StreamCompaction { namespace Efficient { @@ -12,13 +13,110 @@ namespace StreamCompaction { return timer; } + + #define blockSize 128 + int* obuffer; + int* ibuffer; + + + __global__ void upSweep(int n, int* idata, int layer) { + int index = threadIdx.x + (blockIdx.x * blockDim.x); + + int skip = 1 << (layer + 1); // powf(2, layer + 1); + + int i = (index) * skip; + if (i + skip - 1 >= n) { + return; + } + + + idata[int(i + skip - 1)] += idata[int(i + (skip >> 1) - 1)]; + + + } + + + __global__ void downSweep(int n, int* idata, int layer) { + int index = threadIdx.x + (blockIdx.x * blockDim.x); + + + int skip = 1 << (layer + 1); // powf(2, layer + 1); + int i = index * skip; + if (i + skip - 1 >= n) { + return; + } + + + + int t = idata[int(i + (skip >> 1) - 1)]; + + + idata[int(i + (skip >> 1) - 1)] = idata[int(i + skip - 1)]; + idata[int(i + skip - 1)] += t; + } + /** * Performs prefix-sum (aka scan) on idata, storing the result into odata. */ void scan(int n, int *odata, const int *idata) { timer().startGpuTimer(); - // TODO + + + int size = 1 << ilog2ceil(n); + + //Init Buffers + int* obuffer; + int* ibuffer; + cudaMalloc((void**)&obuffer, size * sizeof(int)); + cudaMalloc((void**)&ibuffer, size * sizeof(int)); + cudaMemset(ibuffer, 0, (size) * sizeof(int)); + + cudaMemcpy(obuffer, odata, size * sizeof(int), cudaMemcpyHostToDevice); + cudaMemcpy(ibuffer, idata, size * sizeof(int), cudaMemcpyHostToDevice); + + cudaDeviceSynchronize(); + + // up sweep + int numBlocks; + + + for (int layer = 0; layer <= ilog2ceil(size) - 1; layer++) { + int numThreads = size / int(powf(2, layer + 1)); + numBlocks = (numThreads + blockSize - 1) / blockSize; + upSweep<<>>(size, ibuffer, layer); + cudaDeviceSynchronize(); + + + } + + //Set ibuffer[n - 1] = 0 + cudaMemset(ibuffer + size - 1, 0, sizeof(int)); + + + + // down sweep + for (int layer = ilog2ceil(size) - 1; layer >= 0; layer--) { + int numThreads = size / int(powf(2, layer + 1)); + + + numBlocks = (numThreads + blockSize - 1) / blockSize; + downSweep<<>>(size, ibuffer, layer); + cudaDeviceSynchronize(); + + + } timer().endGpuTimer(); + + + + + cudaMemcpy(odata, ibuffer, n * sizeof(int), cudaMemcpyDeviceToHost); + + + + cudaFree(ibuffer); + cudaFree(obuffer); + } /** @@ -31,10 +129,61 @@ namespace StreamCompaction { * @returns The number of elements remaining after compaction. */ int compact(int n, int *odata, const int *idata) { - timer().startGpuTimer(); - // TODO - timer().endGpuTimer(); - return -1; + + // timer().startGpuTimer(); + //Init Buffers + int* boolBufferDevice; + int* obuffer; + int* ibuffer; + int* scanBuffer = new int[n]; + int* scanBufferDevice; + + cudaMalloc((void**)&scanBufferDevice, n * sizeof(int)); + cudaMalloc((void**)&boolBufferDevice, n * sizeof(int)); + cudaMalloc((void**)&ibuffer, n * sizeof(int)); + cudaMalloc((void**)&obuffer, n * sizeof(int)); + + cudaMemcpy(ibuffer, idata, n * sizeof(int), cudaMemcpyHostToDevice); + + int numBlocks = (n + blockSize - 1) / blockSize; + dim3 fullBlocksPerGrid(numBlocks); + //Map to boolean + Common::kernMapToBoolean<<>>(n, boolBufferDevice, ibuffer); + + int* boolBuffer = new int[n]; + cudaMemcpy(boolBuffer, boolBufferDevice, n * sizeof(int), cudaMemcpyDeviceToHost); + + + //Scan boolean array + scan(n, scanBuffer, boolBuffer); + + + cudaMemcpy(scanBufferDevice, scanBuffer, n * sizeof(int), cudaMemcpyHostToDevice); + + Common::kernScatter <<> > (n, obuffer, ibuffer, boolBufferDevice, scanBufferDevice); + ////Scatter results + + + // timer().endGpuTimer(); + + //Get count by looking at last element of scan + last element of bool + //Last element of scan holds number of elements up to n-1, so + //add bool[n-1] to get the full count + int count; + cudaMemcpy(&count, scanBufferDevice + n - 1, sizeof(int), cudaMemcpyDeviceToHost); + int lastBool; + cudaMemcpy(&lastBool, boolBufferDevice + n - 1, sizeof(int), cudaMemcpyDeviceToHost); + count += lastBool; + cudaFree(boolBuffer); + cudaFree(scanBuffer); + cudaMemcpy(odata, obuffer, n * sizeof(int), cudaMemcpyDeviceToHost); + + + + return count; + + + } } } diff --git a/stream_compaction/naive.cu b/stream_compaction/naive.cu index 43088769..c51b9f82 100644 --- a/stream_compaction/naive.cu +++ b/stream_compaction/naive.cu @@ -2,6 +2,8 @@ #include #include "common.h" #include "naive.h" +#include "cmath" +#include namespace StreamCompaction { namespace Naive { @@ -10,16 +12,98 @@ namespace StreamCompaction { { static PerformanceTimer timer; return timer; - } - // TODO: __global__ + } + + #define blockSize 128 + int* obuffer; + int* ibuffer; + + + //Shift right for exclusive scan + __global__ void shiftRight(int n, int* odata, const int* idata) { + int index = threadIdx.x + (blockIdx.x * blockDim.x); + if (index >= n) { + return; + } + if (index == 0) { + odata[0] = 0; + } + else { + odata[index] = idata[index - 1]; + } + } + + __global__ void kernNaiveScan(int n, int* odata, const int* idata, int offset) { + int index = threadIdx.x + (blockIdx.x * blockDim.x); + + if (index >= n) { + return; + } + + + if (index >= offset) { + odata[index] = idata[(index - offset)] + idata[index]; + + } + else { + odata[index] = idata[index]; + } + + + + + + + } /** * Performs prefix-sum (aka scan) on idata, storing the result into odata. */ void scan(int n, int *odata, const int *idata) { timer().startGpuTimer(); - // TODO + + cudaMalloc((void**)&obuffer, n * sizeof(int)); + cudaMalloc((void**)&ibuffer, n * sizeof(int)); + cudaMemcpy(obuffer, odata, n * sizeof(int), cudaMemcpyHostToDevice); + cudaMemcpy(ibuffer, idata, n * sizeof(int), cudaMemcpyHostToDevice); + + int numBlocks = (n + blockSize - 1) / blockSize; + dim3 fullBlocksPerGrid(numBlocks); + + + + + int count = 0; + for (int layer = 1; layer <= ilog2ceil(n); layer++) { + + int offset = powf(2, layer - 1); + + kernNaiveScan<<>>(n, obuffer, ibuffer, offset); + + cudaDeviceSynchronize(); + count++; + int* temp = ibuffer; + ibuffer = obuffer; + obuffer = temp; + + + + + } + shiftRight << > > (n, obuffer, ibuffer); + cudaDeviceSynchronize(); + + + cudaMemcpy(odata, obuffer, n * sizeof(int), cudaMemcpyDeviceToHost); + + cudaFree(obuffer); + cudaFree(ibuffer); + + odata[0] = 0; + + timer().endGpuTimer(); + } } } diff --git a/stream_compaction/thrust.cu b/stream_compaction/thrust.cu index 1def45e7..ad462ec5 100644 --- a/stream_compaction/thrust.cu +++ b/stream_compaction/thrust.cu @@ -6,6 +6,7 @@ #include "common.h" #include "thrust.h" + namespace StreamCompaction { namespace Thrust { using StreamCompaction::Common::PerformanceTimer; @@ -18,11 +19,24 @@ namespace StreamCompaction { * Performs prefix-sum (aka scan) on idata, storing the result into odata. */ void scan(int n, int *odata, const int *idata) { + int* obuffer; + int* ibuffer; + cudaMalloc(&obuffer, n * sizeof(int)); + cudaMalloc(&ibuffer, n * sizeof(int)); + + cudaMemcpy(ibuffer, idata, n * sizeof(int), cudaMemcpyHostToDevice); + + thrust::device_ptr dev_in = thrust::device_pointer_cast(ibuffer); + thrust::device_ptr dev_out = thrust::device_pointer_cast(obuffer); + timer().startGpuTimer(); - // 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::exclusive_scan(dev_in, dev_in + n, dev_out); timer().endGpuTimer(); + + cudaMemcpy(odata, obuffer, n * sizeof(int), cudaMemcpyDeviceToHost); + + cudaFree(obuffer); + cudaFree(ibuffer); } } }