CIS565-Fall-2021 · asalexan · Sep 16, 2021 · Sep 16, 2021 · Sep 16, 2021 · Sep 17, 2021
diff --git a/README.md b/README.md
@@ -3,12 +3,155 @@ 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)
+* Ashley Alexander-Lee
+  * [LinkedIn](linkedin.com/in/asalexanderlee), [Personal Website](https://asalexanderlee.myportfolio.com/)
+* Tested on: Windows 10, i7-6700 @ 3.40GHz 16GB, Quadro P1000 (Moore 100A Lab)
 
-### (TODO: Your README)
+### Project Description
+This project involved practice writing parallel algorithms in Cuda, chiefly scan, reduction, and compaction. For comparison, I wrote algorithms on the CPU, as well as naive and work-efficient methods that leverage the GPU.
 
-Include analysis, etc. (Remember, this is public, so don't put
-anything here that you don't want to share with the world.)
+#### Parallel Reduction
+Parallel reduction is an algorithm by which you can add all of the numbers in an array, in parallel. You can achieve this by adding pairs of numbers, looping through the algorithm log(n) times, where n is the size of the array.
+
+<p align="left">
+  <img width="420" src="https://github.com/asalexan/Project2-Stream-Compaction/blob/e9dec569b26e1fd189ed6b49b9c8e05a9cee439e/img/reduction_ex2.png">
+</p>
+
+#### Scan (Prefix Sums)
+A scan of an array produces an array where each element represents a sum of the elements that came before it. This algorithm involves adding numbers based on a logarithmic offset, as illustrated below. The last element in the array will represent the reduction. Below is an illustration of a naive parallel scan algorithm, which runs in O(log(n)) time. 
+
+<img width="340" src="https://github.com/asalexan/Project2-Stream-Compaction/blob/087474d36a1086aee0dc950d9e5858fb543d4c87/img/scan_ex2.png">
+
+We can improve upon this algorithm by doing scan in two steps: 1) parallel reduction, and 2) calculation of partial sums to finish the scan. The parallel reduction occurs as described above. The calculation of partial sums works as illustrated below: 
+
+<img width="330" src="https://github.com/asalexan/Project2-Stream-Compaction/blob/087474d36a1086aee0dc950d9e5858fb543d4c87/img/partials_ex2.png">
+
+#### Compaction
+Compaction involves taking an array of values and returning an array filled with only the values that are "truthy". You can do this by performing the following steps:
+
+1. Convert all "truthy" values to 1s and all others to 0s
+2. Perform scan on the resulting array. The result of scan should give you the desired indices in the output compact array for all true values 
+3. Use the indices to write the true values to a compact array
+
+### Performance Analysis
+
+**Array Size vs. Scan Performance**
+
+Below I've charted how each of the scan algorithms perform side-by-side. The CPU algorithm should run in O(n) time, since it must read every element. The naive and work-efficient algorithms should both run in logarithmic time -- the graph mirrors that below. Note that, since the x-axis uses logarithmic scale, a logarithmic function will appear linear. 
+
+![How Array Size Affects Performance](img/arraySizeVPerformance2.png)
+
+The work-efficient and naive algorithms perform better than the thrust implementation. I imagine this is due to two optimizations:
+
+1. For each loop iteration, I halve the number of threads used to call the kernel -- that way, only threads that will do meaningful work are spawned. For example, if there are 8 elements in an array, I spawn 4 threads during the first loop, 2 during the second, and 1 during the last. If we do the math, we see that much fewer threads are spawned:
+
+    | Iteration | Threads (unoptimized) | Threads (optimized) |
+    | --------- | --------------------- | ------------------- |
+    | 1 | 8 | 4 |
+    | 2 | 8 | 2 |
+    | 3 | 8 | 1 |
+    | total | 24 | 7 |
+
+2. I pre-computed the `pow()` operations before calling the kernels, which more than halved the runtime. I show the major results in a later section.
+
+I'll also note that the work-efficient method is more predictable than the thrust method. As you can see from the chart, thrust runtimes are often unpredictable, and can, for example, run in 0.035ms OR 0.115ms given the same number of elements. I imagine the thrust method involves some sort of sort to move the "truthy" values to the front, which could explain the vastly different runtimes, since sort algorithms can have vastly different best and worst-case runtimes based on whether the array is already sorted. 
+
+**Block Size vs. Scan Performance**
+
+![How Block Size Impacts Performance](img/blockSizeVPerformance.png)
+
+I did not notice a significant change in runtime given different block sizes. I imagine this is due to the fact I am using global memory, instead of shared memory. If I was using shared memory, only one block could be used to run a kernel at a time, since all threads access the same shared memory. Therefore, the block size would change how elements it could process in parallel at once. However, since I utilize global memory, the threads don't depend on each other, so block size doesn't impact the performance.
+
+The block size 256 seems to optimize performance best.
+
+### Bottlenecks
+There are a few bottlenecks I noticed, especially in the compaction algorithm. 
+
+1. The major roadblock is memory I/O. I am using global memory, which is very slow, and I have to copy data from device to host and allocate a device array in the compaction algorithm. If I reworked the code so that it used shared memory and cut out the device pointer as a middleman, it would be much faster. 
+2. The threads are regularly accessing memory non-contiguously, leading to a lot of overhead time fetching data.
+3. On the naive algorithm, there is a lot of branching in the kernel, meaning a number of threads are waiting for other threads to finish their execution. 
+
+### Improving Scan Performance
+I wanted to understand why my GPU scan implementations were running so horribly, so I decided to time the reduction and partial sum kernels in the work-efficient algorithm. On an input of 1^8, each kernel was taking about 0.05ms to run. In the interest of improving those times, I performed the `pow()` calculations outside of the kernel and passed them in as parameters, since I realized it was likely a costly operation that I only needed to perform once per kernel call. After removing the `pow()` operations, I more than halved the runtime for both kernels, improving the overall runtime so much that it now regularly outperforms the thrust implementation. I also made this change to my naive implementation, and saw similar gains. You can see some side-by-side comparisons below (tested on array size 1^8).
+
+| Function | Previous Runtime | Improved Runtime|
+| -------- | ---------------- | ---------------- |
+| Reduction (Up-Sweep) | 0.06 | 0.016|
+| Partial Sums (Down-Sweep) | 0.058 | 0.017 |
+| Work-Efficient Scan | 0.1157 | 0.035 |
+| Naive Scan | 0.0563 | 0.0245 |
+
+You can see the huge change in made in performance by comparing the charted performance before the `pow()` removals and after.
+
+**Previous Scan Performance**
+
+![Previous Scan Performance](img/arraySizeVPerformance.png)
+
+**Scan Performance After Improvements**
+
+![Scan Performance After Improvements](img/arraySizeVPerformance2.png)
+
+### Sample Output
+
+```****************
+** SCAN TESTS **
+****************
+    [  31   2  26  37  46  20  20   8  38  37  28  41  20 ...   5   0 ]
+==== cpu scan, power-of-two ====
+   elapsed time: 0.0005ms    (std::chrono Measured)
+    [   0  31  33  59  96 142 162 182 190 228 265 293 334 ... 6434 6439 ]
+==== cpu scan, non-power-of-two ====
+   elapsed time: 0.0004ms    (std::chrono Measured)
+    [   0  31  33  59  96 142 162 182 190 228 265 293 334 ... 6357 6369 ]
+    passed
+==== naive scan, power-of-two ====
+   elapsed time: 0.018432ms    (CUDA Measured)
+    [   0  31  33  59  96 142 162 182 190 228 265 293 334 ... 6434 6439 ]
+    passed
+==== naive scan, non-power-of-two ====
+   elapsed time: 0.018432ms    (CUDA Measured)
+    [   0  31  33  59  96 142 162 182 190 228 265 293 334 ...   0   0 ]
+    passed
+==== work-efficient scan, power-of-two ====
+   elapsed time: 0.039936ms    (CUDA Measured)
+    [   0  31  33  59  96 142 162 182 190 228 265 293 334 ... 6434 6439 ]
+    passed
+==== work-efficient scan, non-power-of-two ====
+   elapsed time: 0.032768ms    (CUDA Measured)
+    [   0  31  33  59  96 142 162 182 190 228 265 293 334 ... 6357 6369 ]
+    passed
+==== thrust scan, power-of-two ====
+   elapsed time: 0.044224ms    (CUDA Measured)
+    [   0  31  33  59  96 142 162 182 190 228 265 293 334 ... 6434 6439 ]
+    passed
+==== thrust scan, non-power-of-two ====
+   elapsed time: 0.032032ms    (CUDA Measured)
+    [   0  31  33  59  96 142 162 182 190 228 265 293 334 ... 6357 6369 ]
+    passed
+
+*****************************
+** STREAM COMPACTION TESTS **
+*****************************
+    [   1   2   0   1   2   0   2   0   0   1   0   1   0 ...   3   0 ]
+==== cpu compact without scan, power-of-two ====
+   elapsed time: 0.0008ms    (std::chrono Measured)
+    [   1   2   1   2   2   1   1   1   2   3   2   1   2 ...   2   3 ]
+    passed
+==== cpu compact without scan, non-power-of-two ====
+   elapsed time: 0.0009ms    (std::chrono Measured)
+    [   1   2   1   2   2   1   1   1   2   3   2   1   2 ...   3   1 ]
+    passed
+==== cpu compact with scan ====
+   elapsed time: 0.0012ms    (std::chrono Measured)
+    [   1   2   1   2   2   1   1   1   2   3   2   1   2 ...   2   3 ]
+    passed
+==== work-efficient compact, power-of-two ====
+   elapsed time: 0.22016ms    (CUDA Measured)
+    [   1   2   1   2   2   1   1   1   2   3   2   1   2 ...   2   3 ]
+    passed
+==== work-efficient compact, non-power-of-two ====
+   elapsed time: 0.257024ms    (CUDA Measured)
+    [   1   2   1   2   2   1   1   1   2   3   2   1   2 ...   3   1 ]
+    passed
+```
 
diff --git a/img/arraySizeVPerformance.png b/img/arraySizeVPerformance.png
diff --git a/img/arraySizeVPerformance2.png b/img/arraySizeVPerformance2.png
diff --git a/img/blockSizeVPerformance.png b/img/blockSizeVPerformance.png
diff --git a/img/naive_scan_ex.JPG b/img/naive_scan_ex.JPG
diff --git a/img/partial_sums_ex.JPG b/img/partial_sums_ex.JPG
diff --git a/img/partials_ex2.png b/img/partials_ex2.png
diff --git a/img/reduction_ex.JPG b/img/reduction_ex.JPG
diff --git a/img/reduction_ex2.png b/img/reduction_ex2.png
diff --git a/img/scan_ex2.png b/img/scan_ex2.png
diff --git a/src/main.cpp b/src/main.cpp
@@ -13,7 +13,7 @@
 #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 << 8;
 const int NPOT = SIZE - 3; // Non-Power-Of-Two
 int *a = new int[SIZE];
 int *b = new int[SIZE];
@@ -38,13 +38,13 @@ int main(int argc, char* argv[]) {
     printDesc("cpu scan, power-of-two");
     StreamCompaction::CPU::scan(SIZE, b, a);
     printElapsedTime(StreamCompaction::CPU::timer().getCpuElapsedTimeForPreviousOperation(), "(std::chrono Measured)");
-    printArray(SIZE, b, true);
+    //printArray(SIZE, b, true);
 
     zeroArray(SIZE, c);
     printDesc("cpu scan, non-power-of-two");
     StreamCompaction::CPU::scan(NPOT, c, a);
     printElapsedTime(StreamCompaction::CPU::timer().getCpuElapsedTimeForPreviousOperation(), "(std::chrono Measured)");
-    printArray(NPOT, b, true);
+    //printArray(NPOT, b, true);
     printCmpResult(NPOT, b, c);
 
     zeroArray(SIZE, c);
@@ -53,48 +53,48 @@ int main(int argc, char* argv[]) {
     printElapsedTime(StreamCompaction::Naive::timer().getGpuElapsedTimeForPreviousOperation(), "(CUDA Measured)");
     //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); */
 
-    zeroArray(SIZE, c);
+    /*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");
     printf("*****************************\n");
     printf("** STREAM COMPACTION TESTS **\n");
@@ -125,28 +125,28 @@ int main(int argc, char* argv[]) {
     expectedNPOT = count;
     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);
     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);
     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);
-
+    
     system("pause"); // stop Win32 console from closing on exit
     delete[] a;
     delete[] b;

diff --git a/stream_compaction/common.cu b/stream_compaction/common.cu
@@ -23,7 +23,8 @@ 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 = (blockDim.x * blockIdx.x) + threadIdx.x;
+            bools[index] = idata[index] ? 1 : 0;
         }
 
         /**
@@ -32,7 +33,11 @@ namespace StreamCompaction {
          */
         __global__ void kernScatter(int n, int *odata,
                 const int *idata, const int *bools, const int *indices) {
-            // TODO
+            int index = (blockDim.x * blockIdx.x) + threadIdx.x;
+            if (bools[index] == 1) {
+                int test = indices[index];
+                odata[indices[index]] = idata[index];
+            }
         }
 
     }

diff --git a/stream_compaction/cpu.cu b/stream_compaction/cpu.cu
@@ -1,6 +1,6 @@
 #include <cstdio>
 #include "cpu.h"
-
+#include <iostream>
 #include "common.h"
 
 namespace StreamCompaction {
@@ -19,7 +19,10 @@ namespace StreamCompaction {
          */
         void scan(int n, int *odata, const int *idata) {
             timer().startCpuTimer();
-            // TODO
+            odata[0] = 0; // explicit starts with 0
+            for (int i = 0; i < n - 1; i++){
+                odata[i+1] = odata[i] + idata[i];   
+            }
             timer().endCpuTimer();
         }
 
@@ -30,9 +33,18 @@ namespace StreamCompaction {
          */
         int compactWithoutScan(int n, int *odata, const int *idata) {
             timer().startCpuTimer();
-            // TODO
+
+            int compactSize = 0;
+            for (int i = 0; i < n; i++){
+                bool elem = idata[i];
+                if (elem){
+                    odata[compactSize] = idata[i]; 
+                    compactSize++;
+                }
+            }
+
             timer().endCpuTimer();
-            return -1;
+            return compactSize;
         }
 
         /**
@@ -42,9 +54,25 @@ namespace StreamCompaction {
          */
         int compactWithScan(int n, int *odata, const int *idata) {
             timer().startCpuTimer();
-            // TODO
+
+            // scan array, replacing true vals with 1s
+            odata[0] = 0;
+            for (int i = 0; i < n; i++) {
+                int elem = idata[i] ? 1 : 0;  
+                odata[i + 1] = odata[i] + elem;
+            }
+
+            // use scatter to produce output
+            int size = odata[n-1];
+            for (int i = 0; i < n; i++){
+                if (idata[i]){ 
+                    int outIndex = odata[i];
+                    odata[outIndex] = idata[i]; 
+                }    
+            }
+
             timer().endCpuTimer();
-            return -1;
+            return size;
         }
     }
 }