Project 2: Yu-Chia Shen

README.md

@@ -3,12 +3,178 @@ 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)
+* Yu-Chia Shen
+  * [LinkedIn](https://www.linkedin.com/in/ycshen0831/)
+* Tested on: Windows 10, i5-11400F @ 4.3GHz 16GB, GTX 3060 12GB (personal)
 
-### (TODO: Your README)
+## Overview
 
-Include analysis, etc. (Remember, this is public, so don't put
-anything here that you don't want to share with the world.)
+This project includes 3 implementations of GPU parallel algorithms, **Scan**, **Stream Compaction**, and **Radix Sort**. Each implementations has several versions for performance comparison.
 
+### Scan  
+Performs an All-Prefix-Sum to an array.
+* **CPU Scan**: Serialized computation over the array. 
+* **GPU Naive Scan**: Parallelized exclusive sum.
+* **GPU Efficient Scan**: Opimized parallelized exclusive sum. Use a balanced binary tree to perform two phase
+  * Up-Sweep (Parallel Reduction)
+  * Down-Sweep - Merge partial sums to build scan in place.
+* **Thrust Implementaion**: A C++ template library for CUDA based on the Standard Template Library (STL).
+
+### Stream Compaction 
+Remove all the 0s from an array. 
+* **CPU Compaction**: Serialized computation over the array. 
+* **CPU Compaction with Scan**: Serialized computation using CPU scan. 
+* **GPU Stream Compaction**: Parallelized computation with 3 steps
+  * Map - Identify every 0s in the array
+  * Scan - Claculate indices for non-0 elements.
+  * Scatter - Write non-0 elements to correct positions.
+
+
+### Radix Sort (Extra Credit)
+Sort the array using partitions that based on one bit.
+* **Thrust Implementation**: A C++ template library that is used as comparison.
+* **GPU Radix Sort**: Parallelized computation with 3 steps for every bits. 
+  * Bit Mapping - From least significant bit to most significant bit, identify a bit mask for current bit.
+  * Scan - Perform scan over the bit mask.
+  * Split - Partition based on current bit.
+
+## Radix Sort (Extra Credit)
+### Usage
+`StreamCompaction::Efficient::radixSort(Size, output, input)`
+
+### Result
+```
+**********************
+** RADIX SORT TESTS **
+**********************
+    [  45  23   3  11  27  32  15  29   1  47  29   3  22 ...  11   0 ]
+==== Thrust sort, power-of-two ====
+   elapsed time: 0.00176ms    (CUDA Measured)
+    [   0   0   0   0   0   0   0   0   0   0   0   0   0 ...  49  49 ]
+==== Radix Sort, power-of-two ====
+   elapsed time: 74.6997ms    (CUDA Measured)
+    [   0   0   0   0   0   0   0   0   0   0   0   0   0 ...  49  49 ]
+    passed
+==== Thrust sort, non-power-of-two ====
+   elapsed time: 0.00192ms    (CUDA Measured)
+    [   0   0   0   0   0   0   0   0   0   0   0   0   0 ...  49  49 ]
+==== Radix Sort, non-power-of-two ====
+   elapsed time: 42.6619ms    (CUDA Measured)
+    [   0   0   0   0   0   0   0   0   0   0   0   0   0 ...  49  49 ]
+    passed
+```
+
+## Performance Analysis
+
+### Execution Time for Scan (Lower is Better)
+![](./img/ScanTime.png)
+_Figure 1: Scan Performance Comparison_
+
+As Figure 1 shows, we have the following conclusions:
+* __CPU Scan__: The worst version of scan. Serialized travel every elements will have a long runtime.
+* __GPU Naive Scan__: This version is slightly better than CPU version, but still has a long runtime.
+* __GPU Efficient Scan__: By using a balance binary tree for optimization, this version reduces almost half of the runtime. However, there is still performance bottlenecks. Due to the frequent access to the global memory, the runtime still increases significantly when the size of array is large. 
+
+  A solution is to use shared memory to replace the global memory. This will greatly improve the performance.
+* __Thrust Implementation__:
+![](./img/nsight.png)
+_Figure 2: Thrust Analysis with NSight_
+
+  As Figure 2 shows, Thrust has the following improvements:
+  * __Memory Allocation__: In the function `DeviceScanKernel`, there are 7696 bytes of Static Shared Memory. Using the shared memory can greatly reduce the memory access time.
+
+  * __Memory Copy Method__: The memory copy method used in Thrust is `cudaMemcpyAsync`. Unlike `cudaMemcpy` will block the host thread, `cudaMemcpyAsync` is non-blocking on the host. Therefore, host can transfer data concurrently, and thus is faster than `cudaMemcpy`. 
+
+### Execution Time for Stream Compaction (Lower is Better)
+![](./img/CompTime.png)
+_Figure 3: Stream Compaction Performance Comparison_
+
+Figure 3 shows the execution time of the 3 versions of Stream Compaction. We can see that the CUDA version still have the best performance. This is because it performs parallelized computation over the array.
+
+Also, we can observe that CPU without scan is better than the version with scan. This is because the stream compaction with scan will iterate the array three times, and thus will have a greater execution time when the array is large.
+
+## Sample Output
+
+```
+Using SIZE: 16 M 
+****************
+** SCAN TESTS **
+****************
+    [  38  44  42   2  16  46  40  33  46  25  24  15  12 ...  39   0 ]
+==== cpu scan, power-of-two ====
+   elapsed time: 27.3899ms    (std::chrono Measured)
+    [   0  38  82 124 126 142 188 228 261 307 332 356 371 ... 410934036 410934075 ]
+==== cpu scan, non-power-of-two ====
+   elapsed time: 27.4172ms    (std::chrono Measured)
+    [   0  38  82 124 126 142 188 228 261 307 332 356 371 ... 410933976 410933979 ]
+    passed
+==== naive scan, power-of-two ====
+   elapsed time: 11.4565ms    (CUDA Measured)
+    [   0  38  82 124 126 142 188 228 261 307 332 356 371 ... 410934036 410934075 ]
+    passed
+==== naive scan, non-power-of-two ====
+   elapsed time: 11.6449ms    (CUDA Measured)
+    [   0  38  82 124 126 142 188 228 261 307 332 356 371 ... 410933976 410933979 ]
+    passed
+==== work-efficient scan, power-of-two ====
+   elapsed time: 4.9359ms    (CUDA Measured)
+    [   0  38  82 124 126 142 188 228 261 307 332 356 371 ... 410934036 410934075 ]
+    passed
+==== work-efficient scan, non-power-of-two ====
+   elapsed time: 4.6543ms    (CUDA Measured)
+    [   0  38  82 124 126 142 188 228 261 307 332 356 371 ... 410933976 410933979 ]
+    passed
+==== thrust scan, power-of-two ====
+   elapsed time: 0.569056ms    (CUDA Measured)
+    [   0  38  82 124 126 142 188 228 261 307 332 356 371 ... 410934036 410934075 ]
+    passed
+==== thrust scan, non-power-of-two ====
+   elapsed time: 0.72192ms    (CUDA Measured)
+    [   0  38  82 124 126 142 188 228 261 307 332 356 371 ... 410933976 410933979 ]
+    passed
+
+*****************************
+** STREAM COMPACTION TESTS **
+*****************************
+    [   3   3   0   0   1   1   1   0   0   1   2   1   1 ...   2   0 ]
+==== cpu compact without scan, power-of-two ====
+   elapsed time: 36.8856ms    (std::chrono Measured)
+    [   3   3   1   1   1   1   2   1   1   2   2   1   2 ...   3   2 ]
+    passed
+==== cpu compact without scan, non-power-of-two ====
+   elapsed time: 36.297ms    (std::chrono Measured)
+    [   3   3   1   1   1   1   2   1   1   2   2   1   2 ...   2   1 ]
+    passed
+==== cpu compact with scan ====
+   elapsed time: 74.6997ms    (std::chrono Measured)
+    [   3   3   1   1   1   1   2   1   1   2   2   1   2 ...   3   2 ]
+    passed
+==== work-efficient compact, power-of-two ====
+   elapsed time: 6.25254ms    (CUDA Measured)
+    [   3   3   1   1   1   1   2   1   1   2   2   1   2 ...   3   2 ]
+    passed
+==== work-efficient compact, non-power-of-two ====
+   elapsed time: 6.69901ms    (CUDA Measured)
+    [   3   3   1   1   1   1   2   1   1   2   2   1   2 ...   2   1 ]
+    passed
+
+**********************
+** RADIX SORT TESTS **
+**********************
+    [  45  23   3  11  27  32  15  29   1  47  29   3  22 ...  11   0 ]
+==== Thrust sort, power-of-two ====
+   elapsed time: 0.00176ms    (CUDA Measured)
+    [   0   0   0   0   0   0   0   0   0   0   0   0   0 ...  49  49 ]
+==== Radix Sort, power-of-two ====
+   elapsed time: 74.6997ms    (CUDA Measured)
+    [   0   0   0   0   0   0   0   0   0   0   0   0   0 ...  49  49 ]
+    passed
+==== Thrust sort, non-power-of-two ====
+   elapsed time: 0.00192ms    (CUDA Measured)
+    [   0   0   0   0   0   0   0   0   0   0   0   0   0 ...  49  49 ]
+==== Radix Sort, non-power-of-two ====
+   elapsed time: 42.6619ms    (CUDA Measured)
+    [   0   0   0   0   0   0   0   0   0   0   0   0   0 ...  49  49 ]
+    passed
+Press any key to continue . . .
+```

img/generate_image.py

@@ -0,0 +1,64 @@
+from cProfile import label
+import matplotlib
+import matplotlib.pyplot as plt
+import numpy as np
+
+if __name__ == "__main__":
+    x = np.array([29, 28, 27, 26, 25, 24, 22])
+    x = list(pow(2, x) / 1e+6)
+
+    ## 2 -> NPOT
+    cpu = [1045.5, 417.7, 215.3, 108.297, 53.2, 27.4, 7.6]
+    cpu2 = [901.8, 420.7, 217.0, 104, 54.2, 27.9, 6.6]
+
+    naive = [482.4, 229.5, 108.56, 50.24, 24.2, 13.6, 2.6]
+    naive2 = [474.8, 225.59, 106.0, 52.4, 24.12, 12.1, 2.9]
+
+    eff = [151.74, 74.14, 37.4, 18.21, 9.64, 5.3, 1.23]
+    eff2 = [151.59, 76.6, 37.06, 18.5, 9.89, 6.4, 1.22]
+
+    thrust = [14.6, 9.385, 3.4, 1.84, 1.024, 0.75, 0.25]
+    thrust2 = [14.3, 7.4, 3.4, 1.88, 1.042, 0.66, 1.83]
+
+    compact_cpu = [1092.2, 577.7, 288, 145.03, 72.21, 37.1, 9.11]
+    compact_cpu_scan = [6129.1, 963, 584, 246.8, 117.26, 72.3, 21.7]
+
+    compact_eff = [213.04, 101.7, 51.7, 26.02, 14.95, 6.2, 1.63]
+
+    sort1 = [6129.18, 963.5, 584.5, 246.8, 117.2, 72.3, 21.7]
+    sort2 = [1340.75, 682.7, 338.9, 164.7, 83.30, 42.5, 12.6]
+
+    thrust_sort = [0.0061, 0.006, 0.003, 0.0026, 0.002, 0.002, 0.0018]
+
+    fig = plt.figure()
+    plt.plot(x, cpu, label="CPU Scan")
+    plt.plot(x, naive, label="CUDA Naive Scan")
+    plt.plot(x, eff, label="CUDA Work-Efficient Scan")
+    plt.plot(x, thrust, label="CUDA Thrust Scan")
+
+    plt.xticks([0, 100, 200, 300, 400, 500], ["0", "100M", "200M", "300M", "400M", "500M"])
+
+    plt.legend(fontsize=15)
+
+    plt.title("Execution Time of Scanning (Lower Is Better)")
+    plt.ylabel('Time (ms)')
+    plt.xlabel("Array Size (Million)")
+
+    plt.grid(linestyle='--')
+
+    fig = plt.figure()
+    plt.plot(x, compact_cpu, label="CPU Compact")
+    plt.plot(x, compact_cpu_scan, label="CPU Compact with Scan")
+    plt.plot(x, compact_eff, label="CUDA Work-Efficient Compact")
+
+    plt.xticks([0, 100, 200, 300, 400, 500], ["0", "100M", "200M", "300M", "400M", "500M"])
+
+    plt.legend(fontsize=15)
+
+    plt.title("Execution Time of Stream Compaction (Lower Is Better)")
+    plt.ylabel('Time (ms)')
+    plt.xlabel("Array Size (Million)")
+
+    plt.grid(linestyle='--')
+
+    plt.show()

src/main.cpp

@@ -13,15 +13,15 @@
 #include <stream_compaction/thrust.h>
 #include "testing_helpers.hpp"
 
-const int SIZE = 1 << 8; // feel free to change the size of array
+const int SIZE = 1 << 26; // 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 main(int argc, char* argv[]) {
     // Scan tests
-
+  printf("Using SIZE: %d M", (SIZE / 1000000));
     printf("\n");
     printf("****************\n");
     printf("** SCAN TESTS **\n");
@@ -51,7 +51,7 @@ 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
@@ -64,35 +64,35 @@ int main(int argc, char* argv[]) {
     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(NPOT, 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");
@@ -137,18 +137,55 @@ 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");
+
+
+    genArray(SIZE - 1, a, 50);  // Leave a 0 at the end to test that edge case
+    a[SIZE - 1] = 0;
+    printArray(SIZE, a, true);
+
+    zeroArray(SIZE, b);
+    printDesc("Thrust sort, power-of-two");
+    StreamCompaction::Thrust::sort(SIZE, b, a);
+    printElapsedTime(StreamCompaction::Thrust::timer().getGpuElapsedTimeForPreviousOperation(), "(CUDA Measured)");
+    printArray(SIZE, b, true);
+
+    zeroArray(SIZE, c);
+    printDesc("Radix Sort, power-of-two");
+    StreamCompaction::Efficient::radixSort(SIZE, c, a);
+    printElapsedTime(StreamCompaction::CPU::timer().getCpuElapsedTimeForPreviousOperation(), "(CUDA Measured)");
+    printArray(SIZE, b, true);
+    printCmpResult(SIZE, b, c);
+
+    zeroArray(SIZE, b);
+    printDesc("Thrust sort, non-power-of-two");
+    StreamCompaction::Thrust::sort(NPOT, b, a);
+    printElapsedTime(StreamCompaction::Thrust::timer().getGpuElapsedTimeForPreviousOperation(), "(CUDA Measured)");
+    printArray(NPOT, b, true);
+
+    zeroArray(SIZE, c);
+    printDesc("Radix Sort, non-power-of-two");
+    StreamCompaction::Efficient::radixSort(NPOT, c, a);
+    printElapsedTime(StreamCompaction::Efficient::timer().getGpuElapsedTimeForPreviousOperation(), "(CUDA Measured)");
+    printArray(NPOT, c, true);
+    printCmpResult(NPOT, b, c);
+
     system("pause"); // stop Win32 console from closing on exit
     delete[] a;
     delete[] b;
     delete[] c;
+
 }