char-fork/llama.cpp
git clone https://git.t4t.associates/char-fork/llama.cpp
c9ca51c1f
char/aliceai-t5-moe
1#version 450 2 3#extension GL_EXT_control_flow_attributes : enable 4#ifdef COOPMAT2 5#extension GL_NV_cooperative_matrix2 : enable 6#extension GL_EXT_shader_explicit_arithmetic_types_float16 : require 7#extension GL_KHR_memory_scope_semantics : enable 8#endif 9 10#ifdef COOPMAT 11#extension GL_KHR_cooperative_matrix : enable 12#extension GL_KHR_shader_subgroup_basic : enable 13#extension GL_EXT_shader_explicit_arithmetic_types_float16 : require 14#extension GL_KHR_memory_scope_semantics : enable 15#endif 16 17#include "types.glsl" 18#include "utils.glsl" 19 20// shape notation: [dim(N), ..., dim(0)] -- stride(dim(j)) >= stride(dim(i)) if i > j 21layout(binding = 0) readonly buffer A { 22 A_TYPE knl_data[]; 23}; // src0 - kernel: [KW, KH, KD, IC*OC] 24 25layout(binding = 1) readonly buffer B { 26 B_TYPE src_data[]; 27}; // src1 - input: [IW, IH, ID, IC*N] -- channel_first format 28 29layout(binding = 2) writeonly buffer D { 30 D_TYPE dst_data[]; 31}; // dst - result: [OW, OH, OD, OC*N] 32 33layout(push_constant) uniform parameter { 34 // I/O channels, batch size 35 uint32_t OC; 36 uint32_t IC; 37 uint32_t N; 38 39 // Tensor spatial sizes: input, output 40 uint32_t IW; 41 uint32_t IH; 42 uint32_t ID; 43 uint32_t OW; 44 uint32_t OH; 45 uint32_t OD; 46 47 // Strides in elements 48 uint32_t nb01; 49 uint32_t nb02; 50 uint32_t nb03; 51 52 uint32_t nb11; 53 uint32_t nb12; 54 uint32_t nb13; 55 56 uint32_t nb1; 57 uint32_t nb2; 58 uint32_t nb3; 59 60 // fastdiv helper values 61 uint32_t OWmp; uint32_t OWL; 62 uint32_t OWOHmp; uint32_t OWOHL; 63 uint32_t OWOHODmp; uint32_t OWOHODL; 64} 65 66p; 67 68layout(local_size_x_id = 0, local_size_y = 1, local_size_z = 1) in; 69// Blocktile sizes 70layout(constant_id = 1) const uint BS_K = 128; 71layout(constant_id = 2) const uint BS_CRS = 16; 72layout(constant_id = 3) const uint BS_NPQ = 128; 73// Thread-tile sizes 74layout(constant_id = 4) const uint TS_K = 8; 75layout(constant_id = 5) const uint SHMEM_PAD = 4; 76// Stride, padding, dilation 77layout(constant_id = 6) const uint s0 = 1; 78layout(constant_id = 7) const uint s1 = 1; 79layout(constant_id = 8) const uint s2 = 1; 80layout(constant_id = 9) const uint p0 = 0; 81layout(constant_id = 10) const uint p1 = 0; 82layout(constant_id = 11) const uint p2 = 0; 83layout(constant_id = 12) const uint d0 = 1; 84layout(constant_id = 13) const uint d1 = 1; 85layout(constant_id = 14) const uint d2 = 1; 86// Kernel spatial sizes 87layout(constant_id = 15) const uint KW = 1; 88layout(constant_id = 16) const uint KH = 1; 89layout(constant_id = 17) const uint KD = 1; 90// when set, skip bounds checks and address clamps (K/CRS/NPQ are tile-aligned) 91layout(constant_id = 18) const uint aligned = 0; 92// stage cm2 result through shmem (Csh) for coalesced stores. cm1 always does this. 93layout(constant_id = 19) const uint csh_store = 0; 94 95#ifdef COOPMAT 96// cm1 subgroup tile: each subgroup computes a WM x WN region as a grid of 97// TM x TN x TK fragments. Requires WM%TM == WN%TN == BS_K%WM == BS_NPQ%WN == 98// BS_CRS%TK == 0, and WG_SIZE == (BS_K/WM) * (BS_NPQ/WN) * subgroup_size. 99layout(constant_id = 20) const uint WM = 32; 100layout(constant_id = 21) const uint WN = 32; 101const uint TM = 16; 102const uint TN = 16; 103const uint TK = 16; 104const uint cms_per_row = WM / TM; 105const uint cms_per_col = WN / TN; 106const uint warps_M = BS_K / WM; 107const uint warps_N = BS_NPQ / WN; 108#endif 109 110// without padding, ID_idx/IH_idx/IW_idx are in bounds by construction 111const bool dhw_in_bounds = (p0 == 0) && (p1 == 0) && (p2 == 0); 112 113uint32_t tid = gl_LocalInvocationID.x; 114const uint32_t WG_SIZE = gl_WorkGroupSize.x; 115 116uint splitWork(uint work_size, uint block_size) { 117 return (block_size + work_size - 1) / block_size; 118} 119 120uint32_t K = p.OC; 121uint32_t CRS = p.IC * KD * KH * KW; 122uint32_t NPQ = p.N * p.OD * p.OH * p.OW; 123 124// Number of blocktiles per input 125uint32_t NB_CRS = splitWork(CRS, BS_CRS); 126 127#if defined(COOPMAT2) || defined(COOPMAT) 128#define SHMEM_TYPE float16_t 129#else 130#define SHMEM_TYPE float 131#endif 132 133const uint32_t Ash_stride = BS_CRS + SHMEM_PAD; 134const uint32_t Bsh_stride = BS_NPQ + SHMEM_PAD; 135 136const uint32_t Ash_len = BS_K * Ash_stride; 137const uint32_t Bsh_len = BS_CRS * Bsh_stride; 138 139shared SHMEM_TYPE Ash[Ash_len]; // K x CRS 140shared SHMEM_TYPE Bsh[Bsh_len]; // CRS x NPQ 141 142#if defined(COOPMAT2) || defined(COOPMAT) 143// stage matC through shmem so global stores are row-major (NPQ-contiguous) 144const uint32_t Csh_stride = BS_NPQ; 145#ifdef COOPMAT 146const uint32_t Csh_len = BS_K * Csh_stride; 147#else 148const uint32_t Csh_len = csh_store != 0 ? BS_K * Csh_stride : 8; // 8 to workaround compiler bug 149#endif 150shared SHMEM_TYPE Csh[Csh_len]; // K x NPQ 151#endif 152 153// Threadtile sizes 154const uint32_t TS_NPQ = BS_K * BS_NPQ / WG_SIZE / TS_K; 155 156// Number of threadtiles per blocktile 157const uint32_t NT_NPQ = BS_NPQ / TS_NPQ; 158 159/* 160Compute 161KxCRS @ CRSxNPQ = K x NPQ 162K=OC 163C=IC 164D,R,S=KD,KH,KW 165Z,P,Q=OD,OH,OW 166*/ 167 168uint32_t B_idx_K = gl_WorkGroupID.x; 169uint32_t B_idx_NPQ = gl_WorkGroupID.y + gl_WorkGroupID.z * 512; 170 171uint32_t T_y = tid / NT_NPQ; 172uint32_t T_x = tid % NT_NPQ; 173 174uint32_t Ar = tid / BS_CRS; 175uint32_t Ac = tid % BS_CRS; 176const uint32_t ArpWg = WG_SIZE / BS_CRS; 177 178uint32_t Br = tid / BS_NPQ; 179uint32_t Bc = tid % BS_NPQ; 180const uint32_t BrpWg = WG_SIZE / BS_NPQ; 181 182void split_crs(uint32_t crs_idx, out uint32_t ic, out uint32_t kd, out uint32_t kh, out uint32_t kw) { 183 const uint32_t KHKW = KH * KW; 184 const uint32_t KDKHKW = KD * KHKW; 185 ic = crs_idx / KDKHKW; 186 uint32_t rem = crs_idx - ic * KDKHKW; 187 kd = rem / KHKW; 188 rem = rem - kd * KHKW; 189 kh = rem / KW; 190 kw = rem - kh * KW; 191} 192 193void split_npq(uint32_t npq_idx, out uint32_t n, out uint32_t od, out uint32_t oh, out uint32_t ow) { 194 const uint32_t OWOH = p.OW * p.OH; 195 n = fastdiv(npq_idx, p.OWOHODmp, p.OWOHODL); 196 uint32_t rem = npq_idx - n * p.OD * OWOH; 197 od = fastdiv(rem, p.OWOHmp, p.OWOHL); 198 rem = rem - od * OWOH; 199 oh = fastdiv(rem, p.OWmp, p.OWL); 200 ow = rem - oh * p.OW; 201} 202 203#ifdef COOPMAT2 204#define ACC_TYPE float16_t 205 206ACC_TYPE perElemOpStore(const in uint32_t r, const in uint32_t c, const in ACC_TYPE elem) 207{ 208 uint32_t K_idx = B_idx_K * BS_K + r; 209 uint32_t NPQ_idx = B_idx_NPQ * BS_NPQ + c; 210 uint32_t N_idx; 211 uint32_t OD_idx; 212 uint32_t OH_idx; 213 uint32_t OW_idx; 214 split_npq(NPQ_idx, N_idx, OD_idx, OH_idx, OW_idx); 215 uint32_t dst_idx = OW_idx + OH_idx * p.nb1 + OD_idx * p.nb2 + (N_idx * p.OC + K_idx) * p.nb3; 216 if (aligned != 0 || (K_idx < K && NPQ_idx < NPQ)) { 217 dst_data[dst_idx] = D_TYPE(elem); 218 } 219 return elem; 220} 221#endif 222 223void main() { 224 if (B_idx_NPQ * BS_NPQ >= NPQ) { 225 return; 226 } 227 228#ifdef COOPMAT2 229 coopmat<ACC_TYPE, gl_ScopeWorkgroup, BS_K, BS_NPQ, gl_MatrixUseAccumulator> matC; 230 matC = coopmat<ACC_TYPE, gl_ScopeWorkgroup, BS_K, BS_NPQ, gl_MatrixUseAccumulator>(0.0); 231#elif defined(COOPMAT) 232 coopmat<float16_t, gl_ScopeSubgroup, TM, TN, gl_MatrixUseAccumulator> sums[cms_per_row * cms_per_col]; 233 [[unroll]] for (uint i = 0; i < cms_per_row * cms_per_col; i++) { 234 sums[i] = coopmat<float16_t, gl_ScopeSubgroup, TM, TN, gl_MatrixUseAccumulator>(0.0); 235 } 236 const uint warp_r = gl_SubgroupID / warps_N; 237 const uint warp_c = gl_SubgroupID % warps_N; 238#else 239 float regC[TS_K][TS_NPQ]; 240 for (uint32_t T_ly = 0; T_ly < TS_K; T_ly++) { 241 for (uint32_t T_lx = 0; T_lx < TS_NPQ; T_lx++) { 242 regC[T_ly][T_lx] = 0.0; 243 } 244 } 245#endif 246 /* Advance block in CRS dim */ 247 [[dont_unroll]] for (uint32_t B_idx_CRS = 0; B_idx_CRS < NB_CRS; B_idx_CRS++) { 248 uint32_t CRS_idx_a = B_idx_CRS * BS_CRS + Ac; 249 uint32_t IC_idx_a; 250 uint32_t KD_idx_a; 251 uint32_t KH_idx_a; 252 uint32_t KW_idx_a; 253 split_crs(CRS_idx_a, IC_idx_a, KD_idx_a, KH_idx_a, KW_idx_a); 254 255 /* Load kernel to A_block: (BS_K x BS_CRS)*/ 256 UNROLL for (uint32_t r_offset = 0; r_offset < BS_K; r_offset += ArpWg) { 257 uint32_t B_ly = r_offset + Ar; 258 uint32_t B_lx = Ac; 259 uint32_t K_idx = B_idx_K * BS_K + B_ly; /* Global K_idx (row index of A)*/ 260 uint32_t knl_idx = KW_idx_a + KH_idx_a * p.nb01 + KD_idx_a * p.nb02 + (K_idx * p.IC + IC_idx_a) * p.nb03; 261 if (aligned == 0) { 262 knl_idx = min(knl_idx, K * CRS - 1); 263 } 264 float val = knl_data[knl_idx]; 265 if (aligned == 0 && (K_idx >= K || CRS_idx_a >= CRS)) { 266 val = 0.0; 267 } 268 Ash[B_ly * Ash_stride + B_lx] = SHMEM_TYPE(val); 269 } 270 /* Load input to B_block: (BS_CRS x BS_NPQ) */ 271 UNROLL for (uint32_t r_offset = 0; r_offset < BS_CRS; r_offset += BrpWg) { 272 uint32_t B_ly = r_offset + Br; /* Row index of B block */ 273 uint32_t B_lx = Bc; 274 uint32_t NPQ_idx = B_idx_NPQ * BS_NPQ + B_lx; /* Global NPQ index (column index of B) */ 275 uint32_t N_idx; 276 uint32_t OD_idx; 277 uint32_t OH_idx; 278 uint32_t OW_idx; 279 split_npq(NPQ_idx, N_idx, OD_idx, OH_idx, OW_idx); 280 281 uint32_t CRS_idx_b = B_idx_CRS * BS_CRS + B_ly; 282 uint32_t IC_idx_b; 283 uint32_t KD_idx_b; 284 uint32_t KH_idx_b; 285 uint32_t KW_idx_b; 286 split_crs(CRS_idx_b, IC_idx_b, KD_idx_b, KH_idx_b, KW_idx_b); 287 288 uint32_t ID_idx = OD_idx * s2 + KD_idx_b * d2 - p2; 289 uint32_t IH_idx = OH_idx * s1 + KH_idx_b * d1 - p1; 290 uint32_t IW_idx = OW_idx * s0 + KW_idx_b * d0 - p0; 291 292 uint32_t src_idx = IW_idx + IH_idx * p.nb11 + ID_idx * p.nb12 + (N_idx * p.IC + IC_idx_b) * p.nb13; 293 // skip clamp when address can't go OOB 294 if (aligned == 0 || !dhw_in_bounds) { 295 src_idx = min(src_idx, p.IC * p.N * p.IW * p.IH * p.ID - 1); 296 } 297 float val = src_data[src_idx]; 298 bool oob = false; 299 if (aligned == 0 && (CRS_idx_b >= CRS || NPQ_idx >= NPQ)) { 300 oob = true; 301 } 302 // also catches lower-bound underflow (idx wraps to 0x80000000+) 303 if (!dhw_in_bounds && (ID_idx >= p.ID || IH_idx >= p.IH || IW_idx >= p.IW)) { 304 oob = true; 305 } 306 if (oob) { 307 val = 0.0; 308 } 309 Bsh[B_ly * Bsh_stride + B_lx] = SHMEM_TYPE(val); 310 } 311 barrier(); 312#ifdef COOPMAT2 313 coopmat<float16_t, gl_ScopeWorkgroup, BS_K, BS_CRS, gl_MatrixUseA> matA; 314 coopmat<float16_t, gl_ScopeWorkgroup, BS_CRS, BS_NPQ, gl_MatrixUseB> matB; 315 316 coopMatLoad(matA, Ash, 0, Ash_stride, gl_CooperativeMatrixLayoutRowMajor); 317 coopMatLoad(matB, Bsh, 0, Bsh_stride, gl_CooperativeMatrixLayoutRowMajor); 318 matC = coopMatMulAdd(matA, matB, matC); 319#elif defined(COOPMAT) 320 // each subgroup multiplies its grid of fragments per TK-sized CRS chunk 321 [[unroll]] for (uint k_step = 0; k_step < BS_CRS / TK; k_step++) { 322 coopmat<float16_t, gl_ScopeSubgroup, TM, TK, gl_MatrixUseA> cache_a[cms_per_row]; 323 [[unroll]] for (uint cm_row = 0; cm_row < cms_per_row; cm_row++) { 324 const uint a_off = (warp_r * WM + cm_row * TM) * Ash_stride + k_step * TK; 325 coopMatLoad(cache_a[cm_row], Ash, a_off, Ash_stride, gl_CooperativeMatrixLayoutRowMajor); 326 } 327 [[unroll]] for (uint cm_col = 0; cm_col < cms_per_col; cm_col++) { 328 coopmat<float16_t, gl_ScopeSubgroup, TK, TN, gl_MatrixUseB> cache_b; 329 const uint b_off = k_step * TK * Bsh_stride + warp_c * WN + cm_col * TN; 330 coopMatLoad(cache_b, Bsh, b_off, Bsh_stride, gl_CooperativeMatrixLayoutRowMajor); 331 [[unroll]] for (uint cm_row = 0; cm_row < cms_per_row; cm_row++) { 332 sums[cm_col * cms_per_row + cm_row] = coopMatMulAdd(cache_a[cm_row], cache_b, sums[cm_col * cms_per_row + cm_row]); 333 } 334 } 335 } 336#else 337 if (T_y * TS_K < K) { 338 UNROLL for (uint32_t CRS_lidx = 0; CRS_lidx < BS_CRS; CRS_lidx++) { 339 float regA[TS_K]; 340 float regB[TS_NPQ]; 341 for (uint32_t T_ly = 0; T_ly < TS_K; T_ly++) { 342 regA[T_ly] = Ash[(T_y * TS_K + T_ly) * Ash_stride + CRS_lidx]; 343 } 344 for (uint32_t T_lx = 0; T_lx < TS_NPQ; T_lx++) { 345 regB[T_lx] = Bsh[CRS_lidx * Bsh_stride + T_x * TS_NPQ + T_lx]; 346 } 347 for (uint32_t T_ly = 0; T_ly < TS_K; T_ly++) { 348 for (uint32_t T_lx = 0; T_lx < TS_NPQ; T_lx++) { 349 regC[T_ly][T_lx] = fma(regA[T_ly], regB[T_lx], regC[T_ly][T_lx]); 350 } 351 } 352 } 353 } 354#endif 355 barrier(); 356 } 357 /* Save C* */ 358#if defined(COOPMAT2) || defined(COOPMAT) 359 // stage matC into Csh, then write to dst with coalesced NPQ-contiguous stores 360#ifdef COOPMAT 361 const bool use_staged_store = true; 362#else 363 const bool use_staged_store = (csh_store != 0); 364#endif 365 if (use_staged_store) { 366#ifdef COOPMAT 367 // cm1: each subgroup stores its fragment grid into its Csh slot 368 [[unroll]] for (uint cm_row = 0; cm_row < cms_per_row; cm_row++) { 369 [[unroll]] for (uint cm_col = 0; cm_col < cms_per_col; cm_col++) { 370 const uint csh_off = (warp_r * WM + cm_row * TM) * Csh_stride + warp_c * WN + cm_col * TN; 371 coopMatStore(sums[cm_col * cms_per_row + cm_row], Csh, csh_off, Csh_stride, gl_CooperativeMatrixLayoutRowMajor); 372 } 373 } 374#else 375 coopMatStore(matC, Csh, 0, Csh_stride, gl_CooperativeMatrixLayoutRowMajor); 376#endif 377 barrier(); 378 379 // cooperative shmem->global: WG threads spread across BS_NPQ (the 380 // contiguous direction of dst), each iter covers store_rows_per_iter K-rows 381 const uint32_t store_rows_per_iter = WG_SIZE / BS_NPQ; 382 const uint32_t store_iters = BS_K / store_rows_per_iter; 383 const uint32_t k_thread_offset = tid / BS_NPQ; 384 const uint32_t npq_thread = tid % BS_NPQ; 385 [[unroll]] for (uint32_t i = 0; i < store_iters; i++) { 386 uint32_t k_local = i * store_rows_per_iter + k_thread_offset; 387 uint32_t K_idx = B_idx_K * BS_K + k_local; 388 uint32_t NPQ_idx = B_idx_NPQ * BS_NPQ + npq_thread; 389 uint32_t N_idx; 390 uint32_t OD_idx; 391 uint32_t OH_idx; 392 uint32_t OW_idx; 393 split_npq(NPQ_idx, N_idx, OD_idx, OH_idx, OW_idx); 394 uint32_t dst_idx = OW_idx + OH_idx * p.nb1 + OD_idx * p.nb2 + (N_idx * p.OC + K_idx) * p.nb3; 395 if (aligned != 0 || (K_idx < K && NPQ_idx < NPQ)) { 396 dst_data[dst_idx] = D_TYPE(Csh[k_local * Csh_stride + npq_thread]); 397 } 398 } 399 } 400#ifdef COOPMAT2 401 else { 402 coopMatPerElementNV(matC, matC, perElemOpStore); 403 } 404#endif 405#else 406 if (T_y * TS_K < K) { 407 for (uint32_t T_ly = 0; T_ly < TS_K; T_ly++) { 408 for (uint32_t T_lx = 0; T_lx < TS_NPQ; T_lx++) { 409 uint32_t K_idx = B_idx_K * BS_K + T_y * TS_K + T_ly; 410 uint32_t NPQ_idx = B_idx_NPQ * BS_NPQ + T_x * TS_NPQ + T_lx; 411 uint32_t N_idx; 412 uint32_t OD_idx; 413 uint32_t OH_idx; 414 uint32_t OW_idx; 415 split_npq(NPQ_idx, N_idx, OD_idx, OH_idx, OW_idx); 416 uint32_t dst_idx = OW_idx + OH_idx * p.nb1 + OD_idx * p.nb2 + (N_idx * p.OC + K_idx) * p.nb3; 417 if (aligned != 0 || (K_idx < K && NPQ_idx < NPQ)) { 418 dst_data[dst_idx] = D_TYPE(regC[T_ly][T_lx]); 419 } 420 } 421 } 422 } 423#endif 424}