4 ms·

Wow, great answer. That deserves more upvotes than it's probably going to get. So what that means for compute shaders is I will need to partition the vectors

Wow, great answer. That deserves more upvotes than it's probably going to get. So what that means for compute shaders is I will need to partition the vectors being multiplied into power of 2 sized partitions and repeatedly call the kernels, but with offsets appropriate to the partition size, and then have the overrun contribute zero as you said. Obviously not as efficient as doing it all GPU side, still faster than doing it on the CPU. Thank you!

ryao 2y ago
You don’t need to repeatedly call the kernels. You just need to loop inside the kernels to have each process every nth multiplication according to their thread id, adding it to the thread’s partial sum and have a branch in the loop skip the addition (or add 0) when you exceed the range. Afterward, you would do a horizontal sum of the partial sums to get the value for the corresponding field in the output vector, and you would have a “grid” do this for all members of the output vector. The original poster’s mathmul kernel does this. See the example under “we can perform a warp sum reduction to combine results from all threads”: https://andrewkchan.dev/posts/yalm.html#section-3.2 https://andrewkchan.dev/posts/yalm.html#section-3.2 It does “if (i >= d) return;“ to handle non-multiples in one dimension and skips the addition in its dot product, matmul_row(), to handle non-multiples in the other dimension. In CUDA, you have threads that are organized into 32 thread warps, that are organized into blocks of warps, that are organized into grids. Threads within a block can share memory, but blocks within a grid cannot. The hardware will effectively loop over blocks in the grid until it finishes. The dimension for the output vector does not require sharing memory, so we can set that to the grid size. The hardware might actually round the grid up, so “if (i >= d) return;” is a check for the case that we are executing past the end of the output vector in case the hardware schedules pass the end of the vector. As for the other dimension, we only have a limited number of threads (in this example 32, which is a single warp), so we need each to do partial sums of the multiplications in the dot product. We essentially break the dot product into thread count (32 here) chunks. Each thread will know which multiplication in the thread count chunk to compute based on its thread id, and we keep adding the thread count (32 here) to that to get the next one until we have finished. The for loop’s loop condition protects us from going past the end of the dot product dimension. Then we have thread count number of partial sums, so we need to do a horizontal sum across the threads in a block to get the final sum. This part is a bit magical and is done via a warp sum reduction. Since the threads share memory, you basically have them start summing their neighbors’ values into their own. Basically, all threads find pairs and one of the pair does the sum, cutting the threads doing the sum in half. Then the ones that did the sums pair up again. This again cuts the number in half. This continues until only 1 thread remains. Then only the remaining thread is allowed to write to “out[i] = rowSum” due to “if (threadIdx.x == 0)”. The warp itself is somewhat special since it can do a really fast horizontal sum via shuffles that cannot be done when using thread blocks ghat utilize multiple warps. That is why the author chose to have the block size be exactly 1 warp rather than a multiple of warps. It helps to think about the threads in the block as a spreadsheet. The columns are the threads and the rows represent instructions. Every thread executes the same instruction, and when a branch says not to do something for only some threads, the ones excluded do noops until execution reaches a point where they all do the same thing. An if else similarly has noops in threads not executing true, followed by noops in the threads not executing false. This is called branch divergence. There are diagrams showing a spreadsheet like view of this, such as the following, which shows a cutdown example: https://image1.slideserve.com/2427671/branch-divergence-l.jpg https://image1.slideserve.com/2427671/branch-divergence-l.jp... I have mentioned a horizontal sum multiple times. That just means summing across the columns inside rows of the “spreadsheet”. It is called horizontal since rows are horizontal. Coincidentally, the partial sum calculation is called a vertical sum since the columns are vertical. Now we need to consider the grid, which is yet another dimension and is the dimension along which we populate the output vector. Basically, you can view the grid as doing multiple spreadsheets in parallel, but not all at once. One comical way to imagine it is to imagine having printed the same spreadsheet hundreds of times an and then feeding them into a paper shredder head first. A normal paper shredder can only do a few at a time, and similarly, the GPU can only do so many at a time too. You can imagine the processing to proceed along the rows in all of the spreadsheets being shredded at the same time.
- upghost 2y ago
  Damn, wow. Good stuff. First of all, I assumed looping within kernels was "always bad" from a performance standpoint (not that I'd done any benchmarks on looping vs repeated kernel invocation). I also assumed it was "always bad" to include conditionals. That certainly opens up a lot of possibilities on algorithms. It also explains why people think I'm crazy when I say I don't know how to multiply matrices that are not in power of 2 dimensions haha. Thank you, not sure why I was limiting myself like that. Not too many people to chat to about low level GPGPU stuff!
  - saagarjha 2y ago
    Loops and conditions within kernels are much better than the alternative, since synchronizing with the CPU is quite expensive. If you unroll appropriately and hide latencies they are very effective.
    - ryao 2y ago
      Sometimes, the compiler will unroll for you. I have noticed that nvcc will unroll this: #define warpSize 32 float warpMax = partial_max; for (int offset = warpSize / 2; offset > 0; offset >>= 1) { float otherValue = __shfl_down_sync(0xFFFFFFFF, warpMax, offset); warpMax = fmaxf(warpMax, otherValue); } Into this: bar.sync 0; shl.b32 %r40, %r2, 2; add.s32 %r8, %r39, %r40; ld.shared.f32 %f24, [%r8]; mov.u32 %r43, 0; mov.b32 %r44, %f24; shfl.sync.down.b32 %r48|%p13, %r44, %r23, %r22, %r24; mov.b32 %f25, %r48; max.ftz.f32 %f26, %f24, %f25; mov.b32 %r49, %f26; shfl.sync.down.b32 %r51|%p14, %r49, %r27, %r22, %r24; mov.b32 %f27, %r51; max.ftz.f32 %f28, %f26, %f27; mov.b32 %r52, %f28; shfl.sync.down.b32 %r54|%p15, %r52, %r30, %r22, %r24; mov.b32 %f29, %r54; max.ftz.f32 %f30, %f28, %f29; mov.b32 %r55, %f30; shfl.sync.down.b32 %r56|%p16, %r55, %r21, %r22, %r24; mov.b32 %f31, %r56; max.ftz.f32 %f32, %f30, %f31; mov.b32 %r57, %f32; shfl.sync.down.b32 %r59|%p17, %r57, %r35, %r22, %r24; mov.b32 %f33, %r59; max.ftz.f32 %f34, %f32, %f33; mov.b32 %r60, %f34; This avoids loop branching, which is unnecessary in this case because we know in advance when the loop terminates. The use of ftz instructions is because I compiled using `nvcc -ftz true -ptx -arch=sm_86 rung.cu -o rung.ptx`. I am not sure if unrolling loops makes sense on a GPU when we don't know when the loop will terminate in advance. If you can suggest any references explaining when that is helpful, I would be happy to read them. :)

ryao 2y ago

You don’t need to repeatedly call the kernels. You just need to loop inside the kernels to have each process every nth multiplication according to their thread id, adding it to the thread’s partial sum and have a branch in the loop skip the addition (or add 0) when you exceed the range. Afterward, you would do a horizontal sum of the partial sums to get the value for the corresponding field in the output vector, and you would have a “grid” do this for all members of the output vector. The original poster’s mathmul kernel does this. See the example under “we can perform a warp sum reduction to combine results from all threads”: https://andrewkchan.dev/posts/yalm.html#section-3.2 https://andrewkchan.dev/posts/yalm.html#section-3.2 It does “if (i >= d) return;“ to handle non-multiples in one dimension and skips the addition in its dot product, matmul_row(), to handle non-multiples in the other dimension. In CUDA, you have threads that are organized into 32 thread warps, that are organized into blocks of warps, that are organized into grids. Threads within a block can share memory, but blocks within a grid cannot. The hardware will effectively loop over blocks in the grid until it finishes. The dimension for the output vector does not require sharing memory, so we can set that to the grid size. The hardware might actually round the grid up, so “if (i >= d) return;” is a check for the case that we are executing past the end of the output vector in case the hardware schedules pass the end of the vector. As for the other dimension, we only have a limited number of threads (in this example 32, which is a single warp), so we need each to do partial sums of the multiplications in the dot product. We essentially break the dot product into thread count (32 here) chunks. Each thread will know which multiplication in the thread count chunk to compute based on its thread id, and we keep adding the thread count (32 here) to that to get the next one until we have finished. The for loop’s loop condition protects us from going past the end of the dot product dimension. Then we have thread count number of partial sums, so we need to do a horizontal sum across the threads in a block to get the final sum. This part is a bit magical and is done via a warp sum reduction. Since the threads share memory, you basically have them start summing their neighbors’ values into their own. Basically, all threads find pairs and one of the pair does the sum, cutting the threads doing the sum in half. Then the ones that did the sums pair up again. This again cuts the number in half. This continues until only 1 thread remains. Then only the remaining thread is allowed to write to “out[i] = rowSum” due to “if (threadIdx.x == 0)”. The warp itself is somewhat special since it can do a really fast horizontal sum via shuffles that cannot be done when using thread blocks ghat utilize multiple warps. That is why the author chose to have the block size be exactly 1 warp rather than a multiple of warps. It helps to think about the threads in the block as a spreadsheet. The columns are the threads and the rows represent instructions. Every thread executes the same instruction, and when a branch says not to do something for only some threads, the ones excluded do noops until execution reaches a point where they all do the same thing. An if else similarly has noops in threads not executing true, followed by noops in the threads not executing false. This is called branch divergence. There are diagrams showing a spreadsheet like view of this, such as the following, which shows a cutdown example: https://image1.slideserve.com/2427671/branch-divergence-l.jpg https://image1.slideserve.com/2427671/branch-divergence-l.jp... I have mentioned a horizontal sum multiple times. That just means summing across the columns inside rows of the “spreadsheet”. It is called horizontal since rows are horizontal. Coincidentally, the partial sum calculation is called a vertical sum since the columns are vertical. Now we need to consider the grid, which is yet another dimension and is the dimension along which we populate the output vector. Basically, you can view the grid as doing multiple spreadsheets in parallel, but not all at once. One comical way to imagine it is to imagine having printed the same spreadsheet hundreds of times an and then feeding them into a paper shredder head first. A normal paper shredder can only do a few at a time, and similarly, the GPU can only do so many at a time too. You can imagine the processing to proceed along the rows in all of the spreadsheets being shredded at the same time.

upghost 2y ago
Damn, wow. Good stuff. First of all, I assumed looping within kernels was "always bad" from a performance standpoint (not that I'd done any benchmarks on looping vs repeated kernel invocation). I also assumed it was "always bad" to include conditionals. That certainly opens up a lot of possibilities on algorithms. It also explains why people think I'm crazy when I say I don't know how to multiply matrices that are not in power of 2 dimensions haha. Thank you, not sure why I was limiting myself like that. Not too many people to chat to about low level GPGPU stuff!
- saagarjha 2y ago
  Loops and conditions within kernels are much better than the alternative, since synchronizing with the CPU is quite expensive. If you unroll appropriately and hide latencies they are very effective.
  - ryao 2y ago
    Sometimes, the compiler will unroll for you. I have noticed that nvcc will unroll this: #define warpSize 32 float warpMax = partial_max; for (int offset = warpSize / 2; offset > 0; offset >>= 1) { float otherValue = __shfl_down_sync(0xFFFFFFFF, warpMax, offset); warpMax = fmaxf(warpMax, otherValue); } Into this: bar.sync 0; shl.b32 %r40, %r2, 2; add.s32 %r8, %r39, %r40; ld.shared.f32 %f24, [%r8]; mov.u32 %r43, 0; mov.b32 %r44, %f24; shfl.sync.down.b32 %r48|%p13, %r44, %r23, %r22, %r24; mov.b32 %f25, %r48; max.ftz.f32 %f26, %f24, %f25; mov.b32 %r49, %f26; shfl.sync.down.b32 %r51|%p14, %r49, %r27, %r22, %r24; mov.b32 %f27, %r51; max.ftz.f32 %f28, %f26, %f27; mov.b32 %r52, %f28; shfl.sync.down.b32 %r54|%p15, %r52, %r30, %r22, %r24; mov.b32 %f29, %r54; max.ftz.f32 %f30, %f28, %f29; mov.b32 %r55, %f30; shfl.sync.down.b32 %r56|%p16, %r55, %r21, %r22, %r24; mov.b32 %f31, %r56; max.ftz.f32 %f32, %f30, %f31; mov.b32 %r57, %f32; shfl.sync.down.b32 %r59|%p17, %r57, %r35, %r22, %r24; mov.b32 %f33, %r59; max.ftz.f32 %f34, %f32, %f33; mov.b32 %r60, %f34; This avoids loop branching, which is unnecessary in this case because we know in advance when the loop terminates. The use of ftz instructions is because I compiled using `nvcc -ftz true -ptx -arch=sm_86 rung.cu -o rung.ptx`. I am not sure if unrolling loops makes sense on a GPU when we don't know when the loop will terminate in advance. If you can suggest any references explaining when that is helpful, I would be happy to read them. :)