diff --git a/bench/benchmark_reduce_sum.py b/bench/benchmark_reduce_sum.py
new file mode 100644
index 0000000..f44ad40
--- /dev/null
+++ b/bench/benchmark_reduce_sum.py
@@ -0,0 +1,47 @@
+import argparse
+
+import torch
+
+from forge_cute_py.ops import reduce_sum
+from forge_cute_py.ref import reduce_sum as ref_reduce_sum
+from forge_cute_py.util.bench import do_bench, estimate_bandwidth, summarize_times
+
+
+def main():
+    parser = argparse.ArgumentParser(description="Benchmark copy/transpose")
+    parser.add_argument("--m", type=int, default=1024)
+    parser.add_argument("--n", type=int, default=1024)
+    parser.add_argument("--dtype", choices=["float16", "bfloat16", "float32"], default="float16")
+    parser.add_argument("--dim", type=int, default=-1)
+    parser.add_argument("--warmup", type=int, default=10)
+    parser.add_argument("--iterations", type=int, default=100)
+    parser.add_argument("--compile-ref", action="store_true")
+    args = parser.parse_args()
+
+    if not torch.cuda.is_available():
+        raise RuntimeError("CUDA required for benchmarking")
+
+    dtype = getattr(torch, args.dtype)
+    x = torch.randn(args.m, args.n, device="cuda", dtype=dtype)
+    dim = args.dim
+
+    def fn():
+        return reduce_sum(x, dim=args.dim)
+
+    times = do_bench(fn, warmup=args.warmup, rep=args.iterations)
+    stats = summarize_times(times)
+    bytes_moved = (x.numel() + x.numel() / x.shape[dim]) * x.element_size()
+    bw = estimate_bandwidth(bytes_moved, stats["p50_ms"])
+    print(f"copy_transpose p50: {stats['p50_ms']:.4f} ms, BW: {bw:.2f} GB/s")
+
+    ref = lambda: ref_reduce_sum(x, dim=dim)
+    if args.compile_ref and hasattr(torch, "compile"):
+        ref = torch.compile(ref, fullgraph=True)
+    ref_times = do_bench(ref, warmup=args.warmup, rep=args.iterations)
+    ref_stats = summarize_times(ref_times)
+    ref_bw = estimate_bandwidth(bytes_moved, ref_stats["p50_ms"])
+    print(f"reference p50: {ref_stats['p50_ms']:.4f} ms, BW: {ref_bw:.2f} GB/s")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/forge_cute_py/kernels/reduce_sum.py b/forge_cute_py/kernels/reduce_sum.py
new file mode 100644
index 0000000..0865009
--- /dev/null
+++ b/forge_cute_py/kernels/reduce_sum.py
@@ -0,0 +1,111 @@
+"""
+Reudction kerne using CuTe DSL
+"""
+
+from typing import Literal, Type
+
+import cutlass
+import cutlass.cute as cute
+
+
+class Reduction:
+    def __init__(
+        self,
+        dtype: Type[cutlass.Numeric],
+        N: int,
+        reduction_dtype: Type[cutlass.Numeric] | None = cutlass.Float32,
+        reduction_op: Literal["sum", "amax", "amin", "prod"] = "sum",
+        dim: int = -1,
+    ):
+        self.dtype = dtype
+        self.N = N
+        self.reduction_dtype = reduction_dtype if reduction_dtype is not None else dtype
+        self.reduction_op = reduction_op
+        self.dim = dim
+
+        if self.dim not in (-1, 0, 1):
+            raise ValueError(f"dim must be either -1, 0 or 1. Got: {self.dim}")
+        if self.reduction_op not in ["sum", "amax", "amin", "prod"]:
+            raise ValueError(
+                f"reduction_op must be either 'sum', 'amax', 'amin', 'prod'. Got: {self.reduction_dtype}"
+            )
+
+        if self.dim not in [-1, 1]:
+            raise NotImplementedError(f"Only support dim=1 or -1, got {self.dim}")
+        if self.reduction_op != "sum":
+            raise NotImplementedError(f"Only support reduction_op=sum, got {self.reduction_op}")
+
+    def _get_tiled_copy(self, vecsize: int = 1):
+        """
+        Adapted from quack's tiles_copy_2d()
+        Reference: https://github.com/Dao-AILab/quack/blob/2e62faaeb6271a780a1360e6c96a003492e47eed/quack/copy_utils.py#L98
+        """
+        threads_per_row = 32
+        num_threads = 128
+        num_blocks_N = cute.ceil_div(self.N // vecsize, threads_per_row)
+        tiler_mn = (num_threads // threads_per_row, vecsize * num_blocks_N * threads_per_row)
+
+        num_copy_bits = vecsize * self.dtype.width
+        copy_op = cute.nvgpu.CopyUniversalOp()
+        copy_atom = cute.make_copy_atom(copy_op, self.dtype, num_bits_per_copy=num_copy_bits)
+        thr_layout = cute.make_ordered_layout(
+            (num_threads // threads_per_row, threads_per_row),
+            order=(1, 0),
+        )
+        val_layout = cute.make_layout((1, vecsize))
+        tiled_copy = cute.make_tiled_copy_tv(copy_atom, thr_layout, val_layout)
+        return tiler_mn, tiled_copy, threads_per_row
+
+    @cute.jit
+    def __call__(self, mX: cute.Tensor, mO: cute.Tensor, stream=None):
+        vecsize = 128 // self.dtype.width
+        tiler_mn, tiled_copy, threads_per_row = self._get_tiled_copy(vecsize=vecsize)
+
+        num_threads = tiled_copy.size
+
+        self.kernel(mX, mO, tiler_mn, tiled_copy, threads_per_row).launch(
+            grid=[cute.ceil_div(mX.shape[0], tiler_mn[0]), 1, 1],
+            block=[num_threads, 1, 1],
+            stream=stream,
+        )
+
+    @cute.kernel
+    def kernel(
+        self,
+        mX: cute.Tensor,
+        mO: cute.Tensor,
+        tiler_mn: cute.Shape,
+        tiled_copy: cute.TiledCopy,
+        threads_per_row: cutlass.Constexpr[int],
+    ):
+        # tv_layout = (thread_layout, value_layout) = ((threads_per_row, num_rows), vec_size)
+        tidx, _, _ = cute.arch.thread_idx()
+        bidx, _, _ = cute.arch.block_idx()
+
+        gX = cute.local_tile(mX, tiler_mn, (bidx, 0))  # (tileM, tileN)
+        # TODO: vectorized store
+        # gO = cute.local_tile(mO, cute.select(tiler_mn, mode=[0]), (bidx,))  # (tileM,)
+
+        thr_copy_X = tiled_copy.get_slice(tidx)
+        # gmem -> rmem
+        tXgX = thr_copy_X.partition_S(gX)
+        tXrX = cute.make_rmem_tensor_like(tXgX)
+        cute.autovec_copy(tXgX, tXrX)
+
+        # reduce with higher precision for numerical stability
+        x = tXrX.load().to(self.reduction_dtype)
+        val = x.reduce(cute.ReductionOp.ADD, init_val=0.0, reduction_profile=0)
+
+        val = cute.arch.warp_reduction_sum(val)
+
+        lane_id = cute.arch.lane_idx()
+        warp_id = cute.arch.warp_idx()
+
+        warps_per_row = threads_per_row // cute.arch.WARP_SIZE
+
+        row_idx = warp_id // warps_per_row
+        col_idx = warp_id % warps_per_row
+
+        # TODO: vetorized store
+        if lane_id == 0 and col_idx == 0:
+            mO[row_idx + tiler_mn[0] * bidx] = val.to(self.dtype)
diff --git a/forge_cute_py/ops/reduce_sum.py b/forge_cute_py/ops/reduce_sum.py
index d501eca..e2f7f03 100644
--- a/forge_cute_py/ops/reduce_sum.py
+++ b/forge_cute_py/ops/reduce_sum.py
@@ -1,4 +1,8 @@
+import cutlass.cute as cute
 import torch
+from cutlass import BFloat16, Float16, Float32
+
+from forge_cute_py.kernels.reduce_sum import Reduction
 
 
 @torch.library.custom_op("forge_cute_py::_reduce_sum", mutates_args={"out"})
@@ -21,12 +25,39 @@ def _reduce_sum(x: torch.Tensor, out: torch.Tensor, dim: int = -1, variant: str
     # Normalize dim to positive index
     dim = dim if dim >= 0 else x.ndim + dim
 
-    # For now, use reference implementation
-    # Future: call kernel implementation based on variant when available
-    from forge_cute_py.ref import reduce_sum as reduce_sum_ref
+    # Map PyTorch dtype to CUTLASS dtype
+    dtype_map = {
+        torch.float16: Float16,
+        torch.float32: Float32,
+        torch.bfloat16: BFloat16,
+    }
+    if x.dtype not in dtype_map:
+        raise ValueError(f"Unsupported dtype: {x.dtype}")
+
+    cute_dtype = dtype_map[x.dtype]
+    compile_key = (cute_dtype, dim, variant, x.shape[dim])
+
+    if compile_key not in _reduce_sum.compile_cache:
+        m = cute.sym_int() if dim != 0 else x.shape[0]
+        n = cute.sym_int() if dim != 1 else x.shape[1]
+        input_shape = (m, n)
+        output_shape = (m,) if dim == 1 else (n,)
+        input_cute = cute.runtime.make_fake_compact_tensor(
+            cute_dtype, input_shape, stride_order=(1, 0)
+        )
+        output_cute = cute.runtime.make_fake_compact_tensor(
+            cute_dtype, output_shape, stride_order=(0)
+        )
+        # Compile and cache the kernel
+        _reduce_sum.compile_cache[compile_key] = cute.compile(
+            Reduction(cute_dtype, n, reduction_op="sum", dim=dim, reduction_dtype=Float32),
+            input_cute,
+            output_cute,
+            cute.runtime.make_fake_stream(use_tvm_ffi_env_stream=True),
+            options="--enable-tvm-ffi",
+        )
 
-    result = reduce_sum_ref(x, dim=dim)
-    out.copy_(result)
+    _reduce_sum.compile_cache[compile_key](x, out)
 
 
 _reduce_sum.compile_cache = {}
@@ -63,3 +94,14 @@ def reduce_sum(x: torch.Tensor, dim: int = -1, variant: str = "shfl") -> torch.T
     out = torch.empty(out_shape, dtype=x.dtype, device=x.device)
     _reduce_sum(x, out, dim, variant)
     return out
+
+
+if __name__ == "__main__":
+    M = 1024
+    N = 1024
+    dtype = torch.float32
+    x = torch.randn(M, N, device="cuda", dtype=dtype)
+    y = reduce_sum(x, dim=-1)
+    ref_y = torch.sum(x, dim=-1)
+
+    torch.testing.assert_close(y, ref_y)
diff --git a/tests/test_reduce_sum.py b/tests/test_reduce_sum.py
index dd773b5..17ae3cc 100644
--- a/tests/test_reduce_sum.py
+++ b/tests/test_reduce_sum.py
@@ -8,8 +8,9 @@
 @pytest.mark.parametrize(
     "shape, dim",
     [
-        ((4, 8), -1),
-        ((8, 4), 0),
+        # ((4, 8), -1),
+        # ((8, 4), 0),
+        ((4096, 1024), 1),
     ],
 )
 @pytest.mark.parametrize(
@@ -25,8 +26,8 @@ def test_reduce_sum_correctness(shape, dim, dtype, atol, rtol, variant):
     x = torch.randn(*shape, device="cuda", dtype=dtype)
     try:
         y = reduce_sum(x, dim=dim, variant=variant)
-    except NotImplementedError:
-        pytest.skip(f"reduce_sum variant {variant} not implemented")
+    except NotImplementedError as e:
+        pytest.skip(f"reduce_sum variant {variant} not implemented\n{e}")
     y_ref = ref_reduce_sum(x, dim=dim)
     torch.testing.assert_close(y, y_ref, atol=atol, rtol=rtol)