I originally posted this to /r/hardware and /r/AMD, and people seemed to like it. I discovered this subreddit, so I'll copy/paste my post. Hopefully someone out there will find it useful!

This seems like a relatively slow subreddit, but I think there are enough beginners here that maybe this post will be useful.

In this topic, I'm just going to stream some ideas about what I know about Vega64. I hope I can inspire some programmers to try to program their GPU! Also, If anyone has more experience programming GPUs (NVidia ones even), please chime in!

For the most part, I assume that the reader is a decent C Programmer who doesn't know anything about GPUs or SIMD.

Vega Introduction

Before going further, I feel like its important to define a few things for AMD's Vega Architecture. I will come back later to better describe some concepts.

• 64 CUs (Compute Units) -- 64 CUs on Vega64. 56 CUs on Vega56.

  • 16kB L1 (Level 1) data-cache per CU

  • 64kB LDS (Local Data Store) per CU

  • 4-vALUs (vector Arithmetic Logic Unit) per CU

    • 16 PE (Processing Elements) per vALU

    • 4 x 256 vGPRs (vector General Purpose Registers) per PE

  • 1-sALU (scalar Arithmetic Logic Unit) per CU

• 8GB of HBM2 RAM

Grand Total: 64 CUs x 4 vALUs x 16 PEs == 4096 "shaders", just as advertised. I'll go into more detail later what a vGPR or sGPR is, but lets first cover the programmer-model.

GPU Programming in a nutshell

Here's some simple C code. Lets assume "x" and "y" are the input to the problem, and "output" is the output:

for(int i=0; i<1000000; i++){ // One Million Items
    output[i] = x[i] + y[i];
}

• "Work Items", (SIMD Threads in CUDA) are the individual units of work that the programmer wishes to accomplish in parallel with each other. Given the example above, a good work item would be "output[i] = x[i] + y[i]". You would have one-million of these commands, and the programmer instinctively knows that all 1-million of these statements could be executed in parallel. OpenCL, CUDA, HCC, and other grossly-parallel languages are designed to help the programmer specify millions of work-items that can be run on a GPU.

• "NDRange" ("Grid" in CUDA) specifies the size of your work items. In the example "for loop" case above, 1000000 would be the NDRange. Aka: there are 1-million things to do. The NDRange or Grid may be 2-dimentional (for 2d images), or 3-dimentional (for videos).

• "Wavefronts" ("Warp" in CUDA) are the smallest group of work-items that a GPU can work on at a time. In the case of Vega, 64-work items constitutes a Wavefront. In the case of the for-loop described earlier, a wave-front would execute between [0, 1, 2, 3... 63] iterations together. A 2nd wave front would execute [64, 65, 66, 67, ... 127] together (and probably in parallel).

• "Workgroups" ("Thread Blocks" in CUDA) are logical groups that the programmer wants to work together. While Wavefronts are what the system actually executes, the Vega system can combine up to 16-wavefronts together and logically work as a single Workgroup. Vega64 supports workgroups of size 1 through 16 Wavefronts, which correlates to 64, 128, ... 1024 WorkItems (1024 == 16 WaveFronts * 64 Threads per Wavefront).

In summary: OpenCL / CUDA programmers setup their code. First, they specify a very large number of work items (or CUDA Threads) which represents parallelism. For example: perhaps you want to calculate something on every pixel of a picture, or calculate individual "Rays" of a Raytracer. The programmer then groups the work items into workgroups. Finally, the GPU itself splits workgroups into Wavefronts (64-threads on Vega).

SIMD Segue

Have you ever tried controlling multiple characters with only one controller? When you hook up one controller, but somehow trick the computer into thinking it is 8-different-controllers? SIMD: Single Instruction Multiple Data, is the GPU-technique for actually executing these thousands-of-threads efficiently.

The chief "innovation" of GPUs is just this multi-control concept, but applied to data instead. Instead of building these huge CPU cores which can execute different threads, you build tiny GPU cores (or shaders) which are forced to play the same program. Instead of 8x wide (like in the GIF I shared), its 64x wide on AMD.

To handle "if" statements or "loops" (which may vary between work-items), there's an additional "execution mask" which the GPU can control. If the execution-mask is "off", an individual thread can be turned off. For example:

if(foo()){
    doA(); // Lets say 10 threads want to do this
} else {
    doB(); // But 54 threads want to do this
}

The 64-threads of the wavefront will be forced to doA() first, with the 10-threads having "execution mask == on", and with the 54-remaining threads having "execution mask == off". Then, doB() will happen next, with 10-threads off, and 54-threads on. This means that any "if-else" statement on a GPU will have BOTH LEGS executed by all threads.

In general, this is called the "thread divergence" problem. The more your threads "split up", the more legs of if-statements (and more generally: loops) have to be executed.

Before I reintroduce Vega's Architecture, keep the multiple-characters / one-controller concept in mind.

Vega Re-Introduction

So here's the crazy part. A single Vega CU doesn't execute just one wavefront at a time. The CU is designed to run upto 40 wavefronts (x 64 threads, so 2560 threads total). These threads don't really all execute simultaneously: the 40-wavefronts are there to give the GPU something to do while waiting for RAM.

Vega's main memory controller can take 350ns or longer to respond. For a 1200MHz system like Vega64, that is 420 cycles of waiting whenever something needs to be fetched from memory. That's a long time to wait! So the overall goal of the system, is to have lots of wavefronts ready to run.

With that out of the way, lets dive back into Vega's architecture. This time focusing on CUs, vALUs, and sALUs.

• 64 CUs (Compute Units) -- 64 CUs on Vega64.

  • 4-vALUs (vector Arithmetic Logic Unit) per CU

    • 16 PE (Processing Elements) per vALU

    • 4 x 256 vGPRs (vector General Purpose Register) per PE

  • 1-sALU (scalar Arithmetic Logic Unit) per CU

The sALU is easiest to explain: sALUs is what handles those "if" statements and "while" statements I talked about in the SIMD section above. sALUs track which threads are "executing" and which aren't. sALUs also handle constants and a couple of other nifty things.

Second order of business: vALUs. The vALUs are where Vega actually gets all of their math power from. While sALUs are needed to build the illusion of wavefronts, vALUs truly execute the wavefront. But how? With only 16-PEs per vALU, how does a wavefront of size 64 actually work?

And btw: your first guess is likely wrong. It is NOT from #vALUs x 16 PEs. Yes, this number is 64, but its an utterly wrong explanation which tripped me up the first time.

The dirty little secret is that each PE repeats itself 4-times in a row, across 4-cycles. This is a hidden fact deep in AMD documentation. In any case, 4-cycles x 16 PE == 64 Workitems per vALU. x4 vALUs == 256 work-items per Compute Unit (every 4 clock cycles).

Why repeat themselves? Because if a simple addition takes 4-clock cycles to operate, then Vega only has to perform 1/4th the math operations while waiting for RAM. (IE: for the 420 cycle wait on a 350ns RAM load... you can "fill" those 420 cycles with only 105 math operations!). Repeating commands over-and-over again helps Vega to hide the memory-latency problem.

Full Occupancy: 4-clocks x 16 PEs x 4 vALUs == 256 Work Items

Full Occupancy, or more like "Occupancy 1", is when each CU (compute unit) has one-work item for each physical thread that could run. Across the 4-clock cycles, 16 PEs, and 4 vALUs per CU, the Compute Unit reaches full occupancy at 256 work items (or 4-Wavefronts).

Alas: RAM is slow. So very, very slow. Even at Occupancy 1 with super-powered HBM2 RAM, Vega would spend too much time waiting for RAM. As such, Vega supports "Occupany 10"... but only IF the programmer can split the limited resources between threads.

In practice, programmers typically reach "Occupancy 4". At occupancy 4, the CU still only executes 256-work items every 4-clock cycles (4-wavefronts), but the 1024 total items (16-wavefronts) give the CU "extra work" to do whenever it notices that one wavefront is waiting for RAM.

Memory hiding problem

Main Memory latency is incredibly slow, but also is variable. RAM may take 350 or more cycles to respond. Even LDS, may respond in a variable amount of time (depending on how many atomic operations are going on, or bank-conflicts).

AMD has two primary mechanisms to hide memory latency.

1. Instruction Level -- AMD's assembly language requires explicit wait-states to hold the pipeline. The "s_waitcnt lgkmcnt(0)" instruction you see in the assembly is just that: wait for local/global/konstant/message counter to be (zero). Careful use of the s_waitcnt instruction can be used to hide latency behind calculations: you can start a memory load to some vGPRs, and then calculate with other vGPRs before waiting.

2. Wavefront Level -- The wavefronts at a system-level allow the CU to find other work, just in case any particular wavefront gets stuck on a s_waitcnt instruction.

While CPUs use out-of-order execution to hide latency and search for instruction-level parallelism... GPUs require the programmer (or compiler) to explicitly put the wait-states in. It is far less flexible, but far cheaper an option to do.

Wavefront level latency hiding is roughly equivalent to a CPU's SMT / Hyperthreading. Except instead of 2-way hyperthreading, the Vega GPU supports 10-way hyperthreads.

Misc. Optimization Notes

• On AMD Systems, 64 is your magic minimum number. Try to have at least 64 threads running at any given time. Ideally, have your workload evenly divisible by 64. For example, 100 threads will be run as 64 thread wavefront + 36 thread wavefront (with 28 wasted vALU states!). 128 threads is more efficient.

• vGPRs (vector General Purpose Registers) are your most precious resource. Each vGPR is a 32-bit of memory that executes at the full speed of Vega (1-operation every 4 clock cycles). Any add, subtract, or multiply in any work-item will have to travel through a vGPR before it can be manipulated. vGPRs roughly correlate to "OpenCL Private Memory", or "CUDA Local Memory".

• At occupancy 1, you can use all 256 vGPRs (1024 bytes). However, "Occupancy 1" is not good enough to keep the GPU busy when its waiting for RAM. The extreme case of "Occupancy 10" gives you only 25 vGPRs to work with (256/10, rounded down). A reasonable occupancy to aim for is Occupancy 4 and above (64 vGPRs at Occupancy 4)

• FP16 Packed Floats will stuff 2x16-bit floats per vGPR. "Pack" things more tightly to save vGPRs and achieve higher occupancy.

• The OpenCL Compiler, as well as HCC, HIP, Vulkan compilers, will overflow OpenCL Private Memory into main-memory (Vega's HBM2) if it doesn't fit into vGPRs. There are compiler flags to tune how many vGPRs the compiler will target. However, your code will be waiting for RAM on an overflow, which is counterproductive. Expect a lot of compiler-tweaking to figure out what the optimal vGPRs for your code will be.

• sGPRs (scalar General Purpose Registers) are similarly precious, but Vega has a lot more of them. I believe Vega has around 800 SGPRs per SIMD unit. That is 4x800 SGPRs per CU. Unfortunately, Vega has an assembly-language limit of 102 SGPRs allocated per wavefront. But an occupancy 8 Vega system should be able to hold 100 sGPRs per wavefront.

• The OpenCL Contant Memory specification is often optimized into sGPRs (but not always). In essence: as long as they are uniform across the 64-item wavefront, an sGPR can be used instead of 64-individual precious vGPRs (and constants most often are uniform. But not always: a vGPR based constant happens if vGPRs index into an array of constants). One example of non-constant use of sGPRs is a uniform for-loop, like "for(int i=0; i<10; i++) {}". Instead of taking up 64 vGPRs (across the 64-work item wavefront), this for loop can be implemented with a single sGPR.

• If you can branch using sGPR registers ("constant" across the whole 64-item wavefront), then you will not need to execute the "else". Effectively, sGPR branching never has a divergence problem. sGPR-based branching and looping has absolutely no penalty on the Vega architecture. (In contrast, vGPR-based branching will cause thread-divergence).

• The sALU can operate on 64-bit integers. sGPRs are of size 32-bits, and so any 64-bit operation will use two sGPRs. There is absolutely no floating-point support on the sALU.

• LDS (Local Data Store) is the 2nd fastest RAM, and is therefore the 2nd most important resource after vGPRs. LDS RAM correlates to "OpenCL Local" and "CUDA Shared". (Yes, "Local" means different things between CUDA and OpenCL. Its very confusing). There is 64kB of LDS per CU.

• LDS can share data between anything within your workgroup. The LDS is the primary reason to use a large 1024-thread workgroup: the workgroup can share the entire LDS space. LDS has full support of atomics (ex: CAS) to provide a basis of thread-safe communications.

• LDS is roughly 32-banks (per CU) of RAM which can be issued every clock tick under ideal circumstances. (/u/Qesa claims ~20 cycles of latency best case). At 1200 MHz (Vega64 base clock), this would give the LDS 153GBps of bandwidth per CU. Across the 64-CUs of Vega64, that's a grand total of 9830.4 GBps bandwidth (and it goes faster as Vega boost-clocks!). Compared to HBM2, which is only 483.8 GBps, you can see why proper use of the LDS can accelerate your code.

• Occupancy will force you to split the LDS. The absolute calculation is harder to formulate, because the LDS is shared by Workgroups (and there can be 1 to 16 wavefronts per workgroup). If you have 40 Workgroups (1-wavefront per workgroup), the 64kB LDS must be split into 1638 bytes between workgroups. However, if there are 5 Workgroups (8-wavefronts aka 512 workitems per workgroup), the 64kB LDS only needs to be split into 13107 chunks between the 5-workgroups, even at max occupancy 10.

• As a rule of thumb: bigger workgroups that share more data will more effectively use the LDS. However, not all workloads allow you to share data easily.

• The minimum workgroup size of 1 wavefront / 64-work items is treated as special. Barriers and synchronization never has to happen! Workgroup size of 1 wavefront (64-work items) by definition executes synchronously with itself. Still, use barrier instructions (and let the compiler figure out that it can turn barriers into no-ops).

• A secondary use of LDS is to use it as a manually managed cache. Don't feel bad if you do this: the LDS is faster than L1 cache.

• L1 vector data cache is 16kB, and slower than even LDS. In general, any program serious about speed will use the LDS explicitly, instead of relying upon the L1 cache. Still, its helpful to know that 16kB of global RAM will be cached for your CU.

• L1 scalar data cache is 16kB, shared between 4 CUs (!!). While this seems smaller than vector L1 Cache, remember that each sALU is running 64-threads / work items. In effect, the 40-wavefronts (x4 == 160 wavefronts max across 4 CUs) represents 10240 threads. But any sALU doesn't store data per-thread... it stores data per wavefront. Despite being small, this L1 scalar data cache can be quite useful in optimized code.

• Profile your code. While the theoretical discussion of this thread may be helpful to understanding why your GPGPU code is slow, you only truly understand performance if you read the hard-data.

• HBM2 Main Memory is very slow (~350 cycles to respond), and relatively low bandwidth ("only" 480 GBps). At Occupancy 1, there will be a total of 16384 workitems (or CUDA Threads) running on your Vega64. The 8GB of HBM2 main-memory can therefore be split up into 512kB.

As Bill Gates used to say, 640kB should be enough for everyone. Unfortunately, GPUs have such huge amounts of parallelism, you really can't even afford to dedicate that much RAM even in an "Occupancy 1" situation. The secret to GPUs is that your work-items will strongly share data with each other.

Yeah yeah yeah, GPUs are "embarassingly parallel", or at least are designed to work that way. But in practice, you MUST share data if you want to get things done. Even with "Occupancy 1", the 512kB of HBM2 RAM per work-item is too small to accomplish most embarassingly parallel tasks.

References

• AMD OpenCL Optimization Guide

• AMD GCN Crash Course

• Advanced Shader Programming on GCN

• GCN Assembly Tutorial -- Seeing the assembly helps understand how sGPR or vGPRs work, and solidify your "wavefront" ideas.

• Vega Assembly Language Manual -- 247 pages of dense, raw, assembly language.

This particular post is due to one curiosity: Raytracers seem to have a highly-divergent pattern: As soon as a "Ray" is emitted, it will either miss all targets, or it will hit a target.

However, even the "hit" targets can become a specular, diffuse, glossy, or subsurface-scattering shader associated with it. Which means that the ray will obviously cause thread-divergence.

Well, I guess I'm getting ahead of myself. Lets first describe raytracing, as well as thread-divergence for the beginners out there.

Raytracing

Raytracing at its core is an incredibly simple idea. In fact, you can implement a raytracer in under 100 lines of code.

A raytracer is a 3d rendering program. From a description of objects ("smallpt" only accepts spheres as its description), it builds a simulated 3d image of those objects. A typical "professional" raytracer would accept a list of triangles to render.

The way light bounces around a scene is modeled as a monte-carlo randomized algorithm. Some light bounced from over here, or over there (at random), to simulate the many paths that light photons can take between light sources and our eyes. Rays start at the camera (usually around 500 Rays per pixel), and randomly bounce around the scene to get an idea of where shadows or objects can be.

for every pixel "p" in the image:
    for i 0 to 500:
        launchRay(starting at pixel p, moving in random direction);

for every ray:
    if (ray hits object){
        bounceOffObject(); 
        // Different for each object: Typically, its color of ray + (launch a new ray, get its color).
        // The new ray is to figure out the color of the object you just hit: you don't know its color
        // because you don't know where the shadows are yet. 
    } else {
        renderBackground(); 
        // HDRi map, or maybe transparent background, or maybe a fog effect
    }

"BounceOffObject" can have many different algorithms: "Glossy" (or Specular) is an algorithm that works well on metals and shiny objects. "Diffuse" is an algorithm that works well on duller objects, like wood. There's subsurface scattering (which looks like gemstones and skin), as well as "Refraction", which works well on water and glass.

The "smallpt" Raytracer only supports Specular, Diffuse, and Refraction effects, as well as a default background of black on miss.

Note: Professional Raytracers get more complex because of the BVH Tree, a data-structure to accelerate the intersection calculation. As well as additional "bounce" algorithms (called BSDFs). Fundamentally though, the Raytracer at its core is incredibly simple. There's just more BSDFs in play, as well as more data-structures for optimization purposes. ("Smallpt" doesn't use a BVH Tree. It just does a many-to-many comparison inefficiently, lol).

Thread Divergence

I've talked about GPU architecture at a beginner level before, and I don't want to fully repeat myself.

But fundamentally, GPUs are SIMD machines. They execute groups of ~32 (NVidia) or 64 (AMD) threads at a time. What this means, is a piece of code like this:

if(foo()){
    doA(); // Lets say 10 threads want to do this
} else {
    doB(); // But 54 threads want to do this
}

Will have both branches executed by all 64 threads. What will happen under the hood is that "doA()" will have 54-threads set to "execution-mask = false", which means those threads will throw away the work.

Similarly, under "doB()", 10 threads will have execution-mask = false, and similarly throw away the work.

Given the issue, it is obvious that the following pattern is inefficient:

switch(ray.getHitType()){
    case Glossy:
        doGlossyReflect();
        break;
    case Diffuse:
        doDiffuseReflect();
        break;
    case SubSurface:
        doSubSurfaceReflect();
        break;
    case MISS:
        doHDRiMap(); // Default: does NOT issue a new ray. Other 3 "reflects" create a new Ray
        break;
}

Note: each of these "doSpecularShading()" statements issues a new ray (aka: thread / work item). So we clearly have a highly-divergent, dynamic work problem! Talk about ultimate thread divergence! None of these 4 cases will execute in parallel, and it is highly likely that all 64 of your rays in a raytracer will decide to do a different branch.

AMD ProRender

I decided to investigate the contradiction above. It is clear that the above "switch/case" statement is inefficient for GPUs. So how does raytracing happen exactly?

AMD has released an open-source GPU raytracer called ProRender, which seems to answer the question: https://github.com/GPUOpen-LibrariesAndSDKs/RadeonProRender-Baikal

In particular, this diagram explains the optimization we're most interested in: https://github.com/GPUOpen-LibrariesAndSDKs/RadeonProRender-Baikal/blob/master/Doc/Images/4.png?raw=true

First things first: there are many GPU Kernels that do little jobs. So instead of one kernel following all the bounces of a ray, there's one kernel that handles all the rays, but only for a part of the bounce calculation. A lot of the code is more CPU-driven than I personally expected. In fact, the entire loop uses maybe 10ish different GPU Kernels.

This has some major benefits:

1. This organization is necessary for the "Raytracer compaction algorithm" I'm going to describe later.

2. Smaller kernels can use more Registers / L1 cache throughout the calculation.

3. You can have more uniformity: All Rays are treated as rays. BSDFs like Glossy vs Diffuse, can have their own kernels to loop through.

In effect, the split-kernel design solves the thread-divergence issue. One kernel handles Glossy. Another kernel handles Diffuse. All of these kernels run in parallel with each other.

Code sample: "adaptive_renderer.cpp". Its a bit hard to tell, but pretty much every subfunction of "RenderTile" becomes an OpenCL kernel call eventually. GeneratePrimaryRays, AccumulateSamples, etc. etc. are all kernels.

"RenderTile" is effectively the entry point. After the first Rays are generated, the "Estimator" class takes over and handles the inner loop.

The overall pseudocode for AMD Prorender is:

AdaptiveRenderer::RenderTile{
    GenerateTileDomain // OpenCL Kernel. 
    GeneratePrimaryRays // Initial Rays
    Estimate() // inner-loop
    AccumulateSamples() // OpenCL Kernel: Adds results of Raytracing to tile
}

// This is the inner-loop of the raytracer
PathTracingEstimator::Estimate (){
    for(i from 0 to "MaxBounces"){
        FillBuffer(); // OpenCL: Initializes the raytracing data-structure
        QueryIntersection(); // OpenCL Kernel, check for intersections
        SampleVolume(); // Volume BSDF: Fog Effects and the like
        ShadeMiss(); // Render the missed rays with the Background
        FilterPathStream();
        Compact(); // FilterPathStream + Compact is the Raytracer Compaction I wanna talk about
        RestorePixelIndices(); // After compaction, move pixel indicies where they should be.
        ShadeSurface(); // Shade hits
        GatherLightSamples();
        Flush(); // Resync data with CPU.
    }
}

I've skipped a ton of steps, but you can see how each individual step is an OpenCL Kernel. A lot of these kernels can run concurrently. The more work you send to the GPU, the more efficiently it can typically perform.

Wait, where's the load-balancing "Raytracer Compaction" description??

Okay, I'm sorry it took this long. But I'm finally ready to talk about it. I think that Raytracer Compaction looks like a broadly-applicable optimization pattern, so I really want to highlight it specifically.

Look back at the ProRender description one more time, and focus on the line "Compact Ray Stream". There's the magic. This is how the programmer keeps 64-work items at 100% utilization on as many wavefronts as possible.

Lets take a wavefront of rays.

R0 | R1 | R2 | R3 | R4 | R5 | .. | R63
A  |A  |A  |A  |A  |A  | A  | .. | A

An "A" means the ray is active. After the "IntersectRays" code however, the many rays will miss the target and become inactive. While others are active.

R0 | R1 | R2 | R3 | R4 | R5 | .. | R63
A  |I  |I  |A  |I  |A  | A  | .. | I

I put an "I" when a ray becomes inactive. For "Real Code" that does this, look at "Intersect_bvh2_lds.cl".

            // Miss here
            hits[index].prim_id = MISS_MARKER;
            hits[index].shape_id = MISS_MARKER;

The "MISS_MARKER" is basically what "I" is. I'm saying R1 or R2 is no longer active. In any case, continuing to loop with this wavefront structure would be inefficient! We've lost half our rays, so we're doing half the work per wavefront!

Compaction itself is handled in CLWParallelPrimitives::Compact. Alas, the code is hard to read and mysterious.

NVidia fortunately had a great paper on the subject: http://www.cse.chalmers.se/~uffe/streamcompaction.pdf

In effect, the "PrefixSum" gives you the address of the compaction. So you run two kernels. The first kernel is a standard PrefixSum (Ex: R0 == 0. R3 ==1. Because R1 and R2 were inactive in the example). The 2nd kernel is a simple "Scatter", where you simply place R0 into 0. Then R3 into 1. Both operations are extremely efficient on a GPU.

As such, the "CLWParallelPrimitives::Compact" is the major lesson to be learned here! It manages to efficiently load-balance 1-dimensional GPU tasks, using the GPU itself! Both NVidia and AMD have used this pattern in various forms, but I figure its worth highlighting it in this post.

Since the CPU has a big loop on the outside of ~10ish GPU kernels, by the time the GPU loops back to the top of the loop, the GPU will handle all Rays in the Raytracer with nearly 100% utilization again. The combination of "split kernel" (~10 kernels in the inner loop, handling different tasks of the Raytracer) + the Stream Compaction algorithm, is what keeps things efficient.

Changelog

• Change 1 -- Added Pseudocode for AdaptiveRenderer::RenderTile and PathTracingEstimator::Estimate, to highlight the split-kernel design.

• Initial Post. Please note any errors, and I'll work to correct them.

Hey everyone,

I'm studying GPUs, but the more I study, the more I realize that this field has a LOT to offer. The SIMD world is small and obscure, but the papers, textbooks, and articles on the subject are often very high quality, with clear expertise in the methodology.

If anyone has a good book, article, or paper to share, please list it here!

My list:

GPU Gems Series

A collection of articles written from NVidia's GPUs covering 2004 to 2011, covering a time just before CUDA was popular to the very start of CUDA.

The articles vary in quality and detail, but overall are pretty good reads.

• GPU Gems 1 (2004)

• GPU Gems 2 (2005)

• GPU Gems 3 (2007)

• GPU Computing Gems Emerald Edition (2011)

• GPU Computing Gems Jade Edition (2011)

Shader X / GPU Pro / GPU Zen

This series of books edited by Wolfgang Engel covers multiple decades worth of SIMD and GPU programming.

• ShaderX covers material from 2002 to 2009.

• GPU Pro covers 2010 to 2016

• GPU Zen is the newest, published in 2017 and 2019 so far.

This huge series of books is listed here: https://www.realtimerendering.com/resources/shaderx/

But I guess I should copy/paste the book titles into this topic for good measure? I won't put all the dates or the full titles.

Like GPU Gems, the quality of each article varies. There's some high-level non-detailed stuff in here, but that's still useful for a quick discussion on some problems. Other articles lead into very in-depth analysis.

• Direct3D ShaderX (2002)

• Shader X2 Introductions and Tutorials with DirectX 9.0 (2003)

• Shader X2 Shader Programming Tips and Tricks with DirectX 9.0

• Shader X3

• Shader X4

• Shader X5

• Shader X6

• Shader X7 (2009)

• GPU Pro: Advanced Rendering Techniques (2010)

• GPU Pro2

• GPU Pro3

• GPU Pro4

• GPU Pro5

• GPU Pro6

• GPU Pro7 (2016)

• GPU Zen (2017)

• GPU Zen 2 (2019)

The "GPU Pro 360" books seem to collect the articles into subjects: one for Lighting, etc. etc. They hold the same information as the GPU Pro books, just by subject instead of by date published.

Vector Models for Data-Parallel Computing by Guy E. Blelloch

Blelloch's PH.D dissertation: https://www.cs.cmu.edu/~guyb/papers/Ble90.pdf

This is a deep dive into prefix-sum operations, using prefix-sum / prefix-max (and other prefix or scan operations) to solve a wide variety of problems.

Though written in 1990, the "Connection Machine" that Blelloch programs for is very similar to modern GPUs. As such, the PH.D Thesis remains surprisingly relevant in today's environment, especially as an introduction to the general power of a prefix-sum.

Technical Manuals

• Vega ISA (https://developer.amd.com/wp-content/resources/Vega_Shader_ISA_28July2017.pdf)

• RDNA ISA (https://developer.amd.com/wp-content/resources/RDNA_Shader_ISA.pdf)

• CUDA PTX (https://docs.nvidia.com/cuda/parallel-thread-execution/index.html)

The above assembly-level (or "near assembly") documents provide the lowest level building block to the modern GPU.

• Cuda Best Practices (https://docs.nvidia.com/cuda/cuda-c-best-practices-guide/index.html)

• AMD OpenCL Optimization guide (http://developer.amd.com/wordpress/media/2013/12/AMD_OpenCL_Programming_Optimization_Guide2.pdf).

Other Books

• Programming Massively Parallel Processors -- CUDA Specific

• EDIT: This webpage has an interesting list: http://www.realtimerendering.com/books.html

• ????

That's all for now. Please post your references!

TLDR: I need guidance for which framework to choose in 2024 (the most promising and vendor agnostic). Most posts related to that in this sub are at least 1 year old. Has something changed since then?

Hi guys, I'm a software engineer interested in HPC and I am completely lost trying to get back to GPGPU. I worked on a research project back in 2017/2018, and I went for OpenCL, as it was very appealing: a cross platform non-vendor specific framework that could run on almost everything. And yeah, it had a good Open Source support, specially from AMD. It sounded promising to me.

I was really excited about newer OpenCL releases, but I moved to other projects in which GPGPU weren't appliacable and lost the track of the framework evolution. Now I'm planning to develop some personal projects and dive deep on GPGPU again, but the ecosystem seems to be screwed up.

OpenCL seems to be diying. No vendor is currently suporting newer versions of the ones they were already supportting in 2017! I researched a bit about SYCL (bought Data Parallel C++ with SYCL book), but again, there is not a wide support or even many projects using SYCL. It also looks like an Intel thing. Vulcan is great, and I might be wrong, but I think it doesn't seem to be suitable for what I want (coding generic algorithms and run it on a GPU), despite it is surely cross platform and open.

It seems now that the only way is to choose a vendor and go for Metal (Apple), CUDA (NVIDIA), HIP (AMD) or SYCL (Intel). So I am basically going to have to write a different backend for every one of those, if I want to be vendor agnostic.

Is there a framework I might be missing? Where would you start in 2024? (considering you are aiming to write code that can run fast on any GPU)

We've just launched the alpha version of tellstory.ai which already is one of the fastest databases in the world. It's GPU accelerated - we use CUDA kernels for query processing. There is an opportunity to join our team at this early stage, if anyone is interested, check out the job ad: https://instarea.com/jobs/cuda-developer/

I'd like to quickly draw your attention to a project I've been working on for the past few months:

hipSYCL is an implementation of SYCL over NVIDIA CUDA/AMD HIP, targeting NVIDIA GPUs and AMD GPUs running ROCm.

It's still work in progress and there are parts of the SYCL specification that are still unimplemented, but it can already be used for many applications.

SYCL is an open standard describing a single-source C++ programming model for heterogeneous systems, originally intended to sit on top of OpenCL. The nice thing about SYCL is that it abstracts away the cumbersome parts (e.g. data migration between host and device and resource management), while still providing access to low-level optimizations such as explicit control over local memory (or shared memory in CUDA).

My SYCL implementation sits on top of HIP/CUDA instead of OpenCL. Since both CUDA (or HIP) and SYCL are single-source programming models based on C++, this allows for an implementation of SYCL as a CUDA library. The SYCL application can then be compiled with the regular CUDA or HIP compilers nvcc (for NVIDIA) and hcc (for AMD). This approach is the general idea behind hipSYCL. In practice, it's a bit more complicated though (there is actually an additional source-to-source transformation step in hipSYCL before code is fed into the NVIDIA/AMD compilers).

There are many advantages to this approach:

• You can write your applications against a vendor-netral, open standard (SYCL) while still being able to use e.g. the latest and greatest CUDA intrinsics or other platform specific optimizations when your SYCL code is compiled for NVIDIA devices. Anything that works in CUDA (or HIP) can in principle also be used with hipSYCL. But please use #ifdefs to remain portable :)

• All debuggers, profilers or other tools for CUDA or HIP also work with hipSYCL applications, since hipSYCL is effectively just another CUDA/HIP library

• Performance is on par with CUDA (or HIP) since the same device compiler is used

• The same code can run on a wide range of devices from CPUs to FPGAs when using the other available SYCL implementations triSYCL and ComputeCpp

• Compared to CUDA, SYCL is much more modern: No more __host__ and __device__ attributes, automatic resource management, out-of-order processing based on implicit task graphs instead of in-order queues and so on.

At the moment, the stage of the project is 'works for me'. If you try hipSYCL, I'd love to have some feedback about what works well, what doesn't work and what features you find most lacking. This helps me to better focus my efforts and to make hipSYCL more robust. Of course, pull requests are also always welcome :)

• Global memory barriers are very slow on GPUs, and can only be executed maybe once per microsecond or so (once every 1000ns). Any global data structure should have "Block local" buffers which only use CUDA-block (or OpenCL Threadgroup) level synchronization, which is far faster instead. In particular, AMD Vega64 seems to compile a global-threadfence into a L1 cache flush.

  • Synchronizing with the CPU (cudaDeviceSynchronize / hipDeviceSynchronize) seems to be only a little bit slower than thread-fences + spinlocks.

• RAM is used up ridiculously quickly. Any GPU will have the ability to run 10,000+ hardware threads. Vega64 should be run with 16384 hardware threads at a minimum for example (and supports up to 163,840 hardware threads at max occupancy). However, 16384 threads will run out of VRAM in just 512kB per thread: you don't even get the traditional "640kB" that should be enough for everyone.

  • Maximize the sharing of data between threads.

  • Because RAM needs to be used with utmost efficiency, you will end up writing your own data-structures rather often. In most cases, you'll use a simple array.

  • array[tail + __popc(writemask & __lanemask_lt())] = someItem; tail+= __prefix_sum(popc(writemask)) is an important pattern. This SIMD-stack paradigm should be your "bread-and-butter" collection due to its simplicity and efficiency. AMD/ROCm users can use __ockl_activelane_u32() to get the current active lane number.

  • SIMD-data structures are commonly "bigger" than the classic data-structures. Each "link" in a linked list should be the same size as the warpSize (32 on NVidia, 64 on AMD cards). Each node in a SIMD-Heap should also be 32+ or 64+ wide to support efficient SIMD-load/store.

• Debugging 10,000+ threads one-at-a-time doesn't really scale. Use the GPU to write GPU-tests per-thread, and then use the CPU to verify data sequentially. Especially if you are hunting threadfence or memory-barrier issues: the only way to catch a memory barrier issue is if you unleash as many threads as possible and run them as long as possible.