diff --git a/README.rst b/README.rst
index ef0a729..dfdb81c 100644
--- a/README.rst
+++ b/README.rst
@@ -1,6 +1,12 @@
 GiMMiK
 ======
-Generator of Matrix Multiplication Kernels - GiMMiK - is a tool for generation of high performance matrix multiplication kernel code for various accelerator platforms. Currently C, CUDA, HIP, ISPC, Metal, and OpenCL are supported.
+Generator of Matrix Multiplication Kernels - GiMMiK - is a tool for generation of high performance matrix multiplication kernel code for various accelerator platforms. Currently C, CUDA, HIP, ISPC, Metal, OpenCL, and SYCL are supported.
+
+The SYCL backend (``platform='sycl'``) emits self-contained launcher functions of
+the form ``sycl::event kname(sycl::queue& q, ...)`` that submit the generated
+kernel to a queue, and can be compiled with any SYCL 2020 compiler (e.g. Intel
+oneAPI DPC++). See ``bench/`` for an OpenCL-vs-SYCL benchmark harness targeting
+Intel GPUs.
 
 What does GiMMiK do?
 --------------------
diff --git a/bench/.gitignore b/bench/.gitignore
new file mode 100644
index 0000000..4ce2f2c
--- /dev/null
+++ b/bench/.gitignore
@@ -0,0 +1,4 @@
+build/
+vbuild/
+*.o
+*.bin
diff --git a/bench/bench_gen.py b/bench/bench_gen.py
new file mode 100644
index 0000000..2b0dd84
--- /dev/null
+++ b/bench/bench_gen.py
@@ -0,0 +1,137 @@
+#!/usr/bin/env python3
+"""Generate OpenCL and SYCL GiMMiK kernels + shared input data for benchmarking.
+
+Produces, under bench/build/:
+  ocl/<name>.cl      - one OpenCL kernel per variant (unique entry name)
+  sycl/<name>.cpp    - one SYCL launcher per variant (unique entry name)
+  B.bin, Cref.bin    - shared fp64 input / reference output (row-major, k x n / m x n)
+  manifest.json      - description of every kernel + problem dims
+"""
+import json
+import os
+import sys
+
+import numpy as np
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+from gimmik import OpenCLMatMul, SYCLMatMul
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+BUILD = os.path.join(HERE, 'build')
+
+
+def make_operator(m, k, sparsity, seed=42):
+    rng = np.random.default_rng(seed)
+    A = rng.standard_normal((m, k))
+    A[rng.random((m, k)) < sparsity] = 0.0
+    # guarantee no fully-zero row so every output column is exercised
+    for j in range(m):
+        if not A[j].any():
+            A[j, rng.integers(k)] = rng.standard_normal()
+    return A
+
+
+def main():
+    m = int(os.environ.get('GMK_M', 32))
+    k = int(os.environ.get('GMK_K', 48))
+    n = int(os.environ.get('GMK_N', 200000))
+    sparsity = float(os.environ.get('GMK_SPARSITY', 0.5))
+
+    ldb = ldc = n
+    A = make_operator(m, k, sparsity)
+    nnz = int((A != 0).sum())
+
+    rng = np.random.default_rng(7)
+    B = rng.standard_normal((k, n))
+    Cref = A @ B
+
+    os.makedirs(os.path.join(BUILD, 'ocl'), exist_ok=True)
+    os.makedirs(os.path.join(BUILD, 'sycl'), exist_ok=True)
+
+    B.astype('<f8').tofile(os.path.join(BUILD, 'B.bin'))
+    Cref.astype('<f8').tofile(os.path.join(BUILD, 'Cref.bin'))
+
+    manifest = {
+        'm': m, 'k': k, 'n': n, 'ldb': ldb, 'ldc': ldc,
+        'nnz': nnz, 'sparsity': sparsity, 'dtype': 'double', 'dsize': 8,
+        'nbix': int(np.count_nonzero(np.any(A != 0, axis=0))),
+        'kernels': [],
+    }
+
+    backends = [('opencl', OpenCLMatMul, 'ocl', 'cl'),
+                ('sycl', SYCLMatMul, 'sycl', 'cpp')]
+
+    for plat, cls, subdir, ext in backends:
+        mm = cls(A, beta=0.0, n=n, ldb=ldb, ldc=ldc)
+        for idx, (src, meta) in enumerate(mm.kernels(np.float64, kname='gimmik_mm')):
+            tpl = meta['tplname']
+            width = meta['width']
+            entry = f'gmk_{plat}_{idx}_{tpl.replace("-", "_")}_w{width}'
+            src = src.replace('gimmik_mm', entry)
+            fname = f'{entry}.{ext}'
+            with open(os.path.join(BUILD, subdir, fname), 'w') as f:
+                f.write(src)
+
+            gws = meta.get('global_work_size')
+            lws = meta.get('local_work_size')
+            manifest['kernels'].append({
+                'platform': plat, 'entry': entry, 'file': f'{subdir}/{fname}',
+                'tpl': tpl, 'width': width,
+                'gws': list(gws) if gws else None,
+                'lws': list(lws) if lws else None,
+                'local_mem_size': meta.get('local_mem_size', 0),
+            })
+
+    with open(os.path.join(BUILD, 'manifest.json'), 'w') as f:
+        json.dump(manifest, f, indent=2)
+
+    # Shared problem dimensions
+    with open(os.path.join(BUILD, 'dims.h'), 'w') as f:
+        f.write('#pragma once\n')
+        f.write(f'#define GMK_M {m}\n#define GMK_K {k}\n#define GMK_N {n}\n')
+        f.write(f'#define GMK_LDB {ldb}\n#define GMK_LDC {ldc}\n')
+        f.write(f'#define GMK_NNZ {nnz}\n#define GMK_NBIX {manifest["nbix"]}\n')
+
+    ocl = [x for x in manifest['kernels'] if x['platform'] == 'opencl']
+    scl = [x for x in manifest['kernels'] if x['platform'] == 'sycl']
+
+    # OpenCL registry (entry name, source file, work sizes)
+    with open(os.path.join(BUILD, 'ocl_registry.h'), 'w') as f:
+        f.write('#pragma once\n')
+        f.write('typedef struct { const char* entry; const char* file; '
+                'const char* tpl; int gdim; size_t g0,g1; size_t l0,l1; } '
+                'OclKernel;\n')
+        f.write('static const OclKernel g_ocl[] = {\n')
+        for x in ocl:
+            g = x['gws']
+            l = x['lws']
+            gdim = len(g)
+            g0 = g[0]
+            g1 = g[1] if gdim > 1 else 1
+            l0 = l[0] if l else 0
+            l1 = l[1] if l else 0
+            f.write(f'  {{"{x["entry"]}", "{x["file"]}", "{x["tpl"]}", '
+                    f'{gdim}, {g0}, {g1}, {l0}, {l1}}},\n')
+        f.write('};\n')
+        f.write(f'static const int g_ocl_n = {len(ocl)};\n')
+
+    # SYCL registry (declarations + function pointer table)
+    with open(os.path.join(BUILD, 'sycl_registry.cpp'), 'w') as f:
+        f.write('#include "sycl_common.hpp"\n')
+        for x in scl:
+            f.write(f'sycl::event {x["entry"]}'
+                    f'(sycl::queue&, const double*, double*);\n')
+        f.write('const SyclKernel g_sycl[] = {\n')
+        for x in scl:
+            f.write(f'  {{"{x["entry"]}", "{x["tpl"]}", &{x["entry"]}}},\n')
+        f.write('};\n')
+        f.write(f'const int g_sycl_n = {len(scl)};\n')
+
+    print(f'A: {m}x{k}, nnz={nnz} ({100*nnz/(m*k):.1f}% dense), '
+          f'used B rows={manifest["nbix"]}')
+    print(f'n={n}  B={B.nbytes/1e6:.1f} MB  C={Cref.nbytes/1e6:.1f} MB')
+    print(f'Generated {len(manifest["kernels"])} kernels into {BUILD}')
+
+
+if __name__ == '__main__':
+    main()
diff --git a/bench/ocl_bench.cpp b/bench/ocl_bench.cpp
new file mode 100644
index 0000000..a70a94e
--- /dev/null
+++ b/bench/ocl_bench.cpp
@@ -0,0 +1,183 @@
+// OpenCL host benchmark for GiMMiK-generated kernels on the Intel GPU.
+#define CL_TARGET_OPENCL_VERSION 300
+#include <CL/cl.h>
+
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <algorithm>
+#include <string>
+#include <vector>
+
+#include "build/dims.h"
+#include "build/ocl_registry.h"
+
+#define CK(x) do { cl_int _e = (x); if (_e != CL_SUCCESS) { \
+    std::fprintf(stderr, "OpenCL error %d at %s:%d (%s)\n", _e, __FILE__, \
+                 __LINE__, #x); std::exit(1); } } while (0)
+
+static std::vector<double> read_bin(const char* path, size_t count) {
+    std::vector<double> v(count);
+    FILE* f = std::fopen(path, "rb");
+    if (!f) { std::fprintf(stderr, "cannot open %s\n", path); std::exit(1); }
+    if (std::fread(v.data(), sizeof(double), count, f) != count) {
+        std::fprintf(stderr, "short read %s\n", path); std::exit(1);
+    }
+    std::fclose(f);
+    return v;
+}
+
+static std::string read_text(const std::string& path) {
+    FILE* f = std::fopen(path.c_str(), "rb");
+    if (!f) { std::fprintf(stderr, "cannot open %s\n", path.c_str()); std::exit(1); }
+    std::fseek(f, 0, SEEK_END);
+    long sz = std::ftell(f);
+    std::fseek(f, 0, SEEK_SET);
+    std::string s(sz, '\0');
+    if (std::fread(&s[0], 1, sz, f) != (size_t)sz) { std::exit(1); }
+    std::fclose(f);
+    return s;
+}
+
+static cl_device_id pick_gpu(cl_platform_id* out_plat) {
+    cl_uint np = 0;
+    CK(clGetPlatformIDs(0, nullptr, &np));
+    std::vector<cl_platform_id> plats(np);
+    CK(clGetPlatformIDs(np, plats.data(), nullptr));
+    for (auto p : plats) {
+        cl_uint nd = 0;
+        if (clGetDeviceIDs(p, CL_DEVICE_TYPE_GPU, 0, nullptr, &nd) != CL_SUCCESS
+            || nd == 0)
+            continue;
+        std::vector<cl_device_id> devs(nd);
+        CK(clGetDeviceIDs(p, CL_DEVICE_TYPE_GPU, nd, devs.data(), nullptr));
+        char vendor[256] = {0};
+        clGetDeviceInfo(devs[0], CL_DEVICE_VENDOR, sizeof(vendor), vendor, nullptr);
+        if (std::strstr(vendor, "Intel") || std::strstr(vendor, "INTEL")) {
+            *out_plat = p;
+            return devs[0];
+        }
+    }
+    std::fprintf(stderr, "No Intel GPU found\n");
+    std::exit(1);
+}
+
+int main(int argc, char** argv) {
+    const int reps = argc > 1 ? std::atoi(argv[1]) : 50;
+    const size_t nB = (size_t)GMK_K * GMK_N;
+    const size_t nC = (size_t)GMK_M * GMK_N;
+
+    auto B = read_bin("bench/build/B.bin", nB);
+    auto Cref = read_bin("bench/build/Cref.bin", nC);
+
+    cl_platform_id plat;
+    cl_device_id dev = pick_gpu(&plat);
+    char name[256] = {0};
+    clGetDeviceInfo(dev, CL_DEVICE_NAME, sizeof(name), name, nullptr);
+    std::printf("# OpenCL device: %s\n", name);
+
+    cl_int err;
+    cl_context ctx = clCreateContext(nullptr, 1, &dev, nullptr, nullptr, &err);
+    CK(err);
+    cl_command_queue_properties qp[] = {CL_QUEUE_PROPERTIES,
+                                        CL_QUEUE_PROFILING_ENABLE, 0};
+    cl_command_queue q = clCreateCommandQueueWithProperties(ctx, dev, qp, &err);
+    CK(err);
+
+    cl_mem dB = clCreateBuffer(ctx, CL_MEM_READ_ONLY, nB * sizeof(double),
+                               nullptr, &err); CK(err);
+    cl_mem dC = clCreateBuffer(ctx, CL_MEM_READ_WRITE, nC * sizeof(double),
+                               nullptr, &err); CK(err);
+    CK(clEnqueueWriteBuffer(q, dB, CL_TRUE, 0, nB * sizeof(double), B.data(),
+                            0, nullptr, nullptr));
+
+    std::vector<double> Cout(nC);
+    const double gflop = 2.0 * GMK_NNZ * (double)GMK_N / 1e9;
+    const double gbyte = (double)(GMK_NBIX + GMK_M) * GMK_N * 8.0 / 1e9;
+
+    std::printf("# %-26s %10s %10s %10s %12s\n",
+                "kernel", "ms", "GFLOP/s", "GB/s", "max_relerr");
+
+    for (int ki = 0; ki < g_ocl_n; ki++) {
+        const OclKernel& K = g_ocl[ki];
+        std::string src = read_text(std::string("bench/build/") + K.file);
+        const char* csrc = src.c_str();
+        size_t slen = src.size();
+
+        cl_program prog = clCreateProgramWithSource(ctx, 1, &csrc, &slen, &err);
+        CK(err);
+        cl_int be = clBuildProgram(prog, 1, &dev, "-cl-std=CL2.0", nullptr, nullptr);
+        if (be != CL_SUCCESS) {
+            size_t ls = 0;
+            clGetProgramBuildInfo(prog, dev, CL_PROGRAM_BUILD_LOG, 0, nullptr, &ls);
+            std::string log(ls, '\0');
+            clGetProgramBuildInfo(prog, dev, CL_PROGRAM_BUILD_LOG, ls, &log[0], nullptr);
+            std::fprintf(stderr, "build failed for %s:\n%s\n", K.entry, log.c_str());
+            std::exit(1);
+        }
+        cl_kernel kern = clCreateKernel(prog, K.entry, &err); CK(err);
+        CK(clSetKernelArg(kern, 0, sizeof(cl_mem), &dB));
+        CK(clSetKernelArg(kern, 1, sizeof(cl_mem), &dC));
+
+        // rounded global work size
+        size_t gws[2] = {(size_t)K.g0, (size_t)K.g1};
+        size_t lws[2] = {(size_t)K.l0, (size_t)K.l1};
+        const size_t* lp = nullptr;
+        if (K.l0 > 0) {
+            lp = lws;
+            gws[0] = ((gws[0] + lws[0] - 1) / lws[0]) * lws[0];
+        }
+
+        // correctness
+        double zero = 0.0;
+        CK(clEnqueueFillBuffer(q, dC, &zero, sizeof(double), 0,
+                               nC * sizeof(double), 0, nullptr, nullptr));
+        CK(clEnqueueNDRangeKernel(q, kern, K.gdim, nullptr, gws, lp, 0,
+                                  nullptr, nullptr));
+        CK(clFinish(q));
+        CK(clEnqueueReadBuffer(q, dC, CL_TRUE, 0, nC * sizeof(double),
+                               Cout.data(), 0, nullptr, nullptr));
+        double maxrel = 0.0;
+        for (size_t i = 0; i < nC; i++) {
+            double d = std::fabs(Cout[i] - Cref[i]);
+            double r = d / (std::fabs(Cref[i]) + 1e-30);
+            maxrel = std::max(maxrel, r);
+        }
+
+        // warmup
+        for (int w = 0; w < 5; w++)
+            CK(clEnqueueNDRangeKernel(q, kern, K.gdim, nullptr, gws, lp, 0,
+                                      nullptr, nullptr));
+        CK(clFinish(q));
+
+        double best_ms = 1e30;
+        for (int r = 0; r < reps; r++) {
+            cl_event ev;
+            CK(clEnqueueNDRangeKernel(q, kern, K.gdim, nullptr, gws, lp, 0,
+                                      nullptr, &ev));
+            CK(clWaitForEvents(1, &ev));
+            cl_ulong t0, t1;
+            clGetEventProfilingInfo(ev, CL_PROFILING_COMMAND_START,
+                                    sizeof(t0), &t0, nullptr);
+            clGetEventProfilingInfo(ev, CL_PROFILING_COMMAND_END,
+                                    sizeof(t1), &t1, nullptr);
+            best_ms = std::min(best_ms, (t1 - t0) * 1e-6);
+            clReleaseEvent(ev);
+        }
+
+        std::printf("%-28s %10.4f %10.1f %10.1f %12.2e  %s\n",
+                    K.entry, best_ms, gflop / (best_ms * 1e-3),
+                    gbyte / (best_ms * 1e-3), maxrel,
+                    maxrel < 1e-9 ? "OK" : "FAIL");
+
+        clReleaseKernel(kern);
+        clReleaseProgram(prog);
+    }
+
+    clReleaseMemObject(dB);
+    clReleaseMemObject(dC);
+    clReleaseCommandQueue(q);
+    clReleaseContext(ctx);
+    return 0;
+}
diff --git a/bench/run.sh b/bench/run.sh
new file mode 100755
index 0000000..9d61bf2
--- /dev/null
+++ b/bench/run.sh
@@ -0,0 +1,31 @@
+#!/usr/bin/env bash
+# Build + run the OpenCL vs SYCL GiMMiK benchmark on the Intel GPU.
+# Run from the GiMMiK repository root.
+set -e
+
+REPS=${1:-50}
+ROOT=$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)
+cd "$ROOT"
+
+source /opt/intel/oneapi/compiler/latest/env/vars.sh >/dev/null 2>&1
+ONEAPI=/opt/intel/oneapi/compiler/latest
+
+echo "== Generating kernels =="
+python3 bench/bench_gen.py
+
+echo "== Building OpenCL host =="
+icpx -O3 -std=c++17 -Ibench -I"$ONEAPI/include" \
+     bench/ocl_bench.cpp -o bench/build/ocl_bench \
+     -L"$ONEAPI/lib" -lOpenCL
+
+echo "== Building SYCL host =="
+icpx -fsycl -O3 -std=c++17 -Ibench \
+     bench/sycl_bench.cpp bench/build/sycl_registry.cpp bench/build/sycl/*.cpp \
+     -o bench/build/sycl_bench
+
+echo
+echo "########## OpenCL ##########"
+bench/build/ocl_bench "$REPS"
+echo
+echo "########## SYCL ##########"
+ONEAPI_DEVICE_SELECTOR='level_zero:gpu' bench/build/sycl_bench "$REPS"
diff --git a/bench/sweep.sh b/bench/sweep.sh
new file mode 100644
index 0000000..23b3f3c
--- /dev/null
+++ b/bench/sweep.sh
@@ -0,0 +1,40 @@
+#!/usr/bin/env bash
+# Sweep several operator shapes and report the best OpenCL vs SYCL kernel each.
+set -e
+ROOT=$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)
+cd "$ROOT"
+REPS=${1:-100}
+
+# m  k   sparsity   label
+CONFIGS=(
+  "32 48 0.5"
+  "64 64 0.6"
+  "96 96 0.7"
+  "24 40 0.3"
+  "125 64 0.8"
+)
+
+printf "%-16s %-22s %10s %10s | %-22s %10s %10s\n" \
+  "shape(m x k)" "OpenCL best" "ms" "GB/s" "SYCL best" "ms" "GB/s"
+printf '%.0s-' {1..112}; echo
+
+for cfg in "${CONFIGS[@]}"; do
+  read m k sp <<< "$cfg"
+  GMK_M=$m GMK_K=$k GMK_N=${GMK_N:-200000} GMK_SPARSITY=$sp \
+    bash bench/run.sh "$REPS" > bench/build/sweep_out.txt 2>&1 || {
+      echo "run failed for $cfg"; cat bench/build/sweep_out.txt; exit 1; }
+
+  # best (min ms) line per platform
+  ocl=$(awk '/OpenCL ####/{p=1;next}/SYCL ####/{p=0}p&&/^gmk_/{print $1,$2,$4}' \
+        bench/build/sweep_out.txt | sort -k2 -n | head -1)
+  scl=$(awk '/SYCL ####/{p=1;next}p&&/^gmk_/{print $1,$2,$4}' \
+        bench/build/sweep_out.txt | sort -k2 -n | head -1)
+
+  oname=$(echo $ocl | awk '{gsub("gmk_opencl_[0-9]+_","",$1);gsub("_w1","",$1);print $1}')
+  oms=$(echo $ocl | awk '{print $2}');  ogb=$(echo $ocl | awk '{print $3}')
+  sname=$(echo $scl | awk '{gsub("gmk_sycl_[0-9]+_","",$1);gsub("_w1","",$1);print $1}')
+  sms=$(echo $scl | awk '{print $2}');  sgb=$(echo $scl | awk '{print $3}')
+
+  printf "%-16s %-22s %10s %10s | %-22s %10s %10s\n" \
+    "${m} x ${k} (${sp})" "$oname" "$oms" "$ogb" "$sname" "$sms" "$sgb"
+done
diff --git a/bench/sycl_bench.cpp b/bench/sycl_bench.cpp
new file mode 100644
index 0000000..89b2a60
--- /dev/null
+++ b/bench/sycl_bench.cpp
@@ -0,0 +1,89 @@
+// SYCL host benchmark for GiMMiK-generated kernels on the Intel GPU.
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <algorithm>
+#include <vector>
+#include <string>
+
+#include <sycl/sycl.hpp>
+#include "sycl_common.hpp"
+#include "build/dims.h"
+
+static std::vector<double> read_bin(const char* path, size_t count) {
+    std::vector<double> v(count);
+    FILE* f = std::fopen(path, "rb");
+    if (!f) { std::fprintf(stderr, "cannot open %s\n", path); std::exit(1); }
+    if (std::fread(v.data(), sizeof(double), count, f) != count) {
+        std::fprintf(stderr, "short read %s\n", path); std::exit(1);
+    }
+    std::fclose(f);
+    return v;
+}
+
+int main(int argc, char** argv) {
+    const int reps = argc > 1 ? std::atoi(argv[1]) : 50;
+    const size_t nB = (size_t)GMK_K * GMK_N;
+    const size_t nC = (size_t)GMK_M * GMK_N;
+
+    auto B = read_bin("bench/build/B.bin", nB);
+    auto Cref = read_bin("bench/build/Cref.bin", nC);
+
+    sycl::queue q{sycl::gpu_selector_v,
+                  sycl::property::queue::enable_profiling()};
+    std::printf("# SYCL device: %s\n",
+                q.get_device().get_info<sycl::info::device::name>().c_str());
+
+    double* dB = sycl::malloc_device<double>(nB, q);
+    double* dC = sycl::malloc_device<double>(nC, q);
+    q.memcpy(dB, B.data(), nB * sizeof(double)).wait();
+
+    std::vector<double> Cout(nC);
+
+    // useful flops = 2*nnz*n ; bytes moved = (used B rows + m)*n*8
+    const double gflop = 2.0 * GMK_NNZ * (double)GMK_N / 1e9;
+    const double gbyte = (double)(GMK_NBIX + GMK_M) * GMK_N * 8.0 / 1e9;
+
+    std::printf("# %-26s %10s %10s %10s %12s\n",
+                "kernel", "ms", "GFLOP/s", "GB/s", "max_relerr");
+
+    for (int ki = 0; ki < g_sycl_n; ki++) {
+        const SyclKernel& K = g_sycl[ki];
+
+        // correctness
+        q.memset(dC, 0, nC * sizeof(double)).wait();
+        K.fn(q, dB, dC).wait();
+        q.memcpy(Cout.data(), dC, nC * sizeof(double)).wait();
+        double maxrel = 0.0;
+        for (size_t i = 0; i < nC; i++) {
+            double d = std::fabs(Cout[i] - Cref[i]);
+            double r = d / (std::fabs(Cref[i]) + 1e-30);
+            maxrel = std::max(maxrel, r);
+        }
+
+        // warmup
+        for (int w = 0; w < 5; w++) K.fn(q, dB, dC);
+        q.wait();
+
+        // timed: use device profiling, take the min
+        double best_ms = 1e30;
+        for (int r = 0; r < reps; r++) {
+            sycl::event e = K.fn(q, dB, dC);
+            e.wait();
+            auto t0 = e.get_profiling_info<
+                sycl::info::event_profiling::command_start>();
+            auto t1 = e.get_profiling_info<
+                sycl::info::event_profiling::command_end>();
+            best_ms = std::min(best_ms, (t1 - t0) * 1e-6);
+        }
+
+        std::printf("%-28s %10.4f %10.1f %10.1f %12.2e  %s\n",
+                    K.name, best_ms, gflop / (best_ms * 1e-3),
+                    gbyte / (best_ms * 1e-3), maxrel,
+                    maxrel < 1e-9 ? "OK" : "FAIL");
+    }
+
+    sycl::free(dB, q);
+    sycl::free(dC, q);
+    return 0;
+}
diff --git a/bench/sycl_common.hpp b/bench/sycl_common.hpp
new file mode 100644
index 0000000..caf5d14
--- /dev/null
+++ b/bench/sycl_common.hpp
@@ -0,0 +1,11 @@
+#pragma once
+#include <sycl/sycl.hpp>
+
+struct SyclKernel {
+    const char* name;
+    const char* tpl;
+    sycl::event (*fn)(sycl::queue&, const double*, double*);
+};
+
+extern const SyclKernel g_sycl[];
+extern const int g_sycl_n;
diff --git a/bench/validate.sh b/bench/validate.sh
new file mode 100644
index 0000000..bbce525
--- /dev/null
+++ b/bench/validate.sh
@@ -0,0 +1,25 @@
+#!/usr/bin/env bash
+# Build + run the correctness-validation suite (beta, fp32/float2, static+dynamic).
+set -e
+ROOT=$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)
+cd "$ROOT"
+source /opt/intel/oneapi/compiler/latest/env/vars.sh >/dev/null 2>&1
+ONEAPI=/opt/intel/oneapi/compiler/latest
+
+echo "== Generating validation kernels =="
+python3 bench/validate_gen.py
+
+echo "== Building OpenCL validator =="
+icpx -O2 -std=c++17 -Ibench -I"$ONEAPI/include" \
+     bench/validate_ocl.cpp -o bench/vbuild/validate_ocl \
+     -L"$ONEAPI/lib" -lOpenCL
+
+echo "== Building SYCL validator =="
+icpx -fsycl -O2 -std=c++17 -Ibench \
+     bench/validate_sycl.cpp bench/vbuild/vsycl_registry.cpp bench/vbuild/sycl/*.cpp \
+     -o bench/vbuild/validate_sycl
+
+echo; echo "########## OpenCL validation ##########"
+bench/vbuild/validate_ocl
+echo; echo "########## SYCL validation ##########"
+ONEAPI_DEVICE_SELECTOR='level_zero:gpu' bench/vbuild/validate_sycl
diff --git a/bench/validate_gen.py b/bench/validate_gen.py
new file mode 100644
index 0000000..798e133
--- /dev/null
+++ b/bench/validate_gen.py
@@ -0,0 +1,169 @@
+#!/usr/bin/env python3
+"""Generate a correctness-validation suite for the OpenCL and SYCL backends.
+
+Covers the axes that the perf benchmark skips:
+  * beta == 0 and beta != 0
+  * fp64 and fp32 (including the fp32 float2 vectorized variants)
+  * static-n and dynamic-n kernel signatures
+
+For every (dtype, beta, mode, aligne) case it emits *all* kernel variants from
+both backends, shared input/reference data, and C/C++ registries so the host
+programs can run+verify each kernel against a NumPy fp64 reference.
+"""
+import json
+import os
+import sys
+
+import numpy as np
+
+sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+from gimmik import OpenCLMatMul, SYCLMatMul
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+BUILD = os.path.join(HERE, 'vbuild')
+
+# (dtype, beta, mode, aligne)
+CASES = [
+    ('float64', 0.0, 'static',  None),
+    ('float64', 0.5, 'static',  None),
+    ('float64', 0.0, 'dynamic', None),
+    ('float64', 0.5, 'dynamic', None),
+    ('float32', 0.0, 'static',  None),
+    ('float32', 0.5, 'static',  None),
+    ('float32', 0.0, 'static',  2),     # enables float2 vector kernels
+    ('float32', 0.0, 'dynamic', None),
+]
+
+M, K, N = 32, 48, 8192   # N even & divisible by 64 for the split/vector kernels
+
+
+def elem_type(dtype, width):
+    if width > 1:
+        return 'sycl::float2'
+    return 'float' if dtype == 'float32' else 'double'
+
+
+def main():
+    os.makedirs(os.path.join(BUILD, 'ocl'), exist_ok=True)
+    os.makedirs(os.path.join(BUILD, 'sycl'), exist_ok=True)
+
+    rng = np.random.default_rng(42)
+    A = rng.standard_normal((M, K))
+    A[rng.random((M, K)) < 0.5] = 0.0
+    for j in range(M):
+        if not A[j].any():
+            A[j, rng.integers(K)] = rng.standard_normal()
+
+    ldb = ldc = N
+    cases_meta = []
+    ocl_kernels = []
+    sycl_kernels = []
+
+    for cid, (dtype, beta, mode, aligne) in enumerate(CASES):
+        npdt = np.dtype(dtype)
+        code = 1 if dtype == 'float32' else 0
+        dsize = npdt.itemsize
+
+        # shared data in the target dtype (exactly what the GPU reads)
+        rB = rng.standard_normal((K, N)).astype(npdt)
+        rCi = rng.standard_normal((M, N)).astype(npdt)
+        Cref = A @ rB.astype(np.float64) + beta * rCi.astype(np.float64)
+
+        rB.tofile(os.path.join(BUILD, f'case{cid}_B.bin'))
+        rCi.tofile(os.path.join(BUILD, f'case{cid}_Ci.bin'))
+        Cref.astype('<f8').tofile(os.path.join(BUILD, f'case{cid}_Cref.bin'))
+
+        cases_meta.append({'id': cid, 'dtype': dtype, 'code': code,
+                           'dsize': dsize, 'beta': beta, 'mode': mode,
+                           'm': M, 'k': K, 'n': N, 'ldb': ldb, 'ldc': ldc})
+
+        if mode == 'static':
+            kw = dict(n=N, ldb=ldb, ldc=ldc)
+        else:
+            kw = dict()
+
+        for plat, cls, sub, ext in [('opencl', OpenCLMatMul, 'ocl', 'cl'),
+                                    ('sycl', SYCLMatMul, 'sycl', 'cpp')]:
+            mm = cls(A, beta=beta, aligne=aligne, **kw)
+            for idx, (src, meta) in enumerate(
+                    mm.kernels(npdt.type, kname='gimmik_mm')):
+                tpl, width = meta['tplname'], meta['width']
+                entry = f'gmk_{plat}_c{cid}_{idx}_{tpl.replace("-", "_")}_w{width}'
+                src = src.replace('gimmik_mm', entry)
+                with open(os.path.join(BUILD, sub, f'{entry}.{ext}'), 'w') as f:
+                    f.write(src)
+
+                lws = meta.get('local_work_size')
+                rec = {'entry': entry, 'file': f'{sub}/{entry}.{ext}',
+                       'tpl': tpl, 'width': width, 'case': cid, 'mode': mode,
+                       'lws': list(lws) if lws else None}
+                if plat == 'opencl':
+                    ocl_kernels.append(rec)
+                else:
+                    rec['etype'] = elem_type(dtype, width)
+                    sycl_kernels.append(rec)
+
+    # ---- cases header ----
+    with open(os.path.join(BUILD, 'vcases.h'), 'w') as f:
+        f.write('#pragma once\n')
+        f.write('typedef struct { int id; int code; int dsize; double beta; '
+                'int is_dynamic; int m,k,n,ldb,ldc; } VCase;\n')
+        f.write('static const VCase g_cases[] = {\n')
+        for c in cases_meta:
+            f.write(f'  {{{c["id"]}, {c["code"]}, {c["dsize"]}, {c["beta"]}, '
+                    f'{1 if c["mode"]=="dynamic" else 0}, {c["m"]}, {c["k"]}, '
+                    f'{c["n"]}, {c["ldb"]}, {c["ldc"]}}},\n')
+        f.write('};\n')
+        f.write(f'static const int g_cases_n = {len(cases_meta)};\n')
+
+    # ---- OpenCL registry ----
+    with open(os.path.join(BUILD, 'vocl_registry.h'), 'w') as f:
+        f.write('#pragma once\n')
+        f.write('typedef struct { const char* entry; const char* file; '
+                'const char* tpl; int cas; int width; int is_dynamic; '
+                'int has_lws; size_t l0,l1; } VOcl;\n')
+        f.write('static const VOcl g_vocl[] = {\n')
+        for r in ocl_kernels:
+            l0 = r['lws'][0] if r['lws'] else 0
+            l1 = r['lws'][1] if r['lws'] else 0
+            f.write(f'  {{"{r["entry"]}", "{r["file"]}", "{r["tpl"]}", '
+                    f'{r["case"]}, {r["width"]}, '
+                    f'{1 if r["mode"]=="dynamic" else 0}, '
+                    f'{1 if r["lws"] else 0}, {l0}, {l1}}},\n')
+        f.write('};\n')
+        f.write(f'static const int g_vocl_n = {len(ocl_kernels)};\n')
+
+    # ---- SYCL registry (typed wrappers -> uniform thunk) ----
+    with open(os.path.join(BUILD, 'vsycl_registry.cpp'), 'w') as f:
+        f.write('#include "vsycl_common.hpp"\n')
+        for r in sycl_kernels:
+            et = r['etype']
+            if r['mode'] == 'static':
+                f.write(f'sycl::event {r["entry"]}(sycl::queue&, '
+                        f'const {et}*, {et}*);\n')
+            else:
+                f.write(f'sycl::event {r["entry"]}(sycl::queue&, int, '
+                        f'const {et}*, int, {et}*, int);\n')
+        for r in sycl_kernels:
+            et = r['entry']
+            pt = r['etype']
+            f.write(f'static sycl::event w_{et}(sycl::queue& q, void* b, '
+                    f'void* c, int n, int ldb, int ldc) {{ return {et}(q, ')
+            if r['mode'] == 'static':
+                f.write(f'(const {pt}*)b, ({pt}*)c); }}\n')
+            else:
+                f.write(f'n, (const {pt}*)b, ldb, ({pt}*)c, ldc); }}\n')
+        f.write('const VSycl g_vsycl[] = {\n')
+        for r in sycl_kernels:
+            f.write(f'  {{"{r["entry"]}", "{r["tpl"]}", {r["case"]}, '
+                    f'&w_{r["entry"]}}},\n')
+        f.write('};\n')
+        f.write(f'const int g_vsycl_n = {len(sycl_kernels)};\n')
+
+    print(f'A: {M}x{K}, nnz={int((A!=0).sum())}')
+    print(f'cases={len(cases_meta)}  opencl_kernels={len(ocl_kernels)}  '
+          f'sycl_kernels={len(sycl_kernels)}')
+
+
+if __name__ == '__main__':
+    main()
diff --git a/bench/validate_ocl.cpp b/bench/validate_ocl.cpp
new file mode 100644
index 0000000..b64fde7
--- /dev/null
+++ b/bench/validate_ocl.cpp
@@ -0,0 +1,135 @@
+// OpenCL correctness validation across dtype/beta/mode cases.
+#define CL_TARGET_OPENCL_VERSION 300
+#include <CL/cl.h>
+
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <string>
+#include <vector>
+
+#include "vbuild/vcases.h"
+#include "vbuild/vocl_registry.h"
+
+#define CK(x) do { cl_int _e=(x); if(_e!=CL_SUCCESS){ \
+    std::fprintf(stderr,"CL err %d @%d (%s)\n",_e,__LINE__,#x); std::exit(1);} } while(0)
+
+static std::vector<char> read_bin(const std::string& p, size_t bytes) {
+    std::vector<char> v(bytes);
+    FILE* f = std::fopen(p.c_str(), "rb");
+    if (!f) { std::fprintf(stderr, "open %s\n", p.c_str()); std::exit(1); }
+    if (std::fread(v.data(), 1, bytes, f) != bytes) { std::exit(1); }
+    std::fclose(f);
+    return v;
+}
+static std::string read_text(const std::string& p) {
+    FILE* f = std::fopen(p.c_str(), "rb");
+    if (!f) { std::fprintf(stderr, "open %s\n", p.c_str()); std::exit(1); }
+    std::fseek(f, 0, SEEK_END); long s = std::ftell(f); std::fseek(f, 0, SEEK_SET);
+    std::string t(s, '\0');
+    if (std::fread(&t[0], 1, s, f) != (size_t)s) std::exit(1);
+    std::fclose(f); return t;
+}
+
+int main() {
+    cl_uint np; CK(clGetPlatformIDs(0, nullptr, &np));
+    std::vector<cl_platform_id> pl(np); CK(clGetPlatformIDs(np, pl.data(), nullptr));
+    cl_device_id dev = nullptr;
+    for (auto p : pl) {
+        cl_uint nd = 0;
+        if (clGetDeviceIDs(p, CL_DEVICE_TYPE_GPU, 0, nullptr, &nd) == CL_SUCCESS && nd) {
+            std::vector<cl_device_id> d(nd);
+            clGetDeviceIDs(p, CL_DEVICE_TYPE_GPU, nd, d.data(), nullptr);
+            char v[128] = {0};
+            clGetDeviceInfo(d[0], CL_DEVICE_VENDOR, sizeof(v), v, nullptr);
+            if (std::strstr(v, "Intel")) { dev = d[0]; break; }
+        }
+    }
+    if (!dev) { std::fprintf(stderr, "no Intel GPU\n"); return 1; }
+    char name[128] = {0};
+    clGetDeviceInfo(dev, CL_DEVICE_NAME, sizeof(name), name, nullptr);
+    std::printf("# OpenCL device: %s\n", name);
+
+    cl_int err;
+    cl_context ctx = clCreateContext(nullptr, 1, &dev, nullptr, nullptr, &err); CK(err);
+    cl_command_queue q = clCreateCommandQueueWithProperties(ctx, dev, nullptr, &err); CK(err);
+
+    std::printf("# %-34s %-9s %-6s %-7s %11s  %s\n",
+                "kernel", "dtype", "beta", "mode", "norm_err", "status");
+
+    int fails = 0;
+    for (int i = 0; i < g_vocl_n; i++) {
+        const VOcl& K = g_vocl[i];
+        const VCase& c = g_cases[K.cas];
+        const size_t nB = (size_t)c.k * c.n, nC = (size_t)c.m * c.n;
+        const bool f32 = (c.code == 1);
+
+        auto B = read_bin("bench/vbuild/case" + std::to_string(c.id) + "_B.bin", nB * c.dsize);
+        auto Ci = read_bin("bench/vbuild/case" + std::to_string(c.id) + "_Ci.bin", nC * c.dsize);
+        auto Cr = read_bin("bench/vbuild/case" + std::to_string(c.id) + "_Cref.bin", nC * sizeof(double));
+        const double* Cref = reinterpret_cast<const double*>(Cr.data());
+
+        cl_mem dB = clCreateBuffer(ctx, CL_MEM_READ_ONLY, nB * c.dsize, nullptr, &err); CK(err);
+        cl_mem dC = clCreateBuffer(ctx, CL_MEM_READ_WRITE, nC * c.dsize, nullptr, &err); CK(err);
+        CK(clEnqueueWriteBuffer(q, dB, CL_TRUE, 0, nB * c.dsize, B.data(), 0, nullptr, nullptr));
+        CK(clEnqueueWriteBuffer(q, dC, CL_TRUE, 0, nC * c.dsize, Ci.data(), 0, nullptr, nullptr));
+
+        std::string src = read_text("bench/vbuild/" + std::string(K.file));
+        const char* cs = src.c_str(); size_t sl = src.size();
+        cl_program pr = clCreateProgramWithSource(ctx, 1, &cs, &sl, &err); CK(err);
+        if (clBuildProgram(pr, 1, &dev, "-cl-std=CL2.0", nullptr, nullptr) != CL_SUCCESS) {
+            size_t ls; clGetProgramBuildInfo(pr, dev, CL_PROGRAM_BUILD_LOG, 0, nullptr, &ls);
+            std::string lg(ls, '\0');
+            clGetProgramBuildInfo(pr, dev, CL_PROGRAM_BUILD_LOG, ls, &lg[0], nullptr);
+            std::fprintf(stderr, "build %s:\n%s\n", K.entry, lg.c_str()); std::exit(1);
+        }
+        cl_kernel kern = clCreateKernel(pr, K.entry, &err); CK(err);
+
+        // args + work sizes depend on static/dynamic
+        size_t gws[2], lws[2] = {(size_t)K.l0, (size_t)K.l1};
+        const size_t* lp = K.has_lws ? lws : nullptr;
+        cl_int nn = c.n, ldb = c.ldb, ldc = c.ldc;
+        size_t cols = (c.n + K.width - 1) / K.width;   // work-items along n
+        if (K.is_dynamic) {
+            CK(clSetKernelArg(kern, 0, sizeof(cl_int), &nn));
+            CK(clSetKernelArg(kern, 1, sizeof(cl_mem), &dB));
+            CK(clSetKernelArg(kern, 2, sizeof(cl_int), &ldb));
+            CK(clSetKernelArg(kern, 3, sizeof(cl_mem), &dC));
+            CK(clSetKernelArg(kern, 4, sizeof(cl_int), &ldc));
+        } else {
+            CK(clSetKernelArg(kern, 0, sizeof(cl_mem), &dB));
+            CK(clSetKernelArg(kern, 1, sizeof(cl_mem), &dC));
+        }
+        gws[0] = cols; gws[1] = K.has_lws ? (size_t)K.l1 : 1;
+        int gdim = K.has_lws ? 2 : 1;
+        if (K.has_lws) gws[0] = ((cols + K.l0 - 1) / K.l0) * K.l0;
+
+        CK(clEnqueueNDRangeKernel(q, kern, gdim, nullptr, gws, lp, 0, nullptr, nullptr));
+        CK(clFinish(q));
+
+        std::vector<char> out(nC * c.dsize);
+        CK(clEnqueueReadBuffer(q, dC, CL_TRUE, 0, nC * c.dsize, out.data(), 0, nullptr, nullptr));
+
+        double maxabs = 0.0, maxerr = 0.0;
+        for (size_t e = 0; e < nC; e++) {
+            double got = f32 ? (double)reinterpret_cast<float*>(out.data())[e]
+                             : reinterpret_cast<double*>(out.data())[e];
+            maxabs = std::max(maxabs, std::fabs(Cref[e]));
+            maxerr = std::max(maxerr, std::fabs(got - Cref[e]));
+        }
+        double norm = maxerr / (maxabs + 1e-300);
+        double tol = f32 ? 1e-5 : 1e-12;
+        bool ok = norm < tol; fails += !ok;
+
+        std::printf("%-36s %-9s %-6.2f %-7s %11.2e  %s\n",
+                    K.entry, f32 ? "float32" : "float64", c.beta,
+                    c.is_dynamic ? "dyn" : "static", norm, ok ? "OK" : "FAIL");
+
+        clReleaseKernel(kern); clReleaseProgram(pr);
+        clReleaseMemObject(dB); clReleaseMemObject(dC);
+    }
+    std::printf("# %d/%d passed\n", g_vocl_n - fails, g_vocl_n);
+    clReleaseCommandQueue(q); clReleaseContext(ctx);
+    return fails ? 1 : 0;
+}
diff --git a/bench/validate_sycl.cpp b/bench/validate_sycl.cpp
new file mode 100644
index 0000000..e8cb5d2
--- /dev/null
+++ b/bench/validate_sycl.cpp
@@ -0,0 +1,75 @@
+// SYCL correctness validation across dtype/beta/mode cases.
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <string>
+#include <vector>
+
+#include <sycl/sycl.hpp>
+#include "vsycl_common.hpp"
+#include "vbuild/vcases.h"
+
+static std::vector<char> read_bin(const std::string& p, size_t bytes) {
+    std::vector<char> v(bytes);
+    FILE* f = std::fopen(p.c_str(), "rb");
+    if (!f) { std::fprintf(stderr, "open %s\n", p.c_str()); std::exit(1); }
+    if (std::fread(v.data(), 1, bytes, f) != bytes) { std::exit(1); }
+    std::fclose(f);
+    return v;
+}
+
+int main() {
+    sycl::queue q{sycl::gpu_selector_v};
+    std::printf("# SYCL device: %s\n",
+                q.get_device().get_info<sycl::info::device::name>().c_str());
+    std::printf("# %-34s %-9s %-6s %-7s %11s  %s\n",
+                "kernel", "dtype", "beta", "mode", "norm_err", "status");
+
+    int fails = 0;
+    for (int i = 0; i < g_vsycl_n; i++) {
+        const VSycl& K = g_vsycl[i];
+        const VCase& c = g_cases[K.cas];
+        const size_t nB = (size_t)c.k * c.n, nC = (size_t)c.m * c.n;
+        const bool f32 = (c.code == 1);
+
+        auto B = read_bin("bench/vbuild/case" + std::to_string(c.id) + "_B.bin",
+                          nB * c.dsize);
+        auto Ci = read_bin("bench/vbuild/case" + std::to_string(c.id) + "_Ci.bin",
+                           nC * c.dsize);
+        auto Cr = read_bin("bench/vbuild/case" + std::to_string(c.id) + "_Cref.bin",
+                           nC * sizeof(double));
+        const double* Cref = reinterpret_cast<const double*>(Cr.data());
+
+        void* dB = sycl::malloc_device(nB * c.dsize, q);
+        void* dC = sycl::malloc_device(nC * c.dsize, q);
+        q.memcpy(dB, B.data(), nB * c.dsize).wait();
+        q.memcpy(dC, Ci.data(), nC * c.dsize).wait();   // beta*C uses initial C
+
+        K.fn(q, dB, dC, c.n, c.ldb, c.ldc).wait();
+
+        std::vector<char> out(nC * c.dsize);
+        q.memcpy(out.data(), dC, nC * c.dsize).wait();
+
+        double maxabs = 0.0, maxerr = 0.0;
+        for (size_t e = 0; e < nC; e++) {
+            double got = f32 ? (double)reinterpret_cast<float*>(out.data())[e]
+                             : reinterpret_cast<double*>(out.data())[e];
+            maxabs = std::max(maxabs, std::fabs(Cref[e]));
+            maxerr = std::max(maxerr, std::fabs(got - Cref[e]));
+        }
+        double norm = maxerr / (maxabs + 1e-300);
+        double tol = f32 ? 1e-5 : 1e-12;
+        bool ok = norm < tol;
+        fails += !ok;
+
+        std::printf("%-36s %-9s %-6.2f %-7s %11.2e  %s\n",
+                    K.name, f32 ? "float32" : "float64", c.beta,
+                    c.is_dynamic ? "dyn" : "static", norm, ok ? "OK" : "FAIL");
+
+        sycl::free(dB, q);
+        sycl::free(dC, q);
+    }
+    std::printf("# %d/%d passed\n", g_vsycl_n - fails, g_vsycl_n);
+    return fails ? 1 : 0;
+}
diff --git a/bench/vsycl_common.hpp b/bench/vsycl_common.hpp
new file mode 100644
index 0000000..cdaa533
--- /dev/null
+++ b/bench/vsycl_common.hpp
@@ -0,0 +1,15 @@
+#pragma once
+#include <sycl/sycl.hpp>
+
+// uniform thunk: (queue, b, c, n, ldb, ldc) -> event
+typedef sycl::event (*VRunFn)(sycl::queue&, void*, void*, int, int, int);
+
+struct VSycl {
+    const char* name;
+    const char* tpl;
+    int cas;
+    VRunFn fn;
+};
+
+extern const VSycl g_vsycl[];
+extern const int g_vsycl_n;
diff --git a/gimmik/__init__.py b/gimmik/__init__.py
index cd21134..3b66583 100644
--- a/gimmik/__init__.py
+++ b/gimmik/__init__.py
@@ -9,6 +9,7 @@
 from gimmik.metal import MetalMatMul
 from gimmik.opencl import OpenCLMatMul
 from gimmik.ptx import PTXMatMul
+from gimmik.sycl import SYCLMatMul
 
 
 def generate_mm(mat, dtype, platform, alpha=1.0, beta=0.0, funcn='gimmik_mm',
@@ -24,7 +25,8 @@ def generate_mm(mat, dtype, platform, alpha=1.0, beta=0.0, funcn='gimmik_mm',
         'ispc': ISPCMatMul,
         'hip': HIPMatMul,
         'opencl': OpenCLMatMul,
-        'ptx': PTXMatMul
+        'ptx': PTXMatMul,
+        'sycl': SYCLMatMul
     }
 
     mm = platmap[platform](alpha*mat, beta, None, n, ldb, ldc)
diff --git a/gimmik/kernels/sycl/bstream-msplit.mako b/gimmik/kernels/sycl/bstream-msplit.mako
new file mode 100644
index 0000000..70864de
--- /dev/null
+++ b/gimmik/kernels/sycl/bstream-msplit.mako
@@ -0,0 +1,102 @@
+<%
+mx = partition(A, into=msplit, by='rows')
+bchunks = chunk(bix, bsz)
+%>\
+#include <sycl/sycl.hpp>
+
+sycl::event
+% if n is None:
+${kname}(sycl::queue& q, int n,
+         const ${dtype}* __restrict b, int ldb,
+         ${dtype}* __restrict c, int ldc)
+{
+    const int nw = (n + ${width} - 1) / ${width};
+  % if width > 1:
+    ldb /= ${width};
+    ldc /= ${width};
+  % endif
+% else:
+${kname}(sycl::queue& q, const ${dtype}* __restrict b, ${dtype}* __restrict c)
+{
+    const int nw = ${-(-n // width)};
+    const ${'long' if k*ldb >= width*2**31 else 'int'} ldb = ${ldb // width};
+    const ${'long' if m*ldc >= width*2**31 else 'int'} ldc = ${ldc // width};
+% endif
+    const int gx = ((nw + ${blockx} - 1) / ${blockx}) * ${blockx};
+    sycl::range<2> global(${msplit}, gx);
+    sycl::range<2> local(${msplit}, ${blockx});
+
+    return q.submit([&](sycl::handler& cgh) {
+        sycl::local_accessor<${dtype}, 1> bsub(
+            sycl::range<1>(${2 * bsz * blockx}), cgh);
+
+        cgh.parallel_for(sycl::nd_range<2>(global, local),
+                         [=](sycl::nd_item<2> it) [[sycl::reqd_work_group_size(${msplit}, ${blockx})]] {
+        const int i = it.get_global_id(1);
+        const int lx = it.get_local_id(1), ly = it.get_local_id(0);
+
+        ${dtype} bv, csub[${-(-m // msplit)}];
+
+## Fill the initial shared memory block
+% for cid in range(msplit):
+        if (i < nw && ly == ${cid})
+        {
+  % for kx in bchunks[0]:
+    % if loop.index % msplit == cid:
+            bsub[${loop.index * blockx} + lx] = b[i + ${kx}*ldb];
+    % endif
+  % endfor
+        }
+% endfor
+        it.barrier(sycl::access::fence_space::local_space);
+
+## Iterate over each row-chunk of B
+% for bb in range(len(bchunks)):
+  ## Iterate over each row-chunk of C
+  % for cid, mcx in enumerate(mx):
+        if (i < nw && ly == ${cid})
+        {
+    ## Start filling the next shared memory block
+    % if not loop.parent.last:
+      % for kx in bchunks[bb + 1]:
+        % if loop.index % msplit == cid:
+            bsub[${((bb + 1) % 2) * bsz * blockx + loop.index * blockx} + lx] = b[i + ${kx}*ldb];
+        % endif
+      % endfor
+    % endif
+    ## Accumulate our dot products
+    % for kx in bchunks[bb]:
+            bv = bsub[${(bb % 2) * bsz * blockx + loop.index * blockx} + lx];
+      % for j, jx in enumerate(A[mcx, kx]):
+        % if jx != 0 and kx == afix[mcx[j]]:
+            csub[${j}] = ${jx}*bv;
+        % elif jx != 0:
+            csub[${j}] += ${jx}*bv;
+        % endif
+        ## If we're done with this dot product then store to global
+        % if kx == alix[mcx[j]] and beta == 0:
+            c[i + ${mcx[j]}*ldc] = csub[${j}];
+        % elif kx == alix[mcx[j]] and beta == 1:
+            c[i + ${mcx[j]}*ldc] += csub[${j}];
+        % elif kx == alix[mcx[j]]:
+            c[i + ${mcx[j]}*ldc] = csub[${j}] + ${beta}*c[i + ${mcx[j]}*ldc];
+        % endif
+      % endfor
+    % endfor
+    ## Handle rows of A which are all zero
+    % if loop.parent.last:
+      % for j, jx in enumerate(afix):
+        % if jx == -1 and j % msplit == cid and beta == 0:
+            c[i + ${j}*ldc] = ${dtype}(0);
+        % elif jx == -1 and j % msplit == cid and beta != 1:
+            c[i + ${j}*ldc] *= ${beta};
+        % endif
+      % endfor
+    % endif
+        }
+  % endfor
+        it.barrier(sycl::access::fence_space::local_space);
+% endfor
+        });
+    });
+}
diff --git a/gimmik/kernels/sycl/bstream.mako b/gimmik/kernels/sycl/bstream.mako
new file mode 100644
index 0000000..386109e
--- /dev/null
+++ b/gimmik/kernels/sycl/bstream.mako
@@ -0,0 +1,64 @@
+#include <sycl/sycl.hpp>
+
+sycl::event
+% if n is None:
+${kname}(sycl::queue& q, int n,
+         const ${dtype}* __restrict bp, int ldb,
+         ${dtype}* __restrict cp, int ldc)
+{
+    const int nw = (n + ${width} - 1) / ${width};
+  % if width > 1:
+    ldb /= ${width};
+    ldc /= ${width};
+  % endif
+% else:
+${kname}(sycl::queue& q, const ${dtype}* __restrict bp, ${dtype}* __restrict cp)
+{
+    const int nw = ${-(-n // width)};
+    const ${'long' if k*ldb >= width*2**31 else 'int'} ldb = ${ldb // width};
+    const ${'long' if m*ldc >= width*2**31 else 'int'} ldc = ${ldc // width};
+% endif
+    const int gx = ((nw + ${blockx} - 1) / ${blockx}) * ${blockx};
+    return q.parallel_for(
+        sycl::nd_range<1>(sycl::range<1>(gx), sycl::range<1>(${blockx})),
+        [=](sycl::nd_item<1> it) [[sycl::reqd_work_group_size(${blockx})]] {
+        const int i = it.get_global_id(0);
+        if (i >= nw)
+            return;
+        // Re-assert non-aliasing: __restrict on the launcher's pointer
+        // parameters is lost once they are captured by value into the kernel
+        // lambda, so restore it on the pointers actually used in the body.
+        const ${dtype}* __restrict b = bp;
+        ${dtype}* __restrict c = cp;
+        ${dtype} bv, csub[${m}];
+
+## Iterate through the used rows of B
+% for kx in bix:
+        bv = b[i + ${kx}*ldb];
+  % for j, jx in enumerate(A[:, kx]):
+    % if jx != 0 and kx == afix[j]:
+        csub[${j}] = ${jx}*bv;
+    % elif jx != 0:
+        csub[${j}] += ${jx}*bv;
+    % endif
+    ##
+    % if kx == alix[j] and beta == 0:
+        c[i + ${j}*ldc] = csub[${j}];
+    % elif kx == alix[j] and beta == 1:
+        c[i + ${j}*ldc] += csub[${j}];
+    % elif kx == alix[j]:
+        c[i + ${j}*ldc] = csub[${j}] + ${beta}*c[i + ${j}*ldc];
+    % endif
+  % endfor
+% endfor
+
+## Handle rows of A which are all zero
+% for j, jx in enumerate(afix):
+  % if jx == -1 and beta == 0:
+        c[i + ${j}*ldc] = ${dtype}(0);
+  % elif jx == -1 and beta != 1:
+        c[i + ${j}*ldc] *= ${beta};
+  % endif
+% endfor
+    });
+}
diff --git a/gimmik/kernels/sycl/cstream-ksplit.mako b/gimmik/kernels/sycl/cstream-ksplit.mako
new file mode 100644
index 0000000..550dc48
--- /dev/null
+++ b/gimmik/kernels/sycl/cstream-ksplit.mako
@@ -0,0 +1,95 @@
+<%
+kparts = partition(A, ksplit, by='cols')
+cchunks = chunk(range(m), csz)
+loaded = set()
+%>\
+#include <sycl/sycl.hpp>
+
+sycl::event
+% if n is None:
+${kname}(sycl::queue& q, int n,
+         const ${dtype}* __restrict b, int ldb,
+         ${dtype}* __restrict c, int ldc)
+{
+    const int nw = (n + ${width} - 1) / ${width};
+  % if width > 1:
+    ldb /= ${width};
+    ldc /= ${width};
+  % endif
+% else:
+${kname}(sycl::queue& q, const ${dtype}* __restrict b, ${dtype}* __restrict c)
+{
+    const int nw = ${-(-n // width)};
+    const ${'long' if k*ldb >= width*2**31 else 'int'} ldb = ${ldb // width};
+    const ${'long' if m*ldc >= width*2**31 else 'int'} ldc = ${ldc // width};
+% endif
+    const int gx = ((nw + ${blockx} - 1) / ${blockx}) * ${blockx};
+    sycl::range<2> global(${ksplit}, gx);
+    sycl::range<2> local(${ksplit}, ${blockx});
+
+    return q.submit([&](sycl::handler& cgh) {
+        sycl::local_accessor<${dtype}, 1> csub(
+            sycl::range<1>(${(ksplit - 1) * csz * blockx}), cgh);
+
+        cgh.parallel_for(sycl::nd_range<2>(global, local),
+                         [=](sycl::nd_item<2> it) [[sycl::reqd_work_group_size(${ksplit}, ${blockx})]] {
+        const int i = it.get_global_id(1);
+        const int lx = it.get_local_id(1), ly = it.get_local_id(0);
+
+        ${dtype} cv[${-(-csz // ksplit)}], bv[${-(-k // ksplit)}], dotp;
+
+## Iterate over the row-partitions of C
+% for cchunk in cchunks:
+  ## Iterate over the row-partitions of B
+  % for bid, kbx in enumerate(kparts):
+        if (i < nw && ly == ${bid})
+        {
+    ## Evaluate our partial dot products
+    % for j in cchunk:
+      ## Load in any missing parts of B
+      % for kx in kbx:
+        % if A[j, kx] != 0 and kx not in loaded:
+            bv[${loop.index}] = b[i + ${kx}*ldb]; <% loaded.add(kx) %>
+        % endif
+      % endfor
+      % if (dotex := dot(lambda kx: f'bv[{kx}]', A[j, kbx])) != '0.0':
+            dotp = ${dotex};
+      % else:
+            dotp = ${dtype}(0);
+      % endif
+      ## Save to a register
+      % if loop.index % ksplit == bid:
+            cv[${loop.index // ksplit}] = dotp;
+      ## Save to shared memory
+      % else:
+            csub[${(bid - (bid > loop.index % ksplit)) * csz * blockx} + ${loop.index * blockx} + lx] = dotp;
+      % endif
+    % endfor
+        }
+  % endfor
+        it.barrier(sycl::access::fence_space::local_space);
+  ## Iterate over the column-partitions of B
+  % for bid, kbx in enumerate(kparts):
+        if (i < nw && ly == ${bid})
+        {
+    ## Sum and output the final set of dot products
+    % for j in cchunk:
+      % if loop.index % ksplit == bid:
+            dotp = cv[${loop.index // ksplit}] + ${' + '.join(f'csub[{i * csz * blockx + loop.index * blockx} + lx]'
+                                                              for i in range(ksplit - 1))};
+        % if beta == 0:
+            c[i + ${j}*ldc] = dotp;
+        % elif beta == 1:
+            c[i + ${j}*ldc] += dotp;
+        % else:
+            c[i + ${j}*ldc] = dotp + ${beta}*c[i + ${j}*ldc];
+        % endif
+      % endif
+    % endfor
+        }
+  % endfor
+        it.barrier(sycl::access::fence_space::local_space);
+% endfor
+        });
+    });
+}
diff --git a/gimmik/kernels/sycl/cstream.mako b/gimmik/kernels/sycl/cstream.mako
new file mode 100644
index 0000000..13e6c3c
--- /dev/null
+++ b/gimmik/kernels/sycl/cstream.mako
@@ -0,0 +1,44 @@
+#include <sycl/sycl.hpp>
+
+sycl::event
+% if n is None:
+${kname}(sycl::queue& q, int n,
+         const ${dtype}* __restrict bp, int ldb,
+         ${dtype}* __restrict cp, int ldc)
+{
+    const int nw = (n + ${width} - 1) / ${width};
+  % if width > 1:
+    ldb /= ${width};
+    ldc /= ${width};
+  % endif
+% else:
+${kname}(sycl::queue& q, const ${dtype}* __restrict bp, ${dtype}* __restrict cp)
+{
+    const int nw = ${-(-n // width)};
+    const ${'long' if k*ldb >= width*2**31 else 'int'} ldb = ${ldb // width};
+    const ${'long' if m*ldc >= width*2**31 else 'int'} ldc = ${ldc // width};
+% endif
+    const int gx = ((nw + ${blockx} - 1) / ${blockx}) * ${blockx};
+    return q.parallel_for(
+        sycl::nd_range<1>(sycl::range<1>(gx), sycl::range<1>(${blockx})),
+        [=](sycl::nd_item<1> it) [[sycl::reqd_work_group_size(${blockx})]] {
+        const int i = it.get_global_id(0);
+        if (i >= nw)
+            return;
+        // Re-assert non-aliasing: __restrict on the launcher's pointer
+        // parameters is lost once they are captured by value into the kernel
+        // lambda, so restore it on the pointers actually used in the body.
+        const ${dtype}* __restrict b = bp;
+        ${dtype}* __restrict c = cp;
+% for j, jx in enumerate(A):
+  % if beta == 0:
+        c[i + ${j}*ldc] = ${dot(lambda kx: f'b[i + {kx}*ldb]', jx)};
+  % elif beta == 1:
+        c[i + ${j}*ldc] += ${dot(lambda kx: f'b[i + {kx}*ldb]', jx)};
+  % else:
+        c[i + ${j}*ldc] = ${dot(lambda kx: f'b[i + {kx}*ldb]', jx)}
+                        + ${beta}*c[i + ${j}*ldc];
+  % endif
+% endfor
+    });
+}
diff --git a/gimmik/sycl.py b/gimmik/sycl.py
new file mode 100644
index 0000000..0854356
--- /dev/null
+++ b/gimmik/sycl.py
@@ -0,0 +1,68 @@
+# -*- coding: utf-8 -*-
+
+from gimmik.base import MatMul
+
+
+class SYCLMatMul(MatMul):
+    platform = 'sycl'
+    basemeta = {'local_work_size': None, 'local_mem_size': 0, 'width': 1}
+
+    def _kernel_generators(self, dtype, dsize, *, local_mem_size=None):
+        max_local_mem = local_mem_size or 1024**3
+
+        # Default 1D work-group size for the non-tiled kernels.  Launching
+        # these via an explicit nd_range (rather than a basic parallel_for)
+        # avoids the SYCL runtime auto-picking a poor work-group size.  256
+        # (four gfx90a wavefronts / eight NVIDIA warps) hides global-memory
+        # latency better than 128 on these bandwidth-bound kernels.
+        sblkx = 256
+
+        # B loading, C streaming kernel
+        yield ('cstream', {'blockx': sblkx}, {'local_work_size': (sblkx,)})
+
+        # B streaming, C accumulation kernel
+        yield ('bstream', {'blockx': sblkx}, {'local_work_size': (sblkx,)})
+
+        # Four-way m-split B streaming, C accumulation kernel
+        ms, bsz, blkx = 4, 16, 64
+        args = {'msplit': ms, 'blockx': blkx, 'bsz': bsz}
+        meta = {'local_work_size': (blkx, ms),
+                'local_mem_size': 2*blkx*bsz*dsize}
+        if meta['local_mem_size'] < max_local_mem:
+            yield ('bstream-msplit', args, meta)
+
+        # Two-way k-split B loading, C streaming kernel
+        ks, csz, blkx = 2, 32, 64
+        args = {'ksplit': ks, 'csz': csz, 'blockx': blkx}
+        meta = {'local_work_size': (blkx, ks),
+                'local_mem_size': (ks - 1)*csz*blkx*dsize}
+        if meta['local_mem_size'] < max_local_mem:
+            yield ('cstream-ksplit', args, meta)
+
+        # At single precision also consider vectorized kernels
+        if (dtype == 'float' and
+                self.aligne is not None and self.aligne % 2 == 0):
+            # Vector B loading, C streaming kernel
+            args = {'dtype': 'sycl::float2', 'width': 2, 'blockx': sblkx}
+            meta = {'width': 2, 'local_work_size': (sblkx,)}
+            yield ('cstream', args, meta)
+
+            # Vector four-way m-split B streaming, C accumulation kernel
+            ms, bsz, blkx = 4, 16, 64
+            args = {'dtype': 'sycl::float2', 'width': 2, 'msplit': ms,
+                    'blockx': blkx, 'bsz': bsz}
+            meta = {'local_work_size': (blkx, ms),
+                    'local_mem_size': 2*blkx*bsz*dsize, 'width': 2}
+            if meta['local_mem_size'] < max_local_mem:
+                yield ('bstream-msplit', args, meta)
+
+    def _process_meta(self, meta):
+        if self.n is not None:
+            lws, width = meta['local_work_size'], meta['width']
+            nx = -(-self.n // width)
+            if lws is None:
+                meta['global_work_size'] = (nx,)
+            elif len(lws) == 1:
+                meta['global_work_size'] = (-(-nx // lws[0]) * lws[0],)
+            else:
+                meta['global_work_size'] = (nx, lws[1])