char-fork/llama.cpp

git clone https://git.t4t.associates/char-fork/llama.cpp

YiChen Lvmetal: per-op source split + parallel compile (#26561)b615f5b4b

char/aliceai-t5-moe
19.5 KiB595 linesraw
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}