|
| 1 | +#version 450 |
| 2 | +#extension GL_EXT_control_flow_attributes : enable |
| 3 | +#extension GL_EXT_debug_printf : enable |
| 4 | +#extension GL_KHR_shader_subgroup_basic : enable |
| 5 | +#extension GL_KHR_shader_subgroup_ballot : enable |
| 6 | +#extension GL_KHR_shader_subgroup_arithmetic : enable |
| 7 | +#extension GL_KHR_shader_subgroup_shuffle : enable |
| 8 | + |
| 9 | +#include "types.glsl" |
| 10 | + |
| 11 | +layout(constant_id = 0) const int BLOCK_SIZE = 1024; |
| 12 | +layout(constant_id = 1) const int SUBGROUP_SIZE = 32; |
| 13 | +layout(constant_id = 2) const int SUBGROUP_SIZE_LOG2 = 5; |
| 14 | + |
| 15 | +layout(local_size_x_id = 0, local_size_y = 1, local_size_z = 1) in; |
| 16 | + |
| 17 | +// Input can either be the source (A) or intermediate values (S). |
| 18 | +// Similarly, output can be either destination (D) or intermediate values (S). |
| 19 | +layout (binding = 0) readonly buffer A {A_TYPE data_a[];}; |
| 20 | +layout (binding = 0) readonly buffer S {ivec2 data_s[];}; |
| 21 | +layout (binding = 1) writeonly buffer D {int data_d[];}; |
| 22 | +layout (binding = 1) writeonly buffer T {ivec2 data_t[];}; |
| 23 | + |
| 24 | +layout (push_constant) uniform parameter { |
| 25 | + uint orig_ncols; |
| 26 | + uint ncols_input; |
| 27 | + uint ncols_output; |
| 28 | + uint nrows; |
| 29 | + uint first_pass; |
| 30 | + uint last_pass; |
| 31 | +} p; |
| 32 | + |
| 33 | +// pairs of (gid, value) |
| 34 | +shared ivec2 dst_row[BLOCK_SIZE]; |
| 35 | + |
| 36 | +shared int counts[SUBGROUP_SIZE]; |
| 37 | +shared int sh_min_idx; |
| 38 | +shared uint sh_total; |
| 39 | +shared uint offset_partials[BLOCK_SIZE / SUBGROUP_SIZE]; |
| 40 | + |
| 41 | +// Map float values to uint such that comparisons still work. |
| 42 | +// Positive values set the high bit, negative values are inverted. |
| 43 | +// +0.0 -> 0x80000000, -0.0 -> 0x7FFFFFFF are in the correct places. |
| 44 | +uint f2ui(float x) { |
| 45 | + uint y = floatBitsToUint(x); |
| 46 | + if ((y & 0x80000000) != 0) { |
| 47 | + y ^= ~0; |
| 48 | + } else { |
| 49 | + y |= 0x80000000; |
| 50 | + } |
| 51 | + return y; |
| 52 | +} |
| 53 | + |
| 54 | +void topk(const uint row) { |
| 55 | + const int tid = int(gl_LocalInvocationID.x); |
| 56 | + |
| 57 | + // initialize indices |
| 58 | + if (gl_GlobalInvocationID.x < p.ncols_input) { |
| 59 | + if (p.first_pass != 0) { |
| 60 | + const uint row_offset = row * p.ncols_input; |
| 61 | + dst_row[tid] = ivec2(gl_GlobalInvocationID.x, floatBitsToInt(data_a[row_offset + gl_GlobalInvocationID.x])); |
| 62 | + } else { |
| 63 | + const uint row_offset = row * p.orig_ncols; |
| 64 | + dst_row[tid] = data_s[row_offset + gl_GlobalInvocationID.x]; |
| 65 | + } |
| 66 | + } else { |
| 67 | + dst_row[tid] = ivec2(p.orig_ncols, 0xFF800000); // -inf |
| 68 | + } |
| 69 | + barrier(); |
| 70 | + |
| 71 | + if (p.ncols_output == 1) { |
| 72 | + // Fast path for single output - just do a max reduction |
| 73 | + [[unroll]] for (int s = BLOCK_SIZE / 2; s >= 1; s /= 2) { |
| 74 | + if (tid < s) { |
| 75 | + ivec2 a = dst_row[tid]; |
| 76 | + ivec2 b = dst_row[tid + s]; |
| 77 | + if (a.x >= p.orig_ncols || |
| 78 | + b.x < p.orig_ncols && b.y > a.y) { |
| 79 | + dst_row[tid] = b; |
| 80 | + } |
| 81 | + } |
| 82 | + barrier(); |
| 83 | + } |
| 84 | + } else { |
| 85 | + // Do an N-ary search to find the K-th largest value. |
| 86 | + // We remap the float values to be comparable as unsigned integers, |
| 87 | + // and split the range into 2^N smaller ranges where N is the |
| 88 | + // subgroup size. Count how many values are in each range, if the K-th |
| 89 | + // largest value is in the middle of one of thee ranges then repeat |
| 90 | + // and split again. |
| 91 | + |
| 92 | + // Mask is the current set of bits we're searching. Shift is the LSB index. |
| 93 | + int shift = 32 - SUBGROUP_SIZE_LOG2; |
| 94 | + uint mask = ((1 << SUBGROUP_SIZE_LOG2) - 1) << shift; |
| 95 | + |
| 96 | + // The current range. |
| 97 | + uint range_min = 0; |
| 98 | + uint range_max = 0xFF800000; |
| 99 | + // How many are above the current range, and how many we need to find. |
| 100 | + uint total = 0; |
| 101 | + uint limit = min(p.ncols_output, p.ncols_input - gl_WorkGroupID.x * BLOCK_SIZE); |
| 102 | + |
| 103 | + while (mask != 0) { |
| 104 | + barrier(); |
| 105 | + // Initialize bucket counts to zero. |
| 106 | + if (tid < SUBGROUP_SIZE) { |
| 107 | + counts[tid] = 0; |
| 108 | + } |
| 109 | + barrier(); |
| 110 | + // Count how many values are in each bucket. |
| 111 | + if (tid < p.ncols_input) { |
| 112 | + float y = intBitsToFloat(dst_row[tid].y); |
| 113 | + uint fy = f2ui(y); |
| 114 | + if (fy >= range_min && fy < range_max) { |
| 115 | + uint bucket = (fy & mask) >> shift; |
| 116 | + atomicAdd(counts[bucket], 1); |
| 117 | + } |
| 118 | + } |
| 119 | + barrier(); |
| 120 | + |
| 121 | + // On the first subgroup, do a scan to count (from the top down) how |
| 122 | + // many elements are in the top N buckets. Find the index of the first |
| 123 | + // that is over the limit. Copy it to the other invocations through |
| 124 | + // shared memory. |
| 125 | + if (tid < SUBGROUP_SIZE) { |
| 126 | + uint partial_sum = counts[SUBGROUP_SIZE - 1 - tid]; |
| 127 | + partial_sum = subgroupInclusiveAdd(partial_sum) + total; |
| 128 | + uint t = subgroupBallotFindLSB(subgroupBallot(partial_sum >= limit)); |
| 129 | + int min_idx = int(SUBGROUP_SIZE - 1 - t); |
| 130 | + total = subgroupShuffle(partial_sum, t); |
| 131 | + sh_min_idx = min_idx; |
| 132 | + sh_total = total; |
| 133 | + } |
| 134 | + barrier(); |
| 135 | + int min_idx = sh_min_idx; |
| 136 | + total = sh_total; |
| 137 | + |
| 138 | + // Update the range, and break if we've found the K-th largest. |
| 139 | + if (min_idx != SUBGROUP_SIZE - 1) { |
| 140 | + range_max = range_min + ((min_idx + 1) << shift); |
| 141 | + } |
| 142 | + range_min = range_min + (min_idx << shift); |
| 143 | + |
| 144 | + if (total == p.ncols_output) { |
| 145 | + break; |
| 146 | + } |
| 147 | + total -= counts[min_idx]; |
| 148 | + mask >>= SUBGROUP_SIZE_LOG2; |
| 149 | + shift -= SUBGROUP_SIZE_LOG2; |
| 150 | + if (shift < 0) { |
| 151 | + shift = 0; |
| 152 | + } |
| 153 | + } |
| 154 | + |
| 155 | + ivec2 v = dst_row[tid]; |
| 156 | + |
| 157 | + // We need to compact these values to the start of the dst_row array. |
| 158 | + // Have each subgroup count how many items it'll store, so other |
| 159 | + // subgroups can compute their base offset. |
| 160 | + bool top = f2ui(intBitsToFloat(v.y)) >= range_min; |
| 161 | + uvec4 b = subgroupBallot(top); |
| 162 | + uint bit_count = subgroupBallotBitCount(b); |
| 163 | + if ((tid % SUBGROUP_SIZE) == 0) { |
| 164 | + offset_partials[tid / SUBGROUP_SIZE] = bit_count; |
| 165 | + } |
| 166 | + barrier(); |
| 167 | + |
| 168 | + uint out_idx = 0; |
| 169 | + [[unroll]] for (int i = 0; i < BLOCK_SIZE / SUBGROUP_SIZE; ++i) { |
| 170 | + if (i < tid / SUBGROUP_SIZE) { |
| 171 | + out_idx += offset_partials[i]; |
| 172 | + } |
| 173 | + } |
| 174 | + |
| 175 | + uint bit_count_ex = subgroupBallotExclusiveBitCount(b); |
| 176 | + if (top) { |
| 177 | + // TODO: Copy directly to the output? |
| 178 | + dst_row[out_idx + bit_count_ex] = v; |
| 179 | + } |
| 180 | + |
| 181 | + barrier(); |
| 182 | + } |
| 183 | + |
| 184 | + if (tid < p.ncols_output && gl_GlobalInvocationID.x < p.orig_ncols) { |
| 185 | + if (p.last_pass != 0) { |
| 186 | + const uint row_offset = row * p.ncols_output; |
| 187 | + data_d[row_offset + tid] = dst_row[tid].x; |
| 188 | + } else { |
| 189 | + const uint row_offset = row * p.orig_ncols + gl_WorkGroupID.x * p.ncols_output; |
| 190 | + data_t[row_offset + tid] = dst_row[tid]; |
| 191 | + } |
| 192 | + } |
| 193 | +} |
| 194 | + |
| 195 | +void main() { |
| 196 | + uint row = gl_WorkGroupID.y; |
| 197 | + while (row < p.nrows) { |
| 198 | + topk(row); |
| 199 | + row += gl_WorkGroupSize.y * gl_NumWorkGroups.y; |
| 200 | + } |
| 201 | +} |
0 commit comments