char-fork/llama.cpp
git clone https://git.t4t.associates/char-fork/llama.cpp
b615f5b4b
char/aliceai-t5-moe
1#include "common.h" 2 3kernel void kernel_argmax_f32( 4 constant ggml_metal_kargs_argmax & args, 5 device const char * src0, 6 device char * dst, 7 threadgroup char * shmem [[threadgroup(0)]], 8 uint tgpig[[threadgroup_position_in_grid]], 9 uint tpitg[[thread_position_in_threadgroup]], 10 uint sgitg[[simdgroup_index_in_threadgroup]], 11 uint tiisg[[thread_index_in_simdgroup]], 12 uint ntg[[threads_per_threadgroup]]) { 13 device const float * x_row = (device const float *) ((device const char *) src0 + tgpig * args.nb01); 14 15 float lmax = -INFINITY; 16 int32_t larg = -1; 17 18 for (int i00 = tpitg; i00 < args.ne00; i00 += ntg) { 19 if (x_row[i00] > lmax) { 20 lmax = x_row[i00]; 21 larg = i00; 22 } 23 } 24 25 // find the argmax value in the block 26 float max_val = simd_max(lmax); 27 int32_t arg_val = simd_max(select(-1, larg, lmax == max_val)); 28 29 device int32_t * dst_i32 = (device int32_t *) dst; 30 31 threadgroup float * shared_maxval = (threadgroup float *) shmem; 32 threadgroup int32_t * shared_argmax = (threadgroup int32_t *) shmem + N_SIMDWIDTH; 33 34 if (ntg > N_SIMDWIDTH) { 35 if (sgitg == 0) { 36 shared_maxval[tiisg] = -INFINITY; 37 shared_argmax[tiisg] = -1; 38 } 39 40 threadgroup_barrier(mem_flags::mem_threadgroup); 41 42 if (tiisg == 0) { 43 shared_maxval[sgitg] = max_val; 44 shared_argmax[sgitg] = arg_val; 45 } 46 47 threadgroup_barrier(mem_flags::mem_threadgroup); 48 49 max_val = shared_maxval[tiisg]; 50 arg_val = shared_argmax[tiisg]; 51 52 float max_val_reduced = simd_max(max_val); 53 int32_t arg_val_reduced = simd_max(select(-1, arg_val, max_val == max_val_reduced)); 54 55 dst_i32[tgpig] = arg_val_reduced; 56 57 return; 58 } 59 60 dst_i32[tgpig] = arg_val; 61} 62 63kernel void kernel_diag_f32( 64 constant ggml_metal_kargs_diag & args, 65 device const char * src0, 66 device char * dst, 67 uint3 tgpig[[threadgroup_position_in_grid]], 68 ushort tiitg[[thread_index_in_threadgroup]]) { 69 constexpr short NW = N_SIMDWIDTH; 70 71 const int32_t i3 = tgpig.z; 72 const int32_t i2 = tgpig.y; 73 const int32_t i1 = tgpig.x; 74 75 device const float * src0_ptr = (device const float *)(src0 + i2*args.nb02 + i3*args.nb03); 76 device float * dst_ptr = (device float *)(dst + i1*args.nb01 + i2*args.nb2 + i3*args.nb3); 77 78 for (int i0 = tiitg; i0 < args.ne0; i0 += NW) { 79 dst_ptr[i0] = i0 == i1 ? src0_ptr[i0] : 0.0f; 80 } 81} 82 83kernel void kernel_roll_f32( 84 constant ggml_metal_kargs_roll & args, 85 device const char * src0, 86 device char * dst, 87 uint3 tgpig[[threadgroup_position_in_grid]], 88 uint3 tpitg[[thread_position_in_threadgroup]], 89 uint3 ntg[[threads_per_threadgroup]]) { 90 91 const int64_t i3 = tgpig.z; 92 const int64_t i2 = tgpig.y; 93 const int64_t i1 = tgpig.x; 94 95 device const float * src0_ptr = (device const float *) src0; 96 device float * dst_ptr = (device float *) dst; 97 98 for (int i0 = tpitg.x; i0 < args.ne0; i0 += ntg.x) { 99 // apply shifts and wrap around 100 int64_t i00 = i0 - args.s0; 101 int64_t i01 = i1 - args.s1; 102 int64_t i02 = i2 - args.s2; 103 int64_t i03 = i3 - args.s3; 104 105 if (i00 < 0) { i00 += args.ne00; } else if (i00 >= args.ne00) { i00 -= args.ne00; } 106 if (i01 < 0) { i01 += args.ne01; } else if (i01 >= args.ne01) { i01 -= args.ne01; } 107 if (i02 < 0) { i02 += args.ne02; } else if (i02 >= args.ne02) { i02 -= args.ne02; } 108 if (i03 < 0) { i03 += args.ne03; } else if (i03 >= args.ne03) { i03 -= args.ne03; } 109 110 int64_t src_idx = i03*args.ne02*args.ne01*args.ne00 + i02*args.ne01*args.ne00 + i01*args.ne00 + i00; 111 int64_t dst_idx = i3 *args.ne2 *args.ne1 *args.ne0 + i2 *args.ne1 *args.ne0 + i1 *args.ne0 + i0; 112 113 dst_ptr[dst_idx] = src0_ptr[src_idx]; 114 } 115} 116 117template <typename T> 118kernel void kernel_pad_impl( 119 constant ggml_metal_kargs_pad & args, 120 device const char * src0, 121 device char * dst, 122 uint3 tgpig[[threadgroup_position_in_grid]], 123 uint3 tpitg[[thread_position_in_threadgroup]], 124 uint3 ntg[[threads_per_threadgroup]]) { 125 const int32_t i3 = tgpig.z; 126 const int32_t i2 = tgpig.y; 127 const int32_t k0 = tgpig.x/args.ne1; 128 const int32_t i1 = tgpig.x - k0*args.ne1; 129 130 const int32_t i03 = i3; 131 const int32_t i02 = i2; 132 const int32_t i01 = i1; 133 134 device const T * src0_ptr = (device const T *) (src0 + i03*args.nb03 + i02*args.nb02 + i01*args.nb01); 135 device T * dst_ptr = (device T *) (dst + i3*args.nb3 + i2*args.nb2 + i1*args.nb1); 136 137 for (int32_t l0 = 0; l0 < 1024; l0 += ntg.x) { 138 const int32_t i0 = k0*1024 + tpitg.x + l0; 139 if (i0 >= args.ne0) { 140 break; 141 } 142 143 if (i0 < args.ne00 && i1 < args.ne01 && i2 < args.ne02 && i3 < args.ne03) { 144 dst_ptr[i0] = src0_ptr[i0]; 145 } else { 146 dst_ptr[i0] = 0.0f; 147 } 148 } 149} 150 151typedef decltype(kernel_pad_impl<float>) kernel_pad_t; 152 153template [[host_name("kernel_pad_f32")]] kernel kernel_pad_t kernel_pad_impl<float>; 154template [[host_name("kernel_pad_f32_4")]] kernel kernel_pad_t kernel_pad_impl<float4>; 155 156// TODO: this is slow - optimize 157kernel void kernel_pad_reflect_1d_f32( 158 constant ggml_metal_kargs_pad_reflect_1d & args, 159 device const char * src0, 160 device char * dst, 161 uint3 tgpig[[threadgroup_position_in_grid]], 162 uint3 tgpg[[threadgroups_per_grid]], 163 uint3 tpitg[[thread_position_in_threadgroup]], 164 uint3 ntg[[threads_per_threadgroup]]) { 165 166 const int64_t i3 = tgpig.z; 167 const int64_t i2 = tgpig.y; 168 const int64_t i1 = tgpig.x; 169 170 const int64_t i03 = i3; 171 const int64_t i02 = i2; 172 const int64_t i01 = i1; 173 174 device const float * src0_ptr = (device const float *) (src0 + i03*args.nb03 + i02*args.nb02 + i01*args.nb01); 175 device float * dst_ptr = (device float *) (dst + i3*args.nb3 + i2*args.nb2 + i1*args.nb1); 176 177 if (i1 < args.ne01 && i2 < args.ne02 && i3 < args.ne03) { 178 for (int i0 = tpitg.x; i0 < args.ne0; i0 += ntg.x) { 179 if (i0 < args.p0) { 180 dst_ptr[i0] = src0_ptr[args.p0 - i0]; 181 } else if (i0 < args.ne0 - args.p1) { 182 dst_ptr[i0] = src0_ptr[i0 - args.p0]; 183 } else { 184 dst_ptr[i0] = src0_ptr[(args.ne0 - args.p1 - args.p0) - (args.p1 + 1 - (args.ne0 - i0)) - 1]; 185 } 186 } 187 } 188} 189 190kernel void kernel_arange_f32( 191 constant ggml_metal_kargs_arange & args, 192 device char * dst, 193 uint3 tgpig[[threadgroup_position_in_grid]], 194 uint3 tpitg[[thread_position_in_threadgroup]], 195 uint3 ntg[[threads_per_threadgroup]]) { 196 197 device float * dst_ptr = (device float *) dst; 198 199 for (int i0 = tpitg.x; i0 < args.ne0; i0 += ntg.x) { 200 dst_ptr[i0] = args.start + args.step * i0; 201 } 202} 203 204kernel void kernel_timestep_embedding_f32( 205 constant ggml_metal_kargs_timestep_embedding & args, 206 device const char * src0, 207 device char * dst, 208 uint3 tgpig[[threadgroup_position_in_grid]], 209 uint3 tpitg[[thread_position_in_threadgroup]], 210 uint3 ntg[[threads_per_threadgroup]]) { 211 212 int i = tgpig.x; 213 device float * embed_data = (device float *)(dst + i*args.nb1); 214 215 int half_ = args.dim / 2; 216 for (int j = tpitg.x; j < half_; j += ntg.x) { 217 float timestep = ((device float *)src0)[i]; 218 float freq = (float)exp(-log((float)args.max_period) * j / half_); 219 float arg = timestep * freq; 220 embed_data[j ] = cos(arg); 221 embed_data[j + half_] = sin(arg); 222 } 223 224 if (args.dim % 2 != 0 && tpitg.x == 0) { 225 embed_data[2 * half_] = 0.f; 226 } 227} 228 229kernel void kernel_opt_step_adamw_f32( 230 constant ggml_metal_kargs_opt_step_adamw & args, 231 device float * x, 232 device const float * g, 233 device float * g_m, 234 device float * g_v, 235 device const float * pars, 236 uint gid[[thread_position_in_grid]]) { 237 238 if (gid >= args.np) { 239 return; 240 } 241 242 const float alpha = pars[0]; 243 const float beta1 = pars[1]; 244 const float beta2 = pars[2]; 245 const float eps = pars[3]; 246 const float wd = pars[4]; 247 const float beta1h = pars[5]; 248 const float beta2h = pars[6]; 249 250 const float gi = g[gid]; 251 const float gmi = g_m[gid] * beta1 + gi * (1.0f - beta1); 252 const float gvi = g_v[gid] * beta2 + gi * gi * (1.0f - beta2); 253 254 g_m[gid] = gmi; 255 g_v[gid] = gvi; 256 257 const float mh = gmi * beta1h; 258 const float vh = sqrt(gvi * beta2h) + eps; 259 260 x[gid] = x[gid] * (1.0f - alpha * wd) - alpha * mh / vh; 261} 262 263kernel void kernel_opt_step_sgd_f32( 264 constant ggml_metal_kargs_opt_step_sgd & args, 265 device float * x, 266 device const float * g, 267 device const float * pars, 268 uint gid[[thread_position_in_grid]]) { 269 270 if (gid >= args.np) { 271 return; 272 } 273 274 x[gid] = x[gid] * (1.0f - pars[0] * pars[1]) - pars[0] * g[gid]; 275} 276 277template<typename T> 278kernel void kernel_memset( 279 constant ggml_metal_kargs_memset & args, 280 device T * dst, 281 uint tpig[[thread_position_in_grid]]) { 282 dst[tpig] = args.val; 283} 284 285typedef decltype(kernel_memset<int64_t>) kernel_memset_t; 286 287template [[host_name("kernel_memset_i64")]] kernel kernel_memset_t kernel_memset<int64_t>; 288 289constant short FC_count_equal_nsg [[function_constant(FC_COUNT_EQUAL + 0)]]; 290 291template<typename T> 292kernel void kernel_count_equal( 293 constant ggml_metal_kargs_count_equal & args, 294 device const char * src0, 295 device const char * src1, 296 device atomic_int * dst, 297 threadgroup int32_t * shmem_i32 [[threadgroup(0)]], 298 uint3 tgpig[[threadgroup_position_in_grid]], 299 ushort3 tpitg[[thread_position_in_threadgroup]], 300 ushort sgitg[[simdgroup_index_in_threadgroup]], 301 ushort tiisg[[thread_index_in_simdgroup]], 302 ushort3 ntg[[threads_per_threadgroup]]) { 303 const short NSG = FC_count_equal_nsg; 304 305 const int i3 = tgpig.z; 306 const int i2 = tgpig.y; 307 const int i1 = tgpig.x; 308 309 if (i3 >= args.ne03 || i2 >= args.ne02 || i1 >= args.ne01) { 310 return; 311 } 312 313 int sum = 0; 314 315 device const char * base0 = src0 + i1*args.nb01 + i2*args.nb02 + i3*args.nb03; 316 device const char * base1 = src1 + i1*args.nb11 + i2*args.nb12 + i3*args.nb13; 317 318 for (int64_t i0 = tpitg.x; i0 < args.ne00; i0 += ntg.x) { 319 const T v0 = *(device const T *)(base0 + i0*args.nb00); 320 const T v1 = *(device const T *)(base1 + i0*args.nb10); 321 sum += (v0 == v1); 322 } 323 324 sum = simd_sum(sum); 325 326 if (tiisg == 0) { 327 shmem_i32[sgitg] = sum; 328 } 329 330 threadgroup_barrier(mem_flags::mem_threadgroup); 331 332 if (sgitg == 0) { 333 float v = 0.0f; 334 if (tpitg.x < NSG) { 335 v = shmem_i32[tpitg.x]; 336 } 337 338 float total = simd_sum(v); 339 if (tpitg.x == 0) { 340 atomic_fetch_add_explicit(dst, (int32_t) total, memory_order_relaxed); 341 } 342 } 343} 344 345typedef decltype(kernel_count_equal<int32_t>) kernel_count_equal_t; 346 347template [[host_name("kernel_count_equal_i32")]] kernel kernel_count_equal_t kernel_count_equal<int32_t>; 348 349template <typename T> 350kernel void kernel_snake( 351 constant ggml_metal_kargs_snake & args, 352 device const T * x, 353 device const float * a, 354 device const float * inv_b, 355 device T * dst, 356 uint tgpig [[threadgroup_position_in_grid]], 357 uint tpitg [[thread_position_in_threadgroup]], 358 uint ntg [[threads_per_threadgroup]]) { 359 360 const int idx = tgpig * ntg + tpitg; 361 if (idx >= args.T * args.C) { 362 return; 363 } 364 365 const int c = idx / args.T; // x is [T, C], a / inv_b collapse to [1, C] 366 const float xi = float(x[idx]); 367 const float si = sin(a[c] * xi); 368 dst[idx] = T(xi + si * si * inv_b[c]); 369} 370 371template [[host_name("kernel_snake_f32")]] kernel void kernel_snake<float>(constant ggml_metal_kargs_snake &, device const float *, device const float *, device const float *, device float *, uint, uint, uint); 372template [[host_name("kernel_snake_f16")]] kernel void kernel_snake<half>(constant ggml_metal_kargs_snake &, device const half *, device const float *, device const float *, device half *, uint, uint, uint); 373#if defined(GGML_METAL_HAS_BF16) 374template [[host_name("kernel_snake_bf16")]] kernel void kernel_snake<bfloat>(constant ggml_metal_kargs_snake &, device const bfloat *, device const float *, device const float *, device bfloat *, uint, uint, uint); 375#endif 376 377template<int N> 378kernel void kernel_fwht_f32( 379 constant ggml_metal_kargs_fwht & args, 380 device const float * src, 381 device float * dst, 382 uint3 tgpig[[threadgroup_position_in_grid]], 383 ushort sgitg[[simdgroup_index_in_threadgroup]], 384 ushort tiisg[[thread_index_in_simdgroup]], 385 ushort3 ntg[[threads_per_threadgroup]]) { 386 387 constexpr int NW = N_SIMDWIDTH; 388 constexpr int NE = N / NW; 389 390 const float scale = 1.0f / sqrt((float) N); 391 392 const int sg_per_tg = ntg.x / NW; 393 const int64_t r = tgpig.x * sg_per_tg + sgitg; 394 if (r >= args.nrows) { 395 return; 396 } 397 398 src += r * N; 399 dst += r * N; 400 401 const int lane = tiisg; 402 403 float reg[NE]; 404 for (int i = 0; i < NE; i++) { 405 reg[i] = src[i*NW + lane]*scale; 406 } 407 for (int i = 1; i < NW; i *= 2) { 408 for (int j = 0; j < NE; j++) { 409 const float val = reg[j]; 410 const float val2 = simd_shuffle_xor(val, i); 411 reg[j] = (lane & i) == 0 ? val2 + val : val2 - val; 412 } 413 } 414 415 for (int i = NW; i < N; i *= 2) { 416 const int step = i / NW; 417 for (int j = 0; j < NE; j += (2 * step)) { 418 for (int k = 0; k < step; k++) { 419 const float x = reg[j + k ]; 420 const float y = reg[j + k + step]; 421 reg[j + k] = x + y; 422 reg[j + k + step] = x - y; 423 } 424 } 425 } 426 427 for (int i = 0; i < NE; i++) { 428 dst[i*NW + lane] = reg[i]; 429 } 430} 431 432typedef decltype(kernel_fwht_f32<64>) kernel_fwht_t; 433 434template [[host_name("kernel_fwht_f32_64")]] kernel kernel_fwht_t kernel_fwht_f32<64>; 435template [[host_name("kernel_fwht_f32_128")]] kernel kernel_fwht_t kernel_fwht_f32<128>; 436template [[host_name("kernel_fwht_f32_256")]] kernel kernel_fwht_t kernel_fwht_f32<256>; 437template [[host_name("kernel_fwht_f32_512")]] kernel kernel_fwht_t kernel_fwht_f32<512>; 438 439kernel void kernel_dsv4_hc_comb_f32( 440 constant ggml_metal_kargs_dsv4_hc_comb & args, 441 device const char * mixes, 442 device const char * scale, 443 device const char * base, 444 device char * dst, 445 uint3 tgpig[[threadgroup_position_in_grid]], 446 ushort tiisg[[thread_index_in_simdgroup]], 447 ushort sgitg[[simdgroup_index_in_threadgroup]], 448 ushort3 ntg[[threads_per_threadgroup]]) { 449 constexpr ushort hc = 4; 450 constexpr ushort comb_offset = 2*hc; 451 452 const int it = tgpig.x*ntg.y + sgitg; 453 if (it >= args.n_tokens) { 454 return; 455 } 456 457 float scale_lane = 0.0f; 458 if (tiisg == 0) { 459 scale_lane = *(device const float *) (scale + 2*args.nb_s0); 460 } 461 const float scale_comb = simd_shuffle(scale_lane, 0); 462 463 float v = 0.0f; 464 if (tiisg < hc*hc) { 465 v = *(device const float *) (mixes + (comb_offset + tiisg)*args.nb_m0 + it*args.nb_m1)*scale_comb 466 + *(device const float *) (base + (comb_offset + tiisg)*args.nb_b0); 467 } 468 469 // Softmax across destinations (the four contiguous lanes for each source). 470 float vmax = max(v, simd_shuffle_xor(v, 1)); 471 vmax = max(vmax, simd_shuffle_xor(vmax, 2)); 472 v = exp(v - vmax); 473 474 float sum = v + simd_shuffle_xor(v, 1); 475 sum += simd_shuffle_xor(sum, 2); 476 v = v/sum + args.eps; 477 478 // Normalize columns: equal destination indices are four lanes apart. 479 sum = v + simd_shuffle_xor(v, 4); 480 sum += simd_shuffle_xor(sum, 8); 481 v /= sum + args.eps; 482 483 for (int i = 1; i < args.n_iter; ++i) { 484 sum = v + simd_shuffle_xor(v, 1); 485 sum += simd_shuffle_xor(sum, 2); 486 v /= sum + args.eps; 487 488 sum = v + simd_shuffle_xor(v, 4); 489 sum += simd_shuffle_xor(sum, 8); 490 v /= sum + args.eps; 491 } 492 493 if (tiisg < hc*hc) { 494 const ushort idst = tiisg & 3; 495 const ushort isrc = tiisg >> 2; 496 *(device float *) (dst + idst*args.nb_d0 + isrc*args.nb_d1 + it*args.nb_d2) = v; 497 } 498} 499 500kernel void kernel_dsv4_hc_pre_f32( 501 constant ggml_metal_kargs_dsv4_hc_pre & args, 502 device const char * x, 503 device const char * weights, 504 device char * dst, 505 uint3 tgpig[[threadgroup_position_in_grid]], 506 ushort tiisg[[thread_index_in_simdgroup]], 507 ushort sgitg[[simdgroup_index_in_threadgroup]], 508 ushort3 ntg[[threads_per_threadgroup]]) { 509 constexpr ushort hc = 4; 510 511 const int it = tgpig.y; 512 const int i0 = ((int) tgpig.x*ntg.y + sgitg)*32 + tiisg; 513 514 float weight_lane = 0.0f; 515 if (tiisg < hc) { 516 weight_lane = *(device const float *) (weights + tiisg*args.nb_w0 + it*args.nb_w1); 517 } 518 519 float w[hc]; 520 FOR_UNROLL (ushort ih = 0; ih < hc; ++ih) { 521 w[ih] = simd_shuffle(weight_lane, ih); 522 } 523 524 if (i0 >= args.n_embd) { 525 return; 526 } 527 528 device const char * xb = x + i0*args.nb_x0 + it*args.nb_x2; 529 float result = 0.0f; 530 FOR_UNROLL (ushort ih = 0; ih < hc; ++ih) { 531 result = fma(*(device const float *) (xb + ih*args.nb_x1), w[ih], result); 532 } 533 534 *(device float *) (dst + i0*args.nb_d0 + it*args.nb_d1) = result; 535} 536 537kernel void kernel_dsv4_hc_post_f32( 538 constant ggml_metal_kargs_dsv4_hc_post & args, 539 device const char * x, 540 device const char * residual, 541 device const char * post, 542 device const char * comb, 543 device char * dst, 544 uint3 tgpig[[threadgroup_position_in_grid]], 545 ushort tiisg[[thread_index_in_simdgroup]], 546 ushort sgitg[[simdgroup_index_in_threadgroup]], 547 ushort3 ntg[[threads_per_threadgroup]]) { 548 constexpr ushort hc = 4; 549 550 const int it = tgpig.y; 551 const int i0 = ((int) tgpig.x*ntg.y + sgitg)*32 + tiisg; 552 553 float coeff_lane = 0.0f; 554 if (tiisg < hc) { 555 coeff_lane = *(device const float *) (post + tiisg*args.nb_p0 + it*args.nb_p1); 556 } else if (tiisg < hc + hc*hc) { 557 const ushort idx = tiisg - hc; 558 const ushort idst = idx & 3; 559 const ushort isrc = idx >> 2; 560 coeff_lane = *(device const float *) (comb + idst*args.nb_c0 + isrc*args.nb_c1 + it*args.nb_c2); 561 } 562 563 float post_reg[hc]; 564 float comb_reg[hc][hc]; 565 FOR_UNROLL (ushort idst = 0; idst < hc; ++idst) { 566 post_reg[idst] = simd_shuffle(coeff_lane, idst); 567 } 568 FOR_UNROLL (ushort isrc = 0; isrc < hc; ++isrc) { 569 FOR_UNROLL (ushort idst = 0; idst < hc; ++idst) { 570 comb_reg[isrc][idst] = simd_shuffle(coeff_lane, hc + idst + hc*isrc); 571 } 572 } 573 574 if (i0 >= args.n_embd) { 575 return; 576 } 577 578 const float xv = *(device const float *) (x + i0*args.nb_x0 + it*args.nb_x1); 579 float result[hc]; 580 FOR_UNROLL (ushort idst = 0; idst < hc; ++idst) { 581 result[idst] = xv*post_reg[idst]; 582 } 583 584 device const char * rb = residual + i0*args.nb_r0 + it*args.nb_r2; 585 FOR_UNROLL (ushort isrc = 0; isrc < hc; ++isrc) { 586 const float rv = *(device const float *) (rb + isrc*args.nb_r1); 587 FOR_UNROLL (ushort idst = 0; idst < hc; ++idst) { 588 result[idst] = fma(rv, comb_reg[isrc][idst], result[idst]); 589 } 590 } 591 592 FOR_UNROLL (ushort idst = 0; idst < hc; ++idst) { 593 *(device float *) (dst + i0*args.nb_d0 + idst*args.nb_d1 + it*args.nb_d2) = result[idst]; 594 } 595}