diff --git a/.vscode/launch.json b/.vscode/launch.json new file mode 100644 index 0000000..f687e19 --- /dev/null +++ b/.vscode/launch.json @@ -0,0 +1,14 @@ +{ + // Use IntelliSense to learn about possible attributes. + // Hover to view descriptions of existing attributes. + // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387 + "version": "0.2.0", + "configurations": [ + { + "name": "CUDA C++: Launch", + "type": "cuda-gdb", + "request": "launch", + "program": "${workspaceFolder}/build/bin/cis5650_stream_compaction_test" + } + ] +} \ No newline at end of file diff --git a/README.md b/README.md index 0e38ddb..b20503f 100644 --- a/README.md +++ b/README.md @@ -3,12 +3,169 @@ CUDA Stream Compaction **University of Pennsylvania, CIS 565: GPU Programming and Architecture, Project 2** -* (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) +* Xinran Tao + - [LinkedIn](https://www.linkedin.com/in/xinran-tao/), [Personal Website](https://www.xinrantao.com/), [GitHub](https://github.com/theBoilingPoint). +* Tested on: + - Ubuntu 22.04, i7-11700K @ 3.60GHz × 16, RAM 32GB, GeForce RTX 3080 Ti 12GB (Personal) -### (TODO: Your README) +# Base Credit +In order to count only the algorithm runtime, GPU timers are put inside the initialisation and memory copy functions. For example: +```cpp +cudaMalloc((void**)&dev_data, n * sizeof(int)); +cudaMemcpy(dev_data, data, n * sizeof(int), cudaMemcpyHostToDevice); -Include analysis, etc. (Remember, this is public, so don't put -anything here that you don't want to share with the world.) +timer().startGpuTimer(); +// The actual algorithm +timer().endGpuTimer(); +cudaFree(dev_data); +``` + +## Block Size Optimisation +A quick block size optimisation is done for the naive and efficient scan. The range of the block size tested is from 32 to 1024 (inclusive), with each step being the power of 2. The size of the input array is fixed at 2^20. The best block size for each algorithm is selected based on the average runtime between array size that is the power of 2 and non-power of 2 (i.e. the power of 2 size - 3). The result is as follows: +| Block Size | Naive | Work Efficient Slow | Work Efficient Fast | +| --------------| --------- | ------------------- | ------------------- | +| 32 | 0.059808 | 0.584002 | 0.602928 | +| 64 | 0.519168 | 0.665728 | 0.35776 | +| 128 | 0.442368 | 0.32856 | 0.322176 | +| 256 | 0.451072 | 0.262784 | 0.199088 | +| 512 | 0.44856 | 0.268496 | 0.255712 | +| 1024 | 0.488448 | 0.389296 | 0.278512 | + +Therefore, the block size chosen for naive scan is **32**, and **256** for both work efficient slow and fast scan. The block sizes are fixed for all subsequent experiments. + +## Analysis of the CPU, Naive, Work Efficient, and Thrust Scan +![](img/analysis/runtime_raw_data.png) + +![](img/analysis/Input%20Array%20Size%20(Power%20of%202)-Runtime%20Graph%20for%20CPU,%20Naive,%20Work%20Efficient,%20and%20Thrust%20Implementation.png) + +![](img/analysis/Input%20Array%20Size%20(Non%20Power%20of%202)-Runtime%20Graph%20for%20CPU,%20Naive,%20Work%20Efficient,%20and%20Thrust%20Implementation.png) + +The **CPU implementation**, both for power-of-two and non-power-of-two sizes, shows a predictable linear increase in runtime as the input size grows. This suggests that the CPU implementation is computation-bound, scaling consistently with the size of the input array. However, for very large inputs, especially beyond size 24, the non-power-of-two runtimes are slightly slower than power-of-two runtimes, which could indicate some inefficiencies in memory handling or alignment issues for irregular sizes. Still, the primary bottleneck here seems to be the inherent computational complexity of the algorithm, as it scales naturally with input size. + +For the **Naive implementation**, the performance bottlenecks are quite clear. The power-of-two version exhibits a steep increase in runtime as the input size grows, becoming particularly inefficient for larger sizes, as seen by the runtime at size 30 (~919 ms). This steep increase suggests that the Naive implementation is computation-bound and inefficiently parallelized, with large overhead in processing each element. Moreover, it might also suffer from poor memory access patterns, as its performance significantly worsens with increased input sizes, indicating that both computation and memory I/O are bottlenecks. The non-power-of-two Naive implementation performs better, which could be attributed to better handling of irregular data sizes, possibly benefiting from more efficient memory management. However, as input sizes increase, this advantage diminishes, and the overall computational inefficiencies of the Naive approach remain evident. + +The **Work Efficient Slow implementation** shows an improvement over the Naive implementation, particularly for larger input sizes. For power-of-two inputs, it performs moderately well, though its runtime increases more slowly than the Naive version, suggesting that it optimizes the computational workload better. However, the slowdown observed at larger sizes implies that memory I/O might still be a bottleneck, as efficient computation becomes less useful if memory access patterns aren't optimized. The non-power-of-two version of Work Efficient Slow performs slightly better, indicating better handling of irregular input sizes, but the performance gap narrows as input sizes grow, suggesting that at larger scales, both computational and memory inefficiencies become more pronounced. + +The **Work Efficient Fast implementation** offers significant performance gains over both the Naive and Slow versions, particularly for larger inputs. For power-of-two sizes, its scaling is far more controlled, and even for larger sizes, it remains competitive with the CPU. This suggests that the Work Efficient Fast implementation has effectively optimized both computation and memory access. It likely uses more efficient memory patterns, reducing the overhead caused by memory I/O, while also optimizing the computational workload, which allows it to remain fast even for larger inputs. The non-power-of-two version of Work Efficient Fast performs similarly, with only marginal degradation for irregular input sizes. This demonstrates that this implementation is relatively well-balanced and doesn't suffer from the same memory I/O bottlenecks that are apparent in the Naive and Slow approaches. + +Finally, the **Thrust implementation** consistently outperforms all other implementations for both power-of-two and non-power-of-two input sizes, particularly for small and medium-sized arrays. It maintains an almost constant runtime for smaller inputs, which indicates that the computational workload is highly optimized and that memory I/O is not a bottleneck, at least for smaller inputs. As input sizes grow, the runtime increases, but it remains the fastest or near the fastest implementation for most input sizes. This suggests that Thrust effectively handles both computation and memory I/O, optimizing for both types of input sizes. However, for the largest input sizes, the Work Efficient Fast implementation surpasses Thrust, indicating that while Thrust is optimized for general efficiency, very large data sets might still cause some memory bottlenecks. + +### Thrust's Implementation +From the Thrust implementation's performance behavior, particularly its near-constant runtime for smaller input sizes and gradual increase as the input size grows, we can infer that Thrust likely has highly optimized memory management and computational routines. The near-constant runtime at smaller sizes suggests minimal overhead for memory allocation and data transfers at this scale, implying that these stages are either negligible or well-pipelined. + +However, for larger input sizes, where the performance begins to degrade, we can speculate that some internal operations—such as repeated memory allocations or device-host memory copies—become more expensive. The increased data transfer costs between the host and device could be a bottleneck as input size increases. Moreover, the parallel computation itself may reach a point where it is bounded by the hardware's memory bandwidth rather than computational power, which Nsight might reveal as periods of memory transfer in the timeline. + +## Test Program Output +Here the input array size is 2^20, and 2^20 - 3 for the non-power of 2 size. +``` +**************** +** SCAN TESTS ** +**************** +==== Input ==== + [ 28 8 45 47 22 24 13 45 29 46 22 48 25 ... 35 0 ] +==== cpu scan, power-of-two ==== + elapsed time: 0.552226ms (std::chrono Measured) + [ 0 28 36 81 128 150 174 187 232 261 307 329 377 ... 25694240 25694275 ] +==== cpu scan, non-power-of-two ==== + elapsed time: 0.5049ms (std::chrono Measured) + [ 0 28 36 81 128 150 174 187 232 261 307 329 377 ... 25694195 25694198 ] + passed +==== naive scan, power-of-two ==== + elapsed time: 0.720896ms (CUDA Measured) + [ 0 28 36 81 128 150 174 187 232 261 307 329 377 ... 25694240 25694275 ] + passed +==== naive scan, non-power-of-two ==== + elapsed time: 0.765952ms (CUDA Measured) + [ 0 28 36 81 128 150 174 187 232 261 307 329 377 ... 0 0 ] + passed +==== work-efficient scan, power-of-two ==== +Using FAST work efficient scan + elapsed time: 0.413472ms (CUDA Measured) + [ 0 28 36 81 128 150 174 187 232 261 307 329 377 ... 25694240 25694275 ] + passed +==== work-efficient scan, non-power-of-two ==== +Using FAST work efficient scan + elapsed time: 0.18592ms (CUDA Measured) + [ 0 28 36 81 128 150 174 187 232 261 307 329 377 ... 25694195 25694198 ] + passed +==== thrust scan, power-of-two ==== + elapsed time: 0.82944ms (CUDA Measured) + [ 0 28 36 81 128 150 174 187 232 261 307 329 377 ... 25694240 25694275 ] + passed +==== thrust scan, non-power-of-two ==== + elapsed time: 0.099328ms (CUDA Measured) + [ 0 28 36 81 128 150 174 187 232 261 307 329 377 ... 25694195 25694198 ] + passed + +***************************** +** STREAM COMPACTION TESTS ** +***************************** +==== Input ==== + [ 1 3 2 1 1 2 3 1 2 3 2 1 1 ... 1 0 ] +==== cpu compact without scan, power-of-two ==== + elapsed time: 2.05752ms (std::chrono Measured) + [ 1 3 2 1 1 2 3 1 2 3 2 1 1 ... 3 1 ] + passed +==== cpu compact without scan, non-power-of-two ==== + elapsed time: 2.08289ms (std::chrono Measured) + [ 1 3 2 1 1 2 3 1 2 3 2 1 1 ... 3 3 ] + passed +==== cpu compact with scan ==== + elapsed time: 4.46571ms (std::chrono Measured) + [ 1 3 2 1 1 2 3 1 2 3 2 1 1 ... 3 1 ] + passed +==== work-efficient compact, power-of-two ==== + elapsed time: 0.3528ms (CUDA Measured) + [ 1 3 2 1 1 2 3 1 2 3 2 1 1 ... 3 1 ] + passed +==== work-efficient compact, non-power-of-two ==== + elapsed time: 0.233184ms (CUDA Measured) + [ 1 3 2 1 1 2 3 1 2 3 2 1 1 ... 3 3 ] + passed +``` + +# Extra Credit + +## GPU Slower than CPU? +As explained in week 4's lecture, the original version of the work efficient scan is launching threads that are not doing any work, as shown below: +![](img/analysis/week4_slide.png) +Only the threads coloured in yellow are doing useful work (in this case the sum operation). And as stated in the slide, the number of threads required is halved at each pass. + +Therefore, to optimise the original work efficient scan, we need to launch a smaller kernel size at each iteration, and read the numbers for adding by offsetting the thread index. The original implementation is the `upsweep_slow` function and the improved version is the `upsweep` function inside `efficient.cu`. To toggle between the two, change the `USE_FAST_UPSWEEP` macro in on top of the file. + +The runtime comparison between the original, the improved, and the cpu version can be found in Section _Analysis of the CPU, Naive, Work Efficient, and Thrust Scan_. It's clear to see that the advantage of using the improved version of the work-efficient scan starts showing when the input array size is large (i.e. above 2^18). + +For small input sizes, the overheads associated with GPU execution—kernel launch, memory transfer, underutilization, and synchronization—can outweigh the benefits of parallel computation. The CPU, in contrast, can handle small workloads more efficiently due to its ability to process data directly in memory with minimal overhead. Therefore, while your Work Efficient Fast implementation is optimized for larger datasets, the GPU's inherent overhead for smaller arrays results in slower performance compared to the CPU. + +## Radix Sort +The `sort.h` and `sort.cu` files are added to the `stream_compaction` directory. It's a simple version that uses the global memory and the fast version of the work-efficient scan. + +## Test Program Output +For the easy case, the input array size is 8. For the other two tests, the input array size is 2^20. The range of values for the normal range test is [0, 50) and [0, 2147483648) for the large range test. The implementation from thrust is used as the reference for correctness. +``` +***************************** +** RADIX SORT TESTS ** +***************************** +==== Input ==== + [ 4 7 2 6 3 5 1 0 ] +==== radix sort thrust, easy case ==== + [ 0 1 2 3 4 5 6 7 ] +==== radix sort with work-efficient scan, easy case ==== + [ 0 1 2 3 4 5 6 7 ] + passed +==== Input ==== + [ 0 0 1 1 0 3 0 1 0 2 1 2 0 ... 3 0 ] +==== radix sort thrust ==== + [ 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 3 3 ] +==== radix sort with work-efficient scan ==== + [ 0 0 0 0 0 0 0 0 0 0 0 0 0 ... 3 3 ] + passed +==== Input ==== + [ 2030340076 367405552 870298961 958594477 1919490608 650582255 502155192 1038962173 812028104 979602990 45151013 1940063754 445481300 ... 1491375275 1238713711 ] +==== radix sort thrust, large range ==== + [ 128 9077 11470 13673 16336 20082 23386 25138 25834 27486 28982 30125 31576 ... 2147481763 2147483586 ] +==== radix sort with work-efficient scan, large range ==== + [ 128 9077 11470 13673 16336 20082 23386 25138 25834 27486 28982 30125 31576 ... 2147481763 2147483586 ] + passed +``` diff --git a/img/analysis/Input Array Size (Non Power of 2)-Runtime Graph for CPU, Naive, Work Efficient, and Thrust Implementation.png b/img/analysis/Input Array Size (Non Power of 2)-Runtime Graph for CPU, Naive, Work Efficient, and Thrust Implementation.png new file mode 100644 index 0000000..f8ddbdf Binary files /dev/null and b/img/analysis/Input Array Size (Non Power of 2)-Runtime Graph for CPU, Naive, Work Efficient, and Thrust Implementation.png differ diff --git a/img/analysis/Input Array Size (Power of 2)-Runtime Graph for CPU, Naive, Work Efficient, and Thrust Implementation.png b/img/analysis/Input Array Size (Power of 2)-Runtime Graph for CPU, Naive, Work Efficient, and Thrust Implementation.png new file mode 100644 index 0000000..3f63f39 Binary files /dev/null and b/img/analysis/Input Array Size (Power of 2)-Runtime Graph for CPU, Naive, Work Efficient, and Thrust Implementation.png differ diff --git a/img/analysis/runtime_raw_data.png b/img/analysis/runtime_raw_data.png new file mode 100644 index 0000000..0869b18 Binary files /dev/null and b/img/analysis/runtime_raw_data.png differ diff --git a/img/analysis/week4_slide.png b/img/analysis/week4_slide.png new file mode 100644 index 0000000..effa729 Binary files /dev/null and b/img/analysis/week4_slide.png differ diff --git a/src/main.cpp b/src/main.cpp index 896ac2b..b182eef 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -11,9 +11,10 @@ #include #include #include +#include #include "testing_helpers.hpp" -const int SIZE = 1 << 8; // feel free to change the size of array +const int SIZE = 1 << 20; // feel free to change the size of array, default is 8 const int NPOT = SIZE - 3; // Non-Power-Of-Two int *a = new int[SIZE]; int *b = new int[SIZE]; @@ -29,6 +30,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; + printDesc("Input"); printArray(SIZE, a, true); // initialize b using StreamCompaction::CPU::scan you implement @@ -51,48 +53,48 @@ int main(int argc, char* argv[]) { printDesc("naive scan, power-of-two"); StreamCompaction::Naive::scan(SIZE, c, a); printElapsedTime(StreamCompaction::Naive::timer().getGpuElapsedTimeForPreviousOperation(), "(CUDA Measured)"); - //printArray(SIZE, c, true); + printArray(SIZE, c, true); printCmpResult(SIZE, b, c); - /* For bug-finding only: Array of 1s to help find bugs in stream compaction or scan - onesArray(SIZE, c); - printDesc("1s array for finding bugs"); - StreamCompaction::Naive::scan(SIZE, c, a); - printArray(SIZE, c, true); */ + /* For bug-finding only: Array of 1s to help find bugs in stream compaction or scan */ + // onesArray(SIZE, c); + // printDesc("1s array for finding bugs"); + // StreamCompaction::Naive::scan(SIZE, c, a); + // printArray(SIZE, c, true); zeroArray(SIZE, c); printDesc("naive scan, non-power-of-two"); StreamCompaction::Naive::scan(NPOT, c, a); printElapsedTime(StreamCompaction::Naive::timer().getGpuElapsedTimeForPreviousOperation(), "(CUDA Measured)"); - //printArray(SIZE, c, true); + printArray(SIZE, c, true); printCmpResult(NPOT, b, c); zeroArray(SIZE, c); printDesc("work-efficient scan, power-of-two"); StreamCompaction::Efficient::scan(SIZE, c, a); printElapsedTime(StreamCompaction::Efficient::timer().getGpuElapsedTimeForPreviousOperation(), "(CUDA Measured)"); - //printArray(SIZE, c, true); + printArray(SIZE, c, true); printCmpResult(SIZE, b, c); zeroArray(SIZE, c); printDesc("work-efficient scan, non-power-of-two"); StreamCompaction::Efficient::scan(NPOT, c, a); printElapsedTime(StreamCompaction::Efficient::timer().getGpuElapsedTimeForPreviousOperation(), "(CUDA Measured)"); - //printArray(NPOT, c, true); + printArray(NPOT, c, true); printCmpResult(NPOT, b, 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,6 +106,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; + printDesc("Input"); printArray(SIZE, a, true); int count, expectedCount, expectedNPOT; @@ -137,18 +140,78 @@ int main(int argc, char* argv[]) { 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); printDesc("work-efficient compact, non-power-of-two"); count = StreamCompaction::Efficient::compact(NPOT, c, a); printElapsedTime(StreamCompaction::Efficient::timer().getGpuElapsedTimeForPreviousOperation(), "(CUDA Measured)"); - //printArray(count, c, true); + printArray(count, c, true); printCmpLenResult(count, expectedNPOT, b, c); + printf("\n"); + printf("*****************************\n"); + printf("** RADIX SORT TESTS **\n"); + printf("*****************************\n"); + + int *d = new int[8]; + + d[0] = 4; // 100 + d[1] = 7; // 111 + d[2] = 2; // 010 + d[3] = 6; // 110 + d[4] = 3; // 011 + d[5] = 5; // 101 + d[6] = 1; // 001 + d[7] = 0; // 000 + + printDesc("Input"); + printArray(8, d, true); + + zeroArray(8, b); + printDesc("radix sort thrust, easy case"); + StreamCompaction::Sort::radix_sort_thrust(8, b, d); + printArray(8, b, false); + + zeroArray(8, c); + printDesc("radix sort with work-efficient scan, easy case"); + StreamCompaction::Sort::radix_sort(8, c, d); + printArray(8, c, false); + printCmpResult(8, b, c); + + printDesc("Input"); + printArray(SIZE, a, true); + + zeroArray(SIZE, b); + printDesc("radix sort thrust"); + StreamCompaction::Sort::radix_sort_thrust(SIZE, b, a); + printArray(SIZE, b, true); + + zeroArray(SIZE, c); + printDesc("radix sort with work-efficient scan"); + StreamCompaction::Sort::radix_sort(SIZE, c, a); + printArray(SIZE, c, true); + printCmpResult(SIZE, b, c); + + genArray(SIZE , a, 2147483648); // 2,147,483,647 is the largest value for int32 + printDesc("Input"); + printArray(SIZE, a, true); + + zeroArray(SIZE, b); + printDesc("radix sort thrust, large range"); + StreamCompaction::Sort::radix_sort_thrust(SIZE, b, a); + printArray(SIZE, b, true); + + zeroArray(SIZE, c); + printDesc("radix sort with work-efficient scan, large range"); + StreamCompaction::Sort::radix_sort(SIZE, c, a); + printArray(SIZE, c, true); + printCmpResult(SIZE, b, c); + system("pause"); // stop Win32 console from closing on exit delete[] a; delete[] b; delete[] c; + delete[] d; } diff --git a/stream_compaction/CMakeLists.txt b/stream_compaction/CMakeLists.txt index 19511ca..bbcb5e3 100644 --- a/stream_compaction/CMakeLists.txt +++ b/stream_compaction/CMakeLists.txt @@ -4,6 +4,7 @@ set(headers "naive.h" "efficient.h" "thrust.h" + "sort.h" ) set(sources @@ -12,6 +13,7 @@ set(sources "naive.cu" "efficient.cu" "thrust.cu" + "sort.cu" ) list(SORT headers) diff --git a/stream_compaction/common.cu b/stream_compaction/common.cu index 2ed6d63..cd9ba00 100644 --- a/stream_compaction/common.cu +++ b/stream_compaction/common.cu @@ -23,7 +23,13 @@ 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 i = threadIdx.x + (blockIdx.x * blockDim.x); + + if (i >= n) { + return; + } + + bools[i] = idata[i] != 0 ? 1 : 0; } /** @@ -32,8 +38,15 @@ namespace StreamCompaction { */ __global__ void kernScatter(int n, int *odata, const int *idata, const int *bools, const int *indices) { - // TODO - } + int i = threadIdx.x + (blockIdx.x * blockDim.x); + if (i >= n) { + return; + } + + if (bools[i] == 1) { + odata[indices[i]] = idata[i]; + } + } } } diff --git a/stream_compaction/cpu.cu b/stream_compaction/cpu.cu index 719fa11..f390053 100644 --- a/stream_compaction/cpu.cu +++ b/stream_compaction/cpu.cu @@ -19,10 +19,23 @@ 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(); } + // Use this in compactWithScan so that we won't have "CPU timer already started" bug + void scan_without_timer(int n, int *odata, const int *idata) { + odata[0] = 0; + for (int i = 1; i < n; i++) { + odata[i] = odata[i - 1] + idata[i - 1]; + } + } + /** * CPU stream compaction without using the scan function. * @@ -30,9 +43,17 @@ namespace StreamCompaction { */ int compactWithoutScan(int n, int *odata, const int *idata) { timer().startCpuTimer(); - // TODO + + int cur_out_idx = 0; + for (int i = 0; i < n; i++) { + if (idata[i] != 0) { + odata[cur_out_idx] = idata[i]; + cur_out_idx++; + } + } + timer().endCpuTimer(); - return -1; + return cur_out_idx; } /** @@ -42,9 +63,28 @@ namespace StreamCompaction { */ int compactWithScan(int n, int *odata, const int *idata) { timer().startCpuTimer(); - // TODO + + // Step 1: Map idata to a 0/1 array + int *bools = new int[n]; // Array to hold the 0s and 1s + for (int i = 0; i < n; i++) { + bools[i] = (idata[i] != 0) ? 1 : 0; + } + + // Step 2: Perform the scan using the existing scan function + int *scanResult = new int[n]; // Array to hold the scan result + scan_without_timer(n, scanResult, bools); // Perform the scan on the bools array + + // Step 3: Scatter non-zero elements into odata + int count = 0; + for (int i = 0; i < n; i++) { + if (bools[i] == 1) { + odata[scanResult[i]] = idata[i]; + count++; + } + } + timer().endCpuTimer(); - return -1; + return count; } } } diff --git a/stream_compaction/efficient.cu b/stream_compaction/efficient.cu index 2db346e..f678b33 100644 --- a/stream_compaction/efficient.cu +++ b/stream_compaction/efficient.cu @@ -3,24 +3,196 @@ #include "common.h" #include "efficient.h" +#define blockSize_slow 256 +#define blockSize_fast 256 +#define USE_FAST_UPSWEEP 1 + namespace StreamCompaction { namespace Efficient { using StreamCompaction::Common::PerformanceTimer; + + // Initialise device variables to use in compact() + int *dev_bools, *dev_idata, *dev_odata, *dev_scanResult; + PerformanceTimer& timer() { static PerformanceTimer timer; return timer; } + // n & (n - 1): This expression clears the lowest set bit of n. + // If n is a power of 2, it has exactly one bit set, and subtracting 1 from it flips + // all the bits after the most significant bit. Performing n & (n - 1) results in + // 0 if and only if n is a power of 2. + // n > 0: This ensures that n is positive, since negative numbers and zero are not powers of 2. + __device__ __host__ bool isPowerOf2(int n) { + return (n > 0) && ((n & (n - 1)) == 0); + } + + __device__ int mod(int a, int b) { + return a - (b * (a / b)); + } + + __global__ void upsweep_slow(int n, int d, int *data) { + int k = threadIdx.x + (blockIdx.x * blockDim.x); + + if (k >= n) { + return; + } + + // is that shifting bits to the left by n positions is equivalent to multiplying the number by 2^n. + // faster than calling pow(2, n) + int two_pow_d_plus_1 = 1 << (d + 1); + + if (mod(k, two_pow_d_plus_1) != 0) { + return; + } + + int two_pow_d = 1 << d; + + data[k + two_pow_d_plus_1 - 1] += data[k + two_pow_d - 1]; + } + + __global__ void upsweep(int n, int d, int *data) { + int k = threadIdx.x + (blockIdx.x * blockDim.x); + + if (k >= n) { + return; + } + + data[((k + 1) << d) - 1] += data[k * (1 << d) + (1 << (d - 1)) - 1]; + } + + __global__ void init_downsweep(int n, int *odata) { + int i = threadIdx.x + (blockIdx.x * blockDim.x); + + if (i == n - 1) { + odata[i] = 0; + } + } + + __global__ void downsweep(int n, int d, int *data) { + int k = threadIdx.x + (blockIdx.x * blockDim.x); + + if (k >= n) { + return; + } + + // is that shifting bits to the left by n positions is equivalent to multiplying the number by 2^n. + // faster than calling pow(2, n) + int two_pow_d_plus_1 = 1 << (d + 1); + + if (mod(k, two_pow_d_plus_1) != 0) { + return; + } + + int two_pow_d = 1 << d; + int t = data[k + two_pow_d - 1]; + data[k + two_pow_d - 1] = data[k + two_pow_d_plus_1 - 1]; + data[k + two_pow_d_plus_1 - 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 - timer().endGpuTimer(); + int *dev_odata_local; + // Your intermediate array sizes will need to be rounded to the next power of two. + int rounded_n = isPowerOf2(n) ? n : 1 << ilog2ceil(n); + + cudaMalloc((void**)&dev_odata_local, rounded_n * sizeof(int)); + + // Copy idata to dev_odata_local first + // Although we might initialise dev_odata_local with more than n elements, idata only contains n elements + cudaMemcpy(dev_odata_local, idata, n * sizeof(int), cudaMemcpyHostToDevice); + + // Set up the grid and block sizes + + int upper_bound = ilog2ceil(rounded_n); + int upper_bound_minus_1 = upper_bound - 1; + + if (USE_FAST_UPSWEEP) { + printf("Using FAST work efficient scan\n"); + + dim3 fullBlocksPerGrid((rounded_n + blockSize_fast - 1) / blockSize_fast); + int gridSize = (rounded_n / 2 + blockSize_fast - 1) / blockSize_fast; + + timer().startGpuTimer(); + + for (int d = 1; d <= upper_bound; d++) { + gridSize = ((rounded_n >> d) + blockSize_fast - 1) / blockSize_fast; + upsweep << > >(rounded_n >> d, d, dev_odata_local); + } + + // downsweep + init_downsweep << > > (rounded_n, dev_odata_local); + + for (int d = upper_bound_minus_1; d >= 0; d--) { + downsweep << > > (rounded_n, d, dev_odata_local); + } + + timer().endGpuTimer(); + } + else { + printf("Using SLOW work efficient scan\n"); + + dim3 fullBlocksPerGrid((rounded_n + blockSize_slow - 1) / blockSize_slow); + + timer().startGpuTimer(); + + // upsweep + for (int d = 0; d < upper_bound; d++) { + upsweep_slow << > > (rounded_n, d, dev_odata_local); + } + + // downsweep + init_downsweep << > > (rounded_n, dev_odata_local); + + for (int d = upper_bound_minus_1; d >= 0; d--) { + downsweep << > > (rounded_n, d, dev_odata_local); + } + + timer().endGpuTimer(); + } + + // Copy the result back to the host + // Note that odata is only supposed to have n elements + cudaMemcpy(odata, dev_odata_local, n * sizeof(int), cudaMemcpyDeviceToHost); + + cudaFree(dev_odata_local); } + /************************************************************************************ + * Define another scan function so that it can be called from compact() + * Here dev_bools has already been initialised in compact() + ************************************************************************************/ + void scan_without_timer(int n, int *odata, const int *idata) { + // Initialise dev_scanResult. dev_bools only has n elements. + cudaMemcpy(odata, idata, n * sizeof(int), cudaMemcpyDeviceToDevice); + + // Your intermediate array sizes will need to be rounded to the next power of two. + int rounded_n = isPowerOf2(n) ? n : 1 << ilog2ceil(n); + int upper_bound = ilog2ceil(rounded_n); + int upper_bound_minus_1 = upper_bound - 1; + + // Set up the grid and block sizes + dim3 fullBlocksPerGrid((rounded_n + blockSize_fast - 1) / blockSize_fast); + int gridSize = (rounded_n / 2 + blockSize_fast - 1) / blockSize_fast; + + // upsweep + for (int d = 1; d <= upper_bound; d++) { + gridSize = ((rounded_n >> d) + blockSize_fast - 1) / blockSize_fast; + upsweep << > >(rounded_n >> d, d, odata); + } + + // downsweep + init_downsweep << > > (rounded_n, odata); + + for (int d = upper_bound_minus_1; d >= 0; d--) { + downsweep << > > (rounded_n, d, odata); + } + } + /** * Performs stream compaction on idata, storing the result into odata. * All zeroes are discarded. @@ -31,10 +203,46 @@ namespace StreamCompaction { * @returns The number of elements remaining after compaction. */ int compact(int n, int *odata, const int *idata) { + // Initialise host variables for returning + int *scanResult = new int[n]; + + // Allocate memory on the device + cudaMalloc((void**)&dev_idata, n * sizeof(int)); + cudaMalloc((void**)&dev_bools, n * sizeof(int)); + cudaMalloc((void**)&dev_scanResult, n * sizeof(int)); + cudaMalloc((void**)&dev_odata, n * sizeof(int)); + + cudaMemcpy(dev_idata, idata, n * sizeof(int), cudaMemcpyHostToDevice); + + // Set up the grid and block sizes + dim3 fullBlocksPerGrid((n + blockSize_fast - 1) / blockSize_fast); + timer().startGpuTimer(); - // TODO + + // Map idata to a 0/1 array + Common::kernMapToBoolean << > > (n, dev_bools, dev_idata); + + // Scan the boolean array + scan_without_timer(n, dev_scanResult, dev_bools); // n will be rounded in the scan function + + // Perform scatter + Common::kernScatter << > > (n, dev_odata, dev_idata, dev_bools, dev_scanResult); + timer().endGpuTimer(); - return -1; + + cudaMemcpy(scanResult, dev_scanResult, n * sizeof(int), cudaMemcpyDeviceToHost); + cudaMemcpy(odata, dev_odata, n * sizeof(int), cudaMemcpyDeviceToHost); + int count = n == 0 ? 0 : scanResult[n - 1] + (idata[n - 1] != 0 ? 1 : 0); + + // Clean up device memory + cudaFree(dev_idata); + cudaFree(dev_bools); + cudaFree(dev_scanResult); + cudaFree(dev_odata); + + delete[] scanResult; + + return count; } } } diff --git a/stream_compaction/efficient.h b/stream_compaction/efficient.h index 803cb4f..ce435f6 100644 --- a/stream_compaction/efficient.h +++ b/stream_compaction/efficient.h @@ -8,6 +8,8 @@ namespace StreamCompaction { void scan(int n, int *odata, const int *idata); + void scan_without_timer(int n, int *odata, const int *idata); + int compact(int n, int *odata, const int *idata); } } diff --git a/stream_compaction/naive.cu b/stream_compaction/naive.cu index 4308876..86be613 100644 --- a/stream_compaction/naive.cu +++ b/stream_compaction/naive.cu @@ -3,6 +3,8 @@ #include "common.h" #include "naive.h" +#define blockSize 32 + namespace StreamCompaction { namespace Naive { using StreamCompaction::Common::PerformanceTimer; @@ -11,15 +13,62 @@ namespace StreamCompaction { static PerformanceTimer timer; return timer; } - // TODO: __global__ + + __global__ void init_buff(int n, int *odata, const int *idata) { + int i = threadIdx.x + (blockIdx.x * blockDim.x); + + if (i == 0) { + odata[i] = 0; + return; + } + + if (i < n) { + odata[i] = idata[i - 1]; + } + } + + __global__ void dev_scan(int pow_d, int *odata, const int *idata) { + int k = threadIdx.x + (blockIdx.x * blockDim.x); + + if (k >= pow_d) { + odata[k] = idata[k - pow_d] + idata[k]; + } + } /** * Performs prefix-sum (aka scan) on idata, storing the result into odata. */ void scan(int n, int *odata, const int *idata) { + int *dev_buff1, *dev_buff2; + + // Allocate memory on the device + cudaMalloc((void**)&dev_buff1, n * sizeof(int)); + cudaMalloc((void**)&dev_buff2, n * sizeof(int)); + + // Copy data to the buffer for initialisation + cudaMemcpy(dev_buff1, idata, n * sizeof(int), cudaMemcpyHostToDevice); + + // Set up the grid and block sizes + dim3 fullBlocksPerGrid((n + blockSize - 1) / blockSize); + + // Because we want the first element to be 0 and the rest filled with the original data, we need to fill the values explicitly + init_buff << > > (n, dev_buff2, dev_buff1); + cudaMemcpy(dev_buff1, dev_buff2, n * sizeof(int), cudaMemcpyDeviceToDevice); + timer().startGpuTimer(); - // TODO + + for (int d = 1; d <= ilog2ceil(n); d++) { + int pow_d = 1 << (d - 1); + dev_scan << > > (pow_d, dev_buff2, dev_buff1); + cudaMemcpy(dev_buff1, dev_buff2, n * sizeof(int), cudaMemcpyDeviceToDevice); + } + timer().endGpuTimer(); + + cudaMemcpy(odata, dev_buff1, n * sizeof(int), cudaMemcpyDeviceToHost); + + cudaFree(dev_buff1); + cudaFree(dev_buff2); } } } diff --git a/stream_compaction/sort.cu b/stream_compaction/sort.cu new file mode 100644 index 0000000..16db8f9 --- /dev/null +++ b/stream_compaction/sort.cu @@ -0,0 +1,132 @@ +#include +#include +#include +#include "thrust/sort.h" +#include "common.h" +#include "efficient.h" +#include "sort.h" + +#define blockSize 128 + +namespace StreamCompaction { + namespace Sort { + using StreamCompaction::Common::PerformanceTimer; + int *dev_idata, *dev_bit_extract_neg, *dev_scan_results, *dev_offset_idx, *dev_final_idx, *dev_odata; + + PerformanceTimer& timer() + { + static PerformanceTimer timer; + return timer; + } + + __global__ void compute_not_of_bits(int n, int bit_idx, int *odata, const int *idata) { + int idx = threadIdx.x + (blockIdx.x * blockDim.x); + + if (idx >= n) { + return; + } + + // Extract the bit at the given index and negate it + odata[idx] = ((idata[idx] >> bit_idx) & 1) ^ 1; + } + + __global__ void compute_offset_indices(int n, int totalFalses, int *odata, const int *idata) { + int idx = threadIdx.x + (blockIdx.x * blockDim.x); + + if (idx >= n) { + return; + } + + odata[idx] = idx - idata[idx] + totalFalses; + } + + __global__ void compute_final_indices(int n, int *odata, const int *b_arr, const int *t_arr, const int *f_arr) { + int idx = threadIdx.x + (blockIdx.x * blockDim.x); + + if (idx >= n) { + return; + } + + odata[idx] = !b_arr[idx] ? t_arr[idx] : f_arr[idx]; + } + + __global__ void scatter(int n, int *odata, const int *idata, const int *indices) { + int idx = threadIdx.x + (blockIdx.x * blockDim.x); + + if (idx >= n) { + return; + } + + odata[indices[idx]] = idata[idx]; + } + + int get_total_falses(int n, const int* dev_scan_results, const int* dev_bit_extract_neg) { + if (n == 0) { + return 0; + } + + // Allocate host memory for just one element (the last one) + int h_last_scan_result; + int h_last_bit_extract_neg; + + // Copy only the last element of dev_scan_results and dev_bit_extract_neg + cudaMemcpy(&h_last_scan_result, &dev_scan_results[n - 1], sizeof(int), cudaMemcpyDeviceToHost); + cudaMemcpy(&h_last_bit_extract_neg, &dev_bit_extract_neg[n - 1], sizeof(int), cudaMemcpyDeviceToHost); + + // Compute totalFalses using only the last elements + int totalFalses = h_last_scan_result + (h_last_bit_extract_neg == 1 ? 1 : 0); + + return totalFalses; + } + + void radix_sort(int n, int *odata, const int *idata) { + cudaMalloc((void**)&dev_idata, n * sizeof(int)); + cudaMalloc((void**)&dev_bit_extract_neg, n * sizeof(int)); + cudaMalloc((void**)&dev_scan_results, n * sizeof(int)); + cudaMalloc((void**)&dev_offset_idx, n * sizeof(int)); + cudaMalloc((void**)&dev_final_idx, n * sizeof(int)); + cudaMalloc((void**)&dev_odata, n * sizeof(int)); + + cudaMemcpy(dev_idata, idata, n * sizeof(int), cudaMemcpyHostToDevice); + + // Set up the grid and block sizes + dim3 fullBlocksPerGrid((n + blockSize - 1) / blockSize); + + // We are assuming 32 bits for each integer + // Lopping from 0 will get us the least significant bit first + for (int i = 0; i < 32; i++) { + // Split function + compute_not_of_bits << > > (n, i, dev_bit_extract_neg, dev_idata); + StreamCompaction::Efficient::scan_without_timer(n, dev_scan_results, dev_bit_extract_neg); + int totalFalses = get_total_falses(n, dev_scan_results, dev_bit_extract_neg); + compute_offset_indices << > > (n, totalFalses, dev_offset_idx, dev_scan_results); + compute_final_indices << > > (n, dev_final_idx, dev_bit_extract_neg, dev_offset_idx, dev_scan_results); + scatter << > > (n, dev_odata, dev_idata, dev_final_idx); + cudaMemcpy(dev_idata, dev_odata, n * sizeof(int), cudaMemcpyDeviceToDevice); + } + + cudaMemcpy(odata, dev_odata, n * sizeof(int), cudaMemcpyDeviceToHost); + + cudaFree(dev_idata); + cudaFree(dev_bit_extract_neg); + cudaFree(dev_scan_results); + cudaFree(dev_offset_idx); + cudaFree(dev_final_idx); + cudaFree(dev_odata); + } + + /** + * Performs radix sort on idata, storing the result into odata. + * This is the reference to the implementation, + * as thrust uses radix sort: https://stackoverflow.com/questions/45267748/which-sorting-algorithm-used-in-thrustsort + */ + void radix_sort_thrust(int n, int *odata, const int *idata) { + thrust::device_vector d_idata(idata, idata + n); // Copy data to device + thrust::device_vector d_odata(n); // Allocate space for output on device + + thrust::sort(d_idata.begin(), d_idata.end()); + + thrust::copy(d_idata.begin(), d_idata.end(), odata); + } + } +} \ No newline at end of file diff --git a/stream_compaction/sort.h b/stream_compaction/sort.h new file mode 100644 index 0000000..35d978b --- /dev/null +++ b/stream_compaction/sort.h @@ -0,0 +1,13 @@ +#pragma once + +#include "common.h" + +namespace StreamCompaction { + namespace Sort { + StreamCompaction::Common::PerformanceTimer& timer(); + + void radix_sort_thrust(int n, int *odata, const int *idata); + + void radix_sort(int n, int *odata, const int *idata); + } +} \ No newline at end of file diff --git a/stream_compaction/thrust.cu b/stream_compaction/thrust.cu index 1def45e..e15f378 100644 --- a/stream_compaction/thrust.cu +++ b/stream_compaction/thrust.cu @@ -14,15 +14,24 @@ namespace StreamCompaction { static PerformanceTimer timer; return timer; } + /** * Performs prefix-sum (aka scan) on idata, storing the result into odata. */ void scan(int n, int *odata, const int *idata) { + // Allocate device memory and copy input data to device + thrust::device_vector d_idata(idata, idata + n); // Copy data to device + thrust::device_vector d_odata(n); // Allocate space for output on device + 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()); + + // Perform exclusive scan on the device + thrust::exclusive_scan(d_idata.begin(), d_idata.end(), d_odata.begin()); + timer().endGpuTimer(); + + // Copy the result back to the host in one step + thrust::copy(d_odata.begin(), d_odata.end(), odata); } } }