@@ -7596,74 +7596,151 @@ void ggml_compute_forward_ssm_conv(
7596
7596
static void ggml_compute_forward_ssm_scan_f32 (
7597
7597
const ggml_compute_params * params,
7598
7598
ggml_tensor * dst) {
7599
- const ggml_tensor * src0 = dst->src [0 ]; // s
7600
- const ggml_tensor * src1 = dst->src [1 ]; // x
7601
- const ggml_tensor * src2 = dst->src [2 ]; // dt
7602
- const ggml_tensor * src3 = dst->src [3 ]; // A
7603
- const ggml_tensor * src4 = dst->src [4 ]; // B
7604
- const ggml_tensor * src5 = dst->src [5 ]; // C
7599
+ const ggml_tensor * src0 = dst->src [0 ]; // s {d_state, dim, n_head, n_seqs+}
7600
+ const ggml_tensor * src1 = dst->src [1 ]; // x {dim, n_head, n_seq_tokens, n_seqs}
7601
+ const ggml_tensor * src2 = dst->src [2 ]; // dt {n_head, n_seq_tokens, n_seqs}
7602
+ const ggml_tensor * src3 = dst->src [3 ]; // A {d_state, n_head} or {1, n_head}
7603
+ const ggml_tensor * src4 = dst->src [4 ]; // B {d_state, n_group, n_seq_tokens, n_seqs}
7604
+ const ggml_tensor * src5 = dst->src [5 ]; // C {d_state, n_group, n_seq_tokens, n_seqs}
7605
+ const ggml_tensor * src6 = dst->src [6 ]; // ids {n_seqs}
7605
7606
7606
7607
const int ith = params->ith ;
7607
7608
const int nth = params->nth ;
7608
7609
7609
- const int64_t nc = src0->ne [0 ]; // d_state
7610
- const int64_t nr = src0->ne [1 ]; // d_inner
7611
- const int64_t n_t = src1->ne [1 ]; // number of tokens per sequence
7612
- const int64_t n_s = src0->ne [2 ]; // number of sequences in the batch
7610
+ const int64_t nc = src0->ne [0 ]; // d_state
7611
+ const int64_t nr = src0->ne [1 ]; // dim
7612
+ const int64_t nh = src1->ne [1 ]; // n_head
7613
+ const int64_t ng = src4->ne [1 ];
7614
+ const int64_t nt = src1->ne [2 ]; // number of tokens per sequence
7615
+ const int64_t ns = src1->ne [3 ]; // number of sequences in the batch
7613
7616
7614
- GGML_ASSERT (ggml_nelements (src1) + ggml_nelements (src0) == ggml_nelements (dst));
7617
+ // can't use ggml_nbytes because src1 is not necessarily contiguous
7618
+ const int64_t s_off = ggml_nelements (src1) * ggml_element_size (src1);
7619
+
7620
+ GGML_ASSERT (ggml_nelements (src1) + nc*nr*nh*ns == ggml_nelements (dst));
7615
7621
GGML_ASSERT (src0->nb [0 ] == sizeof (float ));
7616
7622
GGML_ASSERT (src1->nb [0 ] == sizeof (float ));
7617
7623
GGML_ASSERT (src2->nb [0 ] == sizeof (float ));
7618
7624
GGML_ASSERT (src3->nb [0 ] == sizeof (float ));
7619
7625
GGML_ASSERT (src4->nb [0 ] == sizeof (float ));
7620
7626
GGML_ASSERT (src5->nb [0 ] == sizeof (float ));
7621
- // required for the dot product between s and C
7622
- GGML_ASSERT (src0->nb [1 ] == src0->ne [0 ]*sizeof (float ));
7623
- // required for per-sequence offsets for states
7624
- GGML_ASSERT (src0->nb [2 ] == src0->ne [0 ]*src0->ne [1 ]*sizeof (float ));
7625
- // required to get correct offset for state destination (i.e. src1->nb[3])
7626
- GGML_ASSERT (src1->nb [3 ] == src1->ne [0 ]*src1->ne [1 ]*src1->ne [2 ]*sizeof (float ));
7627
+ GGML_ASSERT (src6->nb [0 ] == sizeof (int32_t ));
7628
+ // allows optimizing the modulo since n_group should be a power of 2
7629
+ GGML_ASSERT ((ng & -ng) == ng);
7630
+
7631
+ // heads per thread
7632
+ const int dh = (nh + nth - 1 )/nth;
7633
+
7634
+ // head range for this thread
7635
+ const int ih0 = dh*ith;
7636
+ const int ih1 = MIN (ih0 + dh, nh);
7637
+
7638
+ const int32_t * ids = (const int32_t *) src6->data ;
7639
+
7640
+ for (int i3 = 0 ; i3 < ns; ++i3) {
7641
+ const float * s0 = (const float *) ((const char *) src0->data + ids[i3]*(src0->nb [3 ])); // {d_state, dim, nh, ns}
7642
+ float * s = ( float *) (( char *) dst->data + i3*(src0->nb [3 ]) + s_off); // {d_state, dim, nh, ns}
7643
+
7644
+ for (int i2 = 0 ; i2 < nt; ++i2) {
7645
+ const float * x = (const float *) ((const char *) src1->data + i2*(src1->nb [2 ]) + i3*(src1->nb [3 ])); // {dim, nh, nt, ns}
7646
+ const float * dt = (const float *) ((const char *) src2->data + i2*(src2->nb [1 ]) + i3*(src2->nb [2 ])); // {nh, nt, ns}
7647
+ const float * A = (const float *) ((const char *) src3->data ); // {d_state, nh} or {1, nh}
7648
+ const float * B = (const float *) ((const char *) src4->data + i2*(src4->nb [2 ]) + i3*(src4->nb [3 ])); // {d_state, ng, nt, ns}
7649
+ const float * C = (const float *) ((const char *) src5->data + i2*(src5->nb [2 ]) + i3*(src5->nb [3 ])); // {d_state, ng, nt, ns}
7650
+ float * y = ( float *) (( char *) dst->data + i2*(nh*nr*sizeof (float )) + i3*(nt*nh*nr*sizeof (float ))); // {dim, nh, nt, ns}
7651
+
7652
+ if (src3->ne [0 ] == 1 ) {
7653
+ // Mamba-2 has a scalar decay factor per head; dA can be outside the state-wise loop
7654
+
7655
+ // n_head
7656
+ for (int h = ih0; h < ih1; ++h) {
7657
+ // ref: https://github.com/state-spaces/mamba/blob/62db608da60f6fc790b8ed9f4b3225e95ca15fde/mamba_ssm/ops/triton/softplus.py#L16
7658
+ const float dt_soft_plus = dt[h] <= 20 .0f ? log1pf (expf (dt[h])) : dt[h];
7659
+ const float dA = expf (dt_soft_plus * A[h]);
7660
+
7661
+ // dim
7662
+ for (int i1 = 0 ; i1 < nr; ++i1) {
7663
+ const int ii = i1 + h*nr;
7664
+ const float x_dt = x[ii] * dt_soft_plus;
7665
+ float sumf = 0 .0f ;
7666
+ #if defined(GGML_SIMD)
7667
+ const int np = (nc & ~(GGML_F32_STEP - 1 ));
7627
7668
7628
- // rows per thread
7629
- const int dr = (nr + nth - 1 )/nth;
7669
+ GGML_F32_VEC sum[GGML_F32_ARR] = { GGML_F32_VEC_ZERO };
7630
7670
7631
- // row range for this thread
7632
- const int ir0 = dr*ith;
7633
- const int ir1 = MIN (ir0 + dr, nr);
7634
- const int ir = ir1 - ir0;
7671
+ GGML_F32_VEC adA = GGML_F32_VEC_SET1 (dA);
7672
+ GGML_F32_VEC axdt = GGML_F32_VEC_SET1 (x_dt);
7635
7673
7636
- for (int i3 = 0 ; i3 < n_s; ++i3) {
7637
- for (int i2 = 0 ; i2 < n_t ; ++i2) {
7638
- const float * s0 = (const float *) ((const char *) src0->data + ir0*(src0->nb [1 ]) + i3*(src0->nb [2 ])); // {d_state, d_inner, n_s}
7639
- const float * x = (const float *) ((const char *) src1->data + ir0*(src1->nb [0 ]) + i2*(src1->nb [1 ]) + i3*(src1->nb [2 ])); // {d_inner, n_t, n_s}
7640
- const float * dt = (const float *) ((const char *) src2->data + ir0*(src2->nb [0 ]) + i2*(src2->nb [1 ]) + i3*(src2->nb [2 ])); // {d_inner, n_t, n_s}
7641
- const float * A = (const float *) ((const char *) src3->data + ir0*(src3->nb [1 ])); // {d_state, d_inner}
7642
- const float * B = (const float *) ((const char *) src4->data + i2*(src4->nb [1 ]) + i3*(src4->nb [2 ])); // {d_state, n_t, n_s}
7643
- const float * C = (const float *) ((const char *) src5->data + i2*(src5->nb [1 ]) + i3*(src5->nb [2 ])); // {d_state, n_t, n_s}
7644
- float * y = ( float *) (( char *) dst->data + ir0*(src1->nb [0 ]) + i2*(src1->nb [1 ]) + i3*(src1->nb [2 ])); // {d_inner, n_t, n_s}
7645
- float * s = ( float *) (( char *) dst->data + ir0*(src0->nb [1 ]) + i3*(src0->nb [2 ]) + src1->nb [3 ]); // {d_state, d_inner, n_s}
7646
-
7647
- // use the output as the source for the next token-wise iterations
7648
- if (i2 > 0 ) { s0 = s; }
7674
+ GGML_F32_VEC ax[GGML_F32_ARR];
7675
+ GGML_F32_VEC ay[GGML_F32_ARR];
7676
+ GGML_F32_VEC az[GGML_F32_ARR];
7649
7677
7650
- // d_inner
7651
- for (int i1 = 0 ; i1 < ir; ++i1) {
7652
- // ref: https://github.com/state-spaces/mamba/blob/34076d664838588a3c97727b263478ab9f621a07/mamba_ssm/ops/triton/selective_state_update.py#L78
7653
- float dt_soft_plus = dt[i1] <= 20 .0f ? log1pf (expf (dt[i1])) : dt[i1];
7654
- float x_dt = x[i1] * dt_soft_plus;
7655
- float sumf = 0 .0f ;
7656
- // d_state
7657
- for (int i0 = 0 ; i0 < nc; ++i0) {
7658
- int i = i0 + i1*nc;
7659
- // state = prev_state * dA + dB * x
7660
- float state = (s0[i] * expf (dt_soft_plus * A[i])) + (B[i0] * x_dt);
7661
- // y = rowwise_dotprod(state, C)
7662
- sumf += state * C[i0];
7663
- s[i] = state;
7678
+ for (int i = 0 ; i < np; i += GGML_F32_STEP) {
7679
+ for (int j = 0 ; j < GGML_F32_ARR; j++) {
7680
+ ax[j] = GGML_F32_VEC_LOAD (s0 + i + j*GGML_F32_EPR + ii*nc);
7681
+ ay[j] = GGML_F32_VEC_LOAD (B + i + j*GGML_F32_EPR + (h & (ng - 1 ))*nc);
7682
+ az[j] = GGML_F32_VEC_LOAD (C + i + j*GGML_F32_EPR + (h & (ng - 1 ))*nc);
7683
+
7684
+ ax[j] = GGML_F32_VEC_MUL (ax[j], adA);
7685
+ ay[j] = GGML_F32_VEC_MUL (ay[j], axdt);
7686
+
7687
+ ax[j] = GGML_F32_VEC_ADD (ax[j], ay[j]);
7688
+
7689
+ sum[j] = GGML_F32_VEC_FMA (sum[j], ax[j], az[j]);
7690
+
7691
+ GGML_F32_VEC_STORE (s + i + j*GGML_F32_EPR + ii*nc, ax[j]);
7692
+ }
7693
+ }
7694
+
7695
+ // reduce sum0..sum3 to sum0
7696
+ GGML_F32_VEC_REDUCE (sumf, sum);
7697
+ #else
7698
+ const int np = 0 ;
7699
+ #endif
7700
+ // d_state
7701
+ for (int i0 = np; i0 < nc; ++i0) {
7702
+ const int i = i0 + ii*nc;
7703
+ const int ig = i0 + (h & (ng - 1 ))*nc;
7704
+ // state = prev_state * dA + dB * x
7705
+ const float state = (s0[i] * dA) + (B[ig] * x_dt);
7706
+ // y = rowwise_dotprod(state, C)
7707
+ sumf += state * C[ig];
7708
+ s[i] = state;
7709
+ }
7710
+ y[ii] = sumf;
7711
+ }
7712
+ }
7713
+ } else {
7714
+ // Mamba-1 has an element-wise decay factor for the states
7715
+
7716
+ // n_head
7717
+ for (int h = ih0; h < ih1; ++h) {
7718
+ // ref: https://github.com/state-spaces/mamba/blob/62db608da60f6fc790b8ed9f4b3225e95ca15fde/mamba_ssm/ops/triton/softplus.py#L16
7719
+ const float dt_soft_plus = dt[h] <= 20 .0f ? log1pf (expf (dt[h])) : dt[h];
7720
+
7721
+ // dim
7722
+ for (int i1 = 0 ; i1 < nr; ++i1) {
7723
+ const int ii = i1 + h*nr;
7724
+ const float x_dt = x[ii] * dt_soft_plus;
7725
+ float sumf = 0 .0f ;
7726
+ // NOTE: can't really use GGML_SIMD here because d_state is usually 16
7727
+ // and also because expf is used within the loop.
7728
+ // d_state
7729
+ for (int i0 = 0 ; i0 < nc; ++i0) {
7730
+ const int i = i0 + ii*nc;
7731
+ const int ig = i0 + (h & (ng - 1 ))*nc;
7732
+ // state = prev_state * dA + dB * x
7733
+ const float state = (s0[i] * expf (dt_soft_plus * A[i0 + h*nc])) + (B[ig] * x_dt);
7734
+ // y = rowwise_dotprod(state, C)
7735
+ sumf += state * C[ig];
7736
+ s[i] = state;
7737
+ }
7738
+ y[ii] = sumf;
7739
+ }
7664
7740
}
7665
- y[i1] = sumf;
7666
7741
}
7742
+ // use the output as the source when it's not the first token-wise iteration
7743
+ s0 = s;
7667
7744
}
7668
7745
}
7669
7746
}
0 commit comments