diff --git a/.github/workflows/build-and-publish.yml b/.github/workflows/build-and-publish.yml
index d1b3a38..2e27fa1 100644
--- a/.github/workflows/build-and-publish.yml
+++ b/.github/workflows/build-and-publish.yml
@@ -64,7 +64,9 @@ jobs:
           CIBW_SKIP: "*-win32 *-manylinux_i686 *-musllinux*"
           CIBW_ARCHS_MACOS: "x86_64 universal2"
           CIBW_ARCHS_LINUX: "auto aarch64"
-
+          CIBW_BEFORE_BUILD_LINUX: "yum install -y openblas-devel || apt-get install -y libopenblas-dev"
+          CIBW_BEFORE_BUILD_WINDOWS: "vcpkg install openblas:x64-windows"
+          CIBW_ENVIRONMENT_WINDOWS: "CMAKE_TOOLCHAIN_FILE=C:/vcpkg/scripts/buildsystems/vcpkg.cmake"
       - name: Check wheels
         shell: bash
         run: |
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000..c391e91
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,29 @@
+cmake_minimum_required(VERSION 3.15)
+project(SparseDiffPy LANGUAGES C)
+
+add_subdirectory(SparseDiffEngine)
+
+find_package(Python3 REQUIRED COMPONENTS Interpreter Development.Module NumPy)
+
+# Linux: CMake find_package(BLAS) doesn't set include dirs;
+# openblas puts cblas.h in /usr/include/openblas/
+if(UNIX AND NOT APPLE)
+    find_path(OPENBLAS_INCLUDE_DIR cblas.h PATHS /usr/include/openblas /usr/include)
+    if(OPENBLAS_INCLUDE_DIR)
+        target_include_directories(dnlp_diff PRIVATE ${OPENBLAS_INCLUDE_DIR})
+    endif()
+endif()
+
+python3_add_library(_sparsediffengine MODULE
+    sparsediffpy/_bindings/bindings.c
+)
+
+target_include_directories(_sparsediffengine PRIVATE
+    ${CMAKE_CURRENT_SOURCE_DIR}/SparseDiffEngine/include
+    ${CMAKE_CURRENT_SOURCE_DIR}/SparseDiffEngine/src
+    ${CMAKE_CURRENT_SOURCE_DIR}/sparsediffpy/_bindings
+)
+
+target_link_libraries(_sparsediffengine PRIVATE dnlp_diff Python3::NumPy)
+
+install(TARGETS _sparsediffengine DESTINATION sparsediffpy)
diff --git a/SparseDiffEngine b/SparseDiffEngine
index 7c84b13..9480159 160000
--- a/SparseDiffEngine
+++ b/SparseDiffEngine
@@ -1 +1 @@
-Subproject commit 7c84b137ea7ac51846eaa5a30fe3ea707bc9d7aa
+Subproject commit 94801591f2f059ed3de6d68ab3db5a370bf3d667
diff --git a/pyproject.toml b/pyproject.toml
index f40d0e7..10465ff 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [build-system]
-requires = ["numpy >= 2.0.0", "setuptools >= 68.1.0", "wheel"]
-build-backend = "setuptools.build_meta"
+requires = ["scikit-build-core >= 0.10", "numpy >= 2.0.0"]
+build-backend = "scikit_build_core.build"
 
 [project]
 name = "sparsediffpy"
@@ -14,8 +14,6 @@ license = "Apache-2.0"
 file = "README.md"
 content-type = "text/markdown"
 
-[tool.setuptools]
-include-package-data = false
-
-[tool.setuptools.packages.find]
-include = ["sparsediffpy*"]
+[tool.scikit-build]
+cmake.minimum-version = "3.15"
+wheel.packages = ["sparsediffpy"]
diff --git a/setup.py b/setup.py
deleted file mode 100644
index 3b572ad..0000000
--- a/setup.py
+++ /dev/null
@@ -1,72 +0,0 @@
-"""
-Copyright 2025, the SparseDiffPy developers
-
-Licensed under the Apache License, Version 2.0 (the "License");
-you may not use this file except in compliance with the License.
-You may obtain a copy of the License at
-
-    http://www.apache.org/licenses/LICENSE-2.0
-
-Unless required by applicable law or agreed to in writing, software
-distributed under the License is distributed on an "AS IS" BASIS,
-WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-See the License for the specific language governing permissions and
-limitations under the License.
-"""
-
-import builtins
-import glob
-import platform
-
-from setuptools import Extension, setup
-from setuptools.command.build_ext import build_ext
-
-
-def not_on_windows(s: str) -> str:
-    return s if platform.system().lower() != "windows" else ""
-
-
-class build_ext_numpy(build_ext):
-    """Custom build_ext that injects NumPy headers."""
-
-    def finalize_options(self) -> None:
-        build_ext.finalize_options(self)
-        builtins.__NUMPY_SETUP__ = False
-        import numpy
-        self.include_dirs.append(numpy.get_include())
-
-
-# Collect all C source files from SparseDiffEngine
-diff_engine_sources = [
-    s for s in glob.glob("SparseDiffEngine/src/**/*.c", recursive=True)
-    if "dnlp_diff_engine" not in s
-] + ["sparsediffpy/_bindings/bindings.c"]
-
-# Define _POSIX_C_SOURCE on Linux for clock_gettime and struct timespec
-defines = []
-if platform.system().lower() == "linux":
-    defines.append(("_POSIX_C_SOURCE", "200809L"))
-
-sparsediffengine = Extension(
-    "sparsediffpy._sparsediffengine",
-    sources=diff_engine_sources,
-    include_dirs=[
-        "SparseDiffEngine/include/",
-        "SparseDiffEngine/src/",
-        "sparsediffpy/_bindings/",
-    ],
-    define_macros=defines,
-    extra_compile_args=[
-        "-O3",
-        "-std=c99",
-        "-Wall",
-        not_on_windows("-Wextra"),
-        '-DDIFF_ENGINE_VERSION="0.1.3"',
-    ],
-    extra_link_args=["-lm"] if platform.system().lower() != "windows" else [],
-)
-
-setup(
-    cmdclass={"build_ext": build_ext_numpy},
-    ext_modules=[sparsediffengine],
-)
diff --git a/sparsediffpy/_bindings/atoms/dense_matmul.h b/sparsediffpy/_bindings/atoms/dense_matmul.h
new file mode 100644
index 0000000..d8aef8f
--- /dev/null
+++ b/sparsediffpy/_bindings/atoms/dense_matmul.h
@@ -0,0 +1,101 @@
+#ifndef ATOM_DENSE_MATMUL_H
+#define ATOM_DENSE_MATMUL_H
+
+#include "bivariate.h"
+#include "common.h"
+
+/* Dense left matrix multiplication: A @ f(x) where A is a dense matrix.
+ *
+ * Python signature:
+ *   make_dense_left_matmul(child, A_data_flat, m, n)
+ *
+ * - child: the child expression capsule f(x).
+ * - A_data_flat: contiguous row-major numpy float64 array of size m*n.
+ * - m, n: dimensions of matrix A. */
+static PyObject *py_make_dense_left_matmul(PyObject *self, PyObject *args)
+{
+    PyObject *child_capsule;
+    PyObject *data_obj;
+    int m, n;
+    if (!PyArg_ParseTuple(args, "OOii", &child_capsule, &data_obj, &m, &n))
+    {
+        return NULL;
+    }
+
+    expr *child = (expr *) PyCapsule_GetPointer(child_capsule, EXPR_CAPSULE_NAME);
+    if (!child)
+    {
+        PyErr_SetString(PyExc_ValueError, "invalid child capsule");
+        return NULL;
+    }
+
+    PyArrayObject *data_array =
+        (PyArrayObject *) PyArray_FROM_OTF(data_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
+    if (!data_array)
+    {
+        return NULL;
+    }
+
+    double *A_data = (double *) PyArray_DATA(data_array);
+
+    expr *node = new_left_matmul_dense(child, m, n, A_data);
+    Py_DECREF(data_array);
+
+    if (!node)
+    {
+        PyErr_SetString(PyExc_RuntimeError,
+                        "failed to create dense_left_matmul node");
+        return NULL;
+    }
+    expr_retain(node);
+    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
+}
+
+/* Dense right matrix multiplication: f(x) @ A where A is a dense matrix.
+ *
+ * Python signature:
+ *   make_dense_right_matmul(child, A_data_flat, m, n)
+ *
+ * - child: the child expression capsule f(x).
+ * - A_data_flat: contiguous row-major numpy float64 array of size m*n.
+ * - m, n: dimensions of matrix A. */
+static PyObject *py_make_dense_right_matmul(PyObject *self, PyObject *args)
+{
+    PyObject *child_capsule;
+    PyObject *data_obj;
+    int m, n;
+    if (!PyArg_ParseTuple(args, "OOii", &child_capsule, &data_obj, &m, &n))
+    {
+        return NULL;
+    }
+
+    expr *child = (expr *) PyCapsule_GetPointer(child_capsule, EXPR_CAPSULE_NAME);
+    if (!child)
+    {
+        PyErr_SetString(PyExc_ValueError, "invalid child capsule");
+        return NULL;
+    }
+
+    PyArrayObject *data_array =
+        (PyArrayObject *) PyArray_FROM_OTF(data_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
+    if (!data_array)
+    {
+        return NULL;
+    }
+
+    double *A_data = (double *) PyArray_DATA(data_array);
+
+    expr *node = new_right_matmul_dense(child, m, n, A_data);
+    Py_DECREF(data_array);
+
+    if (!node)
+    {
+        PyErr_SetString(PyExc_RuntimeError,
+                        "failed to create dense_right_matmul node");
+        return NULL;
+    }
+    expr_retain(node);
+    return PyCapsule_New(node, EXPR_CAPSULE_NAME, expr_capsule_destructor);
+}
+
+#endif /* ATOM_DENSE_MATMUL_H */
\ No newline at end of file
diff --git a/sparsediffpy/_bindings/atoms/left_matmul.h b/sparsediffpy/_bindings/atoms/left_matmul.h
index 27fe3a4..09fa985 100644
--- a/sparsediffpy/_bindings/atoms/left_matmul.h
+++ b/sparsediffpy/_bindings/atoms/left_matmul.h
@@ -5,7 +5,7 @@
 #include "common.h"
 
 /* Left matrix multiplication: A @ f(x) where A is a constant matrix */
-static PyObject *py_make_left_matmul(PyObject *self, PyObject *args)
+static PyObject *py_make_sparse_left_matmul(PyObject *self, PyObject *args)
 {
     PyObject *child_capsule;
     PyObject *data_obj, *indices_obj, *indptr_obj;
diff --git a/sparsediffpy/_bindings/atoms/right_matmul.h b/sparsediffpy/_bindings/atoms/right_matmul.h
index c1c3481..f808058 100644
--- a/sparsediffpy/_bindings/atoms/right_matmul.h
+++ b/sparsediffpy/_bindings/atoms/right_matmul.h
@@ -5,7 +5,7 @@
 #include "common.h"
 
 /* Right matrix multiplication: f(x) @ A where A is a constant matrix */
-static PyObject *py_make_right_matmul(PyObject *self, PyObject *args)
+static PyObject *py_make_sparse_right_matmul(PyObject *self, PyObject *args)
 {
     PyObject *child_capsule;
     PyObject *data_obj, *indices_obj, *indptr_obj;
diff --git a/sparsediffpy/_bindings/bindings.c b/sparsediffpy/_bindings/bindings.c
index 2e51960..9f3c4b0 100644
--- a/sparsediffpy/_bindings/bindings.c
+++ b/sparsediffpy/_bindings/bindings.c
@@ -11,6 +11,7 @@
 #include "atoms/const_vector_mult.h"
 #include "atoms/constant.h"
 #include "atoms/cos.h"
+#include "atoms/dense_matmul.h"
 #include "atoms/diag_vec.h"
 #include "atoms/entr.h"
 #include "atoms/exp.h"
@@ -109,10 +110,14 @@ static PyMethodDef DNLPMethods[] = {
     {"make_entr", py_make_entr, METH_VARARGS, "Create entr node"},
     {"make_logistic", py_make_logistic, METH_VARARGS, "Create logistic node"},
     {"make_xexp", py_make_xexp, METH_VARARGS, "Create xexp node"},
-    {"make_left_matmul", py_make_left_matmul, METH_VARARGS,
-     "Create left matmul node (A @ f(x))"},
-    {"make_right_matmul", py_make_right_matmul, METH_VARARGS,
-     "Create right matmul node (f(x) @ A)"},
+    {"make_sparse_left_matmul", py_make_sparse_left_matmul, METH_VARARGS,
+     "Create sparse left matmul node (A @ f(x))"},
+    {"make_dense_left_matmul", py_make_dense_left_matmul, METH_VARARGS,
+     "Create dense left matmul node (A @ f(x)) where A is dense"},
+    {"make_sparse_right_matmul", py_make_sparse_right_matmul, METH_VARARGS,
+     "Create sparse right matmul node (f(x) @ A)"},
+    {"make_dense_right_matmul", py_make_dense_right_matmul, METH_VARARGS,
+     "Create dense right matmul node (f(x) @ A) where A is dense"},
     {"make_quad_form", py_make_quad_form, METH_VARARGS,
      "Create quadratic form node (x' * Q * x)"},
     {"make_quad_over_lin", py_make_quad_over_lin, METH_VARARGS,