Illustration of ordering dependencies for rendering
This is a course project for CMU 15-418/618: Parallel Computer Architecture and Programming. It was a pair project, so the write-up says "we"; the Learnings at the end are mine.
At a glance
- Course: CMU 15-418/618 Parallel Computer Architecture and Programming, Fall 2022
- Stack: CUDA C++, Thrust
- Scope: a circle renderer, taken through four increasingly parallel implementations to a full-score result
- Team: group project
- Code: not public (course integrity policy)
Introduction
In this project we built a parallel image renderer in CUDA. It takes a scene description file and outputs a rendered image. The scene file specifies the image size, the background color, and a list of circles, each with a center, radius, color and velocity (the circles move around the scene).
For each pixel, the renderer works out which circles overlap it and blends their colors. The input is an array of circles (3D position, velocity, radius, color). The basic sequential algorithm for one frame is:
Clear image
For each circle:
Update position and velocity
For each circle:
Compute screen bounding box
For all pixels in bounding box:
Compute pixel center point
If center point is within the circle:
Compute color of circle at point
Blend contribution of circle into image for this pixel
The circles are semi-transparent, so a pixel's color is not the color of any one circle but the blend of every circle overlapping it (the “blend contribution” step in the pseudocode). The renderer represents the color of a circle via a 4-tuple of red (R), green (G), blue (B), and opacity (alpha, written “α”) values (RGBA). Alpha value α = 1.0 corresponds to a fully opaque circle. Alpha value α = 0.0 corresponds to a fully transparent circle. To draw a semi-transparent circle with color (Cr , Cg , Cb , α) on top of a pixel with color Pr , Pg , Pb , the renderer performs the following computation:
Composition is not commutative (X over Y does not look like Y over X), so circles have to be rendered in the order the application provides them, which is depth order. In Figure 2 the circles are ordered red, green, blue: the left image renders them in that order, the right image reverses it.
Implementation
The next three sections cover our first three implementations, with the nvprof data and the reasoning that pushed us on to the next one. The final implementation gets a fuller description.
Implementation 1: A simple parallel renderer
Our first implementation had two phases. The first set of kernel launches parallelized over blocks of pixels: one thread per 32 by 32 block of pixels, so a CUDA block covered a 32 by 32 grid of such pixel blocks. One CUDA block was launched for every (circle, block of blocks) pair. With many circles that meant several kernel launches, because CUDA limits the z-dimension of the block grid to 65535 (spec).
This first kernel computed which circles overlapped each 32 by 32 block of pixels. It filled a large cudaMalloc'd array with a mapping from (pixel block, circle) to whether that circle touched the block.
The second kernel launched one block per 32 by 32 block of pixels, one thread per pixel. Each thread stepped through every circle in order, checked the mapping from the first kernel, and called shadePixel when the circle was in the block.
This approach had several inefficiencies, as the scores below show. The first kernel launched roughly 1300 times the number of circles in blocks (on 1150 by 1150 images). Shared memory went unused and the threads did not collaborate. And the block-circle mapping stored ints for what was a single bit of information.
------------
Score table:
------------
-------------------------------------------------------------------------
| Scene Name | Target Time | Your Time | Score |
-------------------------------------------------------------------------
| rgb | 0.1850 | 0.1291 | 12 |
| rand10k | 1.5643 | 25.9128 | 2 |
| rand100k | 14.8853 | 209.7083 | 2 |
| pattern | 0.2860 | 3.1447 | 2 |
| snowsingle | 6.2212 | 208.2952 | 2 |
| biglittle | 11.4153 | 53.4674 | 4 |
-------------------------------------------------------------------------
| | Total score: | 24/72 |
----------------------------------------------------------------------
The relevant part of nvprof ./render rand10k -r cuda --bench 0:1 is below. For the rand10k scene, the renderer spent most of its time in kernelCirclesShade. We had no other memory profiling tools, so nvprof and the score tables were our only measurements.
Overall: 0.0455 sec (note units are seconds)
==83822== Profiling application: ./render rand10k -r cuda --bench 0:1
==83822== Profiling result:
Type Time(%) Time Calls Avg Min Max Name
GPU activities: 90.08% 23.218ms 1 23.218ms 23.218ms 23.218ms kernelCirclesShade(int, int, float, float, int, int, char*)
9.54% 2.4582ms 1 2.4582ms 2.4582ms 2.4582ms kernelCirclesInclusion(float, float, int, int, int, int, char*)
0.20% 52.544us 1 52.544us 52.544us 52.544us kernelClearImage(float, float, float, float)
0.17% 44.865us 9 4.9850us 1.1840us 11.841us [CUDA memcpy HtoD]
Implementation 2: Renderer taking advantage of shared memory
This implementation started addressing those issues. To get threads collaborating through shared memory, launch fewer blocks, and touch global memory less, we mapped blocks and threads to work differently. There were still two kernel launches.
The first kernel launched one block of 32 by 32 threads per 32 by 32 block of pixels, and the block's 1024 threads stepped through that pixel block's row of the inclusion mapping. This spawned far fewer blocks than before, when there were far too many relative to the number of CUDA cores. The mapping now stored chars instead of ints, which cut the bytes stored and later loaded from global memory.
The second kernel also launched one block of 32 by 32 threads per pixel block. Its 1024 threads cooperatively loaded 1024 entries of the block's row of the mapping, ran sharedMemExclusiveScan to compute indices for the circles marked as in the block, and wrote those circle indices to the front of a shared array. Then each thread took one pixel and stepped through the selected circles in order, calling shadePixel.
Two problems remained. We still allocated a very large array in global memory with cudaMalloc, and both kernels went to global memory constantly. The harder scenes suffered, as the score table below shows.
------------
Score table:
------------
-------------------------------------------------------------------------
| Scene Name | Target Time | Your Time | Score |
-------------------------------------------------------------------------
| rgb | 0.1991 | 0.1750 | 12 |
| rand10k | 1.9623 | 3.2837 | 9 |
| rand100k | 14.7687 | 32.8902 | 7 |
| pattern | 0.2845 | 0.2918 | 12 |
| snowsingle | 7.7129 | 8.9162 | 12 |
| biglittle | 11.7743 | 33.2722 | 6 |
-------------------------------------------------------------------------
| | Total score: | 58/72 |
-------------------------------------------------------------------------
nvprof ./render rand10k -r cuda --bench 0:1 showed the renderer still spending most of its time in the second kernel, though compared with the previous run the numbers confirm the first kernel got much cheaper.
Overall: 0.0207 sec (note units are seconds)
==84635== Profiling application: ./render rand10k -r cuda --bench 0:1
==84635== Profiling result:
Type Time(%) Time Calls Avg Min Max Name
GPU activities: 91.03% 2.9157ms 1 2.9157ms 2.9157ms 2.9157ms kernelCirclesShade(int, int, float, float, int, int, char*)
5.94% 190.21us 1 190.21us 190.21us 190.21us kernelCirclesInclusion(float, float, int, int, int, char*)
1.63% 52.288us 1 52.288us 52.288us 52.288us kernelClearImage(float, float, float, float)
1.40% 44.832us 9 4.9810us 1.2480us 11.904us [CUDA memcpy HtoD]
Implementation 3: Reducing Global Memory Access
To cut global memory traffic, we merged the two kernels into one, still with one block of 32 by 32 threads per 32 by 32 pixel block. Walking the circle array in chunks of 1024, the block's threads cooperatively loaded the circle data, tested each circle against the block, ran sharedMemExclusiveScan to index the circles that were inside, and wrote those indices to the front of a shared array. Then each thread took one pixel, stepped through the shared array, loaded the circle data from global memory and called shadePixel.
One issue remained: we did not keep a local copy of the pixel data, so every update went through imagePtr in global memory. Even so, the score table below shows a big improvement.
------------
Score table:
------------
-------------------------------------------------------------------------
| Scene Name | Target Time | Your Time | Score |
-------------------------------------------------------------------------
| rgb | 0.1862 | 0.1227 | 12 |
| rand10k | 1.9471 | 2.2303 | 12 |
| rand100k | 15.1151 | 21.5748 | 10 |
| pattern | 0.2704 | 0.2080 | 12 |
| snowsingle | 7.7121 | 4.7711 | 12 |
| biglittle | 13.1609 | 25.3954 | 8 |
-------------------------------------------------------------------------
| | Total score: | 66/72 |
-------------------------------------------------------------------------
The nvprof ./render rand100k -r cuda --bench 0:2 output below is less informative than the earlier ones: with all the work in one kernel, it can only tell us that almost all the time is spent there.
Overall: 0.0288 sec (note units are seconds)
==87385== Profiling application: ./render rand100k -r cuda --bench 0:2
==87385== Warning: 1 records have invalid timestamps due to insufficient device buffer space. You can configure the buffer space using the option --device-buffer-size.
==87385== Profiling result:
Type Time(%) Time Calls Avg Min Max Name
GPU activities: 98.49% 28.671ms 1 28.671ms 28.671ms 28.671ms kernelCirclesShade(int, int, float, float, int, int)
1.15% 334.81us 9 37.201us 1.1840us 99.551us [CUDA memcpy HtoD]
0.36% 104.61us 2 52.304us 52.192us 52.416us kernelClearImage(float, float, float, float)
Implementation 4: Final Implementation
The final improvement was to copy each pixel's data from global memory into a local variable, apply every update there, and write it back once at the end. This removed most of the global memory traffic.
The final renderer works like this. The image is divided into 32 by 32 pixel blocks, and a single kernel launches one block of 32 by 32 threads per pixel block. A thread whose pixel lies inside the image loads that pixel's initial data. The block then walks the circle array in global memory 1024 elements at a time. For each chunk, the 1024 threads load the circle data together, and each thread checks whether its circle overlaps the block, writing a 1 into a shared array if so. A sharedMemExclusiveScan over those 1s and 0s yields an index for every overlapping circle, and the threads use those indices to write the circle indices into a second shared array. Each in-bounds thread then steps through that array, loads the circle data from global memory, and shades its pixel; shadePixel reads and writes the pixel in local memory, not global. The block moves on to the next 1024 circles. When every circle has been processed, each in-bounds thread writes its pixel back to global memory.
Updates are race-free because exactly one thread owns each pixel and applies circles to it one at a time. Circle order is preserved because chunks are processed from the start of the array to the end, and within a chunk the scan gathers the relevant circles without reordering them, so pixels see the circles in their original order.
The score table for our final implementation is included below. We ran this code on GHC 67.
------------
Score table:
------------
-------------------------------------------------------------------------
| Scene Name | Target Time | Your Time | Score |
-------------------------------------------------------------------------
| rgb | 0.1865 | 0.1753 | 12 |
| rand10k | 1.9469 | 1.9401 | 12 |
| rand100k | 18.4898 | 18.3768 | 12 |
| pattern | 0.2724 | 0.2465 | 12 |
| snowsingle | 7.7121 | 6.2476 | 12 |
| biglittle | 14.1915 | 16.7866 | 12 |
-------------------------------------------------------------------------
| | Total score: | 72/72 |
-------------------------------------------------------------------------
nvprof ./render rand100k -r cuda --bench 0:1 confirmed that this version is significantly faster than the previous one.
Overall: 0.0354 sec (note units are seconds)
==88283== Profiling application: ./render rand100k -r cuda --bench 0:1
==88283== Profiling result:
Type Time(%) Time Calls Avg Min Max Name
GPU activities: 97.93% 18.317ms 1 18.317ms 18.317ms 18.317ms kernelCirclesShade(int, int)
1.80% 336.00us 9 37.333us 1.2160us 100.16us [CUDA memcpy HtoD]
0.27% 50.848us 1 50.848us 50.848us 50.848us kernelClearImage(float, float, float, float)
The decomposition, in short: each CUDA block owns one 32 by 32 pixel block, and each thread plays two roles. First it stands for one circle, testing it against the block and taking part in the scan, and if its circle overlaps it writes the circle's index into the shared array. Then it stands for one pixel, stepping through that shared array, loading the circle data and shading. The provided scan code decomposes the scan over warps.
Synchronization happens in two places, both inside the kernel's main loop. The first is the set of __syncthreads() calls inside the provided sharedMemExclusiveScan, after the circle data has been loaded and tested for inclusion. No barrier is needed before the scan, because the code runs in warps and each warp has loaded everything it needs before the warp-level scan; the __syncthreads() after it covers the rest. The second is after the threads have written the overlapping circle indices to shared memory, so that no thread starts shading before the block's full circle list is in place.
We reduced the cost of synchronization in three ways. Merging into one kernel removed a cudaDeviceSynchronize(). Reasoning about sharedMemExclusiveScan as above let us remove redundant __syncthreads() calls. And we removed the barrier at the end of the loop: after the last synchronization point, threads need the last elements of the scan's input and output arrays plus the shared circle-index array, and only the scan's input array can be modified before the next barrier. So we copy that last input value into shared memory before the barrier, and the end-of-loop __syncthreads() becomes unnecessary.
We also eliminated unnecessary global memory accesses. Merging the two kernels meant the intermediate results no longer had to live in global memory: the per-block circle intersection results sit in shared memory, where the shading logic reads them without going back to global memory.
And each thread loads its pixel into local memory (registers) first and writes it back only after every circle's contribution has been applied. This sped the renderer up considerably and brought the score to full marks.
Learnings
The four implementations above are the real story of this project. These are the lessons I'd pull out of them.
- Profile first, then optimize what the profiler points at. Every step from implementation 1 to 4 started with
nvprofoutput, and every time the bottleneck was somewhere slightly different from where we guessed: kernel launches, then global memory, then synchronization. - The biggest win came from re-partitioning the work, not from tuning the kernel. Going from one thread per circle to one block per image tile, with the threads cooperating on the circles that touch it, changed the memory access pattern and the synchronization needs at the same time. Micro-optimizations inside the old decomposition were never going to get there.
- Correctness constraints shape the parallel decomposition. Circles must be composited in input order, which rules out the naive parallel approach outright. Working out which order dependencies are real before choosing a decomposition saved us from optimizing something that could never be correct.
- Shared memory plus a scan is a reusable pattern. Each thread tests a candidate, an exclusive scan compacts the survivors into shared memory, then everyone consumes the compact list. The same shape showed up again in later GPU work, and learning it here with the provided
sharedMemExclusiveScanwas worth more than the assignment itself. - Every
__syncthreads()is a design decision. Removing redundant barriers was measurable, but each removal needed a written argument for why it was safe. Justifying each barrier you keep and each one you drop is a habit I still use.