From f4d4089a8d5b3a9c75796accbaa885225059e053 Mon Sep 17 00:00:00 2001 From: XSquirrelC Date: Fri, 4 Sep 2026 02:30:13 +0000 Subject: [PATCH 1/3] ggml: add GGML_TYPE_I8_S / I2_S and five fused CPU INT8 ops --- CMakeLists.txt | 7 + external/ggml/include/ggml.h | 89 ++++- external/ggml/src/ggml-cpu/ggml-cpu.c | 37 ++ external/ggml/src/ggml-cpu/ops.cpp | 407 ++++++++++++++++++++ external/ggml/src/ggml-cpu/ops.h | 7 + external/ggml/src/ggml-impl.h | 20 + external/ggml/src/ggml-quants.c | 132 +++++++ external/ggml/src/ggml-quants.h | 17 + external/ggml/src/ggml.c | 215 ++++++++++- tests/unittests/test_i8_s_fused_ops.cpp | 487 ++++++++++++++++++++++++ 10 files changed, 1414 insertions(+), 4 deletions(-) create mode 100644 tests/unittests/test_i8_s_fused_ops.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 313f05b7e..8d1537e5c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2254,6 +2254,13 @@ if (ENGINE_BUILD_TESTS) COMMAND conv_lowering_matrix_test ) + add_engine_unittest(i8_s_fused_ops_test tests/unittests/test_i8_s_fused_ops.cpp) + + add_test( + NAME i8_s_fused_ops_test + COMMAND i8_s_fused_ops_test + ) + add_engine_unittest(gguf_tensor_source_test tests/unittests/test_gguf_tensor_source.cpp) target_include_directories(gguf_tensor_source_test PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/tests/unittests) diff --git a/external/ggml/include/ggml.h b/external/ggml/include/ggml.h index 0de79aed5..1626f0105 100644 --- a/external/ggml/include/ggml.h +++ b/external/ggml/include/ggml.h @@ -429,7 +429,15 @@ extern "C" { GGML_TYPE_MXFP4 = 39, // MXFP4 (1 block) GGML_TYPE_NVFP4 = 40, // NVFP4 (4 blocks, E4M3 scale) GGML_TYPE_Q1_0 = 41, - GGML_TYPE_COUNT = 42, + // INT8 / ternary with a single per-tensor scale, stored as one F32 + // immediately after the payload rather than interleaved per block. + // Used by the VibeASR CPU pipeline; see ggml_mul_mat_add(). + // NOTE: 36/37 are deliberately avoided even though VibeASR's own fork + // uses them -- they are retired IQ4_NL_4_4/4_8 slots and reusing an ID + // would silently misread GGUF files that still carry the old type. + GGML_TYPE_I8_S = 42, + GGML_TYPE_I2_S = 43, + GGML_TYPE_COUNT = 44, }; // precision @@ -588,6 +596,14 @@ extern "C" { GGML_OP_GLU, GGML_OP_CONVROT_LINEAR, + + // VibeASR CPU INT8 pipeline. Appended at the tail so every existing + // op keeps its value -- GGML_OP_NAME and GGML_OP_SYMBOL are positional. + GGML_OP_ADD_SCALED, + GGML_OP_RMS_NORM_SCALED, + GGML_OP_MUL_MAT_ADD, + GGML_OP_MUL_MAT_ADD_RELU, + GGML_OP_IM2COL_ASYM, GGML_OP_COUNT, }; @@ -2472,6 +2488,77 @@ extern "C" { struct ggml_tensor * bias, int group_size); + // VibeASR CPU INT8 pipeline (GGML_TYPE_I8_S / GGML_TYPE_I2_S). + // + // Ported from https://github.com/microsoft/VibeASR.cpp, an end-to-end INT8 + // ASR stack (INT8 VAE encoder, ternary-weight language model) built for CPU + // inference on edge devices. + // + // These are CPU-only, mirroring how ggml_convrot_linear and + // ggml_sage_attn2_i8 above are CUDA-only. + // + // Unlike those two, the scale is NOT a separate F32 src tensor. An I8_S + // activation's scale is recomputed from that activation at run time, and a + // ggml node has exactly one output, so the scale has to travel with the + // data: every I8_S/I2_S tensor stores one F32 immediately after its int8 + // payload (see ggml_type_extra_bytes in ggml.c). Weight scales could have + // used the separate-tensor convention, but sharing one representation with + // activations keeps a single kernel per op instead of two. + + // y = a*scale + b, fusing a ConvNeXt LayerScale into its residual add. + // a and b are I8_S and same-shape, scale is F32 per-channel broadcast on + // ne[0]. Output is I8_S and carries a freshly computed per-tensor scale. + GGML_API struct ggml_tensor * ggml_add_scaled( + struct ggml_context * ctx, + struct ggml_tensor * a, + struct ggml_tensor * b, + struct ggml_tensor * scale); + + // y = rms_norm(a) * scale, fused. a is I8_S, scale is F32 per-channel, + // output is I8_S. Equivalent to ggml_mul(ggml_rms_norm(a), scale) but + // avoids materializing the F32 intermediate. + GGML_API struct ggml_tensor * ggml_rms_norm_scaled( + struct ggml_context * ctx, + struct ggml_tensor * a, + struct ggml_tensor * scale, + float eps); + + // y = a*b + bias, with a I8_S weights and b I8_S activations. bias is F32 + // and broadcasts on ne[0]. Output is I8_S. The ternary I2_S weights of the + // language model go through plain ggml_mul_mat, which has no bias to fuse. + // + // a with ne[1] == 1 and ne[2] > 1 selects a depthwise contraction: one + // length-ne[0] filter per channel, output indexed channel-major. + GGML_API struct ggml_tensor * ggml_mul_mat_add( + struct ggml_context * ctx, + struct ggml_tensor * a, + struct ggml_tensor * b, + struct ggml_tensor * bias); + + // As ggml_mul_mat_add, with ReLU folded into the epilogue. + GGML_API struct ggml_tensor * ggml_mul_mat_add_relu( + struct ggml_context * ctx, + struct ggml_tensor * a, + struct ggml_tensor * b, + struct ggml_tensor * bias); + + // im2col with independent left/right padding on the width axis. ggml_im2col + // only takes a single symmetric p0, so causal 1D convolutions otherwise need + // a separate ggml_pad_ext node and a full copy of the activation. + GGML_API struct ggml_tensor * ggml_im2col_asym( + struct ggml_context * ctx, + struct ggml_tensor * a, + struct ggml_tensor * b, + int s0, + int s1, + int lp0, + int rp0, + int p1, + int d0, + int d1, + bool is_2D, + enum ggml_type dst_type); + // MINITTS_FLASH_BIAS_WRAPPER: // Helper for models that already assemble a dense additive attention bias // (for example relative-position scores). The helper expands an optional diff --git a/external/ggml/src/ggml-cpu/ggml-cpu.c b/external/ggml/src/ggml-cpu/ggml-cpu.c index d9ec09939..66d2dce64 100644 --- a/external/ggml/src/ggml-cpu/ggml-cpu.c +++ b/external/ggml/src/ggml-cpu/ggml-cpu.c @@ -2092,6 +2092,26 @@ static void ggml_compute_forward(struct ggml_compute_params * params, struct ggm ggml_compute_forward_opt_step_sgd(params, tensor); } break; + case GGML_OP_ADD_SCALED: + { + ggml_compute_forward_add_scaled(params, tensor); + } break; + case GGML_OP_RMS_NORM_SCALED: + { + ggml_compute_forward_rms_norm_scaled(params, tensor); + } break; + case GGML_OP_MUL_MAT_ADD: + { + ggml_compute_forward_mul_mat_add(params, tensor); + } break; + case GGML_OP_MUL_MAT_ADD_RELU: + { + ggml_compute_forward_mul_mat_add_relu(params, tensor); + } break; + case GGML_OP_IM2COL_ASYM: + { + ggml_compute_forward_im2col_asym(params, tensor); + } break; case GGML_OP_NONE: { // nop @@ -2345,6 +2365,7 @@ static int ggml_get_n_tasks(struct ggml_tensor * node, int n_threads) { } break; case GGML_OP_IM2COL: case GGML_OP_IM2COL_FAST_1D: + case GGML_OP_IM2COL_ASYM: case GGML_OP_IM2COL_BACK: case GGML_OP_IM2COL_3D: case GGML_OP_CONV_2D: @@ -2352,6 +2373,10 @@ static int ggml_get_n_tasks(struct ggml_tensor * node, int n_threads) { case GGML_OP_CONV_2D_DW: case GGML_OP_CONV_TRANSPOSE_1D: case GGML_OP_CONV_TRANSPOSE_2D: + case GGML_OP_ADD_SCALED: + case GGML_OP_RMS_NORM_SCALED: + case GGML_OP_MUL_MAT_ADD: + case GGML_OP_MUL_MAT_ADD_RELU: { n_tasks = n_threads; } break; @@ -2948,6 +2973,18 @@ struct ggml_cplan ggml_graph_plan( { cur = ggml_type_size(node->type)*(n_tasks + node->src[0]->ne[0]*n_tasks); } break; + case GGML_OP_ADD_SCALED: + case GGML_OP_RMS_NORM_SCALED: + case GGML_OP_MUL_MAT_ADD: + case GGML_OP_MUL_MAT_ADD_RELU: + { + // F32 staging for the whole output, plus one absmax per + // thread. The output cannot be quantized in place: its + // scale is only known once every element has been + // computed, so all of them have to be held somewhere + // first. + cur = sizeof(float)*ggml_nelements(node) + sizeof(float)*n_tasks; + } break; case GGML_OP_GATED_DELTA_NET: { const int64_t S_v = node->src[2]->ne[0]; diff --git a/external/ggml/src/ggml-cpu/ops.cpp b/external/ggml/src/ggml-cpu/ops.cpp index 0f0f57399..fc2ae9177 100644 --- a/external/ggml/src/ggml-cpu/ops.cpp +++ b/external/ggml/src/ggml-cpu/ops.cpp @@ -11637,3 +11637,410 @@ void ggml_compute_forward_fwht(const ggml_compute_params * params, ggml_tensor * } } } + +// VibeASR CPU INT8 pipeline +// +// Scalar reference implementations of the five fused I8_S ops. They are written +// to be obviously correct rather than fast: the vectorized kernels live with the +// other SIMD code under arch/, and dispatch to them replaces the inner loops +// here without changing the surrounding structure. +// +// Four of the five produce I8_S, and all four share one shape: +// +// 1. compute the result in F32 into params->wdata, tracking this thread's absmax +// 2. barrier, reduce the per-thread absmaxes to a tensor-wide absmax +// 3. quantize this thread's slice into dst; thread 0 writes the output scale +// +// Step 2 is why these are fused ops at all. The output scale cannot be known +// before the whole output has been computed, so an unfused chain would have to +// materialize each intermediate in F32, rescan it, and write it again. +// +// wdata layout, shared by all four: +// +// [0, n_stage) F32 staging for the result +// [n_stage, n_stage + nth) per-thread absmax +// [n_stage + nth, ...) op-specific scratch (int32 accumulators) + +// Steps 2 and 3. The thread's slice of the output is the rectangle +// {row*row_stride + col : row < n_rows, col0 <= col < col1}; contiguous callers +// pass n_rows = 1, row_stride = 0. Both `stage` and dst use the same indexing. +// +// `relu` clamps after rounding rather than before the absmax scan, so the +// negatives that are about to be zeroed still widen the scale: an output +// spanning [-10, +5] scales by 10, leaving the surviving positives to reach 63 +// of the available 127. That is what the reference implementation does, and +// matching it is what makes a layer-by-layer comparison meaningful, so it is +// reproduced here deliberately rather than improved in passing. +static void ggml_i8_s_requantize( + const ggml_compute_params * params, + ggml_tensor * dst, + const float * stage, + float * thread_max, + int64_t row_stride, + int64_t n_rows, + int64_t col0, + int64_t col1, + float local_absmax, + bool relu) { + + const int ith = params->ith; + const int nth = params->nth; + + thread_max[ith] = local_absmax; + + ggml_barrier(params->threadpool); + + if (ith == 0) { + float global_max = 0.0f; + for (int t = 0; t < nth; t++) { + if (thread_max[t] > global_max) global_max = thread_max[t]; + } + thread_max[0] = global_max; + } + + ggml_barrier(params->threadpool); + + const float global_max = thread_max[0]; + const float id = global_max != 0.0f ? 127.0f / global_max : 0.0f; + + int8_t * dst_data = (int8_t *) dst->data; + + for (int64_t row = 0; row < n_rows; row++) { + for (int64_t col = col0; col < col1; col++) { + const int64_t i = row*row_stride + col; + + float v = stage[i] * id; + v = v < -127.0f ? -127.0f : (v > 127.0f ? 127.0f : v); + + int8_t q = (int8_t) roundf(v); + dst_data[i] = relu && q < 0 ? 0 : q; + } + } + + if (ith == 0) { + // Multiplier, so that dequantizing is q * scale like every other ggml + // quantized type -- see ggml_i8_s_to_float. + *ggml_inband_scale(dst) = global_max != 0.0f ? global_max / 127.0f : 0.0f; + } +} + +// ggml_compute_forward_add_scaled + +void ggml_compute_forward_add_scaled( + const ggml_compute_params * params, + ggml_tensor * dst) { + + const ggml_tensor * src0 = dst->src[0]; // a (I8_S) + const ggml_tensor * src1 = dst->src[1]; // b (I8_S) + const ggml_tensor * src2 = dst->src[2]; // scale (F32, per channel) + + GGML_ASSERT(src0->type == GGML_TYPE_I8_S); + GGML_ASSERT(src1->type == GGML_TYPE_I8_S); + GGML_ASSERT(dst->type == GGML_TYPE_I8_S); + GGML_ASSERT(ggml_are_same_shape(src0, src1)); + GGML_ASSERT(ggml_are_same_shape(src0, dst)); + GGML_ASSERT(ggml_is_contiguous(src0)); + GGML_ASSERT(ggml_is_contiguous(src1)); + + const int ith = params->ith; + const int nth = params->nth; + + const int64_t n = ggml_nelements(src0); + const int64_t ne0 = src0->ne[0]; + + const int8_t * a = (const int8_t *) src0->data; + const int8_t * b = (const int8_t *) src1->data; + + const float a_d = *ggml_inband_scale_const(src0); + const float b_d = *ggml_inband_scale_const(src1); + + // Channel index is i % ne0: the graph is channel-major, so ne[0] is the + // channel axis and one coefficient covers every position of that channel. + const float * gamma = (const float *) src2->data; + + float * stage = (float *) params->wdata; + float * thread_max = stage + n; + + const int64_t dr = (n + nth - 1)/nth; + const int64_t i0 = dr*ith; + const int64_t i1 = MIN(i0 + dr, n); + + float local_absmax = 0.0f; + + for (int64_t i = i0; i < i1; i++) { + const float v = (float) a[i] * a_d * gamma[i % ne0] + (float) b[i] * b_d; + + stage[i] = v; + + const float av = fabsf(v); + if (av > local_absmax) local_absmax = av; + } + + ggml_i8_s_requantize(params, dst, stage, thread_max, 0, 1, i0, i1, local_absmax, false); +} + +// ggml_compute_forward_rms_norm_scaled + +void ggml_compute_forward_rms_norm_scaled( + const ggml_compute_params * params, + ggml_tensor * dst) { + + const ggml_tensor * src0 = dst->src[0]; // a (I8_S) + const ggml_tensor * src1 = dst->src[1]; // scale (F32, per channel) + + GGML_ASSERT(src0->type == GGML_TYPE_I8_S); + GGML_ASSERT(dst->type == GGML_TYPE_I8_S); + GGML_ASSERT(ggml_are_same_shape(src0, dst)); + GGML_ASSERT(ggml_is_contiguous(src0)); + + const int ith = params->ith; + const int nth = params->nth; + + const int64_t ne00 = src0->ne[0]; + const int64_t n = ggml_nelements(src0); + const int64_t nr = n / ne00; + + float eps; + memcpy(&eps, dst->op_params, sizeof(float)); + + const int8_t * x = (const int8_t *) src0->data; + const float * gamma = (const float *) src1->data; + + const float x_d = *ggml_inband_scale_const(src0); + + // The normalization is scale-invariant, so the input scale cancels between + // numerator and denominator and never has to be applied to the int8 values: + // + // real[i] = q[i] * x_d + // out[i] = real[i] / sqrt(mean(real^2) + eps) * gamma[i] + // = q[i] / sqrt(mean(q^2) + eps/x_d^2) * gamma[i] + // + // which keeps the sum of squares in exact int64 arithmetic. Only eps has to + // move into the int8 domain. + const float eps_q = (float) ((double) eps / ((double) x_d * (double) x_d)); + + float * stage = (float *) params->wdata; + float * thread_max = stage + n; + + const int64_t dr = (nr + nth - 1)/nth; + const int64_t ir0 = dr*ith; + const int64_t ir1 = MIN(ir0 + dr, nr); + + float local_absmax = 0.0f; + + for (int64_t ir = ir0; ir < ir1; ir++) { + const int8_t * x_row = x + ir*ne00; + float * stage_row = stage + ir*ne00; + + int64_t sum_sq = 0; + for (int64_t i = 0; i < ne00; i++) { + sum_sq += (int32_t) x_row[i] * (int32_t) x_row[i]; + } + + const float rms_inv = 1.0f/sqrtf((float) sum_sq/(float) ne00 + eps_q); + + for (int64_t i = 0; i < ne00; i++) { + const float v = (float) x_row[i] * rms_inv * gamma[i]; + + stage_row[i] = v; + + const float av = fabsf(v); + if (av > local_absmax) local_absmax = av; + } + } + + ggml_i8_s_requantize(params, dst, stage, thread_max, 0, 1, ir0*ne00, ir1*ne00, local_absmax, false); +} + +// ggml_compute_forward_mul_mat_add + +// Shared by GGML_OP_MUL_MAT_ADD and GGML_OP_MUL_MAT_ADD_RELU. +// +// Two shapes are handled. When src0 is [K, 1, C] the op is depthwise: C +// independent length-K dot products per position, with src1 laid out +// [K, N, C] and the output [C, N] channel-major. Otherwise src0 is an ordinary +// [IC, OC] weight matrix, src1 is [IC, N], and the output is [OC, N]. +static void ggml_compute_forward_mul_mat_add_impl( + const ggml_compute_params * params, + ggml_tensor * dst, + bool relu) { + + const ggml_tensor * src0 = dst->src[0]; // weight (I8_S) + const ggml_tensor * src1 = dst->src[1]; // input (I8_S) + const ggml_tensor * src2 = dst->src[2]; // bias (F32) + + GGML_ASSERT(src0->type == GGML_TYPE_I8_S); + GGML_ASSERT(src1->type == GGML_TYPE_I8_S); + GGML_ASSERT(dst->type == GGML_TYPE_I8_S); + GGML_ASSERT(ggml_is_contiguous(src0)); + GGML_ASSERT(ggml_is_contiguous(src1)); + + const int ith = params->ith; + const int nth = params->nth; + + const int64_t ne00 = src0->ne[0]; // IC, or K in the depthwise case + const int64_t ne01 = src0->ne[1]; // OC, or 1 + const int64_t ne02 = src0->ne[2]; // 1, or C + const int64_t ne11 = src1->ne[1]; // N + + const float * bias = (const float *) src2->data; + + // Both scales are multipliers, so they combine by multiplication and the + // int32 accumulator only has to be scaled once. + const float combined_d = *ggml_inband_scale_const(src0) * *ggml_inband_scale_const(src1); + + const int8_t * w = (const int8_t *) src0->data; + const int8_t * x = (const int8_t *) src1->data; + + const bool depthwise = ne01 == 1 && ne02 > 1; + + // Output channel count, and the number of outputs per column. + const int64_t OC = depthwise ? ne02 : ne01; + + // The loops below walk a single batch. Batched matmul has no caller in the + // VibeASR graphs, so refuse it here rather than silently leaving + // dst->ne[2]/ne[3] beyond the first as uninitialized garbage. In the + // depthwise case ne[2] is the channel axis, which is handled. + GGML_ASSERT(dst->ne[3] == 1); + GGML_ASSERT(depthwise || dst->ne[2] == 1); + GGML_ASSERT(ne11*OC == ggml_nelements(dst)); + + float * stage = (float *) params->wdata; + float * thread_max = stage + ne11*OC; + + // Split over columns: every thread owns whole columns of the output, so the + // int32 accumulators need no cross-thread coordination. + const int64_t dr = (ne11 + nth - 1)/nth; + const int64_t col0 = dr*ith; + const int64_t col1 = MIN(col0 + dr, ne11); + + float local_absmax = 0.0f; + + if (depthwise) { + // stage/dst index is ch*ne11 + col, so each channel is a row of length + // ne11 and this thread owns the [col0, col1) band of every row. + for (int64_t ch = 0; ch < ne02; ch++) { + const int8_t * w_ch = w + ch*ne00; + + for (int64_t col = col0; col < col1; col++) { + // src1 is [K, N, C]: the channel stride is ne11*ne00. + const int8_t * x_col = x + ch*ne11*ne00 + col*ne00; + + int32_t acc = 0; + for (int64_t k = 0; k < ne00; k++) { + acc += (int32_t) w_ch[k] * (int32_t) x_col[k]; + } + + const float v = (float) acc * combined_d + bias[ch]; + + stage[ch*ne11 + col] = v; + + const float av = fabsf(v); + if (av > local_absmax) local_absmax = av; + } + } + + ggml_i8_s_requantize(params, dst, stage, thread_max, ne11, ne02, col0, col1, local_absmax, relu); + } else { + // stage/dst index is col*ne01 + oc, so this thread owns the contiguous + // block [col0*ne01, col1*ne01). + const size_t nb11 = src1->nb[1]; + + for (int64_t col = col0; col < col1; col++) { + const int8_t * x_col = (const int8_t *) ((const char *) src1->data + col*nb11); + + for (int64_t oc = 0; oc < ne01; oc++) { + const int8_t * w_row = w + oc*ne00; + + int32_t acc = 0; + for (int64_t k = 0; k < ne00; k++) { + acc += (int32_t) w_row[k] * (int32_t) x_col[k]; + } + + const float v = (float) acc * combined_d + bias[oc]; + + stage[col*ne01 + oc] = v; + + const float av = fabsf(v); + if (av > local_absmax) local_absmax = av; + } + } + + ggml_i8_s_requantize(params, dst, stage, thread_max, 0, 1, col0*ne01, col1*ne01, local_absmax, relu); + } +} + +void ggml_compute_forward_mul_mat_add( + const ggml_compute_params * params, + ggml_tensor * dst) { + ggml_compute_forward_mul_mat_add_impl(params, dst, false); +} + +void ggml_compute_forward_mul_mat_add_relu( + const ggml_compute_params * params, + ggml_tensor * dst) { + ggml_compute_forward_mul_mat_add_impl(params, dst, true); +} + +// ggml_compute_forward_im2col_asym + +void ggml_compute_forward_im2col_asym( + const ggml_compute_params * params, + ggml_tensor * dst) { + + const ggml_tensor * src0 = dst->src[0]; // kernel, for its shape only + const ggml_tensor * src1 = dst->src[1]; // input (I8_S) + + GGML_ASSERT(src1->type == GGML_TYPE_I8_S); + GGML_ASSERT(dst->type == GGML_TYPE_I8_S); + + const int32_t s0 = ggml_get_op_params_i32(dst, 0); + const int32_t lp0 = ggml_get_op_params_i32(dst, 2); + const int32_t d0 = ggml_get_op_params_i32(dst, 4); + const bool is_2D = ggml_get_op_params_i32(dst, 6) == 1; + + // 1D only. The 2D case has no caller, and guessing at its indexing would + // mean shipping an untested path. + GGML_ASSERT(!is_2D && "ggml_im2col_asym: only the 1D case is implemented"); + + const int ith = params->ith; + const int nth = params->nth; + + const int64_t IC = src0->ne[1]; + const int64_t KW = src0->ne[0]; + const int64_t IW = src1->ne[0]; + const int64_t N = src1->ne[2]; + const int64_t OW = dst->ne[1]; + + const int64_t total = N*OW; + const int64_t per_thread = (total + nth - 1)/nth; + const int64_t start = per_thread*ith; + const int64_t end = MIN(start + per_thread, total); + + int8_t * dst_data = (int8_t *) dst->data; + + for (int64_t idx = start; idx < end; idx++) { + const int64_t in = idx/OW; + const int64_t iow = idx%OW; + + int8_t * dst_col = dst_data + idx*(IC*KW); + + for (int64_t iic = 0; iic < IC; iic++) { + const int8_t * src_ch = (const int8_t *) ((const char *) src1->data + + iic*src1->nb[1] + in*src1->nb[2]); + + for (int64_t ikw = 0; ikw < KW; ikw++) { + const int64_t iiw = iow*s0 + ikw*d0 - lp0; + + // Zero is exact in both directions, so padding needs no + // special handling when the scale is applied later. + dst_col[iic*KW + ikw] = (iiw < 0 || iiw >= IW) ? 0 : src_ch[iiw]; + } + } + } + + if (ith == 0) { + // Pure rearrangement: no value changes, so the scale carries over. + *ggml_inband_scale(dst) = *ggml_inband_scale_const(src1); + } +} diff --git a/external/ggml/src/ggml-cpu/ops.h b/external/ggml/src/ggml-cpu/ops.h index e0b3a7bd2..4783a0134 100644 --- a/external/ggml/src/ggml-cpu/ops.h +++ b/external/ggml/src/ggml-cpu/ops.h @@ -115,6 +115,13 @@ void ggml_compute_forward_opt_step_adamw(const struct ggml_compute_params * para void ggml_compute_forward_mul_mat(const struct ggml_compute_params * params, struct ggml_tensor * dst); void ggml_compute_forward_fwht(const struct ggml_compute_params * params, struct ggml_tensor * dst); void ggml_compute_forward_opt_step_sgd(const struct ggml_compute_params * params, struct ggml_tensor * dst); + +// VibeASR CPU INT8 pipeline +void ggml_compute_forward_add_scaled(const struct ggml_compute_params * params, struct ggml_tensor * dst); +void ggml_compute_forward_rms_norm_scaled(const struct ggml_compute_params * params, struct ggml_tensor * dst); +void ggml_compute_forward_mul_mat_add(const struct ggml_compute_params * params, struct ggml_tensor * dst); +void ggml_compute_forward_mul_mat_add_relu(const struct ggml_compute_params * params, struct ggml_tensor * dst); +void ggml_compute_forward_im2col_asym(const struct ggml_compute_params * params, struct ggml_tensor * dst); #ifdef __cplusplus } #endif diff --git a/external/ggml/src/ggml-impl.h b/external/ggml/src/ggml-impl.h index d4cb6c9b3..2f50da73a 100644 --- a/external/ggml/src/ggml-impl.h +++ b/external/ggml/src/ggml-impl.h @@ -354,6 +354,26 @@ struct ggml_cgraph ggml_graph_view(struct ggml_cgraph * cgraph, int i0, int i1); // ggml-alloc.c: true if the operation can reuse memory from its sources GGML_API bool ggml_op_can_inplace(enum ggml_op op); +// Bytes a type needs past its payload. Non-zero only for GGML_TYPE_I8_S and +// GGML_TYPE_I2_S, which append one F32 per-tensor scale; see ggml.c. +GGML_API size_t ggml_type_extra_bytes(enum ggml_type type); + +// The in-band scale of an I8_S or I2_S tensor, which sits immediately after the +// payload. ggml_nbytes already counts the padded scale, so subtracting it back +// out gives the payload end without duplicating either type's row arithmetic. +// +// The scale is a multiplier in both cases: dequantizing is q * scale. +static inline float * ggml_inband_scale(struct ggml_tensor * tensor) { + const size_t extra = ggml_type_extra_bytes(tensor->type); + GGML_ASSERT(extra > 0 && "type has no in-band scale"); + GGML_ASSERT(ggml_is_contiguous(tensor)); + return (float *) ((char *) tensor->data + ggml_nbytes(tensor) - extra); +} + +static inline const float * ggml_inband_scale_const(const struct ggml_tensor * tensor) { + return ggml_inband_scale((struct ggml_tensor *) tensor); +} + // Memory allocation diff --git a/external/ggml/src/ggml-quants.c b/external/ggml/src/ggml-quants.c index 51f08fa3b..f717d78f8 100644 --- a/external/ggml/src/ggml-quants.c +++ b/external/ggml/src/ggml-quants.c @@ -2407,6 +2407,138 @@ void dequantize_row_tq2_0(const block_tq2_0 * GGML_RESTRICT x, float * GGML_REST } } +// ====================== GGML_TYPE_I8_S / GGML_TYPE_I2_S +// +// Unlike every block-quantized type above, these two carry a single F32 scale +// for the whole tensor rather than one per block, stored immediately after the +// payload (ggml_type_extra_bytes reserves the room). The four functions below +// are therefore whole-tensor: `n` is ggml_nelements(), and the scale lives at +// byte offset ggml_row_size(type, n) rounded down to the payload end. +// +// The scale is a multiplier in both directions -- real = q * scale -- matching +// the `d` of every other ggml quant type. VibeASR's own fork uses a multiplier +// for weights but stores the reciprocal for op-produced activations; the two +// are reconciled by a division deep inside each kernel. Only the multiplier +// convention is used here. Activation scales never reach a file, so this +// changes nothing about how an existing GGUF is read. + +#define I2_S_GROUP 128 // ternary values per packed group +#define I2_S_BYTES 32 // bytes those 128 values occupy + +// Scale slot: the F32 sits directly after the payload. Kept as helpers so the +// offset arithmetic is not repeated in four places. The payload size is the +// only difference between the two types -- I8_S is one byte per value, I2_S +// packs 128 values into 32 bytes. +static inline size_t i8_s_payload_bytes(int64_t n) { return (size_t)n; } +static inline size_t i2_s_payload_bytes(int64_t n) { return (size_t)(n / I2_S_GROUP) * I2_S_BYTES; } + +static inline float * i8_s_scale_ptr (void * x, int64_t n) { return (float *)((char *)x + i8_s_payload_bytes(n)); } +static inline const float * i8_s_scale_ptr_const(const void * x, int64_t n) { return (const float *)((const char *)x + i8_s_payload_bytes(n)); } +static inline float * i2_s_scale_ptr (void * x, int64_t n) { return (float *)((char *)x + i2_s_payload_bytes(n)); } +static inline const float * i2_s_scale_ptr_const(const void * x, int64_t n) { return (const float *)((const char *)x + i2_s_payload_bytes(n)); } + +void ggml_i8_s_to_float(const void * GGML_RESTRICT x, float * GGML_RESTRICT y, int64_t n) { + const int8_t * q = (const int8_t *)x; + const float d = *i8_s_scale_ptr_const(x, n); + + for (int64_t i = 0; i < n; ++i) { + y[i] = (float)q[i] * d; + } +} + +size_t ggml_i8_s_from_float(const float * GGML_RESTRICT x, void * GGML_RESTRICT y, int64_t n) { + int8_t * q = (int8_t *)y; + + float amax = 0.0f; + for (int64_t i = 0; i < n; ++i) { + const float a = fabsf(x[i]); + if (a > amax) amax = a; + } + + // Symmetric range: -127..127 rather than -128..127, so that negating a + // tensor is exactly representable and the kernels can widen to int16 + // without a special case for the one asymmetric value. + const float d = amax / 127.0f; + const float id = d != 0.0f ? 1.0f / d : 0.0f; + + for (int64_t i = 0; i < n; ++i) { + int v = nearest_int(x[i] * id); + if (v > 127) v = 127; + if (v < -127) v = -127; + q[i] = (int8_t)v; + } + + *i8_s_scale_ptr(y, n) = d; + + return i8_s_payload_bytes(n) + ggml_type_extra_bytes(GGML_TYPE_I8_S); +} + +// Bit pair -> ternary value. Codes 1 and 3 both mean zero; the packer only ever +// emits 1, but 3 is accepted so a stray pattern dequantizes to 0 instead of +// reading out of bounds. +static const float i2_s_map[4] = { -1.0f, 0.0f, +1.0f, 0.0f }; + +void ggml_i2_s_to_float(const void * GGML_RESTRICT x, float * GGML_RESTRICT y, int64_t n) { + assert(n % I2_S_GROUP == 0); + + const uint8_t * q = (const uint8_t *)x; + const float d = *i2_s_scale_ptr_const(x, n); + + for (int64_t base = 0; base < n; base += I2_S_GROUP) { + // Byte gp holds the values at base+gp, base+32+gp, base+64+gp and + // base+96+gp, in bit pairs 6, 4, 2, 0 -- strided, not consecutive, so + // that a SIMD load of 32 bytes yields 4 aligned lanes of 32 values. + for (int gp = 0; gp < I2_S_BYTES; ++gp) { + const uint8_t b = q[gp]; + + y[base + 0 + gp] = d * i2_s_map[(b >> 6) & 3]; + y[base + 32 + gp] = d * i2_s_map[(b >> 4) & 3]; + y[base + 64 + gp] = d * i2_s_map[(b >> 2) & 3]; + y[base + 96 + gp] = d * i2_s_map[(b >> 0) & 3]; + } + + q += I2_S_BYTES; + } +} + +size_t ggml_i2_s_from_float(const float * GGML_RESTRICT x, void * GGML_RESTRICT y, int64_t n) { + assert(n % I2_S_GROUP == 0); + + uint8_t * q = (uint8_t *)y; + + float amax = 0.0f; + for (int64_t i = 0; i < n; ++i) { + const float a = fabsf(x[i]); + if (a > amax) amax = a; + } + + // The type stores {-1, 0, +1} * d, so d is the absmax and a ternary input + // round-trips exactly. Non-ternary input is rounded to the nearest of the + // three representable values, which is all this type can express. + const float d = amax; + const float id = d != 0.0f ? 1.0f / d : 0.0f; + + memset(q, 0, i2_s_payload_bytes(n)); + + for (int64_t base = 0; base < n; base += I2_S_GROUP) { + for (int j = 0; j < I2_S_GROUP; ++j) { + int v = nearest_int(x[base + j] * id); + if (v > 1) v = 1; + if (v < -1) v = -1; + + const uint8_t code = (uint8_t)(v + 1); // -1,0,1 -> 0,1,2 + + q[j % I2_S_BYTES] |= (uint8_t)(code << (6 - 2*(j / I2_S_BYTES))); + } + + q += I2_S_BYTES; + } + + *i2_s_scale_ptr(y, n) = d; + + return i2_s_payload_bytes(n) + ggml_type_extra_bytes(GGML_TYPE_I2_S); +} + // ====================== "True" 2-bit (de)-quantization void dequantize_row_iq2_xxs(const block_iq2_xxs * GGML_RESTRICT x, float * GGML_RESTRICT y, int64_t k) { diff --git a/external/ggml/src/ggml-quants.h b/external/ggml/src/ggml-quants.h index 3fb6140b7..66216d6cd 100644 --- a/external/ggml/src/ggml-quants.h +++ b/external/ggml/src/ggml-quants.h @@ -102,6 +102,23 @@ GGML_API size_t quantize_q8_0(const float * GGML_RESTRICT src, void * GGML_RESTR GGML_API size_t quantize_mxfp4(const float * GGML_RESTRICT src, void * GGML_RESTRICT dst, int64_t nrows, int64_t n_per_row, const float * imatrix); GGML_API size_t quantize_nvfp4(const float * GGML_RESTRICT src, void * GGML_RESTRICT dst, int64_t nrows, int64_t n_per_row, const float * imatrix); +// GGML_TYPE_I8_S / GGML_TYPE_I2_S conversions. +// +// These are whole-tensor, not per-row: both types carry a single F32 scale for +// the entire tensor, stored immediately after the payload (see +// ggml_type_extra_bytes), so `n` is ggml_nelements() and a row pointer alone +// cannot locate the scale. That is also why the two traits entries leave +// .to_float / .from_float_ref NULL instead of pointing here, and why +// ggml_quantize_chunk does not list these types -- its +// `result == nrows * row_size` invariant cannot hold for a per-tensor scale. +// +// The from_float direction returns the number of bytes written, payload plus +// the padded scale. +GGML_API void ggml_i8_s_to_float (const void * GGML_RESTRICT x, float * GGML_RESTRICT y, int64_t n); +GGML_API size_t ggml_i8_s_from_float(const float * GGML_RESTRICT x, void * GGML_RESTRICT y, int64_t n); +GGML_API void ggml_i2_s_to_float (const void * GGML_RESTRICT x, float * GGML_RESTRICT y, int64_t n); +GGML_API size_t ggml_i2_s_from_float(const float * GGML_RESTRICT x, void * GGML_RESTRICT y, int64_t n); + GGML_API void iq2xs_init_impl(enum ggml_type type); GGML_API void iq2xs_free_impl(enum ggml_type type); GGML_API void iq3xs_init_impl(int grid_size); diff --git a/external/ggml/src/ggml.c b/external/ggml/src/ggml.c index 0ba560368..869645bdc 100644 --- a/external/ggml/src/ggml.c +++ b/external/ggml/src/ggml.c @@ -927,8 +927,49 @@ static const struct ggml_type_traits type_traits[GGML_TYPE_COUNT] = { .type_size = 0, .is_quantized = false, }, + // .to_float / .from_float_ref are deliberately left NULL for the two types + // below. Both carry a single F32 scale for the whole tensor rather than one + // per block, so a row-scoped callback has no way to find it: the scale is + // past the end of the last row, not the row it is handed. Use + // ggml_i8_s_to_float / ggml_i2_s_to_float (ggml-quants.h), which take an + // element count of ggml_nelements(). For the same reason neither type + // appears in ggml_quantize_chunk -- its `result == nrows * row_size` + // invariant cannot hold when a per-tensor scale is in play. + [GGML_TYPE_I8_S] = { + .type_name = "i8_s", + .blck_size = 1, + .type_size = sizeof(int8_t), + .is_quantized = true, + }, + [GGML_TYPE_I2_S] = { + .type_name = "i2_s", + // 128 ternary values packed into 32 bytes: byte gp holds the values at + // positions gp, 32+gp, 64+gp and 96+gp in bit pairs 6,4,2,0. A partial + // group still consumes a full 32 bytes, so the row size is + // ceil(ne0/128)*32 -- declaring the block as 128/32 makes ggml_row_size + // compute exactly that and assert the ne0 % 128 == 0 the packing needs. + .blck_size = 128, + .type_size = 32, + .is_quantized = true, + }, }; +// I8_S and I2_S append one F32 per-tensor scale after the payload. Padded to 32 +// bytes so the next tensor in an arena stays aligned. +// +// This is the only place the layout is encoded; ggml_nbytes and +// ggml_new_tensor_impl both go through it. Returns 0 for every other type, so +// no existing type changes size. +size_t ggml_type_extra_bytes(enum ggml_type type) { + switch (type) { + case GGML_TYPE_I8_S: + case GGML_TYPE_I2_S: + return 32; + default: + return 0; + } +} + const struct ggml_type_traits * ggml_get_type_traits(enum ggml_type type) { assert(type >= 0); assert(type < GGML_TYPE_COUNT); @@ -1084,9 +1125,15 @@ static const char * GGML_OP_NAME[GGML_OP_COUNT] = { "GLU", "CONVROT_LINEAR", + + "ADD_SCALED", + "RMS_NORM_SCALED", + "MUL_MAT_ADD", + "MUL_MAT_ADD_RELU", + "IM2COL_ASYM", }; -static_assert(GGML_OP_COUNT == 102, "GGML_OP_COUNT != 102"); +static_assert(GGML_OP_COUNT == 107, "GGML_OP_COUNT != 107"); static const char * GGML_OP_SYMBOL[GGML_OP_COUNT] = { "none", @@ -1200,9 +1247,15 @@ static const char * GGML_OP_SYMBOL[GGML_OP_COUNT] = { "glu(x)", "convrot_linear(weight_i8, input, weight_scale, bias)", + + "a*scale+b", + "rms_norm(x)*scale", + "a*b+bias", + "relu(a*b+bias)", + "im2col_asym(x)", }; -static_assert(GGML_OP_COUNT == 102, "GGML_OP_COUNT != 102"); +static_assert(GGML_OP_COUNT == 107, "GGML_OP_COUNT != 107"); static_assert(GGML_OP_POOL_COUNT == 2, "GGML_OP_POOL_COUNT != 2"); @@ -1303,7 +1356,7 @@ size_t ggml_nbytes(const struct ggml_tensor * tensor) { } } - return nbytes; + return nbytes + ggml_type_extra_bytes(tensor->type); } size_t ggml_nbytes_pad(const struct ggml_tensor * tensor) { @@ -1757,6 +1810,7 @@ static struct ggml_tensor * ggml_new_tensor_impl( for (int i = 1; i < n_dims; i++) { data_size *= ne[i]; } + data_size += ggml_type_extra_bytes(type); GGML_ASSERT(view_src == NULL || data_size == 0 || data_size + view_offs <= ggml_nbytes(view_src)); @@ -5594,6 +5648,161 @@ struct ggml_tensor * ggml_convrot_linear( return result; } +// VibeASR CPU INT8 pipeline +// +// Five fused ops that keep an activation chain entirely in GGML_TYPE_I8_S: each +// one consumes int8 with a per-tensor scale, does its arithmetic in int32/F32, +// and re-quantizes its own output to int8 with a freshly measured scale. Fusing +// matters because the requantization needs the output absmax, so an unfused +// chain would have to write, read back and rescan every intermediate. +// +// That per-tensor output scale is why these need GGML_TYPE_I8_S rather than the +// GGML_TYPE_I8-plus-separate-scale-tensor pair used by ggml_convrot_linear +// above: a scale computed from the op's own output has nowhere to live but with +// the data, since a node has exactly one output tensor. + +struct ggml_tensor * ggml_add_scaled( + struct ggml_context * ctx, + struct ggml_tensor * a, + struct ggml_tensor * b, + struct ggml_tensor * scale) { + GGML_ASSERT(ggml_are_same_shape(a, b)); + GGML_ASSERT(scale->type == GGML_TYPE_F32); + // Broadcast along ne[0], the channel axis: one LayerScale coefficient per + // channel, applied to every position. + GGML_ASSERT(scale->ne[0] == a->ne[0]); + + struct ggml_tensor * result = ggml_dup_tensor(ctx, a); + + result->op = GGML_OP_ADD_SCALED; + result->src[0] = a; + result->src[1] = b; + result->src[2] = scale; + + return result; +} + +struct ggml_tensor * ggml_rms_norm_scaled( + struct ggml_context * ctx, + struct ggml_tensor * a, + struct ggml_tensor * scale, + float eps) { + GGML_ASSERT(a->type == GGML_TYPE_I8_S); + GGML_ASSERT(scale->type == GGML_TYPE_F32); + GGML_ASSERT(scale->ne[0] == a->ne[0]); + GGML_ASSERT(eps > 0.0f); + + struct ggml_tensor * result = ggml_dup_tensor(ctx, a); + + ggml_set_op_params(result, &eps, sizeof(eps)); + + result->op = GGML_OP_RMS_NORM_SCALED; + result->src[0] = a; + result->src[1] = scale; + + return result; +} + +// Shared by ggml_mul_mat_add and ggml_mul_mat_add_relu: identical shape rules, +// only the epilogue differs. +static struct ggml_tensor * ggml_mul_mat_add_impl( + struct ggml_context * ctx, + struct ggml_tensor * a, + struct ggml_tensor * b, + struct ggml_tensor * bias, + enum ggml_op op) { + GGML_ASSERT(ggml_can_mul_mat(a, b)); + // I8_S weights only. The ternary I2_S weights of the language model go + // through plain ggml_mul_mat, which has no bias to fuse. + GGML_ASSERT(a->type == GGML_TYPE_I8_S); + GGML_ASSERT(b->type == GGML_TYPE_I8_S); + GGML_ASSERT(bias->type == GGML_TYPE_F32); + + // Bias is per output channel, and which ne holds the output channels + // depends on the shape of a: a [IC, OC] is an ordinary matmul, whereas + // a [K, 1, C] is the depthwise case where the channels sit in ne[2]. + GGML_ASSERT(bias->ne[0] == (a->ne[1] == 1 && a->ne[2] > 1 ? a->ne[2] : a->ne[1])); + + const int64_t ne[4] = { a->ne[1], b->ne[1], b->ne[2], b->ne[3] }; + + // Output is int8 with its own scale, so the chain continues without a + // detour through F32. + struct ggml_tensor * result = ggml_new_tensor(ctx, GGML_TYPE_I8_S, 4, ne); + + result->op = op; + result->src[0] = a; + result->src[1] = b; + result->src[2] = bias; + + return result; +} + +struct ggml_tensor * ggml_mul_mat_add( + struct ggml_context * ctx, + struct ggml_tensor * a, + struct ggml_tensor * b, + struct ggml_tensor * bias) { + return ggml_mul_mat_add_impl(ctx, a, b, bias, GGML_OP_MUL_MAT_ADD); +} + +struct ggml_tensor * ggml_mul_mat_add_relu( + struct ggml_context * ctx, + struct ggml_tensor * a, + struct ggml_tensor * b, + struct ggml_tensor * bias) { + return ggml_mul_mat_add_impl(ctx, a, b, bias, GGML_OP_MUL_MAT_ADD_RELU); +} + +struct ggml_tensor * ggml_im2col_asym( + struct ggml_context * ctx, + struct ggml_tensor * a, + struct ggml_tensor * b, + int s0, + int s1, + int lp0, + int rp0, + int p1, + int d0, + int d1, + bool is_2D, + enum ggml_type dst_type) { + if (is_2D) { + GGML_ASSERT(a->ne[2] == b->ne[2]); + } else { + GGML_ASSERT(a->ne[1] == b->ne[1]); + GGML_ASSERT(b->ne[3] == 1); + } + + const int64_t OH = is_2D ? ggml_calc_conv_output_size(b->ne[1], a->ne[1], s1, p1, d1) : 0; + // Same formula as ggml_calc_conv_output_size but with the two width pads + // counted separately instead of as 2*p0. + const int64_t OW = (b->ne[0] + lp0 + rp0 - d0 * (a->ne[0] - 1) - 1) / s0 + 1; + + GGML_ASSERT((!is_2D || OH > 0) && "b too small compared to a"); + GGML_ASSERT((OW > 0) && "b too small compared to a"); + + const int64_t ne[4] = { + is_2D ? (a->ne[2] * a->ne[1] * a->ne[0]) : a->ne[1] * a->ne[0], + OW, + is_2D ? OH : b->ne[2], + is_2D ? b->ne[3] : 1, + }; + + struct ggml_tensor * result = ggml_new_tensor(ctx, dst_type, 4, ne); + + // Note the order: lp0 and rp0 occupy slots 2 and 3, where GGML_OP_IM2COL + // keeps p0 and p1. The two ops read their own params and never share a + // forward, so the layouts do not have to agree. + int32_t params[] = { s0, s1, lp0, rp0, d0, d1, (is_2D ? 1 : 0), p1 }; + ggml_set_op_params(result, params, sizeof(params)); + + result->op = GGML_OP_IM2COL_ASYM; + result->src[0] = a; + result->src[1] = b; + + return result; +} + struct ggml_tensor * ggml_flash_attn_ext_with_bias_mask( struct ggml_context * ctx, struct ggml_tensor * q, diff --git a/tests/unittests/test_i8_s_fused_ops.cpp b/tests/unittests/test_i8_s_fused_ops.cpp new file mode 100644 index 000000000..ad2148c43 --- /dev/null +++ b/tests/unittests/test_i8_s_fused_ops.cpp @@ -0,0 +1,487 @@ +// Numeric checks for the VibeASR INT8 pipeline additions to ggml: +// GGML_TYPE_I8_S / GGML_TYPE_I2_S and the five fused ops built on them. +// +// Every op is compared against a plain-loop reference computed from the same +// inputs. The references duplicate the requantization ordering deliberately, +// including the parts that lose range -- ggml_mul_mat_add_relu clamps after +// rounding, so the absmax that sets the output scale still counts the negatives +// that are about to become zero. The point of these tests is to pin the +// arithmetic that the SIMD kernels and the model port have to reproduce, not to +// judge whether that arithmetic is optimal. + +#include "test_assert.h" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +// The in-band scale layout is internal to ggml, so the tests reach for the same +// declarations the implementation uses rather than re-deriving the offsets. +extern "C" { +size_t ggml_type_extra_bytes(enum ggml_type type); +void ggml_i8_s_to_float (const void * x, float * y, int64_t n); +size_t ggml_i8_s_from_float(const float * x, void * y, int64_t n); +void ggml_i2_s_to_float (const void * x, float * y, int64_t n); +size_t ggml_i2_s_from_float(const float * x, void * y, int64_t n); +} + +using engine::test::require; +using engine::test::require_close; +using engine::test::require_eq; + +namespace { + +constexpr size_t kCtxBytes = 256u * 1024 * 1024; + +float * inband_scale(ggml_tensor * t) { + return reinterpret_cast( + static_cast(t->data) + ggml_nbytes(t) - ggml_type_extra_bytes(t->type)); +} + +// Quantize into an existing I8_S tensor, returning the scale that was chosen. +float fill_i8_s(ggml_tensor * t, const std::vector & values) { + require_eq(static_cast(values.size()), ggml_nelements(t), "fill_i8_s size"); + ggml_i8_s_from_float(values.data(), t->data, ggml_nelements(t)); + return *inband_scale(t); +} + +std::vector read_i8_s(ggml_tensor * t) { + std::vector out(ggml_nelements(t)); + ggml_i8_s_to_float(t->data, out.data(), ggml_nelements(t)); + return out; +} + +std::vector patterned(size_t n, float phase, float scale) { + std::vector v(n); + for (size_t i = 0; i < n; ++i) { + const float x = static_cast(i); + v[i] = scale * (std::sin(phase + 0.13f * x) + 0.4f * std::cos(0.07f * x - phase)); + } + return v; +} + +// The reference requantizer: absmax over the whole result, then round, then +// clamp. Mirrors ggml_i8_s_requantize. +std::vector requantize_ref(const std::vector & values, bool relu, float * scale_out) { + float amax = 0.0f; + for (float v : values) { + amax = std::max(amax, std::fabs(v)); + } + + const float id = amax != 0.0f ? 127.0f / amax : 0.0f; + + std::vector q(values.size()); + for (size_t i = 0; i < values.size(); ++i) { + float v = values[i] * id; + v = std::max(-127.0f, std::min(127.0f, v)); + const int8_t r = static_cast(std::round(v)); + q[i] = relu && r < 0 ? 0 : r; + } + + *scale_out = amax != 0.0f ? amax / 127.0f : 0.0f; + return q; +} + +void compare_i8(ggml_tensor * got, const std::vector & want, float want_scale, + const std::string & label) { + require_eq(static_cast(want.size()), ggml_nelements(got), label + " size"); + + const auto * q = static_cast(got->data); + for (size_t i = 0; i < want.size(); ++i) { + // Exact: both sides do the same rounding on the same floats. A single + // LSB of drift would mean the scale differs, which is what matters. + require_eq(static_cast(q[i]), static_cast(want[i]), + label + " value at " + std::to_string(i)); + } + + require_close(*inband_scale(got), want_scale, 1e-9f, label + " scale"); +} + +// Runs a single-node graph on the CPU backend with the given thread count. The +// four requantizing ops reduce an absmax across threads, so the result has to be +// independent of nth -- these tests run every case at 1 and 4. +void compute(ggml_context * ctx, ggml_tensor * result, int n_threads) { + ggml_cgraph * gf = ggml_new_graph(ctx); + ggml_build_forward_expand(gf, result); + ggml_graph_compute_with_ctx(ctx, gf, n_threads); +} + +// ---------------------------------------------------------------- round trips + +void test_i8_s_round_trip() { + const int64_t n = 4096; + const auto values = patterned(n, 0.3f, 2.5f); + + std::vector buf(n + ggml_type_extra_bytes(GGML_TYPE_I8_S)); + const size_t written = ggml_i8_s_from_float(values.data(), buf.data(), n); + require_eq(written, buf.size(), "i8_s written bytes"); + + std::vector back(n); + ggml_i8_s_to_float(buf.data(), back.data(), n); + + float amax = 0.0f; + for (float v : values) { + amax = std::max(amax, std::fabs(v)); + } + // One step of the quantization grid is amax/127; rounding puts every value + // within half of that. + const float tol = amax / 127.0f * 0.5f + 1e-6f; + + for (int64_t i = 0; i < n; ++i) { + require_close(back[i], values[i], tol, "i8_s round trip at " + std::to_string(i)); + } + + std::cout << "i8_s round trip: n=" << n << " step=" << amax / 127.0f << " OK\n"; +} + +void test_i2_s_round_trip() { + // Ternary input: the type represents {-d, 0, +d} exactly, so this has to + // round trip bit for bit rather than approximately. + const int64_t n = 128 * 7; + const float d = 0.0731f; + + std::vector values(n); + for (int64_t i = 0; i < n; ++i) { + const int code = static_cast(i * 7 % 3) - 1; // -1, 0, +1 cycling + values[i] = static_cast(code) * d; + } + + const size_t payload = static_cast(n / 128) * 32; + std::vector buf(payload + ggml_type_extra_bytes(GGML_TYPE_I2_S)); + const size_t written = ggml_i2_s_from_float(values.data(), buf.data(), n); + require_eq(written, buf.size(), "i2_s written bytes"); + + std::vector back(n); + ggml_i2_s_to_float(buf.data(), back.data(), n); + + for (int64_t i = 0; i < n; ++i) { + require_close(back[i], values[i], 1e-7f, "i2_s round trip at " + std::to_string(i)); + } + + // Check the packing itself, not just the round trip: byte gp of a group + // holds positions gp, 32+gp, 64+gp, 96+gp in bit pairs 6, 4, 2, 0. A + // self-consistent but differently-ordered packing would pass the round trip + // above while being unreadable by the SIMD kernels. + for (int gp = 0; gp < 32; ++gp) { + const uint8_t b = buf[gp]; + for (int sub = 0; sub < 4; ++sub) { + const int code = (b >> (6 - 2 * sub)) & 3; + const int want = static_cast((sub * 32 + gp) * 7 % 3) - 1; + const float value = static_cast(code - 1); + require_close(value, static_cast(want), 1e-7f, + "i2_s bit layout at byte " + std::to_string(gp) + + " pair " + std::to_string(sub)); + } + } + + std::cout << "i2_s round trip: n=" << n << " exact, bit layout OK\n"; +} + +// ---------------------------------------------------------------------- ops + +void test_add_scaled(int n_threads) { + const int64_t C = 96; // channels, ne[0] + const int64_t L = 37; // positions + + ggml_init_params ip = { kCtxBytes, nullptr, false }; + ggml_context * ctx = ggml_init(ip); + + ggml_tensor * a = ggml_new_tensor_2d(ctx, GGML_TYPE_I8_S, C, L); + ggml_tensor * b = ggml_new_tensor_2d(ctx, GGML_TYPE_I8_S, C, L); + ggml_tensor * gamma = ggml_new_tensor_1d(ctx, GGML_TYPE_F32, C); + + const auto a_f = patterned(C * L, 0.1f, 1.7f); + const auto b_f = patterned(C * L, 1.9f, 0.6f); + + fill_i8_s(a, a_f); + fill_i8_s(b, b_f); + + const auto g = patterned(C, 0.5f, 0.25f); + std::memcpy(gamma->data, g.data(), g.size() * sizeof(float)); + + ggml_tensor * result = ggml_add_scaled(ctx, a, b, gamma); + compute(ctx, result, n_threads); + + // Reference works from the dequantized inputs, so it sees exactly what the + // op sees -- the input quantization error is shared, not compounded. + const auto a_q = read_i8_s(a); + const auto b_q = read_i8_s(b); + + std::vector want(C * L); + for (int64_t i = 0; i < C * L; ++i) { + want[i] = a_q[i] * g[i % C] + b_q[i]; + } + + float want_scale = 0.0f; + const auto want_q = requantize_ref(want, false, &want_scale); + compare_i8(result, want_q, want_scale, "add_scaled nth=" + std::to_string(n_threads)); + + ggml_free(ctx); + std::cout << "add_scaled: C=" << C << " L=" << L << " nth=" << n_threads << " OK\n"; +} + +void test_rms_norm_scaled(int n_threads) { + const int64_t C = 128; + const int64_t L = 29; + const float eps = 1e-5f; + + ggml_init_params ip = { kCtxBytes, nullptr, false }; + ggml_context * ctx = ggml_init(ip); + + ggml_tensor * a = ggml_new_tensor_2d(ctx, GGML_TYPE_I8_S, C, L); + ggml_tensor * gamma = ggml_new_tensor_1d(ctx, GGML_TYPE_F32, C); + + const auto a_f = patterned(C * L, 0.7f, 3.1f); + fill_i8_s(a, a_f); + + const auto g = patterned(C, 1.1f, 0.9f); + std::memcpy(gamma->data, g.data(), g.size() * sizeof(float)); + + ggml_tensor * result = ggml_rms_norm_scaled(ctx, a, gamma, eps); + compute(ctx, result, n_threads); + + // Reference normalizes in the float domain. The op instead cancels the input + // scale and keeps the sum of squares in integers, which is the same value + // computed a different way -- so this also checks that the eps rescaling is + // right, since eps is the one term that does not cancel. + const auto a_q = read_i8_s(a); + + std::vector want(C * L); + for (int64_t row = 0; row < L; ++row) { + double sum_sq = 0.0; + for (int64_t i = 0; i < C; ++i) { + const double v = a_q[row * C + i]; + sum_sq += v * v; + } + const float rms_inv = 1.0f / std::sqrt(static_cast(sum_sq / C) + eps); + for (int64_t i = 0; i < C; ++i) { + want[row * C + i] = a_q[row * C + i] * rms_inv * g[i]; + } + } + + float want_scale = 0.0f; + const auto want_q = requantize_ref(want, false, &want_scale); + + // The two paths reach the same value through different arithmetic, so a + // borderline element can round to either side of a grid step. Compare the + // dequantized results with a tolerance of one step instead of demanding + // identical bytes. + require_close(*inband_scale(result), want_scale, want_scale * 1e-4f, + "rms_norm_scaled scale nth=" + std::to_string(n_threads)); + + const auto * q = static_cast(result->data); + for (int64_t i = 0; i < C * L; ++i) { + require_close(static_cast(q[i]), static_cast(want_q[i]), 1.0f, + "rms_norm_scaled value at " + std::to_string(i)); + } + + ggml_free(ctx); + std::cout << "rms_norm_scaled: C=" << C << " L=" << L << " nth=" << n_threads << " OK\n"; +} + +void test_mul_mat_add(bool relu, int n_threads) { + const int64_t IC = 64; + const int64_t OC = 48; + const int64_t N = 23; + + ggml_init_params ip = { kCtxBytes, nullptr, false }; + ggml_context * ctx = ggml_init(ip); + + ggml_tensor * w = ggml_new_tensor_2d(ctx, GGML_TYPE_I8_S, IC, OC); + ggml_tensor * x = ggml_new_tensor_2d(ctx, GGML_TYPE_I8_S, IC, N); + ggml_tensor * bias = ggml_new_tensor_1d(ctx, GGML_TYPE_F32, OC); + + fill_i8_s(w, patterned(IC * OC, 0.2f, 0.8f)); + fill_i8_s(x, patterned(IC * N, 1.3f, 2.2f)); + + const auto b = patterned(OC, 0.9f, 0.15f); + std::memcpy(bias->data, b.data(), b.size() * sizeof(float)); + + ggml_tensor * result = relu ? ggml_mul_mat_add_relu(ctx, w, x, bias) + : ggml_mul_mat_add(ctx, w, x, bias); + compute(ctx, result, n_threads); + + require_eq(result->ne[0], OC, "mul_mat_add ne0"); + require_eq(result->ne[1], N, "mul_mat_add ne1"); + require_eq(static_cast(result->type), static_cast(GGML_TYPE_I8_S), + "mul_mat_add type"); + + const auto w_q = read_i8_s(w); + const auto x_q = read_i8_s(x); + + std::vector want(OC * N); + for (int64_t col = 0; col < N; ++col) { + for (int64_t oc = 0; oc < OC; ++oc) { + double acc = 0.0; + for (int64_t k = 0; k < IC; ++k) { + acc += static_cast(w_q[oc * IC + k]) * x_q[col * IC + k]; + } + want[col * OC + oc] = static_cast(acc) + b[oc]; + } + } + + float want_scale = 0.0f; + const auto want_q = requantize_ref(want, relu, &want_scale); + + const std::string label = std::string(relu ? "mul_mat_add_relu" : "mul_mat_add") + + " nth=" + std::to_string(n_threads); + + require_close(*inband_scale(result), want_scale, want_scale * 1e-4f, label + " scale"); + + const auto * q = static_cast(result->data); + bool saw_negative_input = false; + for (int64_t i = 0; i < OC * N; ++i) { + require_close(static_cast(q[i]), static_cast(want_q[i]), 1.0f, + label + " value at " + std::to_string(i)); + if (want[i] < 0.0f) saw_negative_input = true; + } + + if (relu) { + // Without negatives in the pre-activation the clamp would be untested. + require(saw_negative_input, label + ": inputs never went negative"); + for (int64_t i = 0; i < OC * N; ++i) { + require(q[i] >= 0, label + " left a negative at " + std::to_string(i)); + } + } + + ggml_free(ctx); + std::cout << label << ": IC=" << IC << " OC=" << OC << " N=" << N << " OK\n"; +} + +void test_mul_mat_add_depthwise(int n_threads) { + const int64_t K = 7; // kernel width + const int64_t C = 40; // channels + const int64_t N = 19; // positions + + ggml_init_params ip = { kCtxBytes, nullptr, false }; + ggml_context * ctx = ggml_init(ip); + + // ne[1] == 1 with ne[2] > 1 is what selects the depthwise path. + ggml_tensor * w = ggml_new_tensor_3d(ctx, GGML_TYPE_I8_S, K, 1, C); + ggml_tensor * x = ggml_new_tensor_3d(ctx, GGML_TYPE_I8_S, K, N, C); + ggml_tensor * bias = ggml_new_tensor_1d(ctx, GGML_TYPE_F32, C); + + fill_i8_s(w, patterned(K * C, 0.4f, 1.1f)); + fill_i8_s(x, patterned(K * N * C, 1.7f, 1.9f)); + + const auto b = patterned(C, 0.2f, 0.3f); + std::memcpy(bias->data, b.data(), b.size() * sizeof(float)); + + ggml_tensor * result = ggml_mul_mat_add(ctx, w, x, bias); + compute(ctx, result, n_threads); + + const auto w_q = read_i8_s(w); + const auto x_q = read_i8_s(x); + + // Output index is ch*N + col: channel-major, one row of N per channel. + std::vector want(C * N); + for (int64_t ch = 0; ch < C; ++ch) { + for (int64_t col = 0; col < N; ++col) { + double acc = 0.0; + for (int64_t k = 0; k < K; ++k) { + acc += static_cast(w_q[ch * K + k]) * x_q[ch * N * K + col * K + k]; + } + want[ch * N + col] = static_cast(acc) + b[ch]; + } + } + + float want_scale = 0.0f; + const auto want_q = requantize_ref(want, false, &want_scale); + + const std::string label = "mul_mat_add depthwise nth=" + std::to_string(n_threads); + require_close(*inband_scale(result), want_scale, want_scale * 1e-4f, label + " scale"); + + const auto * q = static_cast(result->data); + for (int64_t i = 0; i < C * N; ++i) { + require_close(static_cast(q[i]), static_cast(want_q[i]), 1.0f, + label + " value at " + std::to_string(i)); + } + + ggml_free(ctx); + std::cout << label << ": K=" << K << " C=" << C << " N=" << N << " OK\n"; +} + +void test_im2col_asym(int n_threads) { + const int64_t IW = 33; + const int64_t IC = 5; + const int64_t KW = 4; + const int s0 = 2; + const int lp0 = 3; // asymmetric on purpose: the whole reason this op + const int rp0 = 0; // exists is that ggml_im2col cannot express it + const int d0 = 1; + + ggml_init_params ip = { kCtxBytes, nullptr, false }; + ggml_context * ctx = ggml_init(ip); + + ggml_tensor * kernel = ggml_new_tensor_3d(ctx, GGML_TYPE_I8_S, KW, IC, 1); + ggml_tensor * input = ggml_new_tensor_3d(ctx, GGML_TYPE_I8_S, IW, IC, 1); + + fill_i8_s(kernel, patterned(KW * IC, 0.1f, 1.0f)); + const float in_scale = fill_i8_s(input, patterned(IW * IC, 0.8f, 1.4f)); + + ggml_tensor * result = ggml_im2col_asym(ctx, kernel, input, s0, 0, lp0, rp0, 0, + d0, 0, false, GGML_TYPE_I8_S); + + const int64_t OW = (IW + lp0 + rp0 - d0 * (KW - 1) - 1) / s0 + 1; + require_eq(result->ne[0], IC * KW, "im2col_asym ne0"); + require_eq(result->ne[1], OW, "im2col_asym ne1"); + + compute(ctx, result, n_threads); + + const auto * in = static_cast(input->data); + const auto * out = static_cast(result->data); + + for (int64_t iow = 0; iow < OW; ++iow) { + for (int64_t iic = 0; iic < IC; ++iic) { + for (int64_t ikw = 0; ikw < KW; ++ikw) { + const int64_t iiw = iow * s0 + ikw * d0 - lp0; + const int8_t want = (iiw < 0 || iiw >= IW) ? 0 : in[iic * IW + iiw]; + require_eq(static_cast(out[iow * (IC * KW) + iic * KW + ikw]), + static_cast(want), + "im2col_asym at ow=" + std::to_string(iow) + + " ic=" + std::to_string(iic) + " kw=" + std::to_string(ikw)); + } + } + } + + // Rearrangement only, so the scale must pass through untouched. + require_close(*inband_scale(result), in_scale, 0.0f, "im2col_asym scale passthrough"); + + ggml_free(ctx); + std::cout << "im2col_asym: IW=" << IW << " KW=" << KW << " lp0=" << lp0 + << " OW=" << OW << " nth=" << n_threads << " OK\n"; +} + +} // namespace + +int main() { + try { + test_i8_s_round_trip(); + test_i2_s_round_trip(); + + // Every requantizing op reduces an absmax across threads, so each is run + // single-threaded and multi-threaded; a missing barrier shows up as a + // scale that depends on the thread count. + for (int nth : {1, 4}) { + test_add_scaled(nth); + test_rms_norm_scaled(nth); + test_mul_mat_add(false, nth); + test_mul_mat_add(true, nth); + test_mul_mat_add_depthwise(nth); + test_im2col_asym(nth); + } + } catch (const std::exception & e) { + std::cerr << "FAILED: " << e.what() << "\n"; + return 1; + } + + std::cout << "all i8_s fused op tests passed\n"; + return 0; +} From 42dd4a67c8f0d03a7cb6f1fd2d6af5f18099f9d0 Mon Sep 17 00:00:00 2001 From: XSquirrelC Date: Fri, 4 Sep 2026 02:30:13 +0000 Subject: [PATCH 2/3] ggml-cpu: SIMD kernels for the I8_S fused ops --- external/ggml/src/ggml-cpu/ops.cpp | 205 ++++++++++++++++++++---- external/ggml/src/ggml-cpu/vec.cpp | 109 +++++++++++++ external/ggml/src/ggml-cpu/vec.h | 7 + tests/unittests/test_i8_s_fused_ops.cpp | 56 +++++-- 4 files changed, 329 insertions(+), 48 deletions(-) diff --git a/external/ggml/src/ggml-cpu/ops.cpp b/external/ggml/src/ggml-cpu/ops.cpp index fc2ae9177..735030a78 100644 --- a/external/ggml/src/ggml-cpu/ops.cpp +++ b/external/ggml/src/ggml-cpu/ops.cpp @@ -11671,6 +11671,71 @@ void ggml_compute_forward_fwht(const ggml_compute_params * params, ggml_tensor * // of the available 127. That is what the reference implementation does, and // matching it is what makes a layer-by-layer comparison meaningful, so it is // reproduced here deliberately rather than improved in passing. +// Scale a contiguous run of staged F32 into I8_S. +// +// relu is folded into the low clamp: clamping to 0 before rounding and rounding +// before clamping to 0 agree on every negative input, and the float clamp is +// free here. The absmax that produced id was taken before the clamp, so a +// channel about to be zeroed still counts towards the scale. +// +// roundf() would be a PLT call per element. rintf() inlines to one instruction +// and rounds ties to even, which is what _mm256_cvtps_epi32 does as well, so the +// vector body and the scalar tail agree - and it matches ggml's own +// nearest_int(). VibeASR rounds ties away from zero in its scalar tail but to +// even in its vector body, so no tie convention reproduces it exactly. +static inline void ggml_i8_s_quantize_range( + int8_t * GGML_RESTRICT dst, + const float * GGML_RESTRICT src, + int64_t n, + float id, + bool relu) { + + int64_t i = 0; + + const float lo = relu ? 0.0f : -127.0f; + +#if defined(__AVX2__) + const __m256 v_id = _mm256_set1_ps(id); + const __m256 v_lo = _mm256_set1_ps(lo); + const __m256 v_hi = _mm256_set1_ps(127.0f); + + for (; i + 8 <= n; i += 8) { + __m256 vf = _mm256_mul_ps(_mm256_loadu_ps(src + i), v_id); + vf = _mm256_min_ps(_mm256_max_ps(vf, v_lo), v_hi); + + const __m256i vi32 = _mm256_cvtps_epi32(vf); + + __m256i vi16 = _mm256_permute4x64_epi64(_mm256_packs_epi32(vi32, vi32), 0xD8); + __m256i vi8 = _mm256_permute4x64_epi64(_mm256_packs_epi16(vi16, vi16), 0xD8); + + _mm_storel_epi64((__m128i *)(dst + i), _mm256_castsi256_si128(vi8)); + } +#elif defined(__ARM_NEON) && defined(__aarch64__) + const float32x4_t v_id = vdupq_n_f32(id); + const float32x4_t v_lo = vdupq_n_f32(lo); + const float32x4_t v_hi = vdupq_n_f32(127.0f); + + for (; i + 8 <= n; i += 8) { + float32x4_t f0 = vmulq_f32(vld1q_f32(src + i ), v_id); + float32x4_t f1 = vmulq_f32(vld1q_f32(src + i + 4), v_id); + + f0 = vminq_f32(vmaxq_f32(f0, v_lo), v_hi); + f1 = vminq_f32(vmaxq_f32(f1, v_lo), v_hi); + + const int16x8_t vi16 = vcombine_s16(vqmovn_s32(vcvtnq_s32_f32(f0)), + vqmovn_s32(vcvtnq_s32_f32(f1))); + + vst1_s8(dst + i, vqmovn_s16(vi16)); + } +#endif + + for (; i < n; ++i) { + float v = src[i] * id; + v = v < lo ? lo : (v > 127.0f ? 127.0f : v); + dst[i] = (int8_t) rintf(v); + } +} + static void ggml_i8_s_requantize( const ggml_compute_params * params, ggml_tensor * dst, @@ -11705,16 +11770,12 @@ static void ggml_i8_s_requantize( int8_t * dst_data = (int8_t *) dst->data; + // Each row's [col0, col1) slice is contiguous, so the vectorized helper runs + // once per row regardless of how the rectangle is strided. for (int64_t row = 0; row < n_rows; row++) { - for (int64_t col = col0; col < col1; col++) { - const int64_t i = row*row_stride + col; - - float v = stage[i] * id; - v = v < -127.0f ? -127.0f : (v > 127.0f ? 127.0f : v); + const int64_t i = row*row_stride + col0; - int8_t q = (int8_t) roundf(v); - dst_data[i] = relu && q < 0 ? 0 : q; - } + ggml_i8_s_quantize_range(dst_data + i, stage + i, col1 - col0, id, relu); } if (ith == 0) { @@ -11854,6 +11915,81 @@ void ggml_compute_forward_rms_norm_scaled( // ggml_compute_forward_mul_mat_add +#define GGML_MUL_MAT_ADD_OC_CHUNK 64 + +// stage[j] = acc[j]*d + bias[j], returning max(|stage[j]|) folded into amax. +// +// This is the second pass over every output, so leaving it scalar costs about as +// much as the dot products do at the small contraction lengths the conv layers +// use. Vectorizing it is exact: max over floats is order-independent, and the +// multiply and add are kept separate so the tail cannot disagree with the body. +static inline float ggml_i8_s_scale_bias_absmax( + float * GGML_RESTRICT stage, + const int32_t * GGML_RESTRICT acc, + const float * GGML_RESTRICT bias, + int64_t n, + float d, + float amax) { + + int64_t i = 0; + +#if defined(__AVX2__) + const __m256 v_d = _mm256_set1_ps(d); + const __m256 v_abs = _mm256_castsi256_ps(_mm256_set1_epi32(0x7fffffff)); + + __m256 v_amax = _mm256_setzero_ps(); + + for (; i + 8 <= n; i += 8) { + const __m256 v = _mm256_add_ps( + _mm256_mul_ps(_mm256_cvtepi32_ps(_mm256_loadu_si256((const __m256i *)(acc + i))), v_d), + _mm256_loadu_ps(bias + i)); + + _mm256_storeu_ps(stage + i, v); + + v_amax = _mm256_max_ps(v_amax, _mm256_and_ps(v, v_abs)); + } + + if (i > 0) { + __m128 r = _mm_max_ps(_mm256_castps256_ps128(v_amax), _mm256_extractf128_ps(v_amax, 1)); + r = _mm_max_ps(r, _mm_movehl_ps(r, r)); + r = _mm_max_ss(r, _mm_shuffle_ps(r, r, 1)); + + const float v = _mm_cvtss_f32(r); + if (v > amax) amax = v; + } +#elif defined(__ARM_NEON) && defined(__aarch64__) + const float32x4_t v_d = vdupq_n_f32(d); + + float32x4_t v_amax = vdupq_n_f32(0.0f); + + for (; i + 4 <= n; i += 4) { + const float32x4_t v = vaddq_f32( + vmulq_f32(vcvtq_f32_s32(vld1q_s32(acc + i)), v_d), + vld1q_f32(bias + i)); + + vst1q_f32(stage + i, v); + + v_amax = vmaxq_f32(v_amax, vabsq_f32(v)); + } + + if (i > 0) { + const float v = vmaxvq_f32(v_amax); + if (v > amax) amax = v; + } +#endif + + for (; i < n; ++i) { + const float v = (float) acc[i]*d + bias[i]; + + stage[i] = v; + + const float av = fabsf(v); + if (av > amax) amax = av; + } + + return amax; +} + // Shared by GGML_OP_MUL_MAT_ADD and GGML_OP_MUL_MAT_ADD_RELU. // // Two shapes are handled. When src0 is [K, 1, C] the op is depthwise: C @@ -11919,24 +12055,32 @@ static void ggml_compute_forward_mul_mat_add_impl( if (depthwise) { // stage/dst index is ch*ne11 + col, so each channel is a row of length // ne11 and this thread owns the [col0, col1) band of every row. + // + // The batched roles are inverted here relative to the matmul path: the + // filter is the single row and the positions are the nrc rows, since it is + // the positions that are strided by ne00 and the filter that is reused. + // That also makes the results contiguous, so the scale/bias pass vectorizes. + int32_t acc[GGML_MUL_MAT_ADD_OC_CHUNK]; + float bias_v[GGML_MUL_MAT_ADD_OC_CHUNK]; + for (int64_t ch = 0; ch < ne02; ch++) { const int8_t * w_ch = w + ch*ne00; - for (int64_t col = col0; col < col1; col++) { - // src1 is [K, N, C]: the channel stride is ne11*ne00. - const int8_t * x_col = x + ch*ne11*ne00 + col*ne00; + // Broadcast once per channel, not per position. + for (int64_t j = 0; j < GGML_MUL_MAT_ADD_OC_CHUNK; j++) { + bias_v[j] = bias[ch]; + } - int32_t acc = 0; - for (int64_t k = 0; k < ne00; k++) { - acc += (int32_t) w_ch[k] * (int32_t) x_col[k]; - } + for (int64_t col = col0; col < col1; col += GGML_MUL_MAT_ADD_OC_CHUNK) { + const int64_t ncol = MIN((int64_t) GGML_MUL_MAT_ADD_OC_CHUNK, col1 - col); - const float v = (float) acc * combined_d + bias[ch]; + // src1 is [K, N, C]: the channel stride is ne11*ne00. + const int8_t * x_col = x + ch*ne11*ne00 + col*ne00; - stage[ch*ne11 + col] = v; + ggml_vec_dot_i8_i8(ne00, acc, 1, x_col, ne00, w_ch, ncol); - const float av = fabsf(v); - if (av > local_absmax) local_absmax = av; + local_absmax = ggml_i8_s_scale_bias_absmax( + stage + ch*ne11 + col, acc, bias_v, ncol, combined_d, local_absmax); } } @@ -11946,23 +12090,22 @@ static void ggml_compute_forward_mul_mat_add_impl( // block [col0*ne01, col1*ne01). const size_t nb11 = src1->nb[1]; + // Output channels are handled a chunk at a time so the int32 + // accumulators fit a fixed stack buffer and params->wdata does not have + // to grow. Within a chunk the weight rows are contiguous and x_col stays + // hot, which is the whole reason for batching them into one call. + int32_t acc[GGML_MUL_MAT_ADD_OC_CHUNK]; + for (int64_t col = col0; col < col1; col++) { const int8_t * x_col = (const int8_t *) ((const char *) src1->data + col*nb11); - for (int64_t oc = 0; oc < ne01; oc++) { - const int8_t * w_row = w + oc*ne00; - - int32_t acc = 0; - for (int64_t k = 0; k < ne00; k++) { - acc += (int32_t) w_row[k] * (int32_t) x_col[k]; - } - - const float v = (float) acc * combined_d + bias[oc]; + for (int64_t oc0 = 0; oc0 < ne01; oc0 += GGML_MUL_MAT_ADD_OC_CHUNK) { + const int64_t noc = MIN((int64_t) GGML_MUL_MAT_ADD_OC_CHUNK, ne01 - oc0); - stage[col*ne01 + oc] = v; + ggml_vec_dot_i8_i8(ne00, acc, 1, w + oc0*ne00, ne00, x_col, noc); - const float av = fabsf(v); - if (av > local_absmax) local_absmax = av; + local_absmax = ggml_i8_s_scale_bias_absmax( + stage + col*ne01 + oc0, acc, bias + oc0, noc, combined_d, local_absmax); } } diff --git a/external/ggml/src/ggml-cpu/vec.cpp b/external/ggml/src/ggml-cpu/vec.cpp index 7e8135617..624090e1e 100644 --- a/external/ggml/src/ggml-cpu/vec.cpp +++ b/external/ggml/src/ggml-cpu/vec.cpp @@ -393,6 +393,115 @@ void ggml_vec_dot_f16(int n, float * GGML_RESTRICT s, size_t bs, ggml_fp16_t * G *s = sumf; } +#if defined(__AVX2__) +static inline int ggml_i8_hsum_i32_4(const __m128i a) { + const __m128i hi64 = _mm_unpackhi_epi64(a, a); + const __m128i sum64 = _mm_add_epi32(hi64, a); + const __m128i hi32 = _mm_shuffle_epi32(sum64, _MM_SHUFFLE(2, 3, 0, 1)); + return _mm_cvtsi128_si32(_mm_add_epi32(sum64, hi32)); +} +#endif + +void ggml_vec_dot_i8_i8(int n, int32_t * GGML_RESTRICT s, size_t bs, const int8_t * GGML_RESTRICT x, size_t bx, const int8_t * GGML_RESTRICT y, int nrc) { + for (int row = 0; row < nrc; ++row) { + const int8_t * xr = x + (size_t)row*bx; + + int i = 0; + int32_t sumi = 0; + +#if defined(__AVX2__) + // The sign trick: maddubs wants the left operand unsigned, so move the + // sign of x onto y and take |x|. Each maddubs lane holds a sum of two + // products, at most 2*128*127 = 32512, which still fits int16 - but a + // second int16 accumulation could not, so widen to int32 every block + // instead of batching in int16. + // + // Everything reduces into one 128-bit accumulator so that short rows + // (the depthwise filters are 4 or 8 long) neither run the 256-bit loop + // nor pay for reducing a register that stayed zero. + __m128i acc = _mm_setzero_si128(); + + if (i + 32 <= n) { + const __m256i one16 = _mm256_set1_epi16(1); + + __m256i acc256 = _mm256_setzero_si256(); + for (; i + 32 <= n; i += 32) { + const __m256i xq = _mm256_loadu_si256((const __m256i *)(xr + i)); + const __m256i yq = _mm256_loadu_si256((const __m256i *)(y + i)); + + const __m256i ax = _mm256_sign_epi8(xq, xq); + const __m256i sy = _mm256_sign_epi8(yq, xq); + const __m256i dot = _mm256_maddubs_epi16(ax, sy); + + acc256 = _mm256_add_epi32(acc256, _mm256_madd_epi16(dot, one16)); + } + acc = _mm_add_epi32(_mm256_castsi256_si128(acc256), + _mm256_extracti128_si256(acc256, 1)); + } + + if (i + 8 <= n) { + const __m128i one16 = _mm_set1_epi16(1); + + for (; i + 16 <= n; i += 16) { + const __m128i xq = _mm_loadu_si128((const __m128i *)(xr + i)); + const __m128i yq = _mm_loadu_si128((const __m128i *)(y + i)); + + const __m128i ax = _mm_sign_epi8(xq, xq); + const __m128i sy = _mm_sign_epi8(yq, xq); + const __m128i dot = _mm_maddubs_epi16(ax, sy); + + acc = _mm_add_epi32(acc, _mm_madd_epi16(dot, one16)); + } + for (; i + 8 <= n; i += 8) { + const __m128i xq = _mm_loadl_epi64((const __m128i *)(xr + i)); + const __m128i yq = _mm_loadl_epi64((const __m128i *)(y + i)); + + const __m128i ax = _mm_sign_epi8(xq, xq); + const __m128i sy = _mm_sign_epi8(yq, xq); + const __m128i dot = _mm_maddubs_epi16(ax, sy); + + acc = _mm_add_epi32(acc, _mm_madd_epi16(dot, one16)); + } + } + + // i > 0 exactly when some vector block ran. Reducing a register that + // stayed zero costs four instructions per row, which is not free when the + // row is a 4-tap depthwise filter and there is one row per output. + if (i > 0) { + sumi = ggml_i8_hsum_i32_4(acc); + } +#elif defined(__ARM_NEON) && defined(__aarch64__) + int32x4_t acc = vdupq_n_s32(0); + for (; i + 16 <= n; i += 16) { + const int8x16_t xv = vld1q_s8(xr + i); + const int8x16_t yv = vld1q_s8(y + i); + #if defined(__ARM_FEATURE_DOTPROD) + acc = vdotq_s32(acc, xv, yv); + #else + // vmull_s8 tops out at 128*128 = 16384, so the int16 products are + // safe; vpadalq_s16 folds them pairwise straight into int32. + acc = vpadalq_s16(acc, vmull_s8(vget_low_s8 (xv), vget_low_s8 (yv))); + acc = vpadalq_s16(acc, vmull_s8(vget_high_s8(xv), vget_high_s8(yv))); + #endif + } + for (; i + 8 <= n; i += 8) { + const int8x8_t xv = vld1_s8(xr + i); + const int8x8_t yv = vld1_s8(y + i); + acc = vpadalq_s16(acc, vmull_s8(xv, yv)); + } + if (i > 0) { + sumi = vaddvq_s32(acc); + } +#endif + + for (; i < n; ++i) { + sumi += (int32_t)xr[i] * (int32_t)y[i]; + } + + s[(size_t)row*bs] = sumi; + } +} + void ggml_vec_silu_f32(const int n, float * y, const float * x) { int i = 0; #if defined(__AVX512F__) && defined(__AVX512DQ__) diff --git a/external/ggml/src/ggml-cpu/vec.h b/external/ggml/src/ggml-cpu/vec.h index 741902479..bf425d76e 100644 --- a/external/ggml/src/ggml-cpu/vec.h +++ b/external/ggml/src/ggml-cpu/vec.h @@ -43,6 +43,13 @@ void ggml_vec_dot_f32(int n, float * GGML_RESTRICT s, size_t bs, const float * G void ggml_vec_dot_bf16(int n, float * GGML_RESTRICT s, size_t bs, ggml_bf16_t * GGML_RESTRICT x, size_t bx, ggml_bf16_t * GGML_RESTRICT y, size_t by, int nrc); void ggml_vec_dot_f16(int n, float * GGML_RESTRICT s, size_t bs, ggml_fp16_t * GGML_RESTRICT x, size_t bx, ggml_fp16_t * GGML_RESTRICT y, size_t by, int nrc); +// int8 x int8 dot products accumulating in int32, for the GGML_TYPE_I8_S ops. +// The scale is per-tensor and cancels out of the contraction, so these stay +// integral and the caller applies it once. nrc rows of x, each bx bytes apart, +// are contracted against the single row y; result row stride is bs int32s. +// n is arbitrary: whatever the vector width does not cover is done scalar. +void ggml_vec_dot_i8_i8(int n, int32_t * GGML_RESTRICT s, size_t bs, const int8_t * GGML_RESTRICT x, size_t bx, const int8_t * GGML_RESTRICT y, int nrc); + void ggml_vec_silu_f32(const int n, float * y, const float * x); ggml_float ggml_vec_cvar_f32(const int n, float * y, const float * x, const float mean); //it will also center y ( y = y - mean ) ggml_float ggml_vec_soft_max_f32(const int n, float * y, const float * x, float max); diff --git a/tests/unittests/test_i8_s_fused_ops.cpp b/tests/unittests/test_i8_s_fused_ops.cpp index ad2148c43..25c73ea90 100644 --- a/tests/unittests/test_i8_s_fused_ops.cpp +++ b/tests/unittests/test_i8_s_fused_ops.cpp @@ -67,8 +67,13 @@ std::vector patterned(size_t n, float phase, float scale) { return v; } -// The reference requantizer: absmax over the whole result, then round, then -// clamp. Mirrors ggml_i8_s_requantize. +// The reference requantizer, mirroring ggml_i8_s_requantize: absmax over the +// whole result -- taken before the relu clamp, so a value about to be zeroed +// still widens the scale -- then clamp, then round ties to even. +// +// std::rint, not std::round: the op rounds ties to even so that its AVX2 and +// NEON bodies agree with their scalar tails. std::round would round ties away +// from zero and make the byte-exact comparisons below fail on any tie. std::vector requantize_ref(const std::vector & values, bool relu, float * scale_out) { float amax = 0.0f; for (float v : values) { @@ -76,13 +81,13 @@ std::vector requantize_ref(const std::vector & values, bool relu, } const float id = amax != 0.0f ? 127.0f / amax : 0.0f; + const float lo = relu ? 0.0f : -127.0f; std::vector q(values.size()); for (size_t i = 0; i < values.size(); ++i) { float v = values[i] * id; - v = std::max(-127.0f, std::min(127.0f, v)); - const int8_t r = static_cast(std::round(v)); - q[i] = relu && r < 0 ? 0 : r; + v = std::max(lo, std::min(127.0f, v)); + q[i] = static_cast(std::rint(v)); } *scale_out = amax != 0.0f ? amax / 127.0f : 0.0f; @@ -286,9 +291,10 @@ void test_rms_norm_scaled(int n_threads) { std::cout << "rms_norm_scaled: C=" << C << " L=" << L << " nth=" << n_threads << " OK\n"; } -void test_mul_mat_add(bool relu, int n_threads) { - const int64_t IC = 64; - const int64_t OC = 48; +// IC and OC are varied by the caller rather than fixed: IC decides how much of +// the contraction the SIMD path covers and how much falls to the scalar tail, +// and OC decides whether the output-channel chunking loop wraps around. +void test_mul_mat_add(bool relu, int n_threads, int64_t IC, int64_t OC) { const int64_t N = 23; ggml_init_params ip = { kCtxBytes, nullptr, false }; @@ -331,6 +337,7 @@ void test_mul_mat_add(bool relu, int n_threads) { const auto want_q = requantize_ref(want, relu, &want_scale); const std::string label = std::string(relu ? "mul_mat_add_relu" : "mul_mat_add") + + " IC=" + std::to_string(IC) + " OC=" + std::to_string(OC) + " nth=" + std::to_string(n_threads); require_close(*inband_scale(result), want_scale, want_scale * 1e-4f, label + " scale"); @@ -352,13 +359,14 @@ void test_mul_mat_add(bool relu, int n_threads) { } ggml_free(ctx); - std::cout << label << ": IC=" << IC << " OC=" << OC << " N=" << N << " OK\n"; + std::cout << label << " N=" << N << " OK\n"; } -void test_mul_mat_add_depthwise(int n_threads) { - const int64_t K = 7; // kernel width +// K is the contraction length, so it decides how much of the dot product the SIMD +// path covers. N is the number of positions, which is what the depthwise path +// batches -- so it decides whether the chunking loop wraps around. +void test_mul_mat_add_depthwise(int n_threads, int64_t K, int64_t N) { const int64_t C = 40; // channels - const int64_t N = 19; // positions ggml_init_params ip = { kCtxBytes, nullptr, false }; ggml_context * ctx = ggml_init(ip); @@ -395,7 +403,9 @@ void test_mul_mat_add_depthwise(int n_threads) { float want_scale = 0.0f; const auto want_q = requantize_ref(want, false, &want_scale); - const std::string label = "mul_mat_add depthwise nth=" + std::to_string(n_threads); + const std::string label = "mul_mat_add depthwise K=" + std::to_string(K) + + " N=" + std::to_string(N) + + " nth=" + std::to_string(n_threads); require_close(*inband_scale(result), want_scale, want_scale * 1e-4f, label + " scale"); const auto * q = static_cast(result->data); @@ -405,7 +415,7 @@ void test_mul_mat_add_depthwise(int n_threads) { } ggml_free(ctx); - std::cout << label << ": K=" << K << " C=" << C << " N=" << N << " OK\n"; + std::cout << label << ": C=" << C << " OK\n"; } void test_im2col_asym(int n_threads) { @@ -472,9 +482,21 @@ int main() { for (int nth : {1, 4}) { test_add_scaled(nth); test_rms_norm_scaled(nth); - test_mul_mat_add(false, nth); - test_mul_mat_add(true, nth); - test_mul_mat_add_depthwise(nth); + // IC=64: SIMD only, no tail. IC=67: three elements land in the + // scalar tail, which a truncating block count would silently drop. + // IC=13: shorter than one 32-byte step, so entirely scalar. + // OC=48 stays inside one output-channel chunk, OC=100 spans two. + test_mul_mat_add(false, nth, 64, 48); + test_mul_mat_add(true, nth, 64, 48); + test_mul_mat_add(false, nth, 67, 48); + test_mul_mat_add(false, nth, 13, 48); + test_mul_mat_add(false, nth, 67, 100); + // K=7 is shorter than any vector step, so the dot is all scalar; + // K=36 runs the 32-byte body plus a 4-element tail. N=19 stays inside + // one position chunk, N=150 spans three. + test_mul_mat_add_depthwise(nth, 7, 19); + test_mul_mat_add_depthwise(nth, 36, 19); + test_mul_mat_add_depthwise(nth, 7, 150); test_im2col_asym(nth); } } catch (const std::exception & e) { From e2dd409e78d1eea8cef6bb057452888b295576d0 Mon Sep 17 00:00:00 2001 From: XSquirrelC Date: Fri, 4 Sep 2026 03:30:05 +0000 Subject: [PATCH 3/3] ggml-cpu: carry the I8_S in-band scale through dup/cont --- external/ggml/src/ggml-cpu/ops.cpp | 10 ++++++ tests/unittests/test_i8_s_fused_ops.cpp | 41 +++++++++++++++++++++++++ 2 files changed, 51 insertions(+) diff --git a/external/ggml/src/ggml-cpu/ops.cpp b/external/ggml/src/ggml-cpu/ops.cpp index 735030a78..b07bb7823 100644 --- a/external/ggml/src/ggml-cpu/ops.cpp +++ b/external/ggml/src/ggml-cpu/ops.cpp @@ -531,6 +531,16 @@ void ggml_compute_forward_dup( if (src0->type == dst->type) { ggml_compute_forward_dup_bytes(params, dst); + // I8_S keeps one scale for the whole tensor, stored past the last + // element: a byte copy moves the payload but leaves that float + // uninitialized, so a cont/cpy of an I8_S view has to carry it over. + if (dst->type == GGML_TYPE_I8_S && params->ith == 0) { + // Reading through to the parent keeps this working for the permuted + // views that ggml_cont() is usually handed, which have no scale of + // their own. + const ggml_tensor * scale_src = src0->view_src ? src0->view_src : src0; + *ggml_inband_scale(dst) = *ggml_inband_scale_const(scale_src); + } return; } diff --git a/tests/unittests/test_i8_s_fused_ops.cpp b/tests/unittests/test_i8_s_fused_ops.cpp index 25c73ea90..8f58519e6 100644 --- a/tests/unittests/test_i8_s_fused_ops.cpp +++ b/tests/unittests/test_i8_s_fused_ops.cpp @@ -469,6 +469,46 @@ void test_im2col_asym(int n_threads) { << " OW=" << OW << " nth=" << n_threads << " OK\n"; } +// The VibeASR encoder flips its activations between channel-major and +// length-major with ggml_cont(ggml_permute(...)), which lands in +// ggml_compute_forward_dup. A byte copy moves the payload but not the in-band +// scale, and the permuted view it is handed has no scale of its own, so the +// scale has to be read through view_src. Getting this wrong leaves the copy +// holding whatever was in the buffer -- values stay right while everything +// downstream is off by an arbitrary factor. +void test_i8_s_cont_permute(int n_threads) { + const int64_t C = 5; + const int64_t L = 33; + + ggml_init_params ip = { kCtxBytes, nullptr, false }; + ggml_context * ctx = ggml_init(ip); + + ggml_tensor * src = ggml_new_tensor_2d(ctx, GGML_TYPE_I8_S, C, L); + const float src_scale = fill_i8_s(src, patterned(C * L, 0.6f, 1.9f)); + + ggml_tensor * result = ggml_cont(ctx, ggml_permute(ctx, src, 1, 0, 2, 3)); + require_eq(result->ne[0], L, "cont(permute) ne0"); + require_eq(result->ne[1], C, "cont(permute) ne1"); + + compute(ctx, result, n_threads); + + const auto * in = static_cast(src->data); + const auto * out = static_cast(result->data); + for (int64_t ic = 0; ic < C; ++ic) { + for (int64_t il = 0; il < L; ++il) { + require_eq(static_cast(out[ic * L + il]), static_cast(in[il * C + ic]), + "cont(permute) value at c=" + std::to_string(ic) + " l=" + std::to_string(il)); + } + } + + // Rearrangement only: the same scale has to come out the other side. + require_close(*inband_scale(result), src_scale, 0.0f, "cont(permute) scale carried over"); + + ggml_free(ctx); + std::cout << "cont(permute) i8_s: C=" << C << " L=" << L + << " nth=" << n_threads << " OK\n"; +} + } // namespace int main() { @@ -498,6 +538,7 @@ int main() { test_mul_mat_add_depthwise(nth, 36, 19); test_mul_mat_add_depthwise(nth, 7, 150); test_im2col_asym(nth); + test_i8_s_cont_permute(nth); } } catch (const std::exception & e) { std::cerr << "FAILED: " << e.what() << "\n";