diff --git a/conversion/__init__.py b/conversion/__init__.py
index 72960be1a..5bf472889 100644
--- a/conversion/__init__.py
+++ b/conversion/__init__.py
@@ -255,6 +255,7 @@ TEXT_MODEL_MAP: dict[str, str] = {
"Qwen3_5ForConditionalGeneration": "qwen",
"Qwen3_5MoeForCausalLM": "qwen",
"Qwen3_5MoeForConditionalGeneration": "qwen",
+ "Qwen3_5TextModel": "qwen",
"Qwen4ExpForCausalLM": "qwen4exp",
"Qwen4ExpForConditionalGeneration": "qwen4exp",
"RND1": "qwen",
diff --git a/conversion/base.py b/conversion/base.py
index 306a9f0f2..7a93855c8 100644
--- a/conversion/base.py
+++ b/conversion/base.py
@@ -2336,12 +2336,26 @@ class TextModel(ModelBase):
else:
raise NotImplementedError("Only MEAN, CLS, and LAST pooling types supported")
self.gguf_writer.add_pooling_type(pooling_type)
+ else:
+ embedding_config_path = self.dir_model / "embedding_config.json"
+ if embedding_config_path.is_file():
+ with open(embedding_config_path, encoding="utf-8") as f:
+ embedding_config = json.load(f)
+ pooling = embedding_config.get("pooling")
+ if pooling == "last_token":
+ self.gguf_writer.add_pooling_type(gguf.PoolingType.LAST)
+ elif pooling is not None:
+ raise NotImplementedError(f"unsupported embedding_config.json pooling {pooling!r}")
# pooling before a classification head (e.g. ModernBertForSequenceClassification)
if (classifier_pooling := self.hparams.get("classifier_pooling")) is not None:
if classifier_pooling not in ("cls", "mean"):
raise NotImplementedError(f"Unsupported classifier_pooling: {classifier_pooling}")
self.gguf_writer.add_classifier_pooling_type(mode_mapping[classifier_pooling])
+ if (classifier_activation := self.hparams.get("classifier_activation")) is not None:
+ if classifier_activation not in ("gelu", "silu", "tanh"):
+ raise NotImplementedError(f"Unsupported classifier_activation: {classifier_activation}")
+ self.gguf_writer.add_classifier_activation(classifier_activation)
def _set_vocab_glmedge(self):
from transformers import AutoTokenizer
diff --git a/conversion/qwen.py b/conversion/qwen.py
index 0b3918733..3af7a035c 100644
--- a/conversion/qwen.py
+++ b/conversion/qwen.py
@@ -650,11 +650,24 @@ class _Qwen35MRopeMixin:
self.gguf_writer.add_rope_dimension_sections(self._QWEN35_DEFAULT_MROPE_SECTION)
-@ModelBase.register("Qwen3_5ForConditionalGeneration", "Qwen3_5ForCausalLM")
+@ModelBase.register("Qwen3_5ForConditionalGeneration", "Qwen3_5ForCausalLM", "Qwen3_5TextModel")
@ModelBase.example("Qwen/Qwen3.5-9B")
class Qwen3_5TextModel(_Qwen35MRopeMixin, _LinearAttentionVReorderBase):
model_arch = gguf.MODEL_ARCH.QWEN35
+ def __init__(self, dir_model, *args, **kwargs):
+ # Inner TextModel does not own mtp.*. Set no_mtp before mixin bumps block_count.
+ hparams = kwargs.pop("hparams", None)
+ if hparams is None:
+ hparams = ModelBase.load_hparams(dir_model, self.is_mistral_format)
+ if get_model_architecture(hparams, ModelType.TEXT) == "Qwen3_5TextModel":
+ self.no_mtp = True
+ super().__init__(dir_model, *args, hparams=hparams, **kwargs)
+
+ def set_gguf_parameters(self):
+ super().set_gguf_parameters()
+ self._try_set_pooling_type()
+
def _is_openjev_checkpoint(dir_model: Path) -> bool:
return (dir_model / "helper" / "shim.py").is_file() and (dir_model / "config.json").is_file()
diff --git a/docs/backend/SYCL.md b/docs/backend/SYCL.md
index 70f5f2eff..19dc69851 100644
--- a/docs/backend/SYCL.md
+++ b/docs/backend/SYCL.md
@@ -803,6 +803,7 @@ User can use the device management in [docs/multi-gpu.md](https://github.com/ggm
| GGML_SYCL_ENABLE_GRAPH | 0 (default) or 1 | Enable running computations through SYCL Graphs feature. Disabled by default because SYCL Graph is still on development, no better performance. |
| GGML_SYCL_ENABLE_HOST_PINNED_MEM | 0 or 1 (default) | Enable host pinned memory to speed up copy data from host to device. When disable it, host memory will common malloc() on CPU. Disable it when use `--load-model mlock`.|
| GGML_SYCL_HOST_PINNED_MEM_2G | 0 (default) or 1 | Limit the max memory allocation to be no more than 2GB when enable host pinned memory. USM allocations above 2 GiB take the relaxed/large-allocation path, which serializes H2D copies with compute and prevents copy/compute overlap. It will impact the startup time. Need more test. Depend on `GGML_SYCL_ENABLE_HOST_PINNED_MEM=1`.|
+| GGML_SYCL_UPLOAD_STAGING_SLOTS | 4 (default) or non-negative integer | Number of 8 MiB pinned host slots used to stage tensor uploads (model loading), so the host copy of one slot overlaps the transfer of the previous one. Set to 0 to use the old path: a malloc'd bounce buffer and a blocking copy per tensor. |
| GGML_SYCL_GET_MEM_API | 0 (default) or 1 | Set to get memory info (free, total) by Level Zero or SYCL API:
0 - Level Zero API: support more GPUs, only run on Level Zero running time. When there is an error, fallback to call SYCL API. Depend on GGML_SYCL_SUPPORT_LEVEL_ZERO_API.
1 - SYCL API: legacy, support more running time, it can't get the free size of some GPUs (like Arc770). In such case, return the free size as value of total size.|
| GGML_SYCL_USE_LEVEL_ZERO_API | 1 (default) or 0 | Use Level Zero API for device memory allocation instead of SYCL. Reduces system RAM usage on Intel dGPUs by avoiding DMA-buf/TTM host memory staging. Requires GGML_SYCL_SUPPORT_LEVEL_ZERO_API=ON at build time. SYCL backend always runs on Level Zero running time even if it's set as OFF (The SYCL api will be usage for memory allocation).|
| GGML_SYCL_ENABLE_DNN | 0 or 1 (default)| Enable running computations through oneDNN and always use oneMKL. |
@@ -816,6 +817,10 @@ User can use the device management in [docs/multi-gpu.md](https://github.com/ggm
| GGML_SYCL_MKL_FA_DIAG | 0 (default) or 1 | Enable output fingerprinting for MKL flash attention. Dumps the first 64 float output values for the first 6 FA calls with n_kv ≥ 1024, labeled with kernel type (MKL/TILE/VEC) for cross-kernel comparison. |
| GGML_SYCL_ENABLE_FUSION | 0 or 1 (default) | Enable fused-kernel dispatch in graph compute. Unsupported types and layouts fall back to the standalone op kernels. See `ggml_sycl_can_fuse()`. |
| GGML_SYCL_ENABLE_ESIMD | 0 or 1 (default)| Enable ESIMD kernels when available. |
+| GGML_SYCL_XMX_GATHER_TYPES | decimal bitmask, all bits set (default) | Weight formats that may use the XMX dequant-GEMM paths, which dequantize weights straight into the XMX tiles. This speeds up prompt processing of MoE models on GPUs with XMX units (Arc A- and B-series, Arc Pro, Data Center GPU Max), for example pp512 of Qwen3-30B-A3B UD-IQ3_XXS by about 50% on an Arc Pro B60. Bits:
* 1: IQ4_NL, 2: IQ3_S, 4: IQ4_XS, 8: IQ3_XXS, 16: IQ2_XXS, 32: IQ2_XS, 64: IQ2_S, 128: IQ1_S, 256: IQ1_M
* 512: Q8_0, 1024: Q4_K, 2048: Q5_K, 4096: Q6_K (MoE `MUL_MAT_ID` only)
Add values to combine them, for example `3` for IQ4_NL and IQ3_S; `0` disables the paths. A set bit does not force the path: batches of more than 64 tokens per expert or row lengths that are not a multiple of 256 (32 for IQ4_NL and Q8_0) use the library GEMM. |
+| GGML_SYCL_XMX_GATHER_SHAPES | decimal bitmask, 255 (default) | XMX `joint_matrix` combinations the paths of `GGML_SYCL_XMX_GATHER_TYPES` may use; the operand type comes from `GGML_SYCL_DYNAMIC_PRECISION` and the best supported combination is picked automatically (logged as `fg_pick_combo`). Bits:
* Xe2, Xe3, Xe-HPC: 1: f16 8x16x16, 2: f16 16x16x16, 4: f16 32x64x16, 8: f16 32x64x32, 32: tf32 8x16x8, 64: bf16 8x16x16
* Xe-HPG (Arc A770, ARL-H): 16: f16 8x8x16, 128: bf16 8x8x16
Clear a bit to exclude a combination, or set a single bit to force one for testing. |
+| GGML_SYCL_DYNAMIC_PRECISION | `F16` (default with `GGML_SYCL_F16=ON`), `BF16`, `TF32` or `F32` (default otherwise) | Operand type of the XMX dequant-GEMM paths (`GGML_SYCL_XMX_GATHER_TYPES`); accumulation is always f32. `F16` is the fastest, but activations above 65504 overflow. `BF16` keeps the f32 range at a 7-bit mantissa, `TF32` keeps the range and the f16 mantissa but is about 30% slower and needs Xe2, Xe3 or Xe-HPC, and `F32` turns the XMX paths off. Ops that request a higher src1 precision ([TAG_GGML_PREC]) get it regardless of this setting. |
+| GGML_SYCL_DYNAMIC_REQUIRED_PRECISION | `F32` (default), `TF32`, `BF16` or `F16` | Lowest type the XMX paths may use for an op that requests an F32 src1, such as Mistral 4 `ffn_down_exps`. The default runs such ops on the library f32 GEMM; `TF32` or `BF16` trade mantissa for speed while keeping the f32 range. `F16` ignores the request and can overflow; it is meant for testing only. |
| GGML_SYCL_MMVQ_WIDE | 0 or 1 (default) | Use the wide-load variant of the reordered Q8_0 mat-vec kernel, which reads four contiguous dwords per operand instead of one value at a time. Set to 0 to fall back to the per-value loads. Only affects Q8_0 weights in the reordered layout. |
| GGML_SYCL_SPARSE_FA | 0 (default) or 1 | Enable Sparse Flash-attention.|
| GGML_SYCL_SPARSE_FA_DEBUG | 0 (default) or 1 | Enable to debug for Sparse Flash-attention.|
diff --git a/ggml/src/ggml-cuda/gated_delta_net.cu b/ggml/src/ggml-cuda/gated_delta_net.cu
index 6498c6a14..fac12f234 100644
--- a/ggml/src/ggml-cuda/gated_delta_net.cu
+++ b/ggml/src/ggml-cuda/gated_delta_net.cu
@@ -1,7 +1,9 @@
#include "gated_delta_net.cuh"
#include "ggml-cuda/common.cuh"
-template
+constexpr int gdn_cols_per_warp = 4;
+
+template
__global__ void __launch_bounds__((ggml_cuda_get_physical_warp_size() < S_v ? ggml_cuda_get_physical_warp_size() : S_v) * 4, 2)
gated_delta_net_cuda(const float * q,
const float * k,
@@ -33,9 +35,19 @@ gated_delta_net_cuda(const float * q,
int txn_slots) {
const uint32_t h_idx = blockIdx.x;
const uint32_t sequence = blockIdx.y;
- // each warp owns one column, using warp-level primitives to reduce across rows
- const int lane = threadIdx.x;
- const int col = blockIdx.z * blockDim.y + threadIdx.y;
+
+ constexpr int warp_size = ggml_cuda_get_physical_warp_size() < S_v ? ggml_cuda_get_physical_warp_size() : S_v;
+ static_assert(S_v % warp_size == 0, "S_v must be a multiple of warp_size");
+ // the warp is split into cols_per_warp segments of lanes_per_col lanes; each segment owns
+ // one state column and reduces within itself
+ constexpr int lanes_per_col = warp_size / cols_per_warp;
+ constexpr int rows_per_lane = S_v / lanes_per_col;
+ static_assert(S_v % lanes_per_col == 0, "S_v must be a multiple of lanes_per_col");
+
+ const int lane = threadIdx.x;
+ const int col_in_warp = lane / lanes_per_col; // column slot within the warp
+ const int lane_in_col = lane - col_in_warp * lanes_per_col; // lane within the column's reduction segment
+ const int col = (blockIdx.z * blockDim.y + threadIdx.y) * cols_per_warp + col_in_warp;
const uint32_t iq1 = fastmodulo(h_idx, neqk1_magic);
const uint32_t iq3 = fastdiv(sequence, rq3_magic);
@@ -50,16 +62,13 @@ gated_delta_net_cuda(const float * q,
curr_state += state_in_offset + col * S_v;
attn_data += (sequence * n_tokens * H + h_idx) * S_v;
- constexpr int warp_size = ggml_cuda_get_physical_warp_size() < S_v ? ggml_cuda_get_physical_warp_size() : S_v;
- static_assert(S_v % warp_size == 0, "S_v must be a multiple of warp_size");
- constexpr int rows_per_lane = (S_v + warp_size - 1) / warp_size;
float s_shard[rows_per_lane];
// state is stored transposed: M[col][i] = S[i][col], row col is contiguous
ggml_cuda_pdl_sync();
#pragma unroll
for (int r = 0; r < rows_per_lane; r++) {
- const int i = r * warp_size + lane;
+ const int i = r * lanes_per_col + lane_in_col;
s_shard[r] = curr_state[i];
}
@@ -86,7 +95,7 @@ gated_delta_net_cuda(const float * q,
float q_reg[rows_per_lane];
#pragma unroll
for (int r = 0; r < rows_per_lane; r++) {
- const int i = r * warp_size + lane;
+ const int i = r * lanes_per_col + lane_in_col;
k_reg[r] = k_t[i];
q_reg[r] = q_t[i];
}
@@ -100,7 +109,7 @@ gated_delta_net_cuda(const float * q,
for (int r = 0; r < rows_per_lane; r++) {
kv_shard += s_shard[r] * k_reg[r];
}
- float kv_col = warp_reduce_sum(kv_shard);
+ float kv_col = warp_reduce_sum(kv_shard);
// delta[col] = (v[col] - g * kv[col]) * beta
float delta_col = (v_t[col] - g_val * kv_col) * beta_val;
@@ -130,9 +139,9 @@ gated_delta_net_cuda(const float * q,
attn_partial += s_shard[r] * q_reg[r];
}
- float attn_col = warp_reduce_sum(attn_partial);
+ float attn_col = warp_reduce_sum(attn_partial);
- if (lane == 0) {
+ if (lane_in_col == 0) {
attn_data[col] = attn_col * scale;
}
} else {
@@ -140,11 +149,11 @@ gated_delta_net_cuda(const float * q,
float kv_shard = 0.0f;
#pragma unroll
for (int r = 0; r < rows_per_lane; r++) {
- const int i = r * warp_size + lane;
+ const int i = r * lanes_per_col + lane_in_col;
kv_shard += expf(g_t[i]) * s_shard[r] * k_reg[r];
}
- float kv_col = warp_reduce_sum(kv_shard);
+ float kv_col = warp_reduce_sum(kv_shard);
// delta[col] = (v[col] - kv[col]) * beta
float delta_col = (v_t[col] - kv_col) * beta_val;
@@ -154,14 +163,14 @@ gated_delta_net_cuda(const float * q,
float attn_partial = 0.0f;
#pragma unroll
for (int r = 0; r < rows_per_lane; r++) {
- const int i = r * warp_size + lane;
+ const int i = r * lanes_per_col + lane_in_col;
s_shard[r] = expf(g_t[i]) * s_shard[r] + k_reg[r] * delta_col;
attn_partial += s_shard[r] * q_reg[r];
}
- float attn_col = warp_reduce_sum(attn_partial);
+ float attn_col = warp_reduce_sum(attn_partial);
- if (lane == 0) {
+ if (lane_in_col == 0) {
attn_data[col] = attn_col * scale;
}
}
@@ -176,7 +185,7 @@ gated_delta_net_cuda(const float * q,
float * curr_state = state + target_slot * state_slot_stride;
#pragma unroll
for (int r = 0; r < rows_per_lane; r++) {
- const int i = r * warp_size + lane;
+ const int i = r * lanes_per_col + lane_in_col;
curr_state[col * S_v + i] = s_shard[r];
}
}
@@ -186,7 +195,7 @@ gated_delta_net_cuda(const float * q,
if constexpr (!keep_rs_t) {
#pragma unroll
for (int r = 0; r < rows_per_lane; r++) {
- const int i = r * warp_size + lane;
+ const int i = r * lanes_per_col + lane_in_col;
state[col * S_v + i] = s_shard[r];
}
}
@@ -206,8 +215,16 @@ static void launch_gated_delta_net(
float * txn_log = nullptr, float * txn_prior = nullptr, int txn_slots = 0) {
//TODO: Add chunked kernel for even faster pre-fill
const int warp_size = ggml_cuda_info().devices[ggml_cuda_get_device()].warp_size;
- const int num_warps = 4;
- dim3 grid_dims(H, n_seqs, (S_v + num_warps - 1) / num_warps);
+ // four columns per warp (see the kernel); shrink the CTA when the wider CTA would leave
+ // SMs without a CTA, so small head counts keep the device filled
+ const int nsm = ggml_cuda_info().devices[ggml_cuda_get_device()].nsm;
+ const int cols_per_warp = gdn_cols_per_warp;
+ int num_warps = 4;
+ while (num_warps > 1 && H*n_seqs*(S_v / (cols_per_warp * num_warps)) < nsm) {
+ num_warps /= 2;
+ }
+ // one CTA covers cols_per_warp*num_warps columns (see the kernel)
+ dim3 grid_dims(H, n_seqs, (S_v + cols_per_warp * num_warps - 1) / (cols_per_warp * num_warps));
dim3 block_dims(warp_size <= S_v ? warp_size : S_v, num_warps, 1);
const uint3 neqk1_magic = init_fastdiv_values(neqk1);
diff --git a/ggml/src/ggml-cuda/ggml-cuda.cu b/ggml/src/ggml-cuda/ggml-cuda.cu
index d85bfe3be..203a41fa1 100644
--- a/ggml/src/ggml-cuda/ggml-cuda.cu
+++ b/ggml/src/ggml-cuda/ggml-cuda.cu
@@ -5427,9 +5427,10 @@ static bool ggml_backend_cuda_device_supports_op(ggml_backend_dev_t dev, const g
case GGML_UNARY_OP_CEIL:
case GGML_UNARY_OP_ROUND:
case GGML_UNARY_OP_TRUNC:
- // TODO: should become:
- //return ggml_is_contiguous_rows(op->src[0]);
- return ggml_is_contiguous(op->src[0]);
+ if (op->src[0]->type == GGML_TYPE_BF16 && ggml_get_unary_op(op) == GGML_UNARY_OP_XIELU) {
+ return false;
+ }
+ return op->src[0]->type == GGML_TYPE_F16 || op->src[0]->type == GGML_TYPE_F32 || op->src[0]->type == GGML_TYPE_BF16;
default:
return false;
}
@@ -5771,7 +5772,7 @@ static bool ggml_backend_cuda_device_supports_op(ggml_backend_dev_t dev, const g
return max_bias == 0.0f;
}
case GGML_OP_ROLL:
- if(op->src[0]->type == GGML_TYPE_F32 && ggml_is_contiguous(op->src[0])) {
+ if(op->src[0]->type == GGML_TYPE_F32) {
return true;
}
return false;
diff --git a/ggml/src/ggml-cuda/norm.cu b/ggml/src/ggml-cuda/norm.cu
index 5543307be..d804b933b 100644
--- a/ggml/src/ggml-cuda/norm.cu
+++ b/ggml/src/ggml-cuda/norm.cu
@@ -3,38 +3,46 @@
template
static __global__ void norm_f32(
- const float * x, float * dst, const int ncols, const int64_t stride_row, const int64_t stride_channel,
- const int64_t stride_sample, const float eps) {
- const int nrows = gridDim.x;
- const int nchannels = gridDim.y;
+ const float * x, float * dst, const int ncols, const int nchannels, const int nsamples,
+ const int64_t stride_row, const int64_t stride_channel, const int64_t stride_sample, const float eps) {
+ const int nrows = gridDim.x;
+ const int row = blockIdx.x;
+ const int tid = threadIdx.x;
- const int row = blockIdx.x;
- const int channel = blockIdx.y;
- const int sample = blockIdx.z;
- const int tid = threadIdx.x;
-
- x += sample*stride_sample + channel*stride_channel + row*stride_row;
- dst += ((sample*nchannels + channel)*nrows + row)*ncols;
-
- float2 mean_var = make_float2(0.0f, 0.0f);
+ extern __shared__ float2 s_sum2[];
ggml_cuda_pdl_sync();
- for (int col = tid; col < ncols; col += block_size) {
- const float xi = x[col];
- mean_var.x += xi;
- mean_var.y += xi * xi;
- }
- // sum up partial sums
- extern __shared__ float2 s_sum2[];
- mean_var = block_reduce(mean_var, s_sum2);
+ // grid.y and grid.z are clamped to the CUDA limit, iterate over the excess channels/samples
+ for (int sample = blockIdx.z; sample < nsamples; sample += gridDim.z) {
+ for (int channel = blockIdx.y; channel < nchannels; channel += gridDim.y) {
+ const float * xc = x + sample*stride_sample + channel*stride_channel + row*stride_row;
+ float * dstc = dst + ((sample*nchannels + channel)*nrows + row)*ncols;
- const float mean = mean_var.x / ncols;
- const float var = mean_var.y / ncols - mean * mean;
- const float inv_std = rsqrtf(var + eps);
+ float2 mean_var = make_float2(0.0f, 0.0f);
- for (int col = tid; col < ncols; col += block_size) {
- dst[col] = (x[col] - mean) * inv_std;
+ for (int col = tid; col < ncols; col += block_size) {
+ const float xi = xc[col];
+ mean_var.x += xi;
+ mean_var.y += xi * xi;
+ }
+
+ // sum up partial sums
+ mean_var = block_reduce(mean_var, s_sum2);
+
+ const float mean = mean_var.x / ncols;
+ const float var = mean_var.y / ncols - mean * mean;
+ const float inv_std = rsqrtf(var + eps);
+
+ for (int col = tid; col < ncols; col += block_size) {
+ dstc[col] = (xc[col] - mean) * inv_std;
+ }
+
+ if constexpr (block_size > WARP_SIZE) {
+ // sync is needed as we reuse s_sum2 across block_reduce invocations, see #26385
+ __syncthreads();
+ }
+ }
}
}
@@ -77,6 +85,8 @@ template (tmp, s_sum);
+ // grid.y and grid.z are clamped to the CUDA limit, iterate over the excess channels/samples
+ for (int sample = blockIdx.z; sample < nsamples; sample += gridDim.z) {
+ for (int channel = blockIdx.y; channel < nchannels; channel += gridDim.y) {
+ const float * xc = x + sample*stride_sample + channel*stride_channel + row*stride_row;
+ float * dstc = dst + ((sample*nchannels + channel)*nrows + row)*ncols;
- const float mean = tmp / ncols;
- const float scale = rsqrtf(mean + eps);
+ [[maybe_unused]] const float * mulc = nullptr;
+ if constexpr (do_multiply) {
+ const uint32_t mul_row = fastmodulo(row, mul_nrows_packed);
+ const uint32_t mul_channel = fastmodulo(channel, mul_nchannels_packed);
+ const uint32_t mul_sample = fastmodulo(sample, mul_nsamples_packed);
+ mulc = mul + mul_sample * mul_stride_sample + mul_channel * mul_stride_channel + mul_row * mul_stride_row;
+ }
- for (int col = tid; col < ncols; col += block_size) {
- if constexpr (do_multiply && do_add) {
- const int mul_col = fastmodulo(col, mul_ncols_packed);
- const int add_col = fastmodulo(col, add_ncols_packed);
- dst[col] = scale * x[col] * mul[mul_col] + add[add_col];
- } else if constexpr (do_multiply) {
- const int mul_col = fastmodulo(col, mul_ncols_packed);
- dst[col] = scale * x[col] * mul[mul_col];
- } else if constexpr (do_scale) {
- dst[col] = scale_out * (scale * x[col]);
- } else {
- dst[col] = scale * x[col];
+ [[maybe_unused]] const float * addc = nullptr;
+ if constexpr (do_add) {
+ const int add_row = fastmodulo(row, add_nrows_packed);
+ const int add_channel = fastmodulo(channel, add_nchannels_packed);
+ const int add_sample = fastmodulo(sample, add_nsamples_packed);
+ addc = add + add_sample * add_stride_sample + add_channel * add_stride_channel + add_row * add_stride_row;
+ }
+
+ float tmp = 0.0f; // partial sum for thread in warp
+
+ for (int col = tid; col < ncols; col += block_size) {
+ const float xi = xc[col];
+ tmp += xi * xi;
+ }
+
+ // sum up partial sums
+ tmp = block_reduce(tmp, s_sum);
+
+ const float mean = tmp / ncols;
+ const float scale = rsqrtf(mean + eps);
+
+ for (int col = tid; col < ncols; col += block_size) {
+ if constexpr (do_multiply && do_add) {
+ const int mul_col = fastmodulo(col, mul_ncols_packed);
+ const int add_col = fastmodulo(col, add_ncols_packed);
+ dstc[col] = scale * xc[col] * mulc[mul_col] + addc[add_col];
+ } else if constexpr (do_multiply) {
+ const int mul_col = fastmodulo(col, mul_ncols_packed);
+ dstc[col] = scale * xc[col] * mulc[mul_col];
+ } else if constexpr (do_scale) {
+ dstc[col] = scale_out * (scale * xc[col]);
+ } else {
+ dstc[col] = scale * xc[col];
+ }
+ }
+
+ if constexpr (block_size > WARP_SIZE) {
+ // sync is needed as we reuse s_sum across block_reduce invocations, see #26385
+ __syncthreads();
+ }
}
}
}
@@ -247,50 +267,57 @@ static __global__ void rms_norm_back_f32(
template
static __global__ void l2_norm_f32(
- const float * x, float * dst, const int ncols, const int64_t stride_row, const int64_t stride_channel,
- const int64_t stride_sample, const float eps) {
- const int nrows = gridDim.x;
- const int nchannels = gridDim.y;
+ const float * x, float * dst, const int ncols, const int nchannels, const int nsamples,
+ const int64_t stride_row, const int64_t stride_channel, const int64_t stride_sample, const float eps) {
+ const int nrows = gridDim.x;
+ const int row = blockIdx.x;
+ const int tid = threadIdx.x;
- const int row = blockIdx.x;
- const int channel = blockIdx.y;
- const int sample = blockIdx.z;
- const int tid = threadIdx.x;
-
- x += sample*stride_sample + channel*stride_channel + row*stride_row;
- dst += ((sample*nchannels + channel)*nrows + row)*ncols;
-
- float tmp = 0.0f; // partial sum for thread in warp
+ extern __shared__ float s_sum[];
ggml_cuda_pdl_sync();
- for (int col = tid; col < ncols; col += block_size) {
- const float xi = x[col];
- tmp += xi * xi;
- }
- // sum up partial sums
- extern __shared__ float s_sum[];
- tmp = block_reduce(tmp, s_sum);
- ggml_cuda_pdl_lc();
+ // grid.y and grid.z are clamped to the CUDA limit, iterate over the excess channels/samples
+ for (int sample = blockIdx.z; sample < nsamples; sample += gridDim.z) {
+ for (int channel = blockIdx.y; channel < nchannels; channel += gridDim.y) {
+ const float * xc = x + sample*stride_sample + channel*stride_channel + row*stride_row;
+ float * dstc = dst + ((sample*nchannels + channel)*nrows + row)*ncols;
- // from https://pytorch.org/docs/stable/generated/torch.nn.functional.normalize.html
- const float scale = rsqrtf(fmaxf(tmp, eps * eps));
+ float tmp = 0.0f; // partial sum for thread in warp
- for (int col = tid; col < ncols; col += block_size) {
- dst[col] = scale * x[col];
+ for (int col = tid; col < ncols; col += block_size) {
+ const float xi = xc[col];
+ tmp += xi * xi;
+ }
+
+ // sum up partial sums
+ tmp = block_reduce(tmp, s_sum);
+
+ // from https://pytorch.org/docs/stable/generated/torch.nn.functional.normalize.html
+ const float scale = rsqrtf(fmaxf(tmp, eps * eps));
+
+ for (int col = tid; col < ncols; col += block_size) {
+ dstc[col] = scale * xc[col];
+ }
+
+ if constexpr (block_size > WARP_SIZE) {
+ // sync is needed as we reuse s_sum across block_reduce invocations, see #26385
+ __syncthreads();
+ }
+ }
}
}
static void norm_f32_cuda(
const float * x, float * dst, const int ncols, const int nrows, const int nchannels, const int nsamples,
const int64_t stride_row, const int64_t stride_channel, const int64_t stride_sample, const float eps, cudaStream_t stream) {
- const dim3 blocks_num(nrows, nchannels, nsamples);
+ const dim3 blocks_num(nrows, MIN(nchannels, UINT16_MAX), MIN(nsamples, UINT16_MAX));
if (ncols < 1024) {
const dim3 block_dims(WARP_SIZE, 1, 1);
- norm_f32<<>>(x, dst, ncols, stride_row, stride_channel, stride_sample, eps);
+ norm_f32<<>>(x, dst, ncols, nchannels, nsamples, stride_row, stride_channel, stride_sample, eps);
} else {
const dim3 block_dims(1024, 1, 1);
- norm_f32<1024><< WARP_SIZE ? 32 * sizeof(float2): 0, stream>>>(x, dst, ncols, stride_row, stride_channel, stride_sample, eps);
+ norm_f32<1024><< WARP_SIZE ? 32 * sizeof(float2): 0, stream>>>(x, dst, ncols, nchannels, nsamples, stride_row, stride_channel, stride_sample, eps);
}
}
@@ -310,19 +337,19 @@ static void rms_norm_f32_cuda(
const float * x, float * dst, const int ncols, const int nrows, const int nchannels, const int nsamples,
const int64_t stride_row, const int64_t stride_channel, const int64_t stride_sample, const float eps, cudaStream_t stream,
const float scale_out = 1.0f) {
- const dim3 blocks_num(nrows, nchannels, nsamples);
+ const dim3 blocks_num(nrows, MIN(nchannels, UINT16_MAX), MIN(nsamples, UINT16_MAX));
if (ncols < 1024) {
const dim3 block_dims(256, 1, 1);
const ggml_cuda_kernel_launch_params launch_params = {blocks_num, block_dims, block_dims.x > WARP_SIZE ? 32 * sizeof(float): 0, stream};
ggml_cuda_kernel_launch(rms_norm_f32<256, false, false, do_scale>, launch_params,
- x, dst, ncols, stride_row, stride_channel, stride_sample, eps,
+ x, dst, ncols, nchannels, nsamples, stride_row, stride_channel, stride_sample, eps,
// underlying cudaLaunchKernelEx does not support default params
nullptr, 0, 0, 0, make_uint3(0, 0, 0), make_uint3(0, 0, 0), make_uint3(0, 0, 0), make_uint3(0, 0, 0),
nullptr, 0, 0, 0, make_uint3(0, 0, 0), make_uint3(0, 0, 0), make_uint3(0, 0, 0), make_uint3(0, 0, 0), scale_out);
} else {
const dim3 block_dims(1024, 1, 1);
const ggml_cuda_kernel_launch_params launch_params = ggml_cuda_kernel_launch_params{blocks_num, block_dims, block_dims.x > WARP_SIZE ? 32 * sizeof(float): 0, stream};
- ggml_cuda_kernel_launch(rms_norm_f32<1024, false, false, do_scale>, launch_params, x, dst, ncols, stride_row, stride_channel, stride_sample, eps,
+ ggml_cuda_kernel_launch(rms_norm_f32<1024, false, false, do_scale>, launch_params, x, dst, ncols, nchannels, nsamples, stride_row, stride_channel, stride_sample, eps,
// underlying cudaLaunchKernelEx does not support default params
nullptr, 0, 0, 0, make_uint3(0, 0, 0), make_uint3(0, 0, 0), make_uint3(0, 0, 0), make_uint3(0, 0, 0),
nullptr, 0, 0, 0, make_uint3(0, 0, 0), make_uint3(0, 0, 0), make_uint3(0, 0, 0), make_uint3(0, 0, 0), scale_out);
@@ -356,7 +383,7 @@ static void rms_norm_mul_f32_cuda(const float * x,
const uint32_t add_nsamples,
const float eps,
cudaStream_t stream) {
- const dim3 blocks_num(nrows, nchannels, nsamples);
+ const dim3 blocks_num(nrows, MIN(nchannels, UINT16_MAX), MIN(nsamples, UINT16_MAX));
if (mul == nullptr) {
rms_norm_f32_cuda(x, dst, ncols, nrows, nchannels, nsamples, stride_row, stride_channel, stride_sample, eps, stream);
return;
@@ -370,7 +397,7 @@ static void rms_norm_mul_f32_cuda(const float * x,
const dim3 block_dims(256, 1, 1);
const ggml_cuda_kernel_launch_params launch_params = ggml_cuda_kernel_launch_params{blocks_num, block_dims, block_dims.x > WARP_SIZE ? 32 * sizeof(float): 0, stream};
ggml_cuda_kernel_launch(rms_norm_f32<256, true>, launch_params,
- x, dst, ncols, stride_row, stride_channel, stride_sample, eps, mul, mul_stride_row, mul_stride_channel,
+ x, dst, ncols, nchannels, nsamples, stride_row, stride_channel, stride_sample, eps, mul, mul_stride_row, mul_stride_channel,
mul_stride_sample, mul_ncols_packed, mul_nrows_packed, mul_nchannels_packed, mul_nsamples_packed,
// underlying cudaLaunchKernelEx does not support default params
nullptr, 0, 0, 0, make_uint3(0, 0, 0), make_uint3(0, 0, 0), make_uint3(0, 0, 0), make_uint3(0, 0, 0), 1.0f);
@@ -378,7 +405,7 @@ static void rms_norm_mul_f32_cuda(const float * x,
const dim3 block_dims(1024, 1, 1);
const ggml_cuda_kernel_launch_params launch_params = ggml_cuda_kernel_launch_params{blocks_num, block_dims, block_dims.x > WARP_SIZE ? 32 * sizeof(float): 0, stream};
ggml_cuda_kernel_launch(rms_norm_f32<1024, true>, launch_params,
- x, dst, ncols, stride_row, stride_channel, stride_sample, eps, mul, mul_stride_row, mul_stride_channel,
+ x, dst, ncols, nchannels, nsamples, stride_row, stride_channel, stride_sample, eps, mul, mul_stride_row, mul_stride_channel,
mul_stride_sample, mul_ncols_packed, mul_nrows_packed, mul_nchannels_packed, mul_nsamples_packed,
// underlying cudaLaunchKernelEx does not support default params
nullptr, 0, 0, 0, make_uint3(0, 0, 0), make_uint3(0, 0, 0), make_uint3(0, 0, 0), make_uint3(0, 0, 0), 1.0f);
@@ -397,7 +424,7 @@ static void rms_norm_mul_f32_cuda(const float * x,
const dim3 block_dims(256, 1, 1);
const ggml_cuda_kernel_launch_params launch_params = ggml_cuda_kernel_launch_params{blocks_num, block_dims,block_dims.x > WARP_SIZE ? 32 * sizeof(float): 0, stream};
ggml_cuda_kernel_launch(rms_norm_f32<256, true, true>, launch_params,
- x, dst, ncols, stride_row, stride_channel, stride_sample, eps, mul, mul_stride_row, mul_stride_channel,
+ x, dst, ncols, nchannels, nsamples, stride_row, stride_channel, stride_sample, eps, mul, mul_stride_row, mul_stride_channel,
mul_stride_sample, mul_ncols_packed, mul_nrows_packed, mul_nchannels_packed, mul_nsamples_packed, add,
add_stride_row, add_stride_channel, add_stride_sample, add_ncols_packed, add_nrows_packed,
add_nchannels_packed, add_nsamples_packed, 1.0f);
@@ -405,7 +432,7 @@ static void rms_norm_mul_f32_cuda(const float * x,
const dim3 block_dims(1024, 1, 1);
const ggml_cuda_kernel_launch_params launch_params = ggml_cuda_kernel_launch_params{blocks_num, block_dims, block_dims.x > WARP_SIZE ? 32 * sizeof(float): 0, stream};
ggml_cuda_kernel_launch(rms_norm_f32<1024, true, true>, launch_params,
- x, dst, ncols, stride_row, stride_channel, stride_sample, eps, mul, mul_stride_row, mul_stride_channel,
+ x, dst, ncols, nchannels, nsamples, stride_row, stride_channel, stride_sample, eps, mul, mul_stride_row, mul_stride_channel,
mul_stride_sample, mul_ncols_packed, mul_nrows_packed, mul_nchannels_packed, mul_nsamples_packed, add,
add_stride_row, add_stride_channel, add_stride_sample, add_ncols_packed, add_nrows_packed,
add_nchannels_packed, add_nsamples_packed, 1.0f);
@@ -426,15 +453,15 @@ static void rms_norm_back_f32_cuda(const float * grad, const float * xf, float *
static void l2_norm_f32_cuda(
const float * x, float * dst, const int ncols, const int nrows, const int nchannels, const int nsamples,
const int64_t stride_row, const int64_t stride_channel, const int64_t stride_sample, const float eps, cudaStream_t stream) {
- const dim3 blocks_num(nrows, nchannels, nsamples);
+ const dim3 blocks_num(nrows, MIN(nchannels, UINT16_MAX), MIN(nsamples, UINT16_MAX));
if (ncols < 1024) {
const dim3 block_dims(WARP_SIZE, 1, 1);
const ggml_cuda_kernel_launch_params launch_params = ggml_cuda_kernel_launch_params{blocks_num, block_dims, 0, stream};
- ggml_cuda_kernel_launch(l2_norm_f32, launch_params, x, dst, ncols, stride_row, stride_channel, stride_sample, eps);
+ ggml_cuda_kernel_launch(l2_norm_f32, launch_params, x, dst, ncols, nchannels, nsamples, stride_row, stride_channel, stride_sample, eps);
} else {
const dim3 block_dims(1024, 1, 1);
const ggml_cuda_kernel_launch_params launch_params = ggml_cuda_kernel_launch_params{blocks_num, block_dims, block_dims.x > WARP_SIZE ? 32 * sizeof(float): 0, stream};
- ggml_cuda_kernel_launch(l2_norm_f32<1024>, launch_params, x, dst, ncols, stride_row, stride_channel, stride_sample, eps);
+ ggml_cuda_kernel_launch(l2_norm_f32<1024>, launch_params, x, dst, ncols, nchannels, nsamples, stride_row, stride_channel, stride_sample, eps);
}
}
diff --git a/ggml/src/ggml-cuda/roll.cu b/ggml/src/ggml-cuda/roll.cu
index a339dfc1a..c9bbe76f0 100644
--- a/ggml/src/ggml-cuda/roll.cu
+++ b/ggml/src/ggml-cuda/roll.cu
@@ -17,6 +17,10 @@ static __global__ void roll_f32_cuda(const float * __restrict__ src,
const int64_t ne01,
const int64_t ne02,
const int64_t ne03,
+ const int64_t nb00,
+ const int64_t nb01,
+ const int64_t nb02,
+ const int64_t nb03,
const int s0,
const int s1,
const int s2,
@@ -39,7 +43,7 @@ static __global__ void roll_f32_cuda(const float * __restrict__ src,
const int64_t d3 = wrap_index(i3 - s3, ne03);
dst[i3 * (ne00 * ne01 * ne02) + i2 * (ne01 * ne00) + i1 * ne00 + i0] =
- src[d3 * (ne00 * ne01 * ne02) + d2 * (ne01 * ne00) + d1 * ne00 + d0];
+ src[(d3 * nb03 + d2 * nb02 + d1 * nb01 + d0 * nb00) / sizeof(float)];
}
void ggml_cuda_op_roll(ggml_backend_cuda_context & ctx, ggml_tensor * dst) {
@@ -63,5 +67,5 @@ void ggml_cuda_op_roll(ggml_backend_cuda_context & ctx, ggml_tensor * dst) {
int64_t num_blocks = (sz + CUDA_ROLL_BLOCK_SIZE - 1) / CUDA_ROLL_BLOCK_SIZE;
roll_f32_cuda<<>>(
- src0_d, dst_d, ne00, ne01, ne02, ne03, s0, s1, s2, s3);
+ src0_d, dst_d, ne00, ne01, ne02, ne03, nb00, nb01, nb02, nb03, s0, s1, s2, s3);
}
diff --git a/ggml/src/ggml-cuda/rope.cu b/ggml/src/ggml-cuda/rope.cu
index e546fb655..b4fca778f 100644
--- a/ggml/src/ggml-cuda/rope.cu
+++ b/ggml/src/ggml-cuda/rope.cu
@@ -709,7 +709,7 @@ void ggml_cuda_op_rope_fused(ggml_backend_cuda_context & ctx, ggml_tensor * rope
// one block per row: block_reduce gives the norm scale, then each thread applies mul and rope to the elements it owns
template
static __global__ void rms_norm_mul_rope_f32(
- const float * x, D * dst, const int ncols,
+ const float * x, D * dst, const int ncols, const int nchannels, const int nsamples,
const int64_t s01, const int64_t s02, const int64_t s03,
const int64_t s1, const int64_t s2, const int64_t s3,
const float eps,
@@ -724,66 +724,76 @@ static __global__ void rms_norm_mul_rope_f32(
const int64_t * row_indices, const int set_rows_stride,
const bool is_neox) {
ggml_cuda_pdl_lc();
- const int row = blockIdx.x;
- const int channel = blockIdx.y;
- const int sample = blockIdx.z;
- const int tid = threadIdx.x;
-
- x += sample*s03 + channel*s02 + row*s01;
-
- const uint32_t mul_row = fastmodulo(row, mul_nrows_packed);
- const uint32_t mul_channel = fastmodulo(channel, mul_nchannels_packed);
- const uint32_t mul_sample = fastmodulo(sample, mul_nsamples_packed);
- mul += mul_sample*mul_s03 + mul_channel*mul_s02 + mul_row*mul_s01;
-
- float tmp = 0.0f;
-
- ggml_cuda_pdl_sync();
- for (int col = tid; col < ncols; col += block_size) {
- const float xi = x[col];
- tmp += xi * xi;
- }
+ const int row = blockIdx.x;
+ const int tid = threadIdx.x;
extern __shared__ float s_sum[];
- tmp = block_reduce(tmp, s_sum);
- const float scale = rsqrtf(tmp/ncols + eps);
+ ggml_cuda_pdl_sync();
- int64_t idst = sample*s3 + channel*s2 + row*s1;
- if (set_rows_stride != 0) {
- idst = row*s1 + row_indices[channel]*set_rows_stride;
- }
- dst += idst;
+ // grid.y and grid.z are clamped to the CUDA limit, iterate over the excess channels/samples
+ for (int sample = blockIdx.z; sample < nsamples; sample += gridDim.z) {
+ for (int channel = blockIdx.y; channel < nchannels; channel += gridDim.y) {
+ const float * xc = x + sample*s03 + channel*s02 + row*s01;
- for (int i0 = 2*tid; i0 < ncols; i0 += 2*block_size) {
- int ix0;
- int ix1;
- if (is_neox && i0 < n_dims) {
- ix0 = i0/2;
- ix1 = i0/2 + n_dims/2;
- } else {
- ix0 = i0 + 0;
- ix1 = i0 + 1;
+ const uint32_t mul_row = fastmodulo(row, mul_nrows_packed);
+ const uint32_t mul_channel = fastmodulo(channel, mul_nchannels_packed);
+ const uint32_t mul_sample = fastmodulo(sample, mul_nsamples_packed);
+ const float * mulc = mul + mul_sample*mul_s03 + mul_channel*mul_s02 + mul_row*mul_s01;
+
+ float tmp = 0.0f;
+
+ for (int col = tid; col < ncols; col += block_size) {
+ const float xi = xc[col];
+ tmp += xi * xi;
+ }
+
+ tmp = block_reduce(tmp, s_sum);
+
+ const float scale = rsqrtf(tmp/ncols + eps);
+
+ int64_t idst = sample*s3 + channel*s2 + row*s1;
+ if (set_rows_stride != 0) {
+ idst = row*s1 + row_indices[channel]*set_rows_stride;
+ }
+ D * dstc = dst + idst;
+
+ for (int i0 = 2*tid; i0 < ncols; i0 += 2*block_size) {
+ int ix0;
+ int ix1;
+ if (is_neox && i0 < n_dims) {
+ ix0 = i0/2;
+ ix1 = i0/2 + n_dims/2;
+ } else {
+ ix0 = i0 + 0;
+ ix1 = i0 + 1;
+ }
+
+ const float x0 = scale * xc[ix0] * mulc[fastmodulo(ix0, mul_ncols_packed)];
+ const float x1 = scale * xc[ix1] * mulc[fastmodulo(ix1, mul_ncols_packed)];
+
+ if (i0 >= n_dims) {
+ dstc[ix0] = ggml_cuda_cast(x0);
+ dstc[ix1] = ggml_cuda_cast(x1);
+ continue;
+ }
+
+ const float theta_base = pos[channel]*powf(theta_scale, i0/2.0f);
+ const float freq_factor = has_ff ? freq_factors[i0/2] : 1.0f;
+
+ float cos_theta;
+ float sin_theta;
+ rope_yarn(theta_base/freq_factor, freq_scale, corr_dims, i0, ext_factor, attn_factor, cos_theta, sin_theta);
+
+ dstc[ix0] = ggml_cuda_cast(x0*cos_theta - x1*sin_theta);
+ dstc[ix1] = ggml_cuda_cast(x0*sin_theta + x1*cos_theta);
+ }
+
+ if constexpr (block_size > WARP_SIZE) {
+ // sync is needed as we reuse s_sum across block_reduce invocations, see #26385
+ __syncthreads();
+ }
}
-
- const float x0 = scale * x[ix0] * mul[fastmodulo(ix0, mul_ncols_packed)];
- const float x1 = scale * x[ix1] * mul[fastmodulo(ix1, mul_ncols_packed)];
-
- if (i0 >= n_dims) {
- dst[ix0] = ggml_cuda_cast(x0);
- dst[ix1] = ggml_cuda_cast(x1);
- continue;
- }
-
- const float theta_base = pos[channel]*powf(theta_scale, i0/2.0f);
- const float freq_factor = has_ff ? freq_factors[i0/2] : 1.0f;
-
- float cos_theta;
- float sin_theta;
- rope_yarn(theta_base/freq_factor, freq_scale, corr_dims, i0, ext_factor, attn_factor, cos_theta, sin_theta);
-
- dst[ix0] = ggml_cuda_cast(x0*cos_theta - x1*sin_theta);
- dst[ix1] = ggml_cuda_cast(x0*sin_theta + x1*cos_theta);
}
}
@@ -806,7 +816,7 @@ static void rms_norm_mul_rope_cuda(
const bool is_neox, cudaStream_t stream) {
GGML_ASSERT(ncols % 2 == 0);
- const dim3 blocks_num(nrows, nchannels, nsamples);
+ const dim3 blocks_num(nrows, MIN(nchannels, UINT16_MAX), MIN(nsamples, UINT16_MAX));
const float theta_scale = powf(freq_base, -2.0f/n_dims);
@@ -820,13 +830,13 @@ static void rms_norm_mul_rope_cuda(
const ggml_cuda_kernel_launch_params launch_params = {blocks_num, block_dims, 32*sizeof(float), stream};
if (freq_factors == nullptr) {
ggml_cuda_kernel_launch(rms_norm_mul_rope_f32<256, false, D>, launch_params,
- x, dst, ncols, s01, s02, s03, s1, s2, s3, eps, mul, mul_s01, mul_s02, mul_s03,
+ x, dst, ncols, nchannels, nsamples, s01, s02, s03, s1, s2, s3, eps, mul, mul_s01, mul_s02, mul_s03,
mul_ncols_packed, mul_nrows_packed, mul_nchannels_packed, mul_nsamples_packed,
n_dims, pos, freq_scale, ext_factor, attn_factor, corr_dims, theta_scale,
freq_factors, row_indices, set_rows_stride, is_neox);
} else {
ggml_cuda_kernel_launch(rms_norm_mul_rope_f32<256, true, D>, launch_params,
- x, dst, ncols, s01, s02, s03, s1, s2, s3, eps, mul, mul_s01, mul_s02, mul_s03,
+ x, dst, ncols, nchannels, nsamples, s01, s02, s03, s1, s2, s3, eps, mul, mul_s01, mul_s02, mul_s03,
mul_ncols_packed, mul_nrows_packed, mul_nchannels_packed, mul_nsamples_packed,
n_dims, pos, freq_scale, ext_factor, attn_factor, corr_dims, theta_scale,
freq_factors, row_indices, set_rows_stride, is_neox);
@@ -836,13 +846,13 @@ static void rms_norm_mul_rope_cuda(
const ggml_cuda_kernel_launch_params launch_params = {blocks_num, block_dims, 32*sizeof(float), stream};
if (freq_factors == nullptr) {
ggml_cuda_kernel_launch(rms_norm_mul_rope_f32<1024, false, D>, launch_params,
- x, dst, ncols, s01, s02, s03, s1, s2, s3, eps, mul, mul_s01, mul_s02, mul_s03,
+ x, dst, ncols, nchannels, nsamples, s01, s02, s03, s1, s2, s3, eps, mul, mul_s01, mul_s02, mul_s03,
mul_ncols_packed, mul_nrows_packed, mul_nchannels_packed, mul_nsamples_packed,
n_dims, pos, freq_scale, ext_factor, attn_factor, corr_dims, theta_scale,
freq_factors, row_indices, set_rows_stride, is_neox);
} else {
ggml_cuda_kernel_launch(rms_norm_mul_rope_f32<1024, true, D>, launch_params,
- x, dst, ncols, s01, s02, s03, s1, s2, s3, eps, mul, mul_s01, mul_s02, mul_s03,
+ x, dst, ncols, nchannels, nsamples, s01, s02, s03, s1, s2, s3, eps, mul, mul_s01, mul_s02, mul_s03,
mul_ncols_packed, mul_nrows_packed, mul_nchannels_packed, mul_nsamples_packed,
n_dims, pos, freq_scale, ext_factor, attn_factor, corr_dims, theta_scale,
freq_factors, row_indices, set_rows_stride, is_neox);
diff --git a/ggml/src/ggml-cuda/unary.cu b/ggml/src/ggml-cuda/unary.cu
index c9a8cf6b0..436fdfffb 100644
--- a/ggml/src/ggml-cuda/unary.cu
+++ b/ggml/src/ggml-cuda/unary.cu
@@ -134,24 +134,67 @@ static void unary_cuda(const T * x, T * dst, const int k, cudaStream_t stream) {
ggml_cuda_kernel_launch(unary_op_kernel, launch_params, x, dst, k);
}
+template
+static __global__ void unary_op_kernel_strided(const T * x, T * dst, const int k,const int64_t ne00,const int64_t ne01,const int64_t ne02,const size_t nb00,const size_t nb01,const size_t nb02,const size_t nb03) {
+ ggml_cuda_pdl_lc();
+ const int i = blockDim.x*blockIdx.x + threadIdx.x;
+
+ if (i >= k) {
+ return;
+ }
+
+ int64_t rem = i;
+ const int64_t i0 = rem % ne00; rem /= ne00;
+ const int64_t i1 = rem % ne01; rem /= ne01;
+ const int64_t i2 = rem % ne02;
+ const int64_t i3 = rem / ne02;
+ const size_t src_byte_offset = i0 * nb00 + i1 * nb01 + i2 * nb02 + i3 * nb03;
+ const T * src_ptr = (const T *)((const char *)x + src_byte_offset);
+
+ ggml_cuda_pdl_sync();
+ dst[i] = ggml_cuda_cast(op(ggml_cuda_cast(*src_ptr)));
+}
+
+template
+static void unary_cuda_strided(const T * x, T * dst, const int k,const int64_t ne00,const int64_t ne01,const int64_t ne02,const size_t nb00,const size_t nb01,const size_t nb02,const size_t nb03, cudaStream_t stream) {
+ const int num_blocks = (k + CUDA_NEG_BLOCK_SIZE - 1) / CUDA_NEG_BLOCK_SIZE;
+ const ggml_cuda_kernel_launch_params launch_params = ggml_cuda_kernel_launch_params((dim3)num_blocks, CUDA_NEG_BLOCK_SIZE, 0, stream);
+ ggml_cuda_kernel_launch(unary_op_kernel_strided, launch_params, x, dst, k, ne00,ne01,ne02,nb00,nb01,nb02,nb03);
+}
+
template
void ggml_cuda_op_unary(ggml_backend_cuda_context & ctx, ggml_tensor * dst) {
const ggml_tensor * src0 = dst->src[0];
const void * src0_d = src0->data;
void * dst_d = dst->data;
- cudaStream_t stream = ctx.stream();
- GGML_ASSERT(ggml_is_contiguous(src0));
+ cudaStream_t stream = ctx.stream();
GGML_ASSERT(src0->type == GGML_TYPE_F32 || src0->type == GGML_TYPE_F16 || src0->type == GGML_TYPE_BF16);
GGML_ASSERT(src0->type == dst->type);
- if (src0->type == GGML_TYPE_F16) {
- unary_cuda((const half *)src0_d, (half *)dst_d, ggml_nelements(src0), stream);
- } else if (src0->type == GGML_TYPE_BF16) {
- unary_cuda((const nv_bfloat16 *)src0_d, (nv_bfloat16 *)dst_d, ggml_nelements(src0), stream);
+ if (ggml_is_contiguous(src0)) {
+ if (src0->type == GGML_TYPE_F16) {
+ unary_cuda((const half *)src0_d, (half *)dst_d, ggml_nelements(src0), stream);
+ } else if (src0->type == GGML_TYPE_BF16) {
+ unary_cuda((const nv_bfloat16 *)src0_d, (nv_bfloat16 *)dst_d, ggml_nelements(src0), stream);
+ } else {
+ unary_cuda((const float *)src0_d, (float *)dst_d, ggml_nelements(src0), stream);
+ }
} else {
- unary_cuda((const float *)src0_d, (float *)dst_d, ggml_nelements(src0), stream);
+ if (src0->type == GGML_TYPE_F16) {
+ unary_cuda_strided((const half *)src0_d, (half *)dst_d, ggml_nelements(src0),
+ src0->ne[0], src0->ne[1], src0->ne[2],
+ src0->nb[0], src0->nb[1], src0->nb[2], src0->nb[3], stream);
+ } else if (src0->type == GGML_TYPE_BF16) {
+ unary_cuda_strided((const nv_bfloat16 *)src0_d, (nv_bfloat16 *)dst_d, ggml_nelements(src0),
+ src0->ne[0], src0->ne[1], src0->ne[2],
+ src0->nb[0], src0->nb[1], src0->nb[2], src0->nb[3], stream);
+ } else {
+ unary_cuda_strided((const float *)src0_d, (float *)dst_d, ggml_nelements(src0),
+ src0->ne[0], src0->ne[1], src0->ne[2],
+ src0->nb[0], src0->nb[1], src0->nb[2], src0->nb[3], stream);
+ }
}
}
diff --git a/ggml/src/ggml-hexagon/htp/matmul-ops.c b/ggml/src/ggml-hexagon/htp/matmul-ops.c
index a589e0a20..1a513d160 100644
--- a/ggml/src/ggml-hexagon/htp/matmul-ops.c
+++ b/ggml/src/ggml-hexagon/htp/matmul-ops.c
@@ -425,7 +425,7 @@ static void hvx_mv_2d_repacked_##SUFFIX(unsigned int nth, unsigned int ith, void
if (copy_cnt > 0) { \
htp_trace_event_start(tr, HTP_TRACE_EVT_HVX_COMP, ct_end); \
if (src2) { \
- hvx_add_f32_uaa((uint8_t *) &dst_col[src0_start_row], \
+ hvx_add_f32_uuu((uint8_t *) &dst_col[src0_start_row], \
(const uint8_t *) tmp, \
(const uint8_t *) ((const float *) mmctx->vtcm_bias + src0_start_row), \
copy_cnt); \
@@ -1108,7 +1108,7 @@ static void hvx_mv_2d(unsigned int nth, unsigned int ith, void * data) {
if (copy_cnt > 0) {
htp_trace_event_start(tr, HTP_TRACE_EVT_HVX_COMP, src0_end_row);
if (src2) {
- hvx_add_f32_uaa((uint8_t *) &dst_col[src0_start_row],
+ hvx_add_f32_uuu((uint8_t *) &dst_col[src0_start_row],
(const uint8_t *) tmp,
(const uint8_t *) ((const float *) mmctx->vtcm_bias + src0_start_row),
copy_cnt);
diff --git a/ggml/src/ggml-sycl/CMakeLists.txt b/ggml/src/ggml-sycl/CMakeLists.txt
index d2196f74d..1697ecd95 100644
--- a/ggml/src/ggml-sycl/CMakeLists.txt
+++ b/ggml/src/ggml-sycl/CMakeLists.txt
@@ -221,4 +221,27 @@ if (GGML_SYCL_DEVICE_ARCH)
"SHELL:-Xsycl-target-backend=spir64_gen \"-device ${GGML_SYCL_DEVICE_ARCH}\""
-fsycl-max-parallel-link-jobs=${GGML_SYCL_MAX_PARALLEL_LINK_JOBS}
)
+
+ # The XMX dequant-GEMM tiles need the sub-group size of the target: 8 on Xe-HPG (DG2, ARL-H),
+ # 16 on Xe-HPC and Xe2 or newer. ocloc fails on the other size, so build only the one that fits.
+ # 0 (unknown name, mixed list, or no XMX) builds no XMX tile and the path stays off.
+ set(_ggml_sycl_xmx_sg "")
+ string(TOLOWER "${GGML_SYCL_DEVICE_ARCH}" _ggml_sycl_archs)
+ string(REPLACE "," ";" _ggml_sycl_archs "${_ggml_sycl_archs}")
+ foreach(_arch IN LISTS _ggml_sycl_archs)
+ if (_arch MATCHES "^(dg2|acm|ats-m|arl-h|xe-hpg|12\\.5[567]\\.|12\\.74\\.)")
+ set(_sg 8)
+ elseif (_arch MATCHES "^(pvc|bmg|lnl|ptl|wcl|nvl|cri|xe2|xe3|xe-hpc|12\\.60\\.|20\\.|30\\.)")
+ set(_sg 16)
+ else()
+ set(_sg 0)
+ endif()
+ if (_ggml_sycl_xmx_sg STREQUAL "" OR _ggml_sycl_xmx_sg EQUAL _sg)
+ set(_ggml_sycl_xmx_sg ${_sg})
+ else()
+ set(_ggml_sycl_xmx_sg 0)
+ endif()
+ endforeach()
+ message(STATUS "GGML_SYCL_DEVICE_ARCH: XMX dequant-GEMM sub-group size ${_ggml_sycl_xmx_sg} (0 = off)")
+ target_compile_definitions(ggml-sycl PRIVATE GGML_SYCL_XMX_AOT_SG=${_ggml_sycl_xmx_sg})
endif()
diff --git a/ggml/src/ggml-sycl/common.hpp b/ggml/src/ggml-sycl/common.hpp
index 5e904bd2a..475509372 100644
--- a/ggml/src/ggml-sycl/common.hpp
+++ b/ggml/src/ggml-sycl/common.hpp
@@ -65,6 +65,61 @@ extern int g_ggml_sycl_enable_fusion;
extern int g_ggml_sycl_enable_esimd;
extern int g_ggml_sycl_mmvq_wide;
extern int g_ggml_sycl_prioritize_dmmv;
+
+// Which quantized weight formats may take the XMX dequant-GEMM paths. A bitmask rather than one
+// flag per path, so a format can be enabled or measured on its own and adding a format is one bit.
+enum ggml_sycl_xmx_gather_type {
+ GGML_SYCL_XMX_GATHER_IQ4_NL = 1 << 0,
+ GGML_SYCL_XMX_GATHER_IQ3_S = 1 << 1,
+ GGML_SYCL_XMX_GATHER_IQ4_XS = 1 << 2,
+ GGML_SYCL_XMX_GATHER_IQ3_XXS = 1 << 3,
+ GGML_SYCL_XMX_GATHER_IQ2_XXS = 1 << 4,
+ GGML_SYCL_XMX_GATHER_IQ2_XS = 1 << 5,
+ GGML_SYCL_XMX_GATHER_IQ2_S = 1 << 6,
+ GGML_SYCL_XMX_GATHER_IQ1_S = 1 << 7,
+ GGML_SYCL_XMX_GATHER_IQ1_M = 1 << 8,
+ GGML_SYCL_XMX_GATHER_Q8_0 = 1 << 9,
+ GGML_SYCL_XMX_GATHER_Q4_K = 1 << 10,
+ GGML_SYCL_XMX_GATHER_Q5_K = 1 << 11,
+ GGML_SYCL_XMX_GATHER_Q6_K = 1 << 12,
+};
+static constexpr int GGML_SYCL_XMX_GATHER_TYPES_DEFAULT = ~0;
+extern int g_ggml_sycl_xmx_gather_types;
+// Which joint_matrix combinations the XMX dequant-GEMM paths may use, one bit each (see fused-gemm.cpp).
+// GGML_SYCL_DYNAMIC_PRECISION picks the operand type, this mask the combinations of that type.
+static constexpr int GGML_SYCL_XMX_GATHER_SHAPES_DEFAULT = 0xff;
+extern int g_ggml_sycl_xmx_gather_shapes;
+
+// GGML_SYCL_DYNAMIC_PRECISION: operand type of the XMX dequant-GEMM paths. F32 turns them off and
+// keeps the library GEMM in f32. A src1 precision request of an op [TAG_GGML_PREC] is always met.
+enum ggml_sycl_dynamic_precision {
+ GGML_SYCL_DYNAMIC_PRECISION_F16,
+ GGML_SYCL_DYNAMIC_PRECISION_BF16,
+ GGML_SYCL_DYNAMIC_PRECISION_TF32,
+ GGML_SYCL_DYNAMIC_PRECISION_F32,
+};
+#ifdef GGML_SYCL_F16
+static constexpr int GGML_SYCL_DYNAMIC_PRECISION_DEFAULT = GGML_SYCL_DYNAMIC_PRECISION_F16;
+#else
+static constexpr int GGML_SYCL_DYNAMIC_PRECISION_DEFAULT = GGML_SYCL_DYNAMIC_PRECISION_F32;
+#endif
+extern int g_ggml_sycl_dynamic_precision;
+// GGML_SYCL_DYNAMIC_REQUIRED_PRECISION: the XMX type an F32 src1 request may run on instead of f32
+// (TF32, or BF16 which also allows tf32). F32 (default): none. F16: src1 requests are ignored.
+extern int g_ggml_sycl_dynamic_required_precision;
+
+// [TAG_GGML_PREC] src1 precision request of the MUL_MAT/MUL_MAT_ID op dst
+static inline int32_t ggml_sycl_src1_prec(const ggml_tensor * dst) {
+ return g_ggml_sycl_dynamic_required_precision == GGML_SYCL_DYNAMIC_PRECISION_F16 ? GGML_PREC_UNDEFINED :
+ dst->op_params[3];
+}
+
+// [TAG_GGML_PREC] the library GEMM and dmmv may convert src1 of the MUL_MAT/MUL_MAT_ID op dst to f16
+static inline bool ggml_sycl_src1_f16_ok(const ggml_tensor * dst) {
+ const int32_t src1_prec = ggml_sycl_src1_prec(dst);
+ return g_ggml_sycl_dynamic_precision != GGML_SYCL_DYNAMIC_PRECISION_F32 &&
+ (src1_prec == GGML_PREC_UNDEFINED || src1_prec >= GGML_PREC_F16);
+}
extern int g_ggml_sycl_enable_flash_attention;
extern int g_ggml_sycl_dev2dev_memcpy;
extern int g_ggml_sycl_fa_onednn;
@@ -333,6 +388,12 @@ struct mmid_row_mapping {
int32_t i2;
};
+struct ggml_sycl_gg_tile {
+ int32_t expert;
+ int32_t n0;
+ int32_t n1;
+};
+
namespace sycl_ex = sycl::ext::oneapi::experimental;
struct ggml_backend_sycl_context {
int device;
@@ -410,6 +471,7 @@ struct ggml_backend_sycl_context {
std::unique_ptr host_pools[GGML_SYCL_MAX_DEVICES];
std::vector mmid_row_mapping_host;
+ std::vector mmid_tile_schedule_host;
static std::unique_ptr new_pool_for_device(queue_ptr qptr, int device);
diff --git a/ggml/src/ggml-sycl/element_wise.cpp b/ggml/src/ggml-sycl/element_wise.cpp
index 2e926abea..0d263aa76 100644
--- a/ggml/src/ggml-sycl/element_wise.cpp
+++ b/ggml/src/ggml-sycl/element_wise.cpp
@@ -452,6 +452,46 @@ static void unary_mul_sycl(const T * x, const T * g, T * dst, const int64_t k, c
});
}
+// ADD(bias) + UNARY + MUL(scale) with both broadcast over dim 0, the delta-net alpha gate:
+// dst[i] = op(a[i] + bias[i % ne0]) * scale[i % ne0]. k == ne0 makes that the flat index.
+template
+static void add_unary_mul_flat_kernel(const float * a, const float * bias, const float * scale, float * dst,
+ const int64_t k, const sycl::nd_item<1> &item_ct1, F op) {
+ SYCL_GLOBAL_ID_LOOP(k, item_ct1) {
+ dst[i] = op(a[i] + bias[i]) * scale[i];
+ }
+}
+
+template
+static void add_unary_mul_bcast_kernel(const float * a, const float * bias, const float * scale, float * dst,
+ const int64_t k, const sycl::uint3 ne0_fd, const sycl::nd_item<1> &item_ct1, F op) {
+ SYCL_GLOBAL_ID_LOOP(k, item_ct1) {
+ const uint32_t h = fastmodulo((uint32_t) i, ne0_fd);
+ dst[i] = op(a[i] + bias[h]) * scale[h];
+ }
+}
+
+template
+static void add_unary_mul_sycl(const float * a, const float * bias, const float * scale, float * dst,
+ const int64_t k, const int64_t ne0, queue_ptr main_stream, F op) {
+ const size_t num_blocks = ceil_div((size_t) k, (size_t) SYCL_GLU_BLOCK_SIZE);
+ const sycl::nd_range<1> range(num_blocks * sycl::range<1>(SYCL_GLU_BLOCK_SIZE), sycl::range<1>(SYCL_GLU_BLOCK_SIZE));
+
+ if (k == ne0) {
+ main_stream->parallel_for(range, [=](sycl::nd_item<1> item_ct1) [[sycl::reqd_sub_group_size(WARP_SIZE)]] {
+ add_unary_mul_flat_kernel(a, bias, scale, dst, k, item_ct1, op);
+ });
+ return;
+ }
+
+ // 32-bit fastdiv, exact only below 2^31; ggml_sycl_can_fuse() already declined past that
+ GGML_ASSERT(k < ((int64_t) 1 << 31));
+ const sycl::uint3 ne0_fd = init_fastdiv_values((uint32_t) ne0);
+ main_stream->parallel_for(range, [=](sycl::nd_item<1> item_ct1) [[sycl::reqd_sub_group_size(WARP_SIZE)]] {
+ add_unary_mul_bcast_kernel(a, bias, scale, dst, k, ne0_fd, item_ct1, op);
+ });
+}
+
namespace ggml_sycl_detail {
static void acc_f32_sycl(const char *x, const char *y, float *dst,
const int64_t n_elements,
@@ -995,6 +1035,19 @@ static inline void ggml_sycl_op_swiglu(ggml_backend_sycl_context & ctx, ggml_ten
});
}
+// Hands `launch` the functor for the unary op of a fused unary chain. Anything else
+// ggml_sycl_can_fuse() has already declined, so the default is a dispatcher bug.
+template
+static void dispatch_fused_unary_op(ggml_unary_op uop, F && launch) {
+ switch (uop) {
+ case GGML_UNARY_OP_SILU: launch([](float v) { return op_silu(v); }); break;
+ case GGML_UNARY_OP_SIGMOID: launch([](float v) { return op_sigmoid(v); }); break;
+ case GGML_UNARY_OP_SOFTPLUS: launch([](float v) { return op_softplus(v); }); break;
+ default:
+ GGML_ABORT("fused unary chain: unsupported unary op %s", ggml_unary_op_name(uop));
+ }
+}
+
// dst = op(unary_node->src[0]) * other, written straight to the MUL output, saving the
// standalone unary launch. Preconditions come from ggml_sycl_can_fuse(); re-asserted here.
void ggml_sycl_op_unary_mul_fused(ggml_backend_sycl_context & ctx, ggml_tensor * unary_node, ggml_tensor * mul_node) {
@@ -1032,13 +1085,41 @@ void ggml_sycl_op_unary_mul_fused(ggml_backend_sycl_context & ctx, ggml_tensor *
}
};
- switch (ggml_get_unary_op(unary_node)) {
- case GGML_UNARY_OP_SILU: dispatch_type([](float v) { return op_silu(v); }); break;
- case GGML_UNARY_OP_SIGMOID: dispatch_type([](float v) { return op_sigmoid(v); }); break;
- case GGML_UNARY_OP_SOFTPLUS: dispatch_type([](float v) { return op_softplus(v); }); break;
- default:
- GGML_ABORT("fused unary+mul: unsupported unary op %s", ggml_unary_op_name(ggml_get_unary_op(unary_node)));
- }
+ dispatch_fused_unary_op(ggml_get_unary_op(unary_node), dispatch_type);
+}
+
+// dst = op(a + bias) * scale for an ADD + UNARY + MUL chain whose bias and scale broadcast
+// over dim 0. Preconditions come from ggml_sycl_can_fuse(); re-asserted here.
+void ggml_sycl_op_add_unary_mul_fused(ggml_backend_sycl_context & ctx, ggml_tensor * add_node,
+ ggml_tensor * unary_node, ggml_tensor * mul_node) {
+ // the dst-arity convention the other fusions follow; a and bias live on add_node
+ scope_op_debug_print scope_dbg_print(__func__, mul_node, /*num_src=*/2);
+
+ const ggml_tensor * a = add_node->src[0];
+ const ggml_tensor * bias = add_node->src[1];
+ const ggml_tensor * scale = (mul_node->src[0] == unary_node) ? mul_node->src[1] : mul_node->src[0];
+
+ // scale is picked by elimination; ggml_can_fuse()'s single-use rule rules out MUL(unary, unary)
+ GGML_ASSERT(scale != unary_node);
+ GGML_ASSERT(a->type == GGML_TYPE_F32 && bias->type == GGML_TYPE_F32);
+ GGML_ASSERT(scale->type == GGML_TYPE_F32 && mul_node->type == GGML_TYPE_F32);
+ GGML_ASSERT(ggml_are_same_shape(a, mul_node));
+ // a and dst are indexed flat
+ GGML_ASSERT(ggml_is_contiguous(a) && ggml_is_contiguous(mul_node));
+ // bias and scale are one contiguous ne0-length row each, broadcast over the outer dims
+ GGML_ASSERT(bias->ne[0] == a->ne[0] && scale->ne[0] == a->ne[0]);
+ GGML_ASSERT(ggml_nrows(bias) == 1 && ggml_nrows(scale) == 1);
+ GGML_ASSERT(ggml_is_contiguous(bias) && ggml_is_contiguous(scale));
+
+ queue_ptr main_stream = ctx.stream();
+ SYCL_CHECK(ggml_sycl_set_device(ctx.device));
+
+ const auto dispatch_op = [&](auto op) {
+ add_unary_mul_sycl((const float *) a->data, (const float *) bias->data, (const float *) scale->data,
+ (float *) mul_node->data, ggml_nelements(mul_node), mul_node->ne[0], main_stream, op);
+ };
+
+ dispatch_fused_unary_op(ggml_get_unary_op(unary_node), dispatch_op);
}
__dpct_inline__ float ggml_sycl_op_swiglu_oai_single(float x, float g, float alpha = 1.702f, float limit = 7.0f) {
diff --git a/ggml/src/ggml-sycl/element_wise.hpp b/ggml/src/ggml-sycl/element_wise.hpp
index d280066ef..0a5bfa956 100644
--- a/ggml/src/ggml-sycl/element_wise.hpp
+++ b/ggml/src/ggml-sycl/element_wise.hpp
@@ -132,4 +132,9 @@ void ggml_sycl_arange(ggml_backend_sycl_context & ctx, ggml_tensor * dst);
// fused UNARY(silu|sigmoid|softplus) + MUL; see ggml_sycl_can_fuse() for the accepted shapes
void ggml_sycl_op_unary_mul_fused(ggml_backend_sycl_context & ctx, ggml_tensor * unary_node, ggml_tensor * mul_node);
+// fused f32 ADD + UNARY(silu|sigmoid|softplus) + MUL with the bias and the scale broadcast
+// over dim 0; see ggml_sycl_can_fuse() for the accepted shapes
+void ggml_sycl_op_add_unary_mul_fused(ggml_backend_sycl_context & ctx, ggml_tensor * add_node,
+ ggml_tensor * unary_node, ggml_tensor * mul_node);
+
#endif // GGML_SYCL_ELEMENTWISE_HPP
diff --git a/ggml/src/ggml-sycl/fused-gemm.cpp b/ggml/src/ggml-sycl/fused-gemm.cpp
new file mode 100644
index 000000000..38c69bc6b
--- /dev/null
+++ b/ggml/src/ggml-sycl/fused-gemm.cpp
@@ -0,0 +1,1093 @@
+#include "fused-gemm.hpp"
+
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace mx = sycl::ext::oneapi::experimental::matrix;
+
+// FG_ / fg_ is short for fused GEMM: the weights are dequantized inside the GEMM, into the XMX tiles.
+
+// A k step is one 32-value weight sub-block; iq3_s and the other superblock formats split their
+// superblock into steps of this width. The sub-groups of a work-group each walk their own K range
+// and are summed at the end.
+static constexpr int FG_BK = QK4_NL;
+static constexpr int FG_KSPLIT = 4;
+
+// Element traits of one joint_matrix operand type. The A stage and the B pack compute in f32 and
+// convert once, in registers, when they write the element, so any type costs the same one pass.
+// store: storage in SLM (A) and in the packed B buffer
+// mtype: matrix_type in matrix_combinations
+// mode: GGML_SYCL_DYNAMIC_PRECISION value that selects this type
+// src: ggml type that needs no conversion into this type (GGML_TYPE_COUNT: none)
+// slow: XMX throughput class, 0 is fastest. f16 and bf16 share the DPAS rate; tf32 does half the
+// K per instruction. B60, Qwen3-30B-A3B pp512: f16 1108, bf16 1000, tf32 751 t/s
+template struct fg_elem;
+
+template <> struct fg_elem {
+ using store = sycl::half;
+ using pair = sycl::half2;
+ static constexpr mx::matrix_type mtype = mx::matrix_type::fp16;
+ static constexpr int mode = GGML_SYCL_DYNAMIC_PRECISION_F16;
+ static constexpr ggml_type src = GGML_TYPE_F16;
+ static constexpr int mant = 10;
+ static constexpr int slow = 0;
+ static store cvt(float x) { return (store) x; }
+ static pair make(float x, float y) { return pair((store) x, (store) y); }
+};
+
+// tf32 rounds to nearest even with plain bit ops: round_to_tf32 needs a SPIR-V extension that the
+// DG2 AOT target rejects
+static inline uint32_t fg_round_bits(float x, int drop) {
+ const uint32_t u = sycl::bit_cast(x);
+ if ((u & 0x7f800000u) == 0x7f800000u) {
+ return (u & 0x7fffffu) ? u | (1u << drop) : u; // nan stays nan
+ }
+ return u + ((1u << (drop - 1)) - 1) + ((u >> drop) & 1);
+}
+
+struct alignas(4) fg_bf16x2 {
+ sycl::ext::oneapi::bfloat16 x, y;
+};
+
+template <> struct fg_elem {
+ using store = sycl::ext::oneapi::bfloat16;
+ using pair = fg_bf16x2;
+ static constexpr mx::matrix_type mtype = mx::matrix_type::bf16;
+ static constexpr int mode = GGML_SYCL_DYNAMIC_PRECISION_BF16;
+ static constexpr ggml_type src = GGML_TYPE_BF16;
+ static constexpr int mant = 7;
+ static constexpr int slow = 0;
+ static store cvt(float x) { return store(x); }
+ static pair make(float x, float y) { return { cvt(x), cvt(y) }; }
+};
+
+// tf32 keeps f32 range and f16 mantissa, in f32 storage
+template <> struct fg_elem {
+ using store = float;
+ using pair = sycl::float2;
+ static constexpr mx::matrix_type mtype = mx::matrix_type::tf32;
+ static constexpr int mode = GGML_SYCL_DYNAMIC_PRECISION_TF32;
+ static constexpr ggml_type src = GGML_TYPE_COUNT;
+ static constexpr int mant = 10;
+ static constexpr int slow = 1;
+ static store cvt(float x) { return sycl::bit_cast(fg_round_bits(x, 13) & ~0x1fffu); }
+ static pair make(float x, float y) { return pair(cvt(x), cvt(y)); }
+};
+
+// One joint_matrix combination (A type, B type, TM x TN x TK, sub-group size; C and D are f32) and
+// the tiling built on it. A sub-group owns SG_ROWS rows of A (at least 16) and BN (at least 32)
+// columns of B. A and B may differ: the device lists the pairs it supports.
+template struct fg_combo {
+ using ta = TA;
+ using tb = TB;
+ using EA = fg_elem;
+ using EB = fg_elem;
+ using tsa = typename EA::store;
+ using tsb = typename EB::store;
+ static constexpr int TM = TM_;
+ static constexpr int TN = TN_;
+ static constexpr int TK = TK_;
+ static constexpr int SG = SG_;
+ static constexpr int VNNI = 4 / sizeof(tsb); // K rows of B packed in one 32-bit word
+ static constexpr int SG_ROWS = TM > 16 ? TM : 16;
+ static constexpr int RPL = SG_ROWS / SG; // A rows one lane decodes per k step
+ static constexpr int MT = SG_ROWS / TM;
+ static constexpr int BN = TN > 32 ? TN : 32;
+ static constexpr int NT = BN / TN;
+ static constexpr int WG_SIZE = FG_KSPLIT * SG;
+ static constexpr mx::layout b_layout = VNNI == 1 ? mx::layout::row_major : mx::layout::ext_intel_packed;
+ // a 64-wide N is mostly padding here and a 32x64 f32 accumulator needs 128 registers per lane,
+ // so it spills: 13x slower on B60
+ static constexpr bool efficient = TN <= 32;
+ static_assert(SG_ROWS % SG == 0 && SG_ROWS % TM == 0 && BN % TN == 0 && FG_BK % TK == 0, "bad tile");
+ static_assert(BN <= GGML_SYCL_FG_MAX_N, "header gate must cover the tile width");
+};
+
+using fg_half = sycl::half;
+using fg_bf16 = sycl::ext::oneapi::bfloat16;
+using fg_tf32 = mx::precision::tf32;
+
+// One bit of GGML_SYCL_XMX_GATHER_SHAPES per combination. Only combinations some device lists in
+// matrix_combinations are built (appendix of sycl_ext_oneapi_matrix and the runtime's own list).
+template static void fg_visit_combo(int idx, F && f);
+static constexpr int FG_N_COMBOS = 8;
+
+// A spir64_gen AOT build (GGML_SYCL_XMX_AOT_SG) drops the combinations of the other sub-group size
+// entirely: ocloc rejects even an empty kernel that asks for a sub-group size it lacks.
+template static constexpr bool fg_listed() {
+#if defined(GGML_SYCL_XMX_AOT_SG)
+ return SG == GGML_SYCL_XMX_AOT_SG;
+#else
+ return true;
+#endif
+}
+
+template static void fg_call_combo(F && f) {
+ if constexpr (fg_listed()) {
+ f(S{});
+ }
+}
+
+template static void fg_visit_combo(int idx, F && f) {
+ switch (idx) {
+ case 0: fg_call_combo>(f); break; // Xe2, Xe3, Xe-HPC
+ case 1: fg_call_combo>(f); break; // Xe2, Xe3, Xe-HPC
+ case 2: fg_call_combo>(f); break; // Xe2, Xe3, Xe-HPC
+ case 3: fg_call_combo>(f); break; // Xe2, Xe3, Xe-HPC
+ case 4: fg_call_combo>(f); break; // Xe-HPG (Arc A), ARL-H
+ case 5: fg_call_combo>(f); break; // Xe2, Xe3, Xe-HPC
+ case 6: fg_call_combo>(f); break; // Xe2, Xe3, Xe-HPC
+ case 7: fg_call_combo>(f); break; // Xe-HPG (Arc A), ARL-H
+ default: GGML_ABORT("bad XMX combination %d", idx);
+ }
+}
+
+// AOT with -fsycl-targets=intel_gpu_*: compile each tile body only for targets with its sub-group
+// size, since IGC fails on the other ones. A JIT build keeps them all, but each combination lands in
+// its own device image (joint_matrix is an optional kernel feature) and only a combination the
+// device reports is launched, so the runtime never asks IGC for the others.
+#if defined(__SYCL_DEVICE_ONLY__)
+# if __SYCL_TARGET_INTEL_GPU_ACM_G10__ || __SYCL_TARGET_INTEL_GPU_ACM_G11__ || __SYCL_TARGET_INTEL_GPU_ACM_G12__ || \
+ __SYCL_TARGET_INTEL_GPU_ARL_H__
+# define FG_AOT_SG 8
+# elif __SYCL_TARGET_INTEL_GPU_PVC__ || __SYCL_TARGET_INTEL_GPU_PVC_VG__ || __SYCL_TARGET_INTEL_GPU_BMG_G21__ || \
+ __SYCL_TARGET_INTEL_GPU_BMG_G31__ || __SYCL_TARGET_INTEL_GPU_LNL_M__ || __SYCL_TARGET_INTEL_GPU_PTL_H__ || \
+ __SYCL_TARGET_INTEL_GPU_PTL_U__ || __SYCL_TARGET_INTEL_GPU_WCL__ || __SYCL_TARGET_INTEL_GPU_NVL_S__ || \
+ __SYCL_TARGET_INTEL_GPU_NVL_U__ || __SYCL_TARGET_INTEL_GPU_NVL_P__
+# define FG_AOT_SG 16
+# elif __SYCL_TARGET_INTEL_GPU_TGLLP__ || __SYCL_TARGET_INTEL_GPU_RKL__ || __SYCL_TARGET_INTEL_GPU_ADL_S__ || \
+ __SYCL_TARGET_INTEL_GPU_ADL_P__ || __SYCL_TARGET_INTEL_GPU_ADL_N__ || __SYCL_TARGET_INTEL_GPU_DG1__ || \
+ __SYCL_TARGET_INTEL_GPU_MTL_U__ || __SYCL_TARGET_INTEL_GPU_MTL_H__
+# define FG_AOT_SG 0 // no XMX
+# endif
+#endif
+
+template static constexpr bool fg_built() {
+#if defined(FG_AOT_SG)
+ return SG == FG_AOT_SG;
+#else
+ return true;
+#endif
+}
+
+// Upper bound on the tile count when total_rows rows are routed to n_as experts: the worst case gives
+// each expert one row and fills whole tiles with the rest. The bound depends only on the shape, not
+// on the routing, so the pool reuses one buffer every ubatch instead of keeping one per size seen.
+static constexpr int64_t grouped_gemm_max_tiles(int64_t total_rows, int64_t n_as, int64_t BN) {
+ return total_rows <= n_as ? total_rows : n_as + (total_rows - n_as) / BN;
+}
+// Tiles do not cross experts, so the bound is not ceil(total_rows / BN): 34 rows over 2 experts with
+// BN = 16 split 17 + 17 need 2 + 2 tiles, where the ceil gives 3.
+static_assert(grouped_gemm_max_tiles(34, 2, 16) == 4);
+
+// the device lists S with an f32 accumulator and output
+template static bool fg_device_has_combo(const std::vector & combinations) {
+ for (const auto & c : combinations) {
+ if (c.atype == S::EA::mtype && c.btype == S::EB::mtype && c.ctype == mx::matrix_type::fp32 &&
+ c.dtype == mx::matrix_type::fp32 &&
+ (c.max_msize >= (size_t) S::TM || c.msize == (size_t) S::TM) &&
+ (c.max_nsize >= (size_t) S::TN || c.nsize == (size_t) S::TN) &&
+ (c.max_ksize >= (size_t) S::TK || c.ksize == (size_t) S::TK)) {
+ return true;
+ }
+ }
+ return false;
+}
+
+template static const char * fg_type_name() {
+ return std::is_same_v ? "f16" : std::is_same_v ? "bf16" : "tf32";
+}
+
+static std::string fg_combo_name(int idx) {
+ std::string name;
+ fg_visit_combo(idx, [&](auto s) {
+ using S = decltype(s);
+ name = std::string(fg_type_name()) + "x" + fg_type_name() + " " +
+ std::to_string(S::TM) + "x" + std::to_string(S::TN) + "x" + std::to_string(S::TK) + " sg" +
+ std::to_string(S::SG);
+ });
+ return name;
+}
+
+// Combinations this build has kernels for and the device lists, one bit each. Cached per device:
+// on a mixed box the first caller's verdict is not the others'.
+static int fg_device_combos(const sycl::device & dev) {
+ static std::mutex mtx;
+ static std::unordered_map known;
+ std::lock_guard lock(mtx);
+ const auto it = known.find(dev);
+ if (it != known.end()) {
+ return it->second;
+ }
+ int available = 0;
+ try {
+ const auto combinations = dev.get_info();
+ const auto sg_sizes = dev.get_info();
+ for (int idx = 0; idx < FG_N_COMBOS; ++idx) {
+ fg_visit_combo(idx, [&](auto s) {
+ using S = decltype(s);
+ const bool sg_ok = std::find(sg_sizes.begin(), sg_sizes.end(), (size_t) S::SG) != sg_sizes.end();
+ if (sg_ok && fg_device_has_combo(combinations)) {
+ available |= 1 << idx;
+ }
+ });
+ }
+ } catch (const sycl::exception &) {
+ available = 0;
+ }
+ GGML_LOG_INFO("%s: %s: XMX dequant-GEMM combinations available 0x%x, allowed 0x%x\n", __func__,
+ dev.get_info().c_str(), available, g_ggml_sycl_xmx_gather_shapes);
+ known.emplace(dev, available);
+ return available;
+}
+
+// Rank of combination S for a src1 of type src1_type, lower is better. Order:
+// 1. throughput: a tile that does not spill, then the fastest type class of A and B
+// 2. B type equal to the src1 type, so the pack is a plain copy
+// 3. B at least as precise as f16
+// 4. the device's native DPAS tile (8 x SG x 32 bytes of K), then the largest M x K
+// A costs nothing to convert: the A stage emits any type at the same cost.
+template static int64_t fg_rank(ggml_type src1_type) {
+ const int64_t spills = !S::efficient;
+ const int64_t slow = std::max(S::EA::slow, S::EB::slow);
+ const int64_t convert = S::EB::src != src1_type;
+ const int64_t lossy = S::EB::mant < 10;
+ const int64_t foreign = !(S::TM == 8 && S::TN == S::SG);
+ const int64_t mk = 1024 - S::TM * S::TK;
+ return ((((spills * 2 + slow) * 2 + convert) * 2 + lossy) * 2 + foreign) * 2048 + mk;
+}
+
+// whether XMX operands of type mode (a GGML_SYCL_DYNAMIC_PRECISION value) meet the src1 request
+// [TAG_GGML_PREC]. f16 lacks the f32 range that BF16 and F32 ask for; an F32 request goes only as far
+// down as GGML_SYCL_DYNAMIC_REQUIRED_PRECISION allows.
+static bool fg_mode_meets(int mode, int32_t src1_prec) {
+ if (src1_prec == GGML_PREC_UNDEFINED || src1_prec >= GGML_PREC_F16) {
+ return true;
+ }
+ if (src1_prec >= GGML_PREC_BF16) {
+ return mode != GGML_SYCL_DYNAMIC_PRECISION_F16;
+ }
+ switch (g_ggml_sycl_dynamic_required_precision) {
+ case GGML_SYCL_DYNAMIC_PRECISION_TF32: return mode == GGML_SYCL_DYNAMIC_PRECISION_TF32;
+ case GGML_SYCL_DYNAMIC_PRECISION_BF16: return mode != GGML_SYCL_DYNAMIC_PRECISION_F16;
+ default: return false;
+ }
+}
+
+// Best allowed combination for this call, or -1 if none. The type is GGML_SYCL_DYNAMIC_PRECISION if it
+// meets the src1 request; if not, bf16 then tf32 for a BF16 request (fastest first) and tf32 then bf16
+// for an F32 request (most mantissa first).
+static int fg_pick_combo(dpct::queue_ptr stream, ggml_type src1_type, int32_t src1_prec) {
+ const sycl::device dev = stream->get_device();
+ const int allowed = fg_device_combos(dev) & g_ggml_sycl_xmx_gather_shapes;
+ const bool f32_req = src1_prec != GGML_PREC_UNDEFINED && src1_prec < GGML_PREC_BF16;
+ const int modes[] = {
+ g_ggml_sycl_dynamic_precision,
+ f32_req ? GGML_SYCL_DYNAMIC_PRECISION_TF32 : GGML_SYCL_DYNAMIC_PRECISION_BF16,
+ f32_req ? GGML_SYCL_DYNAMIC_PRECISION_BF16 : GGML_SYCL_DYNAMIC_PRECISION_TF32,
+ };
+ int best = -1;
+ int64_t best_rank = 0;
+ for (int i = 0; i < 3 && best < 0; ++i) {
+ const int mode = modes[i];
+ if (!fg_mode_meets(mode, src1_prec)) {
+ continue;
+ }
+ for (int idx = 0; idx < FG_N_COMBOS; ++idx) {
+ if (!(allowed & (1 << idx))) {
+ continue;
+ }
+ fg_visit_combo(idx, [&](auto s) {
+ using S = decltype(s);
+ if (S::EA::mode != mode || S::EB::mode != mode) {
+ return;
+ }
+ const int64_t rank = fg_rank(src1_type);
+ if (best < 0 || rank < best_rank) {
+ best = idx;
+ best_rank = rank;
+ }
+ });
+ }
+ }
+ // log each distinct decision once
+ static std::mutex mtx;
+ static std::set> seen;
+ std::lock_guard lock(mtx);
+ if (seen.emplace(std::hash{}(dev), (int) src1_type, src1_prec, best).second) {
+ GGML_LOG_INFO("%s: src1 %s, src1 prec %d -> %s\n", __func__, ggml_type_name(src1_type), src1_prec,
+ best >= 0 ? fg_combo_name(best).c_str() : "none (library GEMM)");
+ }
+ return best;
+}
+
+// src1 [N][K] -> packed [K/V][Npad][V] so B tiles load straight from global memory
+template
+static void fused_gemm_pack_b(const T_src * y, typename E::store * packed, int N, int Npad, int K,
+ dpct::queue_ptr stream) {
+ constexpr int V = 4 / sizeof(typename E::store);
+ const int kqs = K / V;
+ stream->parallel_for(sycl::range<1>((size_t) Npad * kqs), [=](sycl::id<1> id) {
+ const int idx = id[0];
+ const int n = idx / kqs;
+ const int kq = idx - n * kqs;
+ typename E::store vals[V] = {};
+ if (n < N) {
+ const T_src * src = y + (size_t) n * K + V * kq;
+#pragma unroll
+ for (int v = 0; v < V; ++v) {
+ vals[v] = E::cvt((float) src[v]);
+ }
+ }
+ typename E::store * out = packed + ((size_t) kq * Npad + n) * V;
+#pragma unroll
+ for (int v = 0; v < V; ++v) {
+ out[v] = vals[v];
+ }
+ });
+}
+
+// A stage: one lane owns one row and decodes FG_BK values of it per k step, with every scale
+// folded into the value so the mad below sees plain A elements. One overload per weight format.
+template
+static __dpct_inline__ void fg_stage_a(const block_iq4_nl * __restrict__ xrow, const int kb, typename E::pair * a) {
+ const block_iq4_nl blk = xrow[kb];
+ const float d = (float) blk.d;
+#pragma unroll
+ for (int j = 0; j < QK4_NL / 2; j += 2) {
+ const uint8_t q0 = blk.qs[j];
+ const uint8_t q1 = blk.qs[j + 1];
+ a[j / 2] = E::make(d * kvalues_iq4nl[q0 & 0xf], d * kvalues_iq4nl[q1 & 0xf]);
+ a[j / 2 + 8] = E::make(d * kvalues_iq4nl[q0 >> 4], d * kvalues_iq4nl[q1 >> 4]);
+ }
+}
+
+// iq3_s: k step kb is sub-block kb % 8 of superblock kb / 8. The superblock is 110 bytes, so read
+// only the fields of that sub-block instead of copying the block. Same decode as
+// dequantize_block_iq3_s: grid entries are taken as dwords and the sign bit is a plain shift.
+template
+static __dpct_inline__ void fg_stage_a(const block_iq3_s * __restrict__ xrow, const int kb, typename E::pair * a) {
+ static_assert(QK_K == 256, "the iq3_s A stage assumes 8 sub-blocks per superblock");
+ const block_iq3_s * blk = xrow + kb / (QK_K / 32);
+ const int ib8 = kb % (QK_K / 32);
+ const uint8_t * qs = blk->qs + 8 * ib8;
+ const int qh = blk->qh[ib8];
+ const float d = (float) blk->d * (1 + 2 * ((blk->scales[ib8 / 2] >> (4 * (ib8 % 2))) & 0xf));
+#pragma unroll
+ for (int il = 0; il < 4; ++il) {
+ const uint32_t grid1 = iq3s_grid[qs[2 * il + 0] | ((qh << (8 - 2 * il)) & 256)];
+ const uint32_t grid2 = iq3s_grid[qs[2 * il + 1] | ((qh << (7 - 2 * il)) & 256)];
+ const int signs = blk->signs[4 * ib8 + il];
+#pragma unroll
+ for (int j = 0; j < 2; ++j) {
+ const float g1a = (float) ((grid1 >> (16 * j + 0)) & 0xff);
+ const float g1b = (float) ((grid1 >> (16 * j + 8)) & 0xff);
+ const float g2a = (float) ((grid2 >> (16 * j + 0)) & 0xff);
+ const float g2b = (float) ((grid2 >> (16 * j + 8)) & 0xff);
+ const int s = 2 * j;
+ a[4 * il + j] = E::make(d * ((signs & (1 << (s + 0))) ? -g1a : g1a),
+ d * ((signs & (1 << (s + 1))) ? -g1b : g1b));
+ a[4 * il + j + 2] = E::make(d * ((signs & (1 << (s + 4))) ? -g2a : g2a),
+ d * ((signs & (1 << (s + 5))) ? -g2b : g2b));
+ }
+ }
+}
+
+// values per stored block, so a row of K values is K/qk blocks
+template struct fg_block_traits;
+template <> struct fg_block_traits { static constexpr int qk = QK4_NL; };
+template <> struct fg_block_traits { static constexpr int qk = QK_K; };
+template <> struct fg_block_traits { static constexpr int qk = QK_K; };
+template <> struct fg_block_traits { static constexpr int qk = QK_K; };
+template <> struct fg_block_traits { static constexpr int qk = QK_K; };
+template <> struct fg_block_traits { static constexpr int qk = QK_K; };
+template <> struct fg_block_traits { static constexpr int qk = QK_K; };
+template <> struct fg_block_traits { static constexpr int qk = QK_K; };
+template <> struct fg_block_traits { static constexpr int qk = QK_K; };
+template <> struct fg_block_traits { static constexpr int qk = QK8_0; };
+template <> struct fg_block_traits { static constexpr int qk = QK_K; };
+template <> struct fg_block_traits { static constexpr int qk = QK_K; };
+template <> struct fg_block_traits { static constexpr int qk = QK_K; };
+
+// The A stages below are the dequantize_block_iq* kernels rewritten for one k step. There a
+// work-item handled one quarter (il) of one 32-wide sub-block (ib); here one lane produces the
+// whole step, so il becomes a loop and ib is kb inside the superblock. Each quarter yields 8
+// consecutive values, i.e. 4 pairs at a[4*il], so nothing larger than 8 floats is ever live.
+#define FG_SUPERBLOCK(T) \
+ static_assert(QK_K == 256, "the " #T " A stage assumes 8 sub-blocks per superblock"); \
+ const T * blk = xrow + kb / (QK_K / 32); \
+ const int ib = kb % (QK_K / 32)
+
+template
+static __dpct_inline__ void fg_pack_quarter(const float * __restrict__ t, typename E::pair * a, int il) {
+#pragma unroll
+ for (int j = 0; j < 4; ++j) {
+ a[4 * il + j] = E::make(t[2 * j], t[2 * j + 1]);
+ }
+}
+
+template
+static __dpct_inline__ void fg_stage_a(const block_iq4_xs * __restrict__ xrow, const int kb, typename E::pair * a) {
+ FG_SUPERBLOCK(block_iq4_xs);
+ // low nibbles fill the first half of the step, high nibbles the second, so the two halves
+ // land at a[0..7] and a[8..15] and no quarter loop is needed
+ const float d = (float) blk->d *
+ ((((blk->scales_l[ib / 2] >> (4 * (ib % 2))) & 0xf) | (((blk->scales_h >> (2 * ib)) & 3) << 4)) - 32);
+ const uint8_t * q4 = blk->qs + 16 * ib;
+#pragma unroll
+ for (int j = 0; j < 8; ++j) {
+ a[j] = E::make(d * kvalues_iq4nl[q4[2 * j] & 0xf], d * kvalues_iq4nl[q4[2 * j + 1] & 0xf]);
+ a[8 + j] = E::make(d * kvalues_iq4nl[q4[2 * j] >> 4], d * kvalues_iq4nl[q4[2 * j + 1] >> 4]);
+ }
+}
+
+template
+static __dpct_inline__ void fg_stage_a(const block_iq3_xxs * __restrict__ xrow, const int kb, typename E::pair * a) {
+ FG_SUPERBLOCK(block_iq3_xxs);
+ const uint8_t * q3 = blk->qs + 8 * ib;
+ const uint16_t * gas = (const uint16_t *) (blk->qs + QK_K / 4) + 2 * ib;
+ const uint32_t aux32 = gas[0] | (gas[1] << 16);
+ const float d = (float) blk->d * (0.5f + (aux32 >> 28)) * 0.5f;
+#pragma unroll
+ for (int il = 0; il < 4; ++il) {
+ const uint8_t * grid1 = (const uint8_t *) (iq3xxs_grid + q3[2 * il + 0]);
+ const uint8_t * grid2 = (const uint8_t *) (iq3xxs_grid + q3[2 * il + 1]);
+ const uint8_t signs = ksigns_iq2xs[(aux32 >> (7 * il)) & 127];
+ float t[8];
+#pragma unroll
+ for (int j = 0; j < 4; ++j) {
+ t[j + 0] = d * grid1[j] * (signs & kmask_iq2xs[j + 0] ? -1.f : 1.f);
+ t[j + 4] = d * grid2[j] * (signs & kmask_iq2xs[j + 4] ? -1.f : 1.f);
+ }
+ fg_pack_quarter(t, a, il);
+ }
+}
+
+template
+static __dpct_inline__ void fg_stage_a(const block_iq2_xxs * __restrict__ xrow, const int kb, typename E::pair * a) {
+ FG_SUPERBLOCK(block_iq2_xxs);
+ const uint16_t * q2 = blk->qs + 4 * ib;
+ const uint8_t * aux8 = (const uint8_t *) q2;
+ const uint32_t aux32 = q2[2] | (q2[3] << 16);
+ const float d = (float) blk->d * (0.5f + (aux32 >> 28)) * 0.25f;
+#pragma unroll
+ for (int il = 0; il < 4; ++il) {
+ const uint8_t * grid = (const uint8_t *) (iq2xxs_grid + aux8[il]);
+ const uint8_t signs = ksigns_iq2xs[(aux32 >> (7 * il)) & 127];
+ float t[8];
+#pragma unroll
+ for (int j = 0; j < 8; ++j) {
+ t[j] = d * grid[j] * (signs & kmask_iq2xs[j] ? -1.f : 1.f);
+ }
+ fg_pack_quarter(t, a, il);
+ }
+}
+
+template
+static __dpct_inline__ void fg_stage_a(const block_iq2_xs * __restrict__ xrow, const int kb, typename E::pair * a) {
+ FG_SUPERBLOCK(block_iq2_xs);
+ const uint16_t * q2 = blk->qs + 4 * ib;
+#pragma unroll
+ for (int il = 0; il < 4; ++il) {
+ const uint8_t * grid = (const uint8_t *) (iq2xs_grid + (q2[il] & 511));
+ const float d = (float) blk->d * (0.5f + ((blk->scales[ib] >> (4 * (il / 2))) & 0xf)) * 0.25f;
+ const uint8_t signs = ksigns_iq2xs[q2[il] >> 9];
+ float t[8];
+#pragma unroll
+ for (int j = 0; j < 8; ++j) {
+ t[j] = d * grid[j] * (signs & kmask_iq2xs[j] ? -1.f : 1.f);
+ }
+ fg_pack_quarter(t, a, il);
+ }
+}
+
+template
+static __dpct_inline__ void fg_stage_a(const block_iq2_s * __restrict__ xrow, const int kb, typename E::pair * a) {
+ FG_SUPERBLOCK(block_iq2_s);
+#pragma unroll
+ for (int il = 0; il < 4; ++il) {
+ const uint8_t * grid =
+ (const uint8_t *) (iq2s_grid + (blk->qs[4 * ib + il] | ((blk->qh[ib] << (8 - 2 * il)) & 0x300)));
+ const float d = (float) blk->d * (0.5f + ((blk->scales[ib] >> (4 * (il / 2))) & 0xf)) * 0.25f;
+ const uint8_t signs = blk->qs[QK_K / 8 + 4 * ib + il];
+ float t[8];
+#pragma unroll
+ for (int j = 0; j < 8; ++j) {
+ t[j] = d * grid[j] * (signs & kmask_iq2xs[j] ? -1.f : 1.f);
+ }
+ fg_pack_quarter(t, a, il);
+ }
+}
+
+template
+static __dpct_inline__ void fg_stage_a(const block_iq1_s * __restrict__ xrow, const int kb, typename E::pair * a) {
+ FG_SUPERBLOCK(block_iq1_s);
+ const float delta = blk->qh[ib] & 0x8000 ? -1 - IQ1S_DELTA : -1 + IQ1S_DELTA;
+ const float d = (float) blk->d * (2 * ((blk->qh[ib] >> 12) & 7) + 1);
+#pragma unroll
+ for (int il = 0; il < 4; ++il) {
+ uint32_t grid32[2];
+ const int8_t * q = (const int8_t *) grid32;
+ grid32[0] = iq1s_grid_gpu[blk->qs[4 * ib + il] | (((blk->qh[ib] >> (3 * il)) & 7) << 8)];
+ grid32[1] = (grid32[0] >> 4) & 0x0f0f0f0f;
+ grid32[0] &= 0x0f0f0f0f;
+ float t[8];
+#pragma unroll
+ for (int j = 0; j < 8; ++j) {
+ t[j] = d * (q[j] + delta);
+ }
+ fg_pack_quarter(t, a, il);
+ }
+}
+
+template
+static __dpct_inline__ void fg_stage_a(const block_iq1_m * __restrict__ xrow, const int kb, typename E::pair * a) {
+ FG_SUPERBLOCK(block_iq1_m);
+ const uint16_t * sc = (const uint16_t *) blk->scales;
+ iq1m_scale_t scale;
+ scale.u16 = (sc[0] >> 12) | ((sc[1] >> 8) & 0x00f0) | ((sc[2] >> 4) & 0x0f00) | (sc[3] & 0xf000);
+#pragma unroll
+ for (int il = 0; il < 4; ++il) {
+ const int ib16 = 2 * ib + il / 2;
+ const float d = (float) scale.f16 * (2 * ((sc[ib16 / 4] >> (3 * (ib16 % 4))) & 0x7) + 1);
+ const float delta = blk->qh[2 * ib + il / 2] & (0x08 << (4 * (il % 2))) ? -1 - IQ1M_DELTA : -1 + IQ1M_DELTA;
+ uint32_t grid32[2];
+ const int8_t * q = (const int8_t *) grid32;
+ grid32[0] = iq1s_grid_gpu[blk->qs[4 * ib + il] |
+ (((blk->qh[2 * ib + il / 2] >> (4 * (il % 2))) & 7) << 8)];
+ grid32[1] = (grid32[0] >> 4) & 0x0f0f0f0f;
+ grid32[0] &= 0x0f0f0f0f;
+ float t[8];
+#pragma unroll
+ for (int j = 0; j < 8; ++j) {
+ t[j] = d * (q[j] + delta);
+ }
+ fg_pack_quarter(t, a, il);
+ }
+}
+
+
+// q8_0 and the k-quants. Each decode takes the fields of one 32-value sub-block, so the canonical
+// layout and the reorder (SoA) layout of reorder_qw() share it and differ only in where the fields
+// live, as in dequantize.hpp.
+template
+static __dpct_inline__ void fg_decode_q8_0(const int8_t * __restrict__ qs, const float d, typename E::pair * a) {
+#pragma unroll
+ for (int j = 0; j < QK8_0 / 2; ++j) {
+ a[j] = E::make(d * qs[2 * j], d * qs[2 * j + 1]);
+ }
+}
+
+// same unpack as get_scale_min_k4() in dequantize.hpp
+static __dpct_inline__ void fg_scale_min_k4(const int j, const uint8_t * __restrict__ q, uint8_t & d, uint8_t & m) {
+ if (j < 4) {
+ d = q[j] & 63;
+ m = q[j + 4] & 63;
+ } else {
+ d = (q[j + 4] & 0xF) | ((q[j - 4] >> 6) << 4);
+ m = (q[j + 4] >> 4) | ((q[j - 0] >> 6) << 4);
+ }
+}
+
+// q4_K sub-block ib (0..7) uses scale/min pair ib and the low (even ib) or high (odd ib) nibbles of
+// qs[32 * (ib / 2) ...], as in dequantize_row_q4_K
+template
+static __dpct_inline__ void fg_decode_q4_K(const uint8_t * __restrict__ qs, const uint8_t * __restrict__ scales,
+ const sycl::half2 dm, const int ib, typename E::pair * a) {
+ uint8_t sc, mb;
+ fg_scale_min_k4(ib, scales, sc, mb);
+ const float d = (float) dm[0] * sc;
+ const float m = (float) dm[1] * mb;
+ const uint8_t * q = qs + 32 * (ib / 2);
+ const int shift = 4 * (ib % 2);
+#pragma unroll
+ for (int j = 0; j < 16; ++j) {
+ a[j] = E::make(d * ((q[2 * j] >> shift) & 0xF) - m, d * ((q[2 * j + 1] >> shift) & 0xF) - m);
+ }
+}
+
+// q5_K: q4_K plus one high bit per value, bit ib of qh[l]
+template
+static __dpct_inline__ void fg_decode_q5_K(const uint8_t * __restrict__ qs, const uint8_t * __restrict__ qh,
+ const uint8_t * __restrict__ scales, const sycl::half2 dm, const int ib,
+ typename E::pair * a) {
+ uint8_t sc, mb;
+ fg_scale_min_k4(ib, scales, sc, mb);
+ const float d = (float) dm[0] * sc;
+ const float m = (float) dm[1] * mb;
+ const uint8_t * q = qs + 32 * (ib / 2);
+ const int shift = 4 * (ib % 2);
+#pragma unroll
+ for (int j = 0; j < 16; ++j) {
+ const int l = 2 * j;
+ a[j] = E::make(d * (((q[l] >> shift) & 0xF) | (((qh[l] >> ib) & 1) << 4)) - m,
+ d * (((q[l + 1] >> shift) & 0xF) | (((qh[l + 1] >> ib) & 1) << 4)) - m);
+ }
+}
+
+// q6_K: sub-block ib is quarter r = ib % 4 of half h = ib / 4. The half selects ql + 64h, qh + 32h
+// and scales + 8h; the quarter selects ql + 32(r & 1), the ql nibble r / 2, the qh bit pair r and
+// scales + 2r, as in dequantize_row_q6_K. The scale changes at value 16 of the sub-block.
+template
+static __dpct_inline__ void fg_decode_q6_K(const uint8_t * __restrict__ ql, const uint8_t * __restrict__ qh,
+ const int8_t * __restrict__ scales, const float d, const int ib,
+ typename E::pair * a) {
+ const int h = ib / 4;
+ const int r = ib % 4;
+ const uint8_t * q = ql + 64 * h + 32 * (r & 1);
+ const uint8_t * hb = qh + 32 * h;
+ const int8_t * sc = scales + 8 * h + 2 * r;
+#pragma unroll
+ for (int j = 0; j < 16; ++j) {
+ const int l = 2 * j;
+ const float dl = d * sc[j / 8];
+ const int q0 = (((q[l] >> (4 * (r / 2))) & 0xF) | (((hb[l] >> (2 * r)) & 3) << 4)) - 32;
+ const int q1 = (((q[l + 1] >> (4 * (r / 2))) & 0xF) | (((hb[l + 1] >> (2 * r)) & 3) << 4)) - 32;
+ a[j] = E::make(dl * q0, dl * q1);
+ }
+}
+
+template
+static __dpct_inline__ void fg_stage_a(const block_q8_0 * __restrict__ xrow, const int kb, typename E::pair * a) {
+ fg_decode_q8_0(xrow[kb].qs, (float) xrow[kb].d, a);
+}
+
+template
+static __dpct_inline__ void fg_stage_a(const block_q4_K * __restrict__ xrow, const int kb, typename E::pair * a) {
+ FG_SUPERBLOCK(block_q4_K);
+ fg_decode_q4_K(blk->qs, blk->scales, blk->dm, ib, a);
+}
+
+template
+static __dpct_inline__ void fg_stage_a(const block_q5_K * __restrict__ xrow, const int kb, typename E::pair * a) {
+ FG_SUPERBLOCK(block_q5_K);
+ fg_decode_q5_K(blk->qs, blk->qh, blk->scales, blk->dm, ib, a);
+}
+
+template
+static __dpct_inline__ void fg_stage_a(const block_q6_K * __restrict__ xrow, const int kb, typename E::pair * a) {
+ FG_SUPERBLOCK(block_q6_K);
+ fg_decode_q6_K(blk->ql, blk->qh, blk->scales, (float) blk->d, ib, a);
+}
+
+// Reorder (SoA) layout: each block field is one stream over the nblocks of the matrix (of the expert
+// slice for MUL_MAT_ID), in the order reorder_qw() writes them. Block ib holds k step kb.
+template struct fg_soa {
+ static constexpr bool supported = false;
+};
+
+template <> struct fg_soa {
+ static constexpr bool supported = true;
+ // [qs][d]
+ template
+ static __dpct_inline__ void stage(const uint8_t * x, const size_t nblocks, const size_t ib, const int,
+ typename E::pair * a) {
+ const float d = (float) ((const sycl::half *) (x + nblocks * QK8_0))[ib];
+ fg_decode_q8_0((const int8_t *) x + ib * QK8_0, d, a);
+ }
+};
+
+template <> struct fg_soa {
+ static constexpr bool supported = true;
+ // [qs][scales][dm]
+ template
+ static __dpct_inline__ void stage(const uint8_t * x, const size_t nblocks, const size_t ib, const int kb,
+ typename E::pair * a) {
+ const uint8_t * scales = x + nblocks * (QK_K / 2);
+ const sycl::half2 dm = ((const sycl::half2 *) (scales + nblocks * K_SCALE_SIZE))[ib];
+ fg_decode_q4_K(x + ib * (QK_K / 2), scales + ib * K_SCALE_SIZE, dm, kb % (QK_K / 32), a);
+ }
+};
+
+template <> struct fg_soa {
+ static constexpr bool supported = true;
+ // [qs][qh][scales][dm]
+ template
+ static __dpct_inline__ void stage(const uint8_t * x, const size_t nblocks, const size_t ib, const int kb,
+ typename E::pair * a) {
+ const uint8_t * qh = x + nblocks * (QK_K / 2);
+ const uint8_t * scales = qh + nblocks * (QK_K / 8);
+ const sycl::half2 dm = ((const sycl::half2 *) (scales + nblocks * K_SCALE_SIZE))[ib];
+ fg_decode_q5_K(x + ib * (QK_K / 2), qh + ib * (QK_K / 8), scales + ib * K_SCALE_SIZE, dm,
+ kb % (QK_K / 32), a);
+ }
+};
+
+template <> struct fg_soa {
+ static constexpr bool supported = true;
+ // [ql][qh][scales][d]
+ template
+ static __dpct_inline__ void stage(const uint8_t * x, const size_t nblocks, const size_t ib, const int kb,
+ typename E::pair * a) {
+ const uint8_t * qh = x + nblocks * (QK_K / 2);
+ const uint8_t * scales = qh + nblocks * (QK_K / 4);
+ const float d = (float) ((const sycl::half *) (scales + nblocks * (QK_K / 16)))[ib];
+ fg_decode_q6_K(x + ib * (QK_K / 2), qh + ib * (QK_K / 4), (const int8_t *) scales + ib * (QK_K / 16), d,
+ kb % (QK_K / 32), a);
+ }
+};
+
+
+// one SG_ROWS x BN output tile: B columns [b0, b0 + BN) of packed_b go to dst columns [n0, n1),
+// n1 - n0 <= BN
+template
+static void fused_dequant_gemm_tile(
+ const block_q_t * __restrict__ x,
+ const typename S::tsb * __restrict__ packed_b,
+ float * __restrict__ dst,
+ const int M, const int Npad, const int K, const int ldd,
+ const int b0, const int n0, const int n1,
+ sycl::local_accessor tile_a,
+ sycl::local_accessor tile_c,
+ const sycl::nd_item<2> & item) {
+ if constexpr (fg_built()) {
+ using EA = typename S::EA;
+ using TA = typename S::ta;
+ using TB = typename S::tb;
+ const auto sg = item.get_sub_group();
+ const int sg_id = sg.get_group_id()[0];
+ const int lane = sg.get_local_id()[0];
+ const int m0 = item.get_group(1) * S::SG_ROWS;
+ const int nstep = K / FG_BK;
+ const int a_base = sg_id * S::SG_ROWS * FG_BK;
+ const int c_base = sg_id * S::SG_ROWS * S::BN;
+
+ mx::joint_matrix acc[S::MT][S::NT];
+#pragma unroll
+ for (int mt = 0; mt < S::MT; ++mt) {
+#pragma unroll
+ for (int nt = 0; nt < S::NT; ++nt) {
+ mx::joint_matrix_fill(sg, acc[mt][nt], 0.0f);
+ }
+ }
+
+ // lane decodes rows lane, lane + SG, ... of the sub-group's SG_ROWS
+ constexpr int KPB = fg_block_traits::qk / FG_BK; // k steps per block
+ const size_t bpr = K / fg_block_traits::qk;
+ const size_t nblocks = (size_t) M * bpr;
+ size_t row_blk[S::RPL];
+ bool row_ok[S::RPL];
+ typename EA::pair * a[S::RPL];
+#pragma unroll
+ for (int r = 0; r < S::RPL; ++r) {
+ const int row = m0 + r * S::SG + lane;
+ row_ok[r] = row < M;
+ row_blk[r] = (size_t) (row_ok[r] ? row : 0) * bpr;
+ a[r] = (typename EA::pair *) &tile_a[a_base + (r * S::SG + lane) * FG_BK];
+ }
+
+ const auto b_ptr = sycl::address_space_cast(packed_b);
+ const int b_stride = Npad * S::VNNI;
+
+ const int kb_begin = (sg_id * nstep) / FG_KSPLIT;
+ const int kb_end = ((sg_id + 1) * nstep) / FG_KSPLIT;
+ for (int kb = kb_begin; kb < kb_end; ++kb) {
+#pragma unroll
+ for (int r = 0; r < S::RPL; ++r) {
+ if (row_ok[r]) {
+ if constexpr (reordered) {
+ fg_soa::template stage((const uint8_t *) x, nblocks, row_blk[r] + kb / KPB, kb,
+ a[r]);
+ } else {
+ fg_stage_a(x + row_blk[r], kb, a[r]);
+ }
+ } else {
+#pragma unroll
+ for (int j = 0; j < FG_BK / 2; ++j) {
+ a[r][j] = EA::make(0.0f, 0.0f);
+ }
+ }
+ }
+ sycl::group_barrier(sg);
+
+#pragma unroll
+ for (int kt = 0; kt < FG_BK / S::TK; ++kt) {
+ const int kq0 = (kb * FG_BK + kt * S::TK) / S::VNNI;
+ mx::joint_matrix sub_b[S::NT];
+#pragma unroll
+ for (int nt = 0; nt < S::NT; ++nt) {
+ mx::joint_matrix_load(sg, sub_b[nt], b_ptr + (size_t) kq0 * b_stride + (b0 + nt * S::TN) * S::VNNI, b_stride);
+ }
+#pragma unroll
+ for (int mt = 0; mt < S::MT; ++mt) {
+ mx::joint_matrix sub_a;
+ mx::joint_matrix_load(sg, sub_a,
+ tile_a.template get_multi_ptr() + a_base + (mt * S::TM) * FG_BK + kt * S::TK,
+ FG_BK);
+#pragma unroll
+ for (int nt = 0; nt < S::NT; ++nt) {
+ mx::joint_matrix_mad(sg, acc[mt][nt], sub_a, sub_b[nt], acc[mt][nt]);
+ }
+ }
+ }
+ // the next step overwrites tile_a
+ sycl::group_barrier(sg);
+ }
+
+#pragma unroll
+ for (int mt = 0; mt < S::MT; ++mt) {
+#pragma unroll
+ for (int nt = 0; nt < S::NT; ++nt) {
+ mx::joint_matrix_store(sg, acc[mt][nt],
+ tile_c.template get_multi_ptr() + c_base + (mt * S::TM) * S::BN + nt * S::TN,
+ S::BN, mx::layout::row_major);
+ }
+ }
+ sycl::group_barrier(item.get_group());
+
+ // sum the K splits; consecutive lanes write consecutive rows of one dst column
+ for (int idx = item.get_local_linear_id(); idx < S::SG_ROWS * S::BN; idx += S::WG_SIZE) {
+ const int r = idx % S::SG_ROWS;
+ const int c = idx / S::SG_ROWS;
+ const int m = m0 + r;
+ const int n = n0 + c;
+ if (m < M && n < n1) {
+ float sum = 0.0f;
+#pragma unroll
+ for (int s = 0; s < FG_KSPLIT; ++s) {
+ sum += tile_c[s * S::SG_ROWS * S::BN + r * S::BN + c];
+ }
+ dst[(size_t) n * ldd + m] = sum;
+ }
+ }
+ }
+}
+
+template
+static void fused_dequant_gemm_launch(const void * src0, const typename S::tsb * packed, float * dst, const int M,
+ const int N, const int Npad, const int K, const int ldd,
+ const int64_t groups_n, const int64_t groups_m, dpct::queue_ptr stream) {
+ stream->submit([&](sycl::handler & cgh) {
+ sycl::local_accessor tile_a(FG_KSPLIT * S::SG_ROWS * FG_BK, cgh);
+ sycl::local_accessor tile_c(FG_KSPLIT * S::SG_ROWS * S::BN, cgh);
+ cgh.parallel_for(
+ sycl::nd_range<2>(sycl::range<2>(groups_n, groups_m * S::WG_SIZE), sycl::range<2>(1, S::WG_SIZE)),
+ [=](sycl::nd_item<2> item) [[sycl::reqd_sub_group_size(S::SG)]] {
+ const int n0 = item.get_group(0) * S::BN;
+ fused_dequant_gemm_tile((const block_q_t *) src0, packed, dst, M, Npad, K, ldd,
+ n0, n0, N, tile_a, tile_c, item);
+ });
+ });
+}
+
+// grouped: work-group (t, mt) is tile t of the schedule; its B columns sit at t * BN
+template
+static void grouped_dequant_gemm_launch(const char * src0_dd, const size_t expert_stride,
+ const ggml_sycl_gg_tile * tiles_ptr, const typename S::tsb * packed, float * dst,
+ const int M, const int Npad, const int K, const int64_t n_tiles,
+ const int64_t groups_m, dpct::queue_ptr stream) {
+ stream->submit([&](sycl::handler & cgh) {
+ sycl::local_accessor tile_a(FG_KSPLIT * S::SG_ROWS * FG_BK, cgh);
+ sycl::local_accessor tile_c(FG_KSPLIT * S::SG_ROWS * S::BN, cgh);
+ cgh.parallel_for(
+ sycl::nd_range<2>(sycl::range<2>(n_tiles, groups_m * S::WG_SIZE), sycl::range<2>(1, S::WG_SIZE)),
+ [=](sycl::nd_item<2> item) [[sycl::reqd_sub_group_size(S::SG)]] {
+ const int t = item.get_group(0);
+ const ggml_sycl_gg_tile tile = tiles_ptr[t];
+ const block_q_t * x = (const block_q_t *) (src0_dd + (size_t) tile.expert * expert_stride);
+ fused_dequant_gemm_tile(x, packed, dst, M, Npad, K, M, t * S::BN, tile.n0,
+ tile.n1, tile_a, tile_c, item);
+ });
+ });
+}
+
+// src1 f32 rows -> packed [K/V][n_tiles*BN][V], tile t holds its rows [n0, n1) at columns t*BN..,
+// zero past n1. The column runs fastest so a sub-group writes one contiguous run.
+template
+static void grouped_gemm_pack_b(const float * y, typename S::tsb * packed, const ggml_sycl_gg_tile * tiles, int Npad,
+ int K, dpct::queue_ptr stream) {
+ using E = typename S::EB;
+ constexpr int V = S::VNNI;
+ const int kqs = K / V;
+ stream->parallel_for(sycl::range<1>((size_t) Npad * kqs), [=](sycl::id<1> id) {
+ const size_t idx = id[0];
+ const int kq = idx / Npad;
+ const int n = idx - (size_t) kq * Npad;
+ const ggml_sycl_gg_tile tile = tiles[n / S::BN];
+ const int row = tile.n0 + n % S::BN;
+ // one guarded load run per work-item, as a per-element select costs ~1.5% prefill
+ typename S::tsb vals[V] = {};
+ if (row < tile.n1) {
+ const float * src = y + (size_t) row * K + V * kq;
+#pragma unroll
+ for (int v = 0; v < V; ++v) {
+ vals[v] = E::cvt(src[v]);
+ }
+ }
+ typename S::tsb * out = packed + ((size_t) kq * Npad + n) * V;
+#pragma unroll
+ for (int v = 0; v < V; ++v) {
+ out[v] = vals[v];
+ }
+ });
+}
+
+// q8_0 and the k-quants take only the grouped path. The plain kernel decodes A again for every BN columns
+// of a dense batch, and for these formats that costs more than the one dequantization of the library GEMM.
+template static constexpr bool fg_plain_ok() {
+ return !std::is_same_v && !std::is_same_v && !std::is_same_v &&
+ !std::is_same_v;
+}
+
+template struct fg_tag {
+ using type = T;
+ static constexpr bool reordered = R;
+};
+
+// calls f(fg_tag{}) for the weight format and layout; false if it has no A stage
+template static bool fg_visit_layout(bool reordered, F && f) {
+ if (!reordered) {
+ f(fg_tag{});
+ return true;
+ }
+ if constexpr (fg_soa::supported) {
+ f(fg_tag{});
+ return true;
+ }
+ return false;
+}
+
+template static bool fg_visit_type(ggml_type type, bool reordered, F && f) {
+ switch (type) {
+ case GGML_TYPE_IQ4_NL: return fg_visit_layout(reordered, f);
+ case GGML_TYPE_IQ3_S: return fg_visit_layout(reordered, f);
+ case GGML_TYPE_IQ4_XS: return fg_visit_layout(reordered, f);
+ case GGML_TYPE_IQ3_XXS: return fg_visit_layout(reordered, f);
+ case GGML_TYPE_IQ2_XXS: return fg_visit_layout(reordered, f);
+ case GGML_TYPE_IQ2_XS: return fg_visit_layout(reordered, f);
+ case GGML_TYPE_IQ2_S: return fg_visit_layout(reordered, f);
+ case GGML_TYPE_IQ1_S: return fg_visit_layout(reordered, f);
+ case GGML_TYPE_IQ1_M: return fg_visit_layout(reordered, f);
+ case GGML_TYPE_Q8_0: return fg_visit_layout(reordered, f);
+ case GGML_TYPE_Q4_K: return fg_visit_layout(reordered, f);
+ case GGML_TYPE_Q5_K: return fg_visit_layout(reordered, f);
+ case GGML_TYPE_Q6_K: return fg_visit_layout(reordered, f);
+ default: return false;
+ }
+}
+
+template
+static bool fg_fused_run(ggml_type src0_type, const void * src0, const void * src1, ggml_type src1_type,
+ float * dst, int64_t M, int64_t N, int64_t K, int64_t ldd, ggml_sycl_pool & pool,
+ dpct::queue_ptr stream) {
+ const int64_t groups_n = (N + S::BN - 1) / S::BN;
+ const int64_t groups_m = (M + S::SG_ROWS - 1) / S::SG_ROWS;
+ const int Npad = (int) (groups_n * S::BN);
+
+ // src1 is read in its own type: one pass, converted in registers only if B differs
+ ggml_sycl_pool_alloc packed_b(pool, (size_t) K * Npad);
+ if (src1_type == GGML_TYPE_F16) {
+ fused_gemm_pack_b((const sycl::half *) src1, packed_b.get(), (int) N, Npad, (int) K, stream);
+ } else if (src1_type == GGML_TYPE_BF16) {
+ fused_gemm_pack_b((const fg_bf16 *) src1, packed_b.get(), (int) N, Npad, (int) K, stream);
+ } else {
+ fused_gemm_pack_b((const float *) src1, packed_b.get(), (int) N, Npad, (int) K, stream);
+ }
+
+ const typename S::tsb * packed = packed_b.get();
+ return fg_visit_type(src0_type, false, [&](auto tag) {
+ using block_q_t = typename decltype(tag)::type;
+ if constexpr (fg_plain_ok()) {
+ fused_dequant_gemm_launch(src0, packed, dst, (int) M, (int) N, Npad, (int) K, (int) ldd,
+ groups_n, groups_m, stream);
+ }
+ });
+}
+
+bool ggml_sycl_fused_dequant_gemm(ggml_type src0_type, const void * src0, const void * src1, ggml_type src1_type,
+ int32_t src1_prec, float * dst, int64_t M, int64_t N, int64_t K, int64_t ldd,
+ ggml_sycl_pool & pool, dpct::queue_ptr stream) {
+ // every BN columns dequantize A again, so wide N is left to the library GEMM
+ bool plain_ok = false;
+ fg_visit_type(src0_type, false, [&](auto tag) { plain_ok = fg_plain_ok(); });
+ if (g_ggml_sycl_dynamic_precision == GGML_SYCL_DYNAMIC_PRECISION_F32 ||
+ !ggml_sycl_xmx_gather_type_enabled(src0_type) || !plain_ok) {
+ return false;
+ }
+ if (src1_type != GGML_TYPE_F32 && src1_type != GGML_TYPE_F16 && src1_type != GGML_TYPE_BF16) {
+ return false;
+ }
+ if (!ggml_sycl_fused_dequant_gemm_shape_ok(src0_type, M, N, K, ldd)) {
+ return false;
+ }
+ const int combo = fg_pick_combo(stream, src1_type, src1_prec);
+ if (combo < 0) {
+ return false;
+ }
+ bool launched = false;
+ fg_visit_combo(combo, [&](auto s) {
+ launched = fg_fused_run(src0_type, src0, src1, src1_type, dst, M, N, K, ldd, pool, stream);
+ });
+ return launched;
+}
+
+template
+static bool fg_grouped_run(ggml_type src0_type, bool reordered, const void * src0_base, size_t expert_stride,
+ const float * src1, float * dst, const int64_t * expert_row_offsets, int64_t n_as, int64_t M,
+ int64_t K, std::vector & tiles, ggml_sycl_pool & pool,
+ dpct::queue_ptr stream) {
+ // the host knows every slice, so it lays out the work-groups: no search on the device
+ tiles.clear();
+ for (int64_t e = 0; e < n_as; ++e) {
+ const int64_t end = expert_row_offsets[e + 1];
+ for (int64_t n0 = expert_row_offsets[e]; n0 < end; n0 += S::BN) {
+ tiles.push_back({ (int32_t) e, (int32_t) n0, (int32_t) std::min(n0 + S::BN, end) });
+ }
+ }
+ const int64_t n_tiles = tiles.size();
+ const int64_t groups_m = (M + S::SG_ROWS - 1) / S::SG_ROWS;
+ const int Npad = (int) (n_tiles * S::BN);
+
+ const int64_t max_tiles = grouped_gemm_max_tiles(expert_row_offsets[n_as], n_as, S::BN);
+ GGML_ASSERT(n_tiles <= max_tiles);
+ ggml_sycl_pool_alloc tiles_dev(pool, max_tiles);
+ SYCL_CHECK(CHECK_TRY_ERROR(stream->memcpy(tiles_dev.get(), tiles.data(), n_tiles * sizeof(ggml_sycl_gg_tile))));
+
+ ggml_sycl_pool_alloc packed_b(pool, (size_t) K * max_tiles * S::BN);
+ grouped_gemm_pack_b(src1, packed_b.get(), tiles_dev.get(), Npad, (int) K, stream);
+
+ const typename S::tsb * packed = packed_b.get();
+ const ggml_sycl_gg_tile * tiles_ptr = tiles_dev.get();
+ const char * src0_dd = (const char *) src0_base;
+ return fg_visit_type(src0_type, reordered, [&](auto tag) {
+ using T = decltype(tag);
+ grouped_dequant_gemm_launch(src0_dd, expert_stride, tiles_ptr, packed, dst,
+ (int) M, Npad, (int) K, n_tiles, groups_m,
+ stream);
+ });
+}
+
+bool ggml_sycl_grouped_dequant_gemm(ggml_type src0_type, bool reordered, const void * src0_base, size_t expert_stride,
+ const float * src1, int32_t src1_prec, float * dst,
+ const int64_t * expert_row_offsets, int64_t n_as, int64_t M, int64_t K,
+ int64_t total_rows, std::vector & tiles,
+ ggml_sycl_pool & pool, dpct::queue_ptr stream) {
+ int64_t n_active = 0;
+ for (int64_t e = 0; e < n_as; ++e) {
+ n_active += expert_row_offsets[e + 1] > expert_row_offsets[e];
+ }
+ if (g_ggml_sycl_dynamic_precision == GGML_SYCL_DYNAMIC_PRECISION_F32 ||
+ !ggml_sycl_xmx_gather_type_enabled(src0_type) || !fg_visit_type(src0_type, reordered, [](auto) {})) {
+ return false;
+ }
+ if (!ggml_sycl_grouped_dequant_gemm_shape_ok(src0_type, M, K, total_rows, n_active)) {
+ return false;
+ }
+ const int combo = fg_pick_combo(stream, GGML_TYPE_F32, src1_prec);
+ if (combo < 0) {
+ return false;
+ }
+ bool launched = false;
+ fg_visit_combo(combo, [&](auto s) {
+ launched = fg_grouped_run(src0_type, reordered, src0_base, expert_stride, src1, dst,
+ expert_row_offsets, n_as, M, K, tiles, pool, stream);
+ });
+ return launched;
+}
diff --git a/ggml/src/ggml-sycl/fused-gemm.hpp b/ggml/src/ggml-sycl/fused-gemm.hpp
new file mode 100644
index 000000000..298989204
--- /dev/null
+++ b/ggml/src/ggml-sycl/fused-gemm.hpp
@@ -0,0 +1,88 @@
+#ifndef GGML_SYCL_FUSED_GEMM_HPP
+#define GGML_SYCL_FUSED_GEMM_HPP
+
+#include "common.hpp"
+
+
+// Shape and type gates for the kernels below. Device capability is separate: it needs a queue to ask.
+static constexpr int GGML_SYCL_FG_MAX_N = 64; // widest N taken; each shape covers it in BN-wide tiles
+
+// weight formats the fused A stage decodes; K must cover whole stored blocks
+constexpr bool ggml_sycl_fused_dequant_gemm_type_ok(ggml_type src0_type, int64_t K) {
+ // iq4_nl and q8_0 store 32 values per block; every other format here is a 256-value superblock
+ // that the A stage walks in steps of 32, so K must cover whole superblocks.
+ if (src0_type == GGML_TYPE_IQ4_NL || src0_type == GGML_TYPE_Q8_0) {
+ return K % 32 == 0;
+ }
+ const bool superblock =
+ src0_type == GGML_TYPE_Q4_K ||
+ src0_type == GGML_TYPE_Q5_K ||
+ src0_type == GGML_TYPE_Q6_K ||
+ src0_type == GGML_TYPE_IQ3_S ||
+ src0_type == GGML_TYPE_IQ4_XS ||
+ src0_type == GGML_TYPE_IQ3_XXS ||
+ src0_type == GGML_TYPE_IQ2_XXS ||
+ src0_type == GGML_TYPE_IQ2_XS ||
+ src0_type == GGML_TYPE_IQ2_S ||
+ src0_type == GGML_TYPE_IQ1_S ||
+ src0_type == GGML_TYPE_IQ1_M;
+ return superblock && QK_K == 256 && K % QK_K == 0;
+}
+
+constexpr bool ggml_sycl_fused_dequant_gemm_shape_ok(ggml_type src0_type, int64_t M, int64_t N, int64_t K,
+ int64_t ldd) {
+ return ggml_sycl_fused_dequant_gemm_type_ok(src0_type, K) && M > 0 && N > 0 && K > 0 &&
+ N <= GGML_SYCL_FG_MAX_N &&
+ M <= INT32_MAX && N <= INT32_MAX && K <= INT32_MAX && ldd <= INT32_MAX;
+}
+
+// grouped variant: the per-expert fused kernel is only worth it while each expert is narrow,
+// so wider average slices are left to the per-expert library GEMM loop
+constexpr bool ggml_sycl_grouped_dequant_gemm_shape_ok(ggml_type src0_type, int64_t M, int64_t K,
+ int64_t total_rows, int64_t n_active) {
+ return ggml_sycl_fused_dequant_gemm_shape_ok(src0_type, M, 1, K, M) && total_rows > 0 &&
+ total_rows <= INT32_MAX && total_rows <= n_active * GGML_SYCL_FG_MAX_N;
+}
+
+// Runtime type gate, kept out of the constexpr predicates above so those stay pure.
+inline bool ggml_sycl_xmx_gather_type_enabled(ggml_type src0_type) {
+ switch (src0_type) {
+ case GGML_TYPE_IQ4_NL: return (g_ggml_sycl_xmx_gather_types & GGML_SYCL_XMX_GATHER_IQ4_NL ) != 0;
+ case GGML_TYPE_IQ3_S: return (g_ggml_sycl_xmx_gather_types & GGML_SYCL_XMX_GATHER_IQ3_S ) != 0;
+ case GGML_TYPE_IQ4_XS: return (g_ggml_sycl_xmx_gather_types & GGML_SYCL_XMX_GATHER_IQ4_XS ) != 0;
+ case GGML_TYPE_IQ3_XXS: return (g_ggml_sycl_xmx_gather_types & GGML_SYCL_XMX_GATHER_IQ3_XXS ) != 0;
+ case GGML_TYPE_IQ2_XXS: return (g_ggml_sycl_xmx_gather_types & GGML_SYCL_XMX_GATHER_IQ2_XXS ) != 0;
+ case GGML_TYPE_IQ2_XS: return (g_ggml_sycl_xmx_gather_types & GGML_SYCL_XMX_GATHER_IQ2_XS ) != 0;
+ case GGML_TYPE_IQ2_S: return (g_ggml_sycl_xmx_gather_types & GGML_SYCL_XMX_GATHER_IQ2_S ) != 0;
+ case GGML_TYPE_IQ1_S: return (g_ggml_sycl_xmx_gather_types & GGML_SYCL_XMX_GATHER_IQ1_S ) != 0;
+ case GGML_TYPE_IQ1_M: return (g_ggml_sycl_xmx_gather_types & GGML_SYCL_XMX_GATHER_IQ1_M ) != 0;
+ case GGML_TYPE_Q8_0: return (g_ggml_sycl_xmx_gather_types & GGML_SYCL_XMX_GATHER_Q8_0 ) != 0;
+ case GGML_TYPE_Q4_K: return (g_ggml_sycl_xmx_gather_types & GGML_SYCL_XMX_GATHER_Q4_K ) != 0;
+ case GGML_TYPE_Q5_K: return (g_ggml_sycl_xmx_gather_types & GGML_SYCL_XMX_GATHER_Q5_K ) != 0;
+ case GGML_TYPE_Q6_K: return (g_ggml_sycl_xmx_gather_types & GGML_SYCL_XMX_GATHER_Q6_K ) != 0;
+ default: return false;
+ }
+}
+
+// dst[n*ldd + m] = sum_k dequant(src0)[m*K + k] * src1[n*K + k], src1 is F32, F16 or BF16.
+// The XMX combination is picked per call from the src1 type and its precision request src1_prec
+// (op_params[3], [TAG_GGML_PREC]); the accumulator is f32, which meets any request.
+// q8_0 and the k-quants are not handled here, only in the grouped path below.
+// Returns false when the case is not handled (type, device, precision, or shape).
+bool ggml_sycl_fused_dequant_gemm(ggml_type src0_type, const void * src0, const void * src1, ggml_type src1_type,
+ int32_t src1_prec, float * dst, int64_t M, int64_t N, int64_t K, int64_t ldd,
+ ggml_sycl_pool & pool, dpct::queue_ptr stream);
+
+// One launch for every expert of a MUL_MAT_ID: rows of src1/dst are grouped by expert, expert e
+// owns rows [expert_row_offsets[e], expert_row_offsets[e+1]) and reads its weights at
+// src0_base + e*expert_stride. tiles is host scratch that must stay alive until the queue drains.
+// reordered: each expert slice is in the reorder (SoA) layout of reorder_qw().
+// dst[n*M + m] = sum_k dequant(src0_e)[m*K + k] * src1[n*K + k]
+// Returns false when the case is not handled (type, layout, device, precision, or shape).
+bool ggml_sycl_grouped_dequant_gemm(ggml_type src0_type, bool reordered, const void * src0_base, size_t expert_stride,
+ const float * src1, int32_t src1_prec, float * dst,
+ const int64_t * expert_row_offsets, int64_t n_as, int64_t M, int64_t K,
+ int64_t total_rows, std::vector & tiles,
+ ggml_sycl_pool & pool, dpct::queue_ptr stream);
+
+#endif // GGML_SYCL_FUSED_GEMM_HPP
diff --git a/ggml/src/ggml-sycl/fusion.cpp b/ggml/src/ggml-sycl/fusion.cpp
index d3e995233..79fe13a1d 100644
--- a/ggml/src/ggml-sycl/fusion.cpp
+++ b/ggml/src/ggml-sycl/fusion.cpp
@@ -64,6 +64,12 @@ static bool ggml_sycl_should_fuse_mul_mat_glu(const ggml_tensor * gate, const gg
return true;
}
+// the unary ops the fused unary chains in element_wise.cpp have a functor for
+static bool ggml_sycl_fused_unary_has_kernel(ggml_unary_op unary_op) {
+ return unary_op == GGML_UNARY_OP_SILU || unary_op == GGML_UNARY_OP_SIGMOID ||
+ unary_op == GGML_UNARY_OP_SOFTPLUS;
+}
+
bool ggml_sycl_can_fuse(const ggml_cgraph * cgraph, int node_idx, std::initializer_list ops,
std::initializer_list unary_ops) {
#ifndef NDEBUG
@@ -184,9 +190,7 @@ bool ggml_sycl_can_fuse(const ggml_cgraph * cgraph, int node_idx, std::initializ
return false;
}
- // the ops ggml_sycl_op_unary_mul_fused() has a kernel for
- if (unary_op != GGML_UNARY_OP_SILU && unary_op != GGML_UNARY_OP_SIGMOID &&
- unary_op != GGML_UNARY_OP_SOFTPLUS) {
+ if (!ggml_sycl_fused_unary_has_kernel(unary_op)) {
return false;
}
@@ -233,6 +237,55 @@ bool ggml_sycl_can_fuse(const ggml_cgraph * cgraph, int node_idx, std::initializ
return true;
}
+ // ADD(bias) + UNARY + MUL(scale): the delta-net alpha gate, softplus(alpha + dt) * a.
+ // The broadcast is what stops the same-shape UNARY + MUL branch above firing past one token.
+ if (ops.size() == 3 && ops.begin()[0] == GGML_OP_ADD && ops.begin()[1] == GGML_OP_UNARY &&
+ ops.begin()[2] == GGML_OP_MUL && unary_ops.size() == 1) {
+ const ggml_tensor * add = cgraph->nodes[node_idx];
+ const ggml_tensor * unary = cgraph->nodes[node_idx + 1];
+ const ggml_tensor * mul = cgraph->nodes[node_idx + 2];
+
+ const ggml_unary_op unary_op = ggml_get_unary_op(unary);
+ if (unary_op != unary_ops.begin()[0]) {
+ return false;
+ }
+
+ if (!ggml_sycl_fused_unary_has_kernel(unary_op)) {
+ return false;
+ }
+
+ // ggml_can_fuse() has already pinned the chain: unary consumes add, mul consumes
+ // unary, add and unary have one use each, and all three have the same shape
+ const ggml_tensor * a = add->src[0];
+ const ggml_tensor * bias = add->src[1];
+ const ggml_tensor * scale = (mul->src[0] == unary) ? mul->src[1] : mul->src[0];
+
+ if (a->type != GGML_TYPE_F32 || bias->type != GGML_TYPE_F32 ||
+ scale->type != GGML_TYPE_F32 || mul->type != GGML_TYPE_F32) {
+ return false;
+ }
+
+ // the activation and the destination are indexed flat
+ if (!ggml_is_contiguous(a) || !ggml_is_contiguous(mul) || !ggml_are_same_shape(a, mul)) {
+ return false;
+ }
+
+ // the kernel reads the bias and the scale as v[col], so each must be a single
+ // contiguous row spanning ne0
+ if (bias->ne[0] != a->ne[0] || scale->ne[0] != a->ne[0] ||
+ ggml_nrows(bias) != 1 || ggml_nrows(scale) != 1 ||
+ !ggml_is_contiguous(bias) || !ggml_is_contiguous(scale)) {
+ return false;
+ }
+
+ // the 32-bit fastdiv is inexact past 2^31; decline, the unfused path handles it
+ if (ggml_nelements(mul) >= ((int64_t) 1 << 31)) {
+ return false;
+ }
+
+ return true;
+ }
+
if (ops.size() == 3 && ops.begin()[0] == GGML_OP_SSM_CONV && ops.begin()[1] == GGML_OP_ADD &&
ops.begin()[2] == GGML_OP_UNARY && unary_ops.size() == 1 && unary_ops.begin()[0] == GGML_UNARY_OP_SILU) {
const ggml_tensor * ssm_conv = cgraph->nodes[node_idx];
diff --git a/ggml/src/ggml-sycl/fwht.cpp b/ggml/src/ggml-sycl/fwht.cpp
index fb48d7fec..5af549c23 100644
--- a/ggml/src/ggml-sycl/fwht.cpp
+++ b/ggml/src/ggml-sycl/fwht.cpp
@@ -46,8 +46,8 @@ static constexpr float H20[20][20] = {
#undef P
#undef N
-template
-static void fwht_kernel(const float * __restrict__ src, float * __restrict__ dst, const int64_t n_rows,
+template
+static void fwht_kernel(const T * __restrict__ src, float * __restrict__ dst, const int64_t n_rows,
const float scale, const sycl::nd_item<2> & item) {
const sycl::sub_group sg = item.get_sub_group();
@@ -67,7 +67,7 @@ static void fwht_kernel(const float * __restrict__ src, float * __restrict__ dst
#pragma unroll
for (int i = 0; i < el_w; ++i) {
- reg[i] = src[i * WARP_SIZE + lane] * scale;
+ reg[i] = static_cast(src[i * WARP_SIZE + lane]) * scale;
}
// Butterflies inside the sub-group. The partner of a lane with bit h clear is the
@@ -107,8 +107,8 @@ static void fwht_kernel(const float * __restrict__ src, float * __restrict__ dst
}
}
-template
-static void launch_fwht(const float * src, float * dst, const int64_t n_rows, const float scale,
+template
+static void launch_fwht(const T * src, float * dst, const int64_t n_rows, const float scale,
dpct::queue_ptr stream) {
constexpr int rows_per_block = 4;
@@ -120,7 +120,7 @@ static void launch_fwht(const float * src, float * dst, const int64_t n_rows, co
stream->parallel_for(sycl::nd_range<2>(global, local),
[=](sycl::nd_item<2> item) [[sycl::reqd_sub_group_size(WARP_SIZE)]] {
- fwht_kernel(src, dst, n_rows, scale, item);
+ fwht_kernel(src, dst, n_rows, scale, item);
});
}
@@ -128,8 +128,8 @@ static void launch_fwht(const float * src, float * dst, const int64_t n_rows, co
// keeps N/NT values rather than N/WARP_SIZE. Butterflies below the sub-group width
// still shuffle; those up to NT go through work-group local memory; the rest stay
// in registers.
-template
-static void fwht_kernel_wide(const float * __restrict__ src,
+template
+static void fwht_kernel_wide(const T * __restrict__ src,
float * __restrict__ dst,
const int64_t n_rows,
const float scale,
@@ -151,7 +151,7 @@ static void fwht_kernel_wide(const float * __restrict__ src,
float reg[el_w];
#pragma unroll
for (int i = 0; i < el_w; ++i) {
- reg[i] = src[i * NT + tid] * scale;
+ reg[i] = static_cast(src[i * NT + tid]) * scale;
}
const sycl::sub_group sg = item.get_sub_group();
@@ -207,8 +207,8 @@ static void fwht_kernel_wide(const float * __restrict__ src,
}
}
-template
-static void launch_fwht_wide(const float * src,
+template
+static void launch_fwht_wide(const T * src,
float * dst,
const int64_t n_rows,
const float scale,
@@ -220,13 +220,13 @@ static void launch_fwht_wide(const float * src,
sycl::local_accessor smem(sycl::range<1>(N), cgh);
cgh.parallel_for(sycl::nd_range<2>(global, local),
[=](sycl::nd_item<2> item) [[sycl::reqd_sub_group_size(WARP_SIZE)]] {
- fwht_kernel_wide(src, dst, n_rows, scale, item, get_pointer(smem));
+ fwht_kernel_wide(src, dst, n_rows, scale, item, get_pointer(smem));
});
});
}
-template
-static void kronecker_kernel(const float * __restrict__ src,
+template