Repository: coreylowman/dfdx (language: rust)

Issue:
Use `gemm` package for cpu matrix multiplication
Recently f16 support was added to gemm, and this is the package that supports `faer`.

Integration would include calling [gemm::gemm](https://docs.rs/gemm/0.14.8/gemm/fn.gemm.html).

From benchmarking f16 is at LEAST 1500x more efficient than the current dfdx implementation


Discussion / hints:
Blocked by https://github.com/sarah-ek/gemm/issues/2. Either that needs to be resolved, or we need to wait until half-rs releases a new version on crates.io and both update to that version
[can this be closed now or is still anything missing?]
Nope, still requires commits use the `gemm` package instead of matrixmultiply

Fix patch:
diff --git a/.github/workflows/cargo-check-features.yml b/.github/workflows/cargo-check-features.yml
index 6752a598d..6e622e25f 100644
--- a/.github/workflows/cargo-check-features.yml
+++ b/.github/workflows/cargo-check-features.yml
@@ -11,9 +11,9 @@ jobs:
       matrix:
         config:
           - toolchain: stable
-            command: cargo hack check --feature-powerset --no-dev-deps --depth 2 --skip default,nightly,cpu-mkl-matmul,cuda,cudnn
+            command: cargo hack check --feature-powerset --no-dev-deps --depth 2 --skip default,nightly,cuda,cudnn
           - toolchain: nightly
-            command: cargo hack check --each-feature --no-dev-deps --features nightly --skip default,cpu-mkl-matmul,cuda,cudnn
+            command: cargo hack check --each-feature --no-dev-deps --features nightly --skip default,cuda,cudnn
 
     steps:
       - uses: actions/checkout@v2
diff --git a/Cargo.toml b/Cargo.toml
index 52c4ca8de..e715a3ff6 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -29,7 +29,6 @@ no-std-compat = { version = "0.4.1", default-features = false, features = [ "all
 spin = { version = "0.9.8", default-features = false, features = ["spin_mutex", "rwlock", "portable_atomic"], optional = true }
 rand = { version = "0.8.5", default-features = false, features = ["std_rng"] }
 rand_distr = { version = "0.4.3", default-features = false }
-matrixmultiply = { version = "0.3.2", default-features = false, optional = true }
 zip = { version = "0.6.2", default-features = false, optional = true }
 cblas-sys = { version = "0.1.4", default-features = false, optional = true }
 libc = { version = "0.2", default-features = false, optional = true }
@@ -38,6 +37,7 @@ num-traits = { version = "0.2.15", default-features = false }
 safetensors = { version = "0.3", default-features = false, optional = true }
 memmap2 = { version = "0.5", default-features = false, optional = true }
 half = { git = "https://github.com/starkat99/half-rs.git", branch = "main", optional = true, features = ["num-traits", "rand_distr"] }
+gemm = { version = "0.15.3", default-features = false, optional = true }
 
 [dev-dependencies]
 tempfile = "3.3.0"
@@ -48,17 +48,14 @@ indicatif = "0.17.3"
 glob = { version = "0.3.1", optional = true }
 
 [features]
-default = ["std", "fast-alloc", "cpu-par-matmul"]
-nightly = ["half?/use-intrinsics"]
+default = ["std", "fast-alloc", "cpu"]
+nightly = ["half?/use-intrinsics", "gemm?/nightly"]
 
-std = ["cudarc?/std", "matrixmultiply?/std", "rand_distr/std_math"]
+std = ["cudarc?/std", "rand_distr/std_math", "gemm?/std"]
 fast-alloc = ["std"]
 no-std = ["no-std-compat", "dep:spin", "cudarc?/no-std", "num-traits/libm"]
 
-cpu-seq-matmul = ["dep:matrixmultiply"]
-cpu-par-matmul = ["std", "dep:matrixmultiply", "matrixmultiply?/threading"]
-cpu-mkl-matmul = ["dep:cblas-sys", "dep:libc"]
-
+cpu = ["dep:gemm"]
 cuda = ["dep:cudarc", "dep:glob"]
 cudnn = ["cuda", "cudarc?/cudnn"]
 
diff --git a/README.md b/README.md
index ace2a83d6..9079bcf1b 100644
--- a/README.md
+++ b/README.md
@@ -44,13 +44,6 @@ allocations when tensors are cloned.
 
 Enable the `cuda` feature to start using the `Cuda` device! Requires the installation of nvidia's cuda toolkit. See [feature flags docs](https://docs.rs/dfdx/latest/dfdx/feature_flags/index.html) for more info.
 
-## BLAS libraries
-
-The [matrixmultiply crate](https://crates.io/crates/matrixmultiply) is the default BLAS library. **You don't need
-to do download/install anything for this to work!**
-
-To link to the `Intel MKL` libraries (assuming you installed it already) use the `cpu-mkl-matmul` feature. See [feature flags docs](https://docs.rs/dfdx/latest/dfdx/feature_flags/index.html) for more info.
-
 ## API Preview
 
 Check [examples/](examples/) for more details.
diff --git a/build.rs b/build.rs
index 1b4bd1a9d..6f19b61e5 100644
--- a/build.rs
+++ b/build.rs
@@ -9,9 +9,6 @@ fn main() {
 
     #[cfg(feature = "cuda")]
     cuda::build_ptx();
-
-    #[cfg(feature = "cpu-mkl-matmul")]
-    intel_mkl::link().unwrap();
 }
 
 fn maybe_enable_nightly() {
@@ -213,212 +210,3 @@ mod cuda {
         }
     }
 }
-
-#[cfg(feature = "cpu-mkl-matmul")]
-mod intel_mkl {
-    //! This script links to Intel MKL when the `cpu-mkl-matmul` feature is enabled.
-    //! The dynamically linked and threaded implementation is chosen as
-    //! a good default, because:
-    //! 1. Dynamic compiles faster and makes smaller binaries
-    //! 2. Threaded implementation is faster for large workloads.
-    //!
-    //! As described [by Intel here](https://www.intel.com/content/www/us/en/developer/articles/technical/intel-math-kernel-library-intel-mkl-and-pkg-config-tool.html).
-    //!
-    //! **NOTE** lp64 vs ilp64 is chosen based on [target_pointer_width](https://doc.rust-lang.org/reference/conditional-compilation.html#target_pointer_width) cfg variable.
-    //!
-    //! # Adding a new system
-    //!
-    //! To add a new target system, the following blocks need to exist:
-    //!
-    //! - `SHARED_LIB_DIRS`: This should contain directories relative to `ONEAPI_ROOT` environment variable that contain the
-    //!     shared libraries. `main()` will check if any of these directories are not on the `PATH` environment variable and crash if not.
-    //! - `LINK_DIRS`: The directory where `.lib` files are. `main()` will output a `cargo:rustc-link-search` for each of these.
-    //!
-    //! # Supported systems
-    //!
-    //! - [x] Windows 32 bit
-    //! - [x] Windows 64 bit
-    //! - [x] Linux 32 bit
-    //! - [x] Linux 64 bit
-    //! - [x] MacOS 64 bit
-    //!
-    //! This script also creates a "nightly" feature if the crate is compiled on a nightly branch
-
-    pub const MKL_VERSION: &str = "2022.1.0";
-    pub const STATIC: bool = false;
-    pub const DYNAMIC: bool = true;
-    pub const SEQUENTIAL: bool = false;
-    pub const THREADED: bool = true;
-
-    pub const LINK_TYPE: &str = if STATIC { "static" } else { "dylib" };
-    pub const LIB_POSTFIX: &str = if cfg!(windows) && DYNAMIC { "_dll" } else { "" };
-    pub const LD_DIR: &str = if cfg!(windows) {
-        "PATH"
-    } else if cfg!(target_os = "linux") {
-        "LD_LIBRARY_PATH"
-    } else if cfg!(target_os = "macos") {
-        "DYLD_LIBRARY_PATH"
-    } else {
-        ""
-    };
-
-    pub const DEFAULT_ONEAPI_ROOT: &str = if cfg!(windows) {
-        "C:/Program Files (x86)/Intel/oneAPI/"
-    } else {
-        "/opt/intel/oneapi/"
-    };
-
-    pub const MKL_CORE: &str = "mkl_core";
-    pub const MKL_THREAD: &str = if SEQUENTIAL {
-        "mkl_sequential"
-    } else {
-        "mkl_intel_thread"
-    };
-    pub const THREADING_LIB: &str = if cfg!(windows) { "libiomp5md" } else { "iomp5" };
-    pub const MKL_INTERFACE: &str = if cfg!(target_pointer_width = "32") {
-        "mkl_intel_ilp64"
-    } else {
-        "mkl_intel_lp64"
-    };
-
-    #[cfg(not(any(target_os = "linux", target_os = "windows", target_os = "macos")))]
-    pub const UNSUPPORTED_OS_ERROR: _ = "Target OS is not supported. Please contact me";
-
-    #[cfg(all(target_os = "windows", target_pointer_width = "32"))]
-    pub const LINK_DIRS: &[&str] = &[
-        "compiler/latest/windows/compiler/lib/ia32_win",
-        "mkl/latest/lib/ia32",
-    ];
-
-    #[cfg(all(target_os = "windows", target_pointer_width = "64"))]
-    pub const LINK_DIRS: &[&str] = &[
-        "compiler/latest/windows/compiler/lib/intel64_win",
-        "mkl/latest/lib/intel64",
-    ];
-
-    #[cfg(all(target_os = "windows", target_pointer_width = "32"))]
-    pub const SHARED_LIB_DIRS: &[&str] = &[
-        "compiler/latest/windows/redist/ia32_win/compiler",
-        "mkl/latest/redist/ia32",
-    ];
-
-    #[cfg(all(target_os = "windows", target_pointer_width = "64"))]
-    pub const SHARED_LIB_DIRS: &[&str] = &[
-        "compiler/latest/windows/redist/intel64_win/compiler",
-        "mkl/latest/redist/intel64",
-    ];
-
-    #[cfg(all(target_os = "linux", target_pointer_width = "32"))]
-    pub const LINK_DIRS: &[&str] = &[
-        "compiler/latest/linux/compiler/lib/ia32_lin",
-        "mkl/latest/lib/ia32",
-    ];
-
-    #[cfg(all(target_os = "linux", target_pointer_width = "64"))]
-    pub const LINK_DIRS: &[&str] = &[
-        "compiler/latest/linux/compiler/lib/intel64_lin",
-        "mkl/latest/lib/intel64",
-    ];
-
-    #[cfg(target_os = "linux")]
-    pub const SHARED_LIB_DIRS: &[&str] = LINK_DIRS;
-
-    #[cfg(target_os = "macos")]
-    const MACOS_COMPILER_PATH: &str = "compiler/latest/mac/compiler/lib";
-
-    #[cfg(target_os = "macos")]
-    pub const LINK_DIRS: &[&str] = &[MACOS_COMPILER_PATH, "mkl/latest/lib"];
-
-    #[cfg(target_os = "macos")]
-    pub const SHARED_LIB_DIRS: &[&str] = &["mkl/latest/lib"];
-
-    #[derive(Debug)]
-    pub enum BuildError {
-        OneAPINotFound(std::path::PathBuf),
-        OneAPINotADir(std::path::PathBuf),
-        PathNotFound(std::env::VarError),
-        AddSharedLibDirToPath(String),
-    }
-
-    pub fn link() -> Result<(), BuildError> {
-        let root = std::env::var("ONEAPI_ROOT").unwrap_or_else(|_| DEFAULT_ONEAPI_ROOT.to_string());
-        println!("Using '{root}' as ONEAPI_ROOT");
-
-        let path = match std::env::var(LD_DIR) {
-            Ok(path) => path,
-            Err(e) => {
-                // On macOS it's unusual to set $DYLD_LIBRARY_PATH, so we want to provide a helpful message
-                println!(
-                    "Library path env variable '{LD_DIR}' was not found. Run `{}`",
-                    suggest_setvars_cmd(&root)
-                );
-                return Err(BuildError::PathNotFound(e));
-            }
-        };
-
-        if DYNAMIC {
-            // check to make sure that things in `SHARED_LIB_DIRS` are in `$LD_DIR`.
-            let path = path.replace('\\', "/");
-            for shared_lib_dir in SHARED_LIB_DIRS {
-                let versioned_dir = shared_lib_dir.replace("latest", MKL_VERSION);
-
-                println!("Checking that '{shared_lib_dir}' or '{versioned_dir}' is in {LD_DIR}");
-                if !path.contains(shared_lib_dir) && !path.contains(&versioned_dir) {
-                    println!(
-                        "'{shared_lib_dir}' not found in library path. Run `{}`",
-                        suggest_setvars_cmd(&root)
-                    );
-                    return Err(BuildError::AddSharedLibDirToPath(
-                        shared_lib_dir.to_string(),
-                    ));
-                }
-            }
-        }
-
-        let root: std::path::PathBuf = root.into();
-
-        if !root.exists() {
-            return Err(BuildError::OneAPINotFound(root));
-        }
-        if !root.is_dir() {
-            return Err(BuildError::OneAPINotADir(root));
-        }
-
-        for rel_lib_dir in LINK_DIRS {
-            let lib_dir = root.join(rel_lib_dir);
-            println!("cargo:rustc-link-search={}", lib_dir.display());
-        }
-
-        println!("cargo:rustc-link-lib={LINK_TYPE}={MKL_INTERFACE}{LIB_POSTFIX}");
-        println!("cargo:rustc-link-lib={LINK_TYPE}={MKL_THREAD}{LIB_POSTFIX}");
-        println!("cargo:rustc-link-lib={LINK_TYPE}={MKL_CORE}{LIB_POSTFIX}");
-        if THREADED {
-            println!("cargo:rustc-link-lib=dylib={THREADING_LIB}");
-        }
-
-        if !cfg!(windows) {
-            println!("cargo:rustc-link-lib=pthread");
-            println!("cargo:rustc-link-lib=m");
-            println!("cargo:rustc-link-lib=dl");
-        }
-
-        #[cfg(target_os = "macos")]
-        {
-            println!(
-                "cargo:rustc-link-arg=-Wl,-rpath,{}/{MACOS_COMPILER_PATH}",
-                root.display(),
-            );
-        }
-
-        Ok(())
-    }
-
-    // This section creates a feature "nightly" enabled if built on a nightly branch
-    fn suggest_setvars_cmd(root: &str) -> String {
-        if cfg!(windows) {
-            format!("{root}/setvars.bat")
-        } else {
-            format!("source {root}/setvars.sh")
-        }
-    }
-}
diff --git a/src/feature_flags.rs b/src/feature_flags.rs
index ba93ab145..6cbf64f56 100644
--- a/src/feature_flags.rs
+++ b/src/feature_flags.rs
@@ -9,12 +9,12 @@
 //!
 //! Cuda (with cudnn):
 //! ```toml
-//! dfdx = { version = "...", default-features = false, features = ["std", "cuda", "cudnn"]}
+//! dfdx = { version = "...", default-features = false, features = ["std", "fast-alloc", "cuda", "cudnn"]}
 //! ```
 //!
 //! Cpu:
 //! ```toml
-//! dfdx = { version = "...", default-features = false, features = ["std", "cpu-par-matmul"]}
+//! dfdx = { version = "...", default-features = false, features = ["std", "fast-alloc", "cpu"]}
 //! ```
 //!
 //! # "std"
@@ -70,38 +70,9 @@
 //! dfdx = { version = "...", features = ["cudnn"] }
 //! ```
 //!
-//! # "cpu-seq-matmul"
+//! # "cpu"
 //!
-//! Used to enable using `matrixmultiply` for matmul operations.
-//!
-//! # "cpu-par-matmul"
-//!
-//! Used to enable the threading feature of `matrixmultiply`. This makes matmuls
-//! substantially faster!
-//!
-//! # "cpu-mkl-matmul"
-//!
-//! Enables using the `Intel MKL` libraries (assuming you installed it already) for matrix multiplication.
-//!
-//! Linking is currently tested & verified on the following platforms:
-//!
-//! - [x] Windows
-//! - [x] Linux
-//! - [x] macOS
-//!
-//! Example:
-//! ```toml
-//! dfdx = { version = "...", features = ["cpu-mkl-matmul"] }
-//! ```
-//!
-//! #### Installing Intel MKL libraries
-//!
-//! It's pretty easy!
-//!
-//! You will need to install Intel MKL on your own from
-//! [this page](https://www.intel.com/content/www/us/en/developer/tools/oneapi/base-toolkit-download.html).
-//!
-//! `build.rs` will fail helpfully if you don't have the correct path/environment variables.
+//! Used to enable using `gemm` for matmul operations.
 //!
 //! # "numpy"
 //!
diff --git a/src/tensor_ops/matmul/cpu_kernel.rs b/src/tensor_ops/matmul/cpu_kernel.rs
index 4fb8522fa..65c29cd9e 100644
--- a/src/tensor_ops/matmul/cpu_kernel.rs
+++ b/src/tensor_ops/matmul/cpu_kernel.rs
@@ -5,18 +5,6 @@ use crate::tensor::{Cpu, Tensor, ZerosTensor};
 
 use std::sync::Arc;
 
-#[cfg(feature = "cpu-mkl-matmul")]
-use cblas_sys::{
-    cblas_dgemm as dgemm, cblas_sgemm as sgemm, CblasColMajor as ColMajor, CblasNoTrans as NoTr,
-    CblasRowMajor as RowMajor, CblasTrans as Tr,
-};
-
-#[cfg(all(
-    any(feature = "cpu-seq-matmul", feature = "cpu-par-matmul"),
-    not(feature = "cpu-mkl-matmul")
-))]
-use matrixmultiply::{dgemm, sgemm};
-
 #[allow(unused)]
 fn naive_gemm<F: num_traits::Float + std::ops::AddAssign, M: Dim, K: Dim, N: Dim>(
     (m, k, n): (M, K, N),
@@ -53,6 +41,48 @@ pub(crate) trait MatMulImpl<E> {
     );
 }
 
+#[cfg(feature = "f16")]
+impl MatMulImpl<half::f16> for Cpu {
+    #[inline]
+    fn matmul<M: Dim, K: Dim, N: Dim>(
+        (m, k, n): (M, K, N),
+        ap: *const half::f16,
+        a_strides: [usize; 2],
+        bp: *const half::f16,
+        b_strides: [usize; 2],
+        cp: *mut half::f16,
+        c_strides: [usize; 2],
+    ) {
+        #[cfg(not(feature = "cpu"))]
+        naive_gemm((m, k, n), ap, a_strides, bp, b_strides, cp, c_strides);
+
+        #[cfg(feature = "cpu")]
+        unsafe {
+            gemm::gemm(
+                m.size(),
+                n.size(),
+                k.size(),
+                cp as *mut gemm::f16,
+                c_strides[1] as isize,
+                c_strides[0] as isize,
+                true,
+                ap as *mut gemm::f16,
+                a_strides[1] as isize,
+                a_strides[0] as isize,
+                bp as *mut gemm::f16,
+                b_strides[1] as isize,
+                b_strides[0] as isize,
+                gemm::f16::ONE,
+                gemm::f16::ONE,
+                false,
+                false,
+                false,
+                gemm::Parallelism::None,
+            )
+        }
+    }
+}
+
 impl MatMulImpl<f32> for Cpu {
     #[inline]
     fn matmul<M: Dim, K: Dim, N: Dim>(
@@ -64,42 +94,33 @@ impl MatMulImpl<f32> for Cpu {
         cp: *mut f32,
         c_strides: [usize; 2],
     ) {
-        let (m, k, n) = (m.size(), k.size(), n.size());
-        #[cfg(feature = "cpu-mkl-matmul")]
-        unsafe {
-            let (lda, a_tr) = super::matrix_strides((m, k), a_strides);
-            let (ldb, b_tr) = super::matrix_strides((k, n), b_strides);
-            let (ldc, c_tr) = super::matrix_strides((m, n), c_strides);
-            let (m, n, k) = (m as libc::c_int, n as libc::c_int, k as libc::c_int);
-            let layout = if c_tr { ColMajor } else { RowMajor };
-            let (a_tr, b_tr) = if c_tr {
-                (if a_tr { NoTr } else { Tr }, if b_tr { NoTr } else { Tr })
-            } else {
-                (if a_tr { Tr } else { NoTr }, if b_tr { Tr } else { NoTr })
-            };
-            sgemm(
-                layout, a_tr, b_tr, m, n, k, 1.0, ap, lda as i32, bp, ldb as i32, 1.0, cp,
-                ldc as i32,
-            );
-        }
+        #[cfg(not(feature = "cpu"))]
+        naive_gemm((m, k, n), ap, a_strides, bp, b_strides, cp, c_strides);
 
-        #[cfg(all(
-            any(feature = "cpu-seq-matmul", feature = "cpu-par-matmul"),
-            not(feature = "cpu-mkl-matmul")
-        ))]
+        #[cfg(feature = "cpu")]
         unsafe {
-            let [ar, ac] = a_strides.map(|x| x as isize);
-            let [br, bc] = b_strides.map(|x| x as isize);
-            let [cr, cc] = c_strides.map(|x| x as isize);
-            sgemm(m, k, n, 1.0, ap, ar, ac, bp, br, bc, 1.0, cp, cr, cc);
+            gemm::gemm(
+                m.size(),
+                n.size(),
+                k.size(),
+                cp,
+                c_strides[1] as isize,
+                c_strides[0] as isize,
+                true,
+                ap,
+                a_strides[1] as isize,
+                a_strides[0] as isize,
+                bp,
+                b_strides[1] as isize,
+                b_strides[0] as isize,
+                1.0,
+                1.0,
+                false,
+                false,
+                false,
+                gemm::Parallelism::None,
+            )
         }
-
-        #[cfg(not(any(
-            feature = "cpu-seq-matmul",
-            feature = "cpu-par-matmul",
-            feature = "cpu-mkl-matmul"
-        )))]
-        naive_gemm((m, k, n), ap, a_strides, bp, b_strides, cp, c_strides);
     }
 }
 
@@ -114,45 +135,33 @@ impl MatMulImpl<f64> for Cpu {
         cp: *mut f64,
         c_strides: [usize; 2],
     ) {
-        let (m, k, n) = (m.size(), k.size(), n.size());
-
-        #[cfg(feature = "cpu-mkl-matmul")]
-        unsafe {
-            let (lda, a_tr) = super::matrix_strides((m, k), a_strides);
-            let (ldb, b_tr) = super::matrix_strides((k, n), b_strides);
-            let (ldc, c_tr) = super::matrix_strides((m, n), c_strides);
-            let (m, n, k) = (m as libc::c_int, n as libc::c_int, k as libc::c_int);
-            let layout = if c_tr { ColMajor } else { RowMajor };
-            let (a_tr, b_tr) = if c_tr {
-                (if a_tr { NoTr } else { Tr }, if b_tr { NoTr } else { Tr })
-            } else {
-                (if a_tr { Tr } else { NoTr }, if b_tr { Tr } else { NoTr })
-            };
-            dgemm(
-                layout, a_tr, b_tr, m, n, k, 1.0, ap, lda as i32, bp, ldb as i32, 1.0, cp,
-                ldc as i32,
-            );
-            return;
-        }
+        #[cfg(not(feature = "cpu"))]
+        naive_gemm((m, k, n), ap, a_strides, bp, b_strides, cp, c_strides);
 
-        #[cfg(all(
-            any(feature = "cpu-seq-matmul", feature = "cpu-par-matmul"),
-            not(feature = "cpu-mkl-matmul")
-        ))]
+        #[cfg(feature = "cpu")]
         unsafe {
-            let [ar, ac] = a_strides.map(|x| x as isize);
-            let [br, bc] = b_strides.map(|x| x as isize);
-            let [cr, cc] = c_strides.map(|x| x as isize);
-            dgemm(m, k, n, 1.0, ap, ar, ac, bp, br, bc, 1.0, cp, cr, cc);
-            return;
+            gemm::gemm(
+                m.size(),
+                n.size(),
+                k.size(),
+                cp,
+                c_strides[1] as isize,
+                c_strides[0] as isize,
+                true,
+                ap,
+                a_strides[1] as isize,
+                a_strides[0] as isize,
+                bp,
+                b_strides[1] as isize,
+                b_strides[0] as isize,
+                1.0,
+                1.0,
+                false,
+                false,
+                false,
+                gemm::Parallelism::None,
+            )
         }
-
-        #[cfg(not(any(
-            feature = "cpu-seq-matmul",
-            feature = "cpu-par-matmul",
-            feature = "cpu-mkl-matmul"
-        )))]
-        naive_gemm((m, k, n), ap, a_strides, bp, b_strides, cp, c_strides);
     }
 }
 
diff --git a/src/tensor_ops/matmul/f16.rs b/src/tensor_ops/matmul/f16.rs
deleted file mode 100644
index a115b591c..000000000
--- a/src/tensor_ops/matmul/f16.rs
+++ /dev/null
@@ -1,88 +0,0 @@
-#![allow(clippy::needless_range_loop)]
-
-use crate::{shapes::Dim, tensor::Cpu};
-
-use super::cpu_kernel::MatMulImpl;
-
-impl MatMulImpl<half::f16> for Cpu {
-    #[inline]
-    fn matmul<M: Dim, K: Dim, N: Dim>(
-        (m, k, n): (M, K, N),
-        ap: *const half::f16,
-        a_strides: [usize; 2],
-        bp: *const half::f16,
-        b_strides: [usize; 2],
-        cp: *mut half::f16,
-        c_strides: [usize; 2],
-    ) {
-        const CHUNK_SIZE: usize = 8;
-
-        let mut a_chunk: [f32; CHUNK_SIZE] = Default::default();
-        let mut b_chunk: [[f32; CHUNK_SIZE]; CHUNK_SIZE] = Default::default();
-        let mut c_chunk: [f32; CHUNK_SIZE] = Default::default();
-        for i_m in 0..m.size() {
-            let a_m = unsafe { ap.add(a_strides[0] * i_m) };
-            let c_m = unsafe { cp.add(c_strides[0] * i_m) };
-            for i_n_base in (0..n.size()).step_by(CHUNK_SIZE) {
-                let n_chunk_size = CHUNK_SIZE.min(n.size() - i_n_base);
-
-                {
-                    // load c into chunk
-                    for i_chunk_n in 0..n_chunk_size {
-                        let c_mn = unsafe { *c_m.add(c_strides[1] * (i_n_base + i_chunk_n)) };
-                        c_chunk[i_chunk_n] = c_mn.into();
-                    }
-                    for i_chunk_n in n_chunk_size..CHUNK_SIZE {
-                        c_chunk[i_chunk_n] = 0.0;
-                    }
-                }
-
-                for i_k_base in (0..k.size()).step_by(CHUNK_SIZE) {
-                    let k_chunk_size = CHUNK_SIZE.min(k.size() - i_k_base);
-
-                    {
-                        // load a & b into chunk
-                        for i_chunk_k in 0..k_chunk_size {
-                            let i_k = i_k_base + i_chunk_k;
-                            let a_mk = unsafe { *a_m.add(a_strides[1] * i_k) };
-                            let b_k = unsafe { bp.add(b_strides[0] * i_k) };
-
-                            // load a
-                            a_chunk[i_chunk_k] = a_mk.into();
-
-                            // load b
-                            for i_chunk_n in 0..n_chunk_size {
-                                let i_n = i_n_base + i_chunk_n;
-                                let b_kn = unsafe { *b_k.add(b_strides[1] * i_n) };
-                                b_chunk[i_chunk_k][i_chunk_n] = b_kn.into();
-                            }
-                            for i_chunk_n in n_chunk_size..CHUNK_SIZE {
-                                b_chunk[i_chunk_k][i_chunk_n] = 0.0;
-                            }
-                        }
-                        for i_chunk_k in k_chunk_size..CHUNK_SIZE {
-                            a_chunk[i_chunk_k] = 0.0;
-                            for i_chunk_n in n_chunk_size..CHUNK_SIZE {
-                                b_chunk[i_chunk_k][i_chunk_n] = 0.0;
-                            }
-                        }
-                    }
-
-                    // do the computation
-                    for i_n_chunk in 0..CHUNK_SIZE {
-                        for i_k_chunk in 0..CHUNK_SIZE {
-                            c_chunk[i_n_chunk] +=
-                                a_chunk[i_k_chunk] * b_chunk[i_k_chunk][i_n_chunk];
-                        }
-                    }
-                }
-
-                // store c chunk in memory
-                for (i_chunk, &c_f32) in c_chunk.iter().enumerate().take(n_chunk_size) {
-                    let c = half::f16::from_f32(c_f32);
-                    unsafe { *c_m.add(c_strides[1] * (i_n_base + i_chunk)) = c };
-                }
-            }
-        }
-    }
-}
diff --git a/src/tensor_ops/matmul/mod.rs b/src/tensor_ops/matmul/mod.rs
index 48943ed6d..d4f4c6977 100644
--- a/src/tensor_ops/matmul/mod.rs
+++ b/src/tensor_ops/matmul/mod.rs
@@ -1,8 +1,5 @@
 #![allow(clippy::type_complexity)]
 
-#[cfg(feature = "f16")]
-mod f16;
-
 pub(super) mod cpu_kernel;
 
 #[cfg(feature = "cuda")]


Test patch: