Improve the ring backend initialization

.circleci/config.yml

@@ -7,9 +7,18 @@ parameters:
   nightly_build:
     type: boolean
     default: false
+  weekly_build:
+    type: boolean
+    default: false
   test_release:
     type: boolean
     default: false
+  linux_release:
+    type: boolean
+    default: false
+  cuda_release:
+    type: boolean
+    default: false
 
 jobs:
   build_documentation:
@@ -64,9 +73,9 @@ jobs:
                  git push -f origin gh-pages
 
   linux_build_and_test:
-    machine:
-      image: ubuntu-2204:current
-      resource_class: large
+    docker:
+      - image: cimg/python:3.9
+
     steps:
       - checkout
       - run:
@@ -78,36 +87,34 @@ jobs:
       - run:
           name: Install dependencies
           command: |
-            export DEBIAN_FRONTEND=noninteractive
-            export NEEDRESTART_MODE=a
+            pip install --upgrade cmake
+            pip install nanobind==2.4.0
+            pip install numpy
             sudo apt-get update
-            sudo apt-get install -y libblas-dev liblapack-dev liblapacke-dev
+            sudo apt-get install libblas-dev liblapack-dev liblapacke-dev
             sudo apt-get install openmpi-bin openmpi-common libopenmpi-dev
-            curl -LsSf https://astral.sh/uv/install.sh | sh
       - run:
           name: Install Python package
           command: |
-            uv venv
-            uv pip install cmake
-            DEBUG=1 CMAKE_ARGS="-DCMAKE_COMPILE_WARNING_AS_ERROR=ON" \
-              uv pip install -e ".[dev]" -v
+            CMAKE_ARGS="-DMLX_BUILD_METAL=OFF" \
+              python3 setup.py build_ext --inplace
+            CMAKE_ARGS="-DMLX_BUILD_METAL=OFF" \
+              python3 setup.py develop
       - run:
           name: Generate package stubs
           command: |
-            uv pip install typing_extensions
-            uv run --no-project setup.py generate_stubs
+            echo "stubs"
+            pip install typing_extensions
+            python setup.py generate_stubs
       - run:
           name: Run Python tests
           command: |
-            source .venv/bin/activate
-            python -m unittest discover python/tests -v
+            python3 -m unittest discover python/tests -v
             mpirun --bind-to none -host localhost:8 -np 8 python python/tests/mpi_test_distributed.py
-            mlx.launch --verbose -n 8 python/tests/ring_test_distributed.py -v 2> >(tee -a stderr.log >&2)
-            if $(grep "\[WARN\]" stderr.log); then echo "Distributed ring test failed"; exit 1; fi
+            mlx.launch --verbose -n 8 python/tests/ring_test_distributed.py
       - run:
           name: Build CPP only
           command: |
-            source .venv/bin/activate
             mkdir -p build && cd build 
             cmake .. -DMLX_BUILD_METAL=OFF -DCMAKE_BUILD_TYPE=DEBUG
             make -j `nproc`
@@ -133,49 +140,50 @@ jobs:
       - run:
           name: Install dependencies
           command: |
-            HOMEBREW_NO_AUTO_UPDATE=1 HOMEBREW_NO_INSTALL_CLEANUP=1 \
-              brew install openmpi uv
+            brew install python@3.9
+            brew install openmpi
+            python3.9 -m venv env
+            source env/bin/activate
+            pip install --upgrade pip
+            pip install --upgrade cmake
+            pip install nanobind==2.4.0
+            pip install numpy
+            pip install torch
+            pip install tensorflow
+            pip install unittest-xml-reporting
       - run:
           name: Install Python package
           command: |
-            uv venv --python 3.9
-            uv pip install \
-              nanobind==2.4.0 \
-              cmake \
-              numpy \
-              torch \
-              tensorflow \
-              unittest-xml-reporting
+            source env/bin/activate
             DEBUG=1 CMAKE_ARGS="-DCMAKE_COMPILE_WARNING_AS_ERROR=ON" \
-              uv pip install -e . -v
+              pip install -e . -v
       - run:
           name: Generate package stubs
           command: |
-            uv pip install typing_extensions
-            uv run --no-project setup.py generate_stubs
+            source env/bin/activate
+            pip install typing_extensions
+            python setup.py generate_stubs
       - run:
           name: Run Python tests
           command: |
-            source .venv/bin/activate
+            source env/bin/activate
             LOW_MEMORY=1 DEVICE=cpu python -m xmlrunner discover -v python/tests -o test-results/cpu
             LOW_MEMORY=1 DEVICE=gpu METAL_DEVICE_WRAPPER_TYPE=1 METAL_DEBUG_ERROR_MODE=0 python -m xmlrunner discover -v python/tests -o test-results/gpu
             mpirun --bind-to none -host localhost:8 -np 8 -x DYLD_LIBRARY_PATH=/opt/homebrew/lib/ python python/tests/mpi_test_distributed.py
-            mlx.launch --verbose -n 8 python/tests/ring_test_distributed.py -v 2> >(tee -a stderr.log >&2)
-            if $(grep "\[WARN\]" stderr.log); then echo "Distributed ring test failed"; exit 1; fi
+            mlx.launch --verbose -n 8 python/tests/ring_test_distributed.py
       - run:
           name: Build example extension
           command: |
-            source .venv/bin/activate
+            source env/bin/activate
             cd examples/extensions
-            uv pip install -r requirements.txt
-            uv run --no-project setup.py build_ext --inplace
-            uv run --no-project python test.py
+            pip install -r requirements.txt
+            python setup.py build_ext -j8
       - store_test_results:
           path: test-results
       - run:
           name: Build CPP only
           command: |
-            source .venv/bin/activate
+            source env/bin/activate
             mkdir -p build && cd build && cmake .. && make -j `sysctl -n hw.ncpu`
       - run:
           name: Run CPP tests
@@ -184,7 +192,7 @@ jobs:
       - run:
           name: Build small binary
           command: |
-            source .venv/bin/activate
+            source env/bin/activate
             cd build/
             cmake .. -DCMAKE_BUILD_TYPE=MinSizeRel \
               -DBUILD_SHARED_LIBS=ON \
@@ -196,61 +204,34 @@ jobs:
       - run:
           name: Run Python tests with JIT
           command: |
+            source env/bin/activate
             CMAKE_ARGS="-DMLX_METAL_JIT=ON" \
-              uv pip install -e .
+              pip install -e . -v
             LOW_MEMORY=1 DEVICE=gpu METAL_DEVICE_WRAPPER_TYPE=1 \
               METAL_DEBUG_ERROR_MODE=0 \
-              uv run --no-project python -m xmlrunner discover \
-                -v python/tests \
-                -o test-results/gpu_jit
+              python -m xmlrunner discover -v python/tests -o test-results/gpu_jit
 
   cuda_build_and_test:
-    parameters:
-      image_date:
-        type: string
-        default: "2023.11.1"
     machine:
-      image: "linux-cuda-12:<< parameters.image_date >>"
+      image: linux-cuda-12:default
       resource_class: gpu.nvidia.small.gen2
     steps:
       - checkout
-      - restore_cache:
-          keys:
-            - cuda-<< parameters.image_date >>-{{ arch }}-
-      - run:
-          name: Install dependencies
-          command: |
-            sudo apt-get update
-            sudo apt-get install libcudnn9-dev-cuda-12
-            sudo apt-get install libblas-dev liblapack-dev liblapacke-dev
-            sudo apt-get install libnccl2 libnccl-dev
-            curl -sL https://github.com/ccache/ccache/releases/download/v4.11.3/ccache-4.11.3-linux-x86_64.tar.xz | tar xJf -
-            sudo mv ccache-4.11.3-linux-x86_64/ccache /usr/bin/ccache
-            rm -rf ccache-4.11.3-linux-x86_64
-            curl -LsSf https://astral.sh/uv/install.sh | sh
       - run:
           name: Install Python package
           command: |
-            uv venv
-            DEBUG=1 CMAKE_ARGS="-DMLX_BUILD_CUDA=ON -DCMAKE_COMPILE_WARNING_AS_ERROR=ON -DCMAKE_CUDA_COMPILER=`which nvcc`" \
-              uv pip install -e ".[dev]" -v
+            sudo apt-get update
+            sudo apt-get install libblas-dev liblapack-dev liblapacke-dev
+            python -m venv env
+            source env/bin/activate
+            CMAKE_ARGS="-DMLX_BUILD_CUDA=ON -DCMAKE_CUDA_COMPILER=`which nvcc`" \
+              pip install -e ".[dev]"
       - run:
           name: Run Python tests
           command: |
-            source .venv/bin/activate
+            source env/bin/activate
             LOW_MEMORY=1 DEVICE=cpu python -m unittest discover python/tests -v
             LOW_MEMORY=1 DEVICE=gpu python -m tests discover python/tests -v
-      - run:
-          name: CCache report
-          command: |
-            ccache --show-stats
-            ccache --zero-stats
-            ccache --max-size 400MB
-            ccache --cleanup
-      - save_cache:
-          key: cuda-<< parameters.image_date >>-{{ arch }}-{{ epoch }}
-          paths:
-            - /home/circleci/.cache/ccache
 
   build_release:
     parameters:
@@ -303,18 +284,7 @@ jobs:
           name: Build Python package
           command: |
             source env/bin/activate
-            python setup.py clean --all
-            << parameters.build_env >> MLX_BUILD_STAGE=1 python -m build -w
-      - when:
-          condition:
-            equal: ["3.9", << parameters.python_version >>]
-          steps:
-            - run:
-                name: Build common package
-                command: |
-                  source env/bin/activate
-                  python setup.py clean --all
-                  << parameters.build_env >> MLX_BUILD_STAGE=2 python -m build -w
+            << parameters.build_env >> python -m build -w
       - when:
           condition: << parameters.build_env >>
           steps:
@@ -331,100 +301,87 @@ jobs:
       python_version:
         type: string
         default: "3.9"
-      build_env:
+      extra_env:
         type: string
-        default: ""
-    machine:
-      image: ubuntu-2204:current
-      resource_class: large
+        default: "DEV_RELEASE=1"
+    docker:
+      - image: ubuntu:20.04
     steps:
       - checkout
       - run:
           name: Build wheel
           command: |
             PYTHON=python<< parameters.python_version >>
-            export DEBIAN_FRONTEND=noninteractive
-            export NEEDRESTART_MODE=a
-            sudo apt-get update
-            TZ=Etc/UTC sudo apt-get -y install tzdata
-            sudo add-apt-repository -y ppa:deadsnakes/ppa
-            sudo apt-get install -y $PYTHON $PYTHON-dev $PYTHON-full
-            sudo apt-get install -y libblas-dev liblapack-dev liblapacke-dev
+            apt-get update
+            apt-get upgrade -y
+            DEBIAN_FRONTEND=noninteractive TZ=Etc/UTC apt-get -y install tzdata
+            apt-get install -y apt-utils
+            apt-get install -y software-properties-common
+            add-apt-repository -y ppa:deadsnakes/ppa
+            apt-get install -y $PYTHON $PYTHON-dev $PYTHON-full
+            apt-get install -y libblas-dev liblapack-dev liblapacke-dev
+            apt-get install -y build-essential git
             $PYTHON -m venv env
             source env/bin/activate
             pip install --upgrade pip
             pip install --upgrade cmake
+            pip install nanobind==2.4.0
+            pip install --upgrade setuptools
+            pip install numpy
             pip install auditwheel
             pip install patchelf
             pip install build
             pip install twine
-            << parameters.build_env >> pip install ".[dev]" -v
+            << parameters.extra_env >> pip install . -v
             pip install typing_extensions
             python setup.py generate_stubs
-            python setup.py clean --all
-            MLX_BUILD_STAGE=1 << parameters.build_env >> python -m build -w
-            bash python/scripts/repair_linux.sh
-      - when:
-          condition:
-            equal: ["3.9", << parameters.python_version >>]
-          steps:
+            << parameters.extra_env >> python -m build --wheel
+            auditwheel show dist/*
+            auditwheel repair dist/* --plat manylinux_2_31_x86_64
       - run:
-                name: Build common package
+          name: Upload package
           command: |
             source env/bin/activate
-                  python setup.py clean --all
-                  << parameters.build_env >> MLX_BUILD_STAGE=2 \
-                    python -m build -w
-                  auditwheel repair dist/mlx_cpu*.whl --plat manylinux_2_35_x86_64
-      - when:
-          condition: << parameters.build_env >>
-          steps:
-            - run:
-                name: Upload packages
-                command: |
-                  source env/bin/activate
-                  twine upload wheelhouse/*.whl
+            twine upload wheelhouse/*
       - store_artifacts:
           path: wheelhouse/
 
   build_cuda_release:
     parameters:
-      build_env:
+      python_version:
         type: string
-        default: ""
+        default: "3.9"
+      extra_env:
+        type: string
+        default: "DEV_RELEASE=1"
     machine:
-      image: ubuntu-2204:current
-      resource_class: large
+      image: linux-cuda-12:default
+      resource_class: gpu.nvidia.small.gen2
     steps:
       - checkout
       - run:
           name: Build wheel
           command: |
-            export DEBIAN_FRONTEND=noninteractive
-            export NEEDRESTART_MODE=a
-            wget https://developer.download.nvidia.com/compute/cuda/repos/ubuntu2404/x86_64/cuda-keyring_1.1-1_all.deb
-            sudo dpkg -i cuda-keyring_1.1-1_all.deb
             sudo apt-get update
-            sudo apt-get install cuda-toolkit-12-9 libcudnn9-dev-cuda-12
-            sudo apt-get install libnccl2 libnccl-dev
             sudo apt-get install libblas-dev liblapack-dev liblapacke-dev
-            sudo apt-get install zip
+            python -m venv env
+            source env/bin/activate
             pip install auditwheel
             pip install patchelf
             pip install build
             pip install twine
-            export PATH=/usr/local/cuda/bin${PATH:+:${PATH}}
-            export LD_LIBRARY_PATH=/usr/local/cuda/lib64${LD_LIBRARY_PATH:+:${LD_LIBRARY_PATH}}
-            << parameters.build_env >> MLX_BUILD_STAGE=2 \
+            << parameters.extra_env >> \
               CMAKE_ARGS="-DMLX_BUILD_CUDA=ON -DCMAKE_CUDA_COMPILER=`which nvcc`" \
-              python -m build -w
+              pip install ".[dev]" -v
+            python setup.py generate_stubs
+            << parameters.extra_env >> \
+              CMAKE_ARGS="-DMLX_BUILD_CUDA=ON -DCMAKE_CUDA_COMPILER=`which nvcc`" \
+              python -m build --wheel
             bash python/scripts/repair_cuda.sh
-      - when:
-          condition: << parameters.build_env >>
-          steps:
       - run:
           name: Upload package
           command: |
+            source env/bin/activate
             twine upload wheelhouse/*.whl
       - store_artifacts:
           path: wheelhouse/
@@ -437,6 +394,7 @@ workflows:
             pattern: "^(?!pull/)[-\\w]+$"
             value: << pipeline.git.branch >>
         - not: << pipeline.parameters.nightly_build >>
+        - not: << pipeline.parameters.weekly_build >>
         - not: << pipeline.parameters.test_release >>
     jobs:
       - mac_build_and_test:
@@ -444,16 +402,14 @@ workflows:
             parameters:
               macosx_deployment_target: ["13.5", "14.0"]
       - linux_build_and_test
-      - cuda_build_and_test:
-          matrix:
-            parameters:
-              image_date: ["2023.11.1", "2025.05.1"]
+      - cuda_build_and_test 
       - build_documentation 
 
   build_pypi_release:
     when:
       and:
         - not: << pipeline.parameters.nightly_build >>
+        - not: << pipeline.parameters.weekly_build >>
         - not: << pipeline.parameters.test_release >>
     jobs:
       - build_release:
@@ -545,16 +501,7 @@ workflows:
           matrix:
             parameters:
               python_version: ["3.9", "3.10", "3.11", "3.12", "3.13"]
-              build_env: ["PYPI_RELEASE=1"]
-      - build_cuda_release:
-          filters:
-            tags:
-              only: /^v.*/
-            branches:
-              ignore: /.*/
-          matrix:
-            parameters:
-              build_env: ["PYPI_RELEASE=1"]
+              extra_env: ["PYPI_RELEASE=1"]
 
   prb:
     when:
@@ -575,9 +522,6 @@ workflows:
           requires: [ hold ]
       - cuda_build_and_test:
           requires: [ hold ]
-          matrix:
-            parameters:
-              image_date: ["2023.11.1", "2025.05.1"]
   nightly_build:
     when:
       and:
@@ -636,17 +580,11 @@ workflows:
               - macosx_deployment_target: "15.0"
                 xcode_version: "15.0.0"
                 python_version: "3.13"
-      - build_linux_release:
-          matrix:
-            parameters:
-              python_version: ["3.9", "3.10", "3.11", "3.12", "3.13"]
-      - build_cuda_release
-
-  build_dev_release:
+  weekly_build:
     when:
       and:
         - equal: [ main, << pipeline.git.branch >> ]
-        - << pipeline.parameters.test_release >>
+        - << pipeline.parameters.weekly_build >>
     jobs:
       - build_release:
           matrix:
@@ -716,12 +654,25 @@ workflows:
                 xcode_version: "15.0.0"
                 python_version: "3.13"
                 build_env: "DEV_RELEASE=1"
+  linux_test_release:
+    when:
+      and:
+        - equal: [ main, << pipeline.git.branch >> ]
+        - << pipeline.parameters.linux_release >>
+    jobs:
       - build_linux_release:
           matrix:
             parameters:
               python_version: ["3.9", "3.10", "3.11", "3.12", "3.13"]
-              build_env: ["DEV_RELEASE=1"]
+              extra_env: ["PYPI_RELEASE=1"]
+  cuda_test_release:
+    when:
+      and:
+        - equal: [ main, << pipeline.git.branch >> ]
+        - << pipeline.parameters.cuda_release >>
+    jobs:
       - build_cuda_release:
           matrix:
             parameters:
-              build_env: ["DEV_RELEASE=1"]
+              python_version: ["3.9", "3.10", "3.11", "3.12", "3.13"]
+              extra_env: ["PYPI_RELEASE=1"]

ACKNOWLEDGMENTS.md

@@ -19,17 +19,11 @@ MLX was developed with contributions from the following individuals:
 - Gleb Pobudzey: Added the `where` primitive, and groups in 1D and 2D convolutions.
 - Paul Paczuski: Improved stability of BCE loss calculation
 - Max-Heinrich Laves: Added `conv_transpose1d`, `conv_transpose2d`, and `conv_transpose3d` ops.
-- Gökdeniz Gülmez: Added the `Muon (MomentUm Orthogonalized by Newton-schulz)` optimizer.
 
 <a href="https://github.com/ml-explore/mlx/graphs/contributors">
   <img class="dark-light" src="https://contrib.rocks/image?repo=ml-explore/mlx&anon=0&columns=20&max=100&r=true" />
 </a>
 
-# Organizations
-
-MLX has received contributions from the following companies:
-- NVIDIA Corporation & Affiliates
-
 # Third-Party Software
 
 MLX leverages several third-party software, listed here together with

CMakeLists.txt

@@ -41,9 +41,7 @@ option(MLX_BUILD_GGUF "Include support for GGUF format" ON)
 option(MLX_BUILD_SAFETENSORS "Include support for safetensors format" ON)
 option(MLX_BUILD_BLAS_FROM_SOURCE "Build OpenBLAS from source code" OFF)
 option(MLX_METAL_JIT "Use JIT compilation for Metal kernels" OFF)
-option(MLX_USE_CCACHE "Use CCache for compilation cache when available" ON)
 option(BUILD_SHARED_LIBS "Build mlx as a shared library" OFF)
-option(USE_SYSTEM_FMT "Use system's provided fmt library" OFF)
 
 # --------------------- Processor tests -------------------------
 message(
@@ -66,17 +64,10 @@ if(${CMAKE_SYSTEM_NAME} MATCHES "Darwin")
       message(WARNING "Building for x86_64 arch is not officially supported.")
     endif()
   endif()
+
 else()
   set(MLX_BUILD_METAL OFF)
-endif()
-
-if(MLX_USE_CCACHE)
-  find_program(CCACHE_PROGRAM ccache)
-  if(CCACHE_PROGRAM)
-    set(CMAKE_C_COMPILER_LAUNCHER "${CCACHE_PROGRAM}")
-    set(CMAKE_CXX_COMPILER_LAUNCHER "${CCACHE_PROGRAM}")
-    set(CMAKE_CUDA_COMPILER_LAUNCHER "${CCACHE_PROGRAM}")
-  endif()
+  message(WARNING "MLX is prioritised for Apple silicon systems using macOS.")
 endif()
 
 # ----------------------------- Lib -----------------------------
@@ -243,16 +234,12 @@ target_include_directories(
 # Do not add mlx_EXPORTS define for shared library.
 set_target_properties(mlx PROPERTIES DEFINE_SYMBOL "")
 
-if(USE_SYSTEM_FMT)
-  find_package(fmt REQUIRED)
-else()
-  FetchContent_Declare(
+FetchContent_Declare(
   fmt
   GIT_REPOSITORY https://github.com/fmtlib/fmt.git
   GIT_TAG 10.2.1
   EXCLUDE_FROM_ALL)
-  FetchContent_MakeAvailable(fmt)
-endif()
+FetchContent_MakeAvailable(fmt)
 target_link_libraries(mlx PRIVATE $<BUILD_INTERFACE:fmt::fmt-header-only>)
 
 if(MLX_BUILD_PYTHON_BINDINGS)

README.md

@@ -68,23 +68,18 @@ in the documentation.
 
 ## Installation
 
-MLX is available on [PyPI](https://pypi.org/project/mlx/). To install MLX on
-macOS, run:
+MLX is available on [PyPI](https://pypi.org/project/mlx/). To install the Python API, run:
 
-```bash
+**With `pip`**:
+
+```
 pip install mlx
 ```
 
-To install the CUDA backend on Linux, run:
+**With `conda`**:
 
-```bash
-pip install mlx[cuda]
 ```
-
-To install a CPU-only Linux package, run:
-
-```bash
-pip install mlx[cpu]
+conda install -c conda-forge mlx
 ```
 
 Checkout the

cmake/FindNCCL.cmake

@@ -1,54 +0,0 @@
-# FindNCCL.cmake This module finds the NVIDIA NCCL library and its include
-# directories.
-
-set(NCCL_ROOT_DIR
-    $ENV{NCCL_ROOT_DIR}
-    CACHE PATH "Folder contains NVIDIA NCCL")
-
-find_path(
-  NCCL_INCLUDE_DIRS
-  NAMES nccl.h
-  HINTS ${NCCL_INCLUDE_DIR} ${NCCL_ROOT_DIR} ${NCCL_ROOT_DIR}/include
-        ${CUDA_TOOLKIT_ROOT_DIR}/include)
-
-if($ENV{USE_STATIC_NCCL})
-  message(
-    STATUS "USE_STATIC_NCCL detected. Linking against static NCCL library")
-  set(NCCL_LIBNAME "libnccl_static.a")
-else()
-  set(NCCL_LIBNAME "nccl")
-endif()
-
-find_library(
-  NCCL_LIBRARIES
-  NAMES ${NCCL_LIBNAME}
-  HINTS ${NCCL_LIB_DIR}
-        ${NCCL_ROOT_DIR}
-        ${NCCL_ROOT_DIR}/lib
-        ${NCCL_ROOT_DIR}/lib/x86_64-linux-gnu
-        ${NCCL_ROOT_DIR}/lib64
-        ${CUDA_TOOLKIT_ROOT_DIR}/lib
-        ${CUDA_TOOLKIT_ROOT_DIR}/lib64)
-
-include(FindPackageHandleStandardArgs)
-find_package_handle_standard_args(NCCL DEFAULT_MSG NCCL_INCLUDE_DIRS
-                                  NCCL_LIBRARIES)
-
-if(NCCL_FOUND)
-  set(NCCL_HEADER_FILE "${NCCL_INCLUDE_DIRS}/nccl.h")
-  message(
-    STATUS "Determining NCCL version from the header file: ${NCCL_HEADER_FILE}")
-  file(
-    STRINGS ${NCCL_HEADER_FILE} NCCL_MAJOR_VERSION_DEFINED
-    REGEX "^[ \t]*#define[ \t]+NCCL_MAJOR[ \t]+[0-9]+.*$"
-    LIMIT_COUNT 1)
-  if(NCCL_MAJOR_VERSION_DEFINED)
-    string(REGEX REPLACE "^[ \t]*#define[ \t]+NCCL_MAJOR[ \t]+" ""
-                         NCCL_MAJOR_VERSION ${NCCL_MAJOR_VERSION_DEFINED})
-    message(STATUS "NCCL_MAJOR_VERSION: ${NCCL_MAJOR_VERSION}")
-  endif()
-  message(
-    STATUS
-      "Found NCCL (include: ${NCCL_INCLUDE_DIRS}, library: ${NCCL_LIBRARIES})")
-  mark_as_advanced(NCCL_ROOT_DIR NCCL_INCLUDE_DIRS NCCL_LIBRARIES)
-endif()

docs/requirements.txt

@@ -1,5 +1,4 @@
 sphinx
 breathe
 sphinx-book-theme
-sphinx-copybutton
 mlx

docs/src/conf.py

@@ -18,7 +18,6 @@ release = version
 # -- General configuration ---------------------------------------------------
 
 extensions = [
-    "sphinx_copybutton",
     "sphinx.ext.autodoc",
     "sphinx.ext.autosummary",
     "sphinx.ext.intersphinx",

docs/src/dev/custom_metal_kernels.rst

@@ -128,7 +128,6 @@ relying on a copy from ``ensure_row_contiguous``:
       input_names=["inp"],
       output_names=["out"],
       source=source
-      ensure_row_contiguous=False,
   )
 
   def exp_elementwise(a: mx.array):
@@ -139,6 +138,7 @@ relying on a copy from ``ensure_row_contiguous``:
           threadgroup=(256, 1, 1),
           output_shapes=[a.shape],
           output_dtypes=[a.dtype],
+          ensure_row_contiguous=False,
       )
       return outputs[0]
 

docs/src/dev/extensions.rst

@@ -138,13 +138,13 @@ more concrete:
         * representing the vectorized computation and the axis which
         * corresponds to the output vectorized dimension.
         */
-        std::pair<std::vector<array>, std::vector<int>> vmap(
+        virtual std::pair<std::vector<array>, std::vector<int>> vmap(
             const std::vector<array>& inputs,
             const std::vector<int>& axes) override;
 
-        /** The name of primitive. */
-        const char* name() const override {
-          return "Axpby";
+        /** Print the primitive. */
+        void print(std::ostream& os) override {
+            os << "Axpby";
         }
 
         /** Equivalence check **/
@@ -394,14 +394,14 @@ below.
         out.set_data(allocator::malloc(out.nbytes()));
 
         // Resolve name of kernel
-        std::stream kname;
-        kname = "axpby_general_" + type_to_name(out);
+        std::ostringstream kname;
+        kname << "axpby_" << "general_" << type_to_name(out);
 
         // Load the metal library
-        auto lib = d.get_library("mlx_ext", current_binary_dir());
+        auto lib = d.get_library("mlx_ext");
 
         // Make a kernel from this metal library
-        auto kernel = d.get_kernel(kname, lib);
+        auto kernel = d.get_kernel(kname.str(), lib);
 
         // Prepare to encode kernel
         auto& compute_encoder = d.get_command_encoder(s.index);

docs/src/index.rst

@@ -70,7 +70,6 @@ are the CPU and GPU.
    python/fft
    python/linalg
    python/metal
-   python/cuda
    python/memory_management
    python/nn
    python/optimizers

docs/src/install.rst

@@ -13,7 +13,7 @@ silicon computer is
 
     pip install mlx
 
-To install from PyPI your system must meet the following requirements:
+To install from PyPI you must meet the following requirements:
 
 - Using an M series chip (Apple silicon)
 - Using a native Python >= 3.9
@@ -23,39 +23,22 @@ To install from PyPI your system must meet the following requirements:
     MLX is only available on devices running macOS >= 13.5
     It is highly recommended to use macOS 14 (Sonoma)
 
+
+MLX is also available on conda-forge. To install MLX with conda do:
+
+.. code-block:: shell
+
+   conda install conda-forge::mlx
+
 CUDA
 ^^^^
 
-MLX has a CUDA backend which you can install with:
+MLX has a CUDA backend which you can use on any Linux platform with CUDA 12
+and SM 7.0 (Volta) and up. To install MLX with CUDA support, run:
 
 .. code-block:: shell
 
-    pip install mlx[cuda]
-
-To install the CUDA package from PyPi your system must meet the following
-requirements:
-
-- Nvidia architecture >= SM 7.0 (Volta)
-- Nvidia driver >= 550.54.14
-- CUDA toolkit >= 12.0
-- Linux distribution with glibc >= 2.35
-- Python >= 3.9
-
-
-CPU-only (Linux)
-^^^^^^^^^^^^^^^^
-
-For a CPU-only version of MLX that runs on Linux use:
-
-.. code-block:: shell
-
-    pip install mlx[cpu]
-
-To install the CPU-only package from PyPi your system must meet the following
-requirements:
-
-- Linux distribution with glibc >= 2.35
-- Python >= 3.9
+    pip install mlx-cuda
 
 
 Troubleshooting
@@ -271,7 +254,7 @@ and the CUDA toolkit. For example on Ubuntu, run the following:
    dpkg -i cuda-keyring_1.1-1_all.deb
    apt-get update -y
    apt-get -y install cuda-toolkit-12-9
-   apt-get install libblas-dev liblapack-dev liblapacke-dev libcudnn9-dev-cuda-12 -y
+   apt-get install libblas-dev liblapack-dev liblapacke-dev -y
 
 
 When building either the Python or C++ APIs make sure to pass the cmake flag

docs/src/python/cuda.rst

@@ -1,9 +0,0 @@
-CUDA
-=====
-
-.. currentmodule:: mlx.core.cuda
-
-.. autosummary::
-  :toctree: _autosummary
-
-  is_available

docs/src/python/fast.rst

@@ -13,4 +13,3 @@ Fast
   rope
   scaled_dot_product_attention
   metal_kernel
-  cuda_kernel

docs/src/python/optimizers.rst

@@ -51,14 +51,14 @@ the saved state. Here's a simple example:
    optimizer.update(model, grads)
 
    # Save the state
-   state = tree_flatten(optimizer.state, destination={})
-   mx.save_safetensors("optimizer.safetensors", state)
+   state = tree_flatten(optimizer.state)
+   mx.save_safetensors("optimizer.safetensors", dict(state))
 
    # Later on, for example when loading from a checkpoint,
    # recreate the optimizer and load the state
    optimizer = optim.Adam(learning_rate=1e-2)
 
-   state = tree_unflatten(mx.load("optimizer.safetensors"))
+   state = tree_unflatten(list(mx.load("optimizer.safetensors").items()))
    optimizer.state = state
 
 Note, not every optimizer configuation parameter is saved in the state. For

docs/src/python/optimizers/common_optimizers.rst

@@ -19,4 +19,3 @@ Common Optimizers
    Adamax
    Lion
    MultiOptimizer
-   Muon

docs/src/usage/compile.rst

@@ -225,7 +225,7 @@ In some cases returning updated state can be pretty inconvenient. Hence,
   def fun(x, y):
       z = x + y
       state.append(z)
-      return mx.exp(z)
+      return mx.exp(z), state
 
   fun(mx.array(1.0), mx.array(2.0))
   # Prints [array(3, dtype=float32)]

docs/src/usage/export.rst

@@ -151,7 +151,7 @@ parameters, pass them as inputs to the ``call`` wrapper:
      model.update(tree_unflatten(list(params.items())))
      return model(x)
  
-   params = tree_flatten(model.parameters(), destination={})
+   params = dict(tree_flatten(model.parameters()))
    mx.export_function("model.mlxfn", call, (mx.zeros(4),), params)
 
 

examples/extensions/axpby/axpby.cpp

@@ -1,6 +1,5 @@
 // Copyright © 2023-2025 Apple Inc.
 
-#include <dlfcn.h>
 #include <iostream>
 #include <sstream>
 
@@ -17,19 +16,6 @@
 
 namespace my_ext {
 
-// A helper function to find the location of the current binary on disk.
-// The Metal library ("mlx_ext.mtllib"), should be in the same directory.
-std::string current_binary_dir() {
-  static std::string binary_dir = []() {
-    Dl_info info;
-    if (!dladdr(reinterpret_cast<void*>(&current_binary_dir), &info)) {
-      throw std::runtime_error("Unable to get current binary dir.");
-    }
-    return std::filesystem::path(info.dli_fname).parent_path().string();
-  }();
-  return binary_dir;
-}
-
 ///////////////////////////////////////////////////////////////////////////////
 // Operation Implementation
 ///////////////////////////////////////////////////////////////////////////////
@@ -181,15 +167,16 @@ void Axpby::eval_gpu(
   }
 
   // Resolve name of kernel (corresponds to axpby.metal)
-  std::string kname = "axpby_";
-  kname += (contiguous_kernel ? "contiguous_" : "general_");
-  kname += type_to_name(out);
+  std::ostringstream kname;
+  kname << "axpby_";
+  kname << (contiguous_kernel ? "contiguous_" : "general_");
+  kname << type_to_name(out);
 
   // Load the metal library
-  auto lib = d.get_library("mlx_ext", current_binary_dir());
+  auto lib = d.get_library("mlx_ext");
 
   // Make a kernel from this metal library
-  auto kernel = d.get_kernel(kname, lib);
+  auto kernel = d.get_kernel(kname.str(), lib);
 
   // Prepare to encode kernel
   auto& compute_encoder = d.get_command_encoder(s.index);

examples/extensions/axpby/axpby.h

@@ -74,9 +74,9 @@ class Axpby : public mx::Primitive {
       const std::vector<mx::array>& inputs,
       const std::vector<int>& axes) override;
 
-  /** The name of primitive. */
-  const char* name() const override {
-    return "Axpby";
+  /** Print the primitive. */
+  void print(std::ostream& os) override {
+    os << "Axpby";
   }
 
   /** Equivalence check **/

examples/extensions/requirements.txt

@@ -1,4 +1,4 @@
 setuptools>=42
 cmake>=3.25
 mlx>=0.21.0
-nanobind==2.4.0
+nanobind==2.2.0

examples/extensions/test.py

@@ -3,10 +3,8 @@ from mlx_sample_extensions import axpby
 
 a = mx.ones((3, 4))
 b = mx.ones((3, 4))
-c_cpu = axpby(a, b, 4.0, 2.0, stream=mx.cpu)
-c_gpu = axpby(a, b, 4.0, 2.0, stream=mx.gpu)
+c = axpby(a, b, 4.0, 2.0, stream=mx.cpu)
 
-print(f"c shape: {c_cpu.shape}")
-print(f"c dtype: {c_cpu.dtype}")
-print(f"c_cpu correct: {mx.all(c_cpu == 6.0).item()}")
-print(f"c_gpu correct: {mx.all(c_gpu == 6.0).item()}")
+print(f"c shape: {c.shape}")
+print(f"c dtype: {c.dtype}")
+print(f"c correct: {mx.all(c == 6.0).item()}")

mlx/array.h

@@ -10,7 +10,6 @@
 #include "mlx/allocator.h"
 #include "mlx/dtype.h"
 #include "mlx/event.h"
-#include "mlx/small_vector.h"
 
 namespace mlx::core {
 
@@ -19,8 +18,8 @@ class Primitive;
 
 using Deleter = std::function<void(allocator::Buffer)>;
 using ShapeElem = int32_t;
-using Shape = SmallVector<ShapeElem>;
-using Strides = SmallVector<int64_t>;
+using Shape = std::vector<ShapeElem>;
+using Strides = std::vector<int64_t>;
 
 class array {
   /* An array is really a node in a graph. It contains a shared ArrayDesc

mlx/backend/common/utils.cpp

@@ -1,20 +1,14 @@
 // Copyright © 2023-2024 Apple Inc.
 
-#include <dlfcn.h>
-
 #include "mlx/backend/common/utils.h"
+#include "mlx/primitives.h"
 
 namespace mlx::core {
 
-std::filesystem::path current_binary_dir() {
-  static std::filesystem::path binary_dir = []() {
-    Dl_info info;
-    if (!dladdr(reinterpret_cast<void*>(&current_binary_dir), &info)) {
-      throw std::runtime_error("Unable to get current binary dir.");
-    }
-    return std::filesystem::path(info.dli_fname).parent_path();
-  }();
-  return binary_dir;
+std::string get_primitive_string(Primitive* primitive) {
+  std::ostringstream op_t;
+  primitive->print(op_t);
+  return op_t.str();
 }
 
 std::tuple<Shape, std::vector<Strides>> collapse_contiguous_dims(

mlx/backend/common/utils.h

@@ -2,7 +2,6 @@
 
 #pragma once
 
-#include <filesystem>
 #include <tuple>
 #include <vector>
 
@@ -10,8 +9,7 @@
 
 namespace mlx::core {
 
-// Return the directory that contains current shared library.
-std::filesystem::path current_binary_dir();
+std::string get_primitive_string(Primitive* primitive);
 
 inline int64_t
 elem_to_loc(int elem, const Shape& shape, const Strides& strides) {
@@ -197,7 +195,7 @@ void shared_buffer_reshape(
     array& out);
 
 template <typename T>
-inline SmallVector<T> remove_index(SmallVector<T> vec, size_t index) {
+inline std::vector<T> remove_index(std::vector<T> vec, size_t index) {
   vec.erase(std::next(vec.begin(), index));
   return vec;
 }

mlx/backend/cpu/cholesky.cpp

@@ -20,7 +20,7 @@ void cholesky_impl(const array& a, array& factor, bool upper, Stream stream) {
 
   // The decomposition is computed in place, so just copy the input to the
   // output.
-  copy_cpu(
+  copy(
       a,
       factor,
       a.flags().row_contiguous ? CopyType::Vector : CopyType::General,

mlx/backend/cpu/compiled.cpp

@@ -157,12 +157,10 @@ inline void build_kernel(
 #endif
 
   // Start the kernel
-  os << "void " << kernel_name
-     << "(int* shape, int64_t** strides, void** args) {" << std::endl;
+  os << "void " << kernel_name << "(void** args) {" << std::endl;
 
   // Add the input arguments
   int cnt = 0;
-  int strides_index = 1;
   for (size_t i = 0; i < inputs.size(); ++i) {
     // Skip constants from the input list
     if (is_constant(i)) {
@@ -177,8 +175,8 @@ inline void build_kernel(
        << "];" << std::endl;
     // Scalars and contiguous need no strides
     if (!is_scalar(x) && !contiguous) {
-      os << "  const int64_t* " << xname << "_strides = strides["
-         << strides_index++ << "];" << std::endl;
+      os << "  const size_t* " << xname << "_strides = (size_t*)args[" << cnt++
+         << "];" << std::endl;
     }
   }
 
@@ -188,8 +186,10 @@ inline void build_kernel(
     os << "  " << tstr << "* " << namer.get_name(x) << " = (" << tstr
        << "*)args[" << cnt++ << "];" << std::endl;
   }
-  // Add output size
-  if (contiguous) {
+  // Add output strides and shape to extract the indices.
+  if (!contiguous) {
+    os << "  const int* shape = (int*)args[" << cnt++ << "];" << std::endl;
+  } else {
     os << "  const size_t size = (size_t)args[" << cnt++ << "];" << std::endl;
   }
 
@@ -231,7 +231,7 @@ inline void build_kernel(
       os << "static_cast<" << get_type_string(x.dtype()) << ">(tmp_"
          << namer.get_name(x.inputs()[0]) << ");" << std::endl;
     } else {
-      os << x.primitive().name();
+      x.primitive().print(os);
       os << "()(";
       for (int i = 0; i < x.inputs().size() - 1; i++) {
         os << "tmp_" << namer.get_name(x.inputs()[i]) << ", ";
@@ -290,6 +290,7 @@ void Compiled::eval_cpu(
 
   // Collect function input arguments.
   std::vector<void*> args;
+  int strides_index = 1;
   for (size_t i = 0; i < inputs.size(); ++i) {
     if (is_constant_(i)) {
       continue;
@@ -297,6 +298,9 @@ void Compiled::eval_cpu(
     const auto& x = inputs[i];
     encoder.set_input_array(x);
     args.push_back((void*)x.data<void>());
+    if (!contiguous && !is_scalar(x)) {
+      args.push_back(strides[strides_index++].data());
+    }
   }
 
   // Get the kernel name from the lib
@@ -331,20 +335,16 @@ void Compiled::eval_cpu(
     args.push_back(x.data<void>());
     encoder.set_output_array(x);
   }
-  if (contiguous) {
+  if (!contiguous) {
+    args.push_back((void*)shape.data());
+  } else {
     args.push_back((void*)outputs[0].data_size());
   }
-  auto fun = reinterpret_cast<void (*)(int*, int64_t**, void**)>(fn_ptr);
+  auto fun = (void (*)(void**))fn_ptr;
   encoder.dispatch([fun,
                     args = std::move(args),
                     strides = std::move(strides),
-                    shape = std::move(shape)]() mutable {
-    SmallVector<int64_t*> strides_ptrs;
-    for (auto& s : strides) {
-      strides_ptrs.push_back(s.data());
-    }
-    fun(shape.data(), strides_ptrs.data(), args.data());
-  });
+                    shape = std::move(shape)]() mutable { fun(args.data()); });
 }
 
 } // namespace mlx::core

mlx/backend/cpu/conv.cpp

@@ -883,7 +883,7 @@ void explicit_gemm_conv_1D_cpu(
   // Fill with zeros
   std::vector<array> temps;
   temps.push_back(array(0, conv_dtype));
-  copy_cpu(temps.back(), in_padded, CopyType::Scalar, stream);
+  copy(temps.back(), in_padded, CopyType::Scalar, stream);
 
   // Pick input slice from padded
   size_t data_offset = padding_lo[0] * in_padded.strides()[1];
@@ -895,7 +895,7 @@ void explicit_gemm_conv_1D_cpu(
       in_padded_slice.size(),
       data_offset);
   // Copy input values into the slice
-  copy_cpu_inplace(in, in_padded_slice, CopyType::GeneralGeneral, stream);
+  copy_inplace(in, in_padded_slice, CopyType::GeneralGeneral, stream);
   temps.push_back(in_padded_slice);
 
   // Make strided view
@@ -920,7 +920,7 @@ void explicit_gemm_conv_1D_cpu(
   // Materialize strided view
   Shape strided_reshape = {N * oH, wH * C};
   array in_strided(strided_reshape, in_strided_view.dtype(), nullptr, {});
-  copy_cpu(in_strided_view, in_strided, CopyType::General, stream);
+  copy(in_strided_view, in_strided, CopyType::General, stream);
   temps.push_back(in_strided);
 
   // Check wt dtype and prepare
@@ -938,13 +938,13 @@ void explicit_gemm_conv_1D_cpu(
         wt.size(),
         0);
     gemm_wt = array(wt_transpose.shape(), float32, nullptr, {});
-    copy_cpu(wt_transpose, gemm_wt, CopyType::General, stream);
+    copy(wt_transpose, gemm_wt, CopyType::General, stream);
     temps.push_back(gemm_wt);
   } else if (wt.dtype() != float32 || !wt.flags().row_contiguous) {
     auto ctype =
         wt.flags().row_contiguous ? CopyType::Vector : CopyType::General;
     gemm_wt = array(wt.shape(), float32, nullptr, {});
-    copy_cpu(wt, gemm_wt, ctype, stream);
+    copy(wt, gemm_wt, ctype, stream);
     temps.push_back(gemm_wt);
   }
 
@@ -991,7 +991,7 @@ void explicit_gemm_conv_1D_cpu(
 
   // Copy results if needed
   if (out.dtype() != float32) {
-    copy_cpu_inplace(gemm_out, out, CopyType::Vector, stream);
+    copy_inplace(gemm_out, out, CopyType::Vector, stream);
   }
   encoder.add_temporaries(std::move(temps));
 }
@@ -1029,7 +1029,7 @@ void explicit_gemm_conv_2D_cpu(
   // Fill with zeros
   std::vector<array> temps;
   temps.push_back(array(0, conv_dtype));
-  copy_cpu(temps.back(), in_padded, CopyType::Scalar, stream);
+  copy(temps.back(), in_padded, CopyType::Scalar, stream);
 
   // Pick input slice from padded
   size_t data_offset = padding_lo[0] * in_padded.strides()[1] +
@@ -1044,7 +1044,7 @@ void explicit_gemm_conv_2D_cpu(
   temps.push_back(in_padded_slice);
 
   // Copy input values into the slice
-  copy_cpu_inplace(in, in_padded_slice, CopyType::GeneralGeneral, stream);
+  copy_inplace(in, in_padded_slice, CopyType::GeneralGeneral, stream);
 
   // Make strided view
   Shape strided_shape = {N, oH, oW, wH, wW, C};
@@ -1065,7 +1065,7 @@ void explicit_gemm_conv_2D_cpu(
   // Materialize strided view
   Shape strided_reshape = {N * oH * oW, wH * wW * C};
   array in_strided(strided_reshape, in_strided_view.dtype(), nullptr, {});
-  copy_cpu(in_strided_view, in_strided, CopyType::General, stream);
+  copy(in_strided_view, in_strided, CopyType::General, stream);
   temps.push_back(in_strided);
 
   // Check wt dtype and prepare
@@ -1076,7 +1076,7 @@ void explicit_gemm_conv_2D_cpu(
     auto ctype =
         wt.flags().row_contiguous ? CopyType::Vector : CopyType::General;
     gemm_wt = array(wt.shape(), float32, nullptr, {});
-    copy_cpu(wt, gemm_wt, ctype, stream);
+    copy(wt, gemm_wt, ctype, stream);
     temps.push_back(gemm_wt);
   }
 
@@ -1116,7 +1116,7 @@ void explicit_gemm_conv_2D_cpu(
 
   // Copy results if needed
   if (out.dtype() != float32) {
-    copy_cpu_inplace(gemm_out, out, CopyType::Vector, stream);
+    copy_inplace(gemm_out, out, CopyType::Vector, stream);
   }
   encoder.add_temporaries(std::move(temps));
 }
@@ -1156,7 +1156,7 @@ void explicit_gemm_conv_ND_cpu(
 
   // Fill with zeros
   std::vector<array> temps = {array(0, conv_dtype)};
-  copy_cpu(temps.back(), in_padded, CopyType::Scalar, stream);
+  copy(temps.back(), in_padded, CopyType::Scalar, stream);
 
   // Pick input slice from padded
   size_t data_offset = 0;
@@ -1173,7 +1173,7 @@ void explicit_gemm_conv_ND_cpu(
       data_offset);
 
   // Copy input values into the slice
-  copy_cpu_inplace(in, in_padded_slice, CopyType::GeneralGeneral, stream);
+  copy_inplace(in, in_padded_slice, CopyType::GeneralGeneral, stream);
   temps.push_back(in_padded_slice);
 
   // Make strided view
@@ -1212,7 +1212,7 @@ void explicit_gemm_conv_ND_cpu(
   }
 
   array in_strided(strided_reshape, in_strided_view.dtype(), nullptr, {});
-  copy_cpu(in_strided_view, in_strided, CopyType::General, stream);
+  copy(in_strided_view, in_strided, CopyType::General, stream);
   temps.push_back(in_strided);
 
   // Check wt dtype and prepare
@@ -1223,13 +1223,13 @@ void explicit_gemm_conv_ND_cpu(
     auto ctype =
         wt.flags().row_contiguous ? CopyType::Vector : CopyType::General;
     gemm_wt = array(wt.shape(), float32, nullptr, {});
-    copy_cpu(wt, gemm_wt, ctype, stream);
+    copy(wt, gemm_wt, ctype, stream);
     temps.push_back(gemm_wt);
   }
 
   if (flip) {
     auto gemm_wt_ = array(gemm_wt.shape(), float32, nullptr, {});
-    copy_cpu(gemm_wt, gemm_wt_, CopyType::Vector, stream);
+    copy(gemm_wt, gemm_wt_, CopyType::Vector, stream);
     temps.push_back(gemm_wt_);
 
     // Calculate the total size of the spatial dimensions
@@ -1284,7 +1284,7 @@ void explicit_gemm_conv_ND_cpu(
 
   // Copy results if needed
   if (out.dtype() != float32) {
-    copy_cpu_inplace(gemm_out, out, CopyType::Vector, stream);
+    copy_inplace(gemm_out, out, CopyType::Vector, stream);
   }
   encoder.add_temporaries(std::move(temps));
 }

mlx/backend/cpu/copy.cpp

@@ -295,11 +295,7 @@ inline void copy_inplace_dispatch(
 
 } // namespace
 
-void copy_cpu_inplace(
-    const array& src,
-    array& dst,
-    CopyType ctype,
-    Stream stream) {
+void copy_inplace(const array& src, array& dst, CopyType ctype, Stream stream) {
   auto& encoder = cpu::get_command_encoder(stream);
   encoder.set_input_array(src);
   encoder.set_output_array(dst);
@@ -309,7 +305,7 @@ void copy_cpu_inplace(
        ctype]() mutable { copy_inplace_dispatch(src, dst, ctype); });
 }
 
-void copy_cpu(const array& src, array& dst, CopyType ctype, Stream stream) {
+void copy(const array& src, array& dst, CopyType ctype, Stream stream) {
   bool donated = set_copy_output_data(src, dst, ctype);
   if (donated && src.dtype() == dst.dtype()) {
     // If the output has the same type as the input then there is nothing to
@@ -319,10 +315,10 @@ void copy_cpu(const array& src, array& dst, CopyType ctype, Stream stream) {
   if (ctype == CopyType::GeneralGeneral) {
     ctype = CopyType::General;
   }
-  copy_cpu_inplace(src, dst, ctype, stream);
+  copy_inplace(src, dst, ctype, stream);
 }
 
-void copy_cpu_inplace(
+void copy_inplace(
     const array& src,
     array& dst,
     const Shape& data_shape,
@@ -377,10 +373,4 @@ void copy_cpu_inplace(
       });
 }
 
-array contiguous_copy_cpu(const array& arr, Stream stream) {
-  array arr_copy(arr.shape(), arr.dtype(), nullptr, {});
-  copy_cpu(arr, arr_copy, CopyType::General, stream);
-  return arr_copy;
-}
-
 } // namespace mlx::core

mlx/backend/cpu/copy.h

@@ -10,14 +10,10 @@
 
 namespace mlx::core {
 
-void copy_cpu(const array& src, array& dst, CopyType ctype, Stream stream);
-void copy_cpu_inplace(
-    const array& src,
-    array& dst,
-    CopyType ctype,
-    Stream stream);
+void copy(const array& src, array& dst, CopyType ctype, Stream stream);
+void copy_inplace(const array& src, array& dst, CopyType ctype, Stream stream);
 
-void copy_cpu_inplace(
+void copy_inplace(
     const array& src,
     array& dst,
     const Shape& data_shape,
@@ -30,7 +26,4 @@ void copy_cpu_inplace(
     const std::optional<array>& dynamic_i_offset = std::nullopt,
     const std::optional<array>& dynamic_o_offset = std::nullopt);
 
-// Return a contiguous array with same shape that copies the data of |arr|.
-array contiguous_copy_cpu(const array& arr, Stream stream);
-
 } // namespace mlx::core

mlx/backend/cpu/distributed.cpp

@@ -13,7 +13,9 @@ std::pair<array, bool> ensure_row_contiguous(const array& arr, Stream stream) {
   if (arr.flags().row_contiguous) {
     return {arr, false};
   } else {
-    return {contiguous_copy_cpu(arr, stream), true};
+    array arr_copy(arr.shape(), arr.dtype(), nullptr, {});
+    copy(arr, arr_copy, CopyType::General, stream);
+    return {arr_copy, true};
   }
 };
 
@@ -32,7 +34,8 @@ void AllReduce::eval_cpu(
       }
       return in;
     } else {
-      array arr_copy = contiguous_copy_cpu(in, s);
+      array arr_copy(in.shape(), in.dtype(), nullptr, {});
+      copy(in, arr_copy, CopyType::General, s);
       out.copy_shared_buffer(arr_copy);
       return arr_copy;
     }

mlx/backend/cpu/eig.cpp

@@ -135,7 +135,7 @@ void Eig::eval_cpu(
       : array(a.shape(), complex64, nullptr, {});
 
   auto a_copy = array(a.shape(), a.dtype(), nullptr, {});
-  copy_cpu(
+  copy(
       a,
       a_copy,
       a.flags().row_contiguous ? CopyType::Vector : CopyType::General,

mlx/backend/cpu/eigh.cpp

@@ -196,7 +196,7 @@ void Eigh::eval_cpu(
 
   values.set_data(allocator::malloc(values.nbytes()));
 
-  copy_cpu(
+  copy(
       a,
       vectors,
       a.flags().row_contiguous ? CopyType::Vector : CopyType::General,

mlx/backend/cpu/hadamard.cpp

@@ -96,7 +96,7 @@ void Hadamard::eval_cpu(const std::vector<array>& inputs, array& out) {
   if (in.flags().row_contiguous && in.is_donatable()) {
     out.copy_shared_buffer(in);
   } else {
-    copy_cpu(
+    copy(
         in,
         out,
         in.flags().row_contiguous ? CopyType::Vector : CopyType::General,

mlx/backend/cpu/indexing.cpp

@@ -517,7 +517,7 @@ void Scatter::eval_cpu(const std::vector<array>& inputs, array& out) {
   // Copy src into out (copy allocates memory for out)
   auto ctype =
       src.flags().row_contiguous ? CopyType::Vector : CopyType::General;
-  copy_cpu(src, out, ctype, stream());
+  copy(src, out, ctype, stream());
 
   auto& encoder = cpu::get_command_encoder(stream());
   std::vector<array> inds;
@@ -686,7 +686,7 @@ void ScatterAxis::eval_cpu(const std::vector<array>& inputs, array& out) {
   // Copy src into out (copy allocates memory for out)
   auto ctype =
       src.flags().row_contiguous ? CopyType::Vector : CopyType::General;
-  copy_cpu(src, out, ctype, stream());
+  copy(src, out, ctype, stream());
 
   auto& encoder = cpu::get_command_encoder(stream());
   encoder.set_input_array(idx);

mlx/backend/cpu/inverse.cpp

@@ -115,7 +115,7 @@ void inverse_impl(
   //   (A⁻¹)ᵀ = (Aᵀ)⁻¹
 
   // The inverse is computed in place, so just copy the input to the output.
-  copy_cpu(
+  copy(
       a,
       inv,
       a.flags().row_contiguous ? CopyType::Vector : CopyType::General,

mlx/backend/cpu/jit_compiler.cpp

@@ -2,7 +2,6 @@
 
 #include "mlx/backend/cpu/jit_compiler.h"
 
-#include <algorithm>
 #include <sstream>
 #include <vector>
 

mlx/backend/cpu/lapack.h

@@ -47,7 +47,7 @@ INSTANTIATE_LAPACK_REAL(orgqr)
 INSTANTIATE_LAPACK_REAL(syevd)
 INSTANTIATE_LAPACK_REAL(geev)
 INSTANTIATE_LAPACK_REAL(potrf)
-INSTANTIATE_LAPACK_REAL(gesdd)
+INSTANTIATE_LAPACK_REAL(gesvdx)
 INSTANTIATE_LAPACK_REAL(getrf)
 INSTANTIATE_LAPACK_REAL(getri)
 INSTANTIATE_LAPACK_REAL(trtri)

mlx/backend/cpu/logsumexp.cpp

@@ -87,7 +87,8 @@ void LogSumExp::eval_cpu(const std::vector<array>& inputs, array& out) {
     if (x.flags().contiguous && x.strides()[x.ndim() - 1] == 1) {
       return x;
     } else {
-      array x_copy = contiguous_copy_cpu(x, s);
+      auto x_copy = array(x.shape(), x.dtype(), nullptr, {});
+      copy(x, x_copy, CopyType::General, s);
       encoder.add_temporary(x_copy);
       return x_copy;
     }

mlx/backend/cpu/luf.cpp

@@ -31,7 +31,7 @@ void luf_impl(
   strides[ndim - 1] = M;
   strides[ndim - 2] = 1;
   lu.set_data(allocator::malloc(lu.nbytes()), lu.nbytes(), strides, flags);
-  copy_cpu_inplace(
+  copy_inplace(
       a,
       lu,
       a.shape(),

mlx/backend/cpu/masked_mm.cpp

@@ -124,20 +124,21 @@ void BlockMaskedMM::eval_cpu(const std::vector<array>& inputs, array& out) {
         if (!expand_all && stx == arr.shape(-1) && sty == 1) {
           if (do_copy) {
             array arr_copy(arr.shape(), arr.dtype(), nullptr, {});
-            copy_cpu(arr, arr_copy, CopyType::Vector, s);
+            copy(arr, arr_copy, CopyType::Vector, s);
             return std::make_tuple(false, stx, arr_copy, true);
           }
           return std::make_tuple(false, stx, arr, false);
         } else if (!expand_all && stx == 1 && sty == arr.shape(-2)) {
           if (do_copy) {
             array arr_copy(arr.shape(), arr.dtype(), nullptr, {});
-            copy_cpu(arr, arr_copy, CopyType::Vector, s);
+            copy(arr, arr_copy, CopyType::Vector, s);
             return std::make_tuple(true, sty, arr_copy, true);
           }
           return std::make_tuple(true, sty, arr, false);
         } else {
+          array arr_copy(arr.shape(), arr.dtype(), nullptr, {});
+          copy(arr, arr_copy, CopyType::General, s);
           int64_t stx = arr.shape(-1);
-          array arr_copy = contiguous_copy_cpu(arr, s);
           return std::make_tuple(false, stx, arr_copy, true);
         }
       };
@@ -385,7 +386,7 @@ void GatherMM::eval_cpu(const std::vector<array>& inputs, array& out) {
       return std::make_tuple(true, sty, arr);
     } else {
       temps.push_back(array(arr.shape(), arr.dtype(), nullptr, {}));
-      copy_cpu(arr, temps.back(), CopyType::General, s);
+      copy(arr, temps.back(), CopyType::General, s);
       int64_t stx = arr.shape(-1);
       return std::make_tuple(false, stx, temps.back());
     }
@@ -503,7 +504,7 @@ void SegmentedMM::eval_cpu(const std::vector<array>& inputs, array& out) {
       return std::make_tuple(true, sty, x);
     } else {
       array xc(x.shape(), x.dtype(), nullptr, {});
-      copy_cpu(x, xc, CopyType::General, s);
+      copy(x, xc, CopyType::General, s);
       encoder.add_temporary(xc);
       int64_t stx = x.shape(-1);
       return std::make_tuple(false, stx, xc);

mlx/backend/cpu/matmul.cpp

@@ -81,7 +81,7 @@ void matmul_general(
       return std::make_tuple(true, sty, arr);
     } else {
       temps.push_back(array(arr.shape(), arr.dtype(), nullptr, {}));
-      copy_cpu(arr, temps.back(), CopyType::General, stream);
+      copy(arr, temps.back(), CopyType::General, stream);
       stx = arr.shape(-1);
       return std::make_tuple(false, stx, temps.back());
     }
@@ -142,7 +142,7 @@ void AddMM::eval_cpu(const std::vector<array>& inputs, array& out) {
   CopyType ctype = c.data_size() == 1
       ? CopyType::Scalar
       : (c.flags().row_contiguous ? CopyType::Vector : CopyType::General);
-  copy_cpu(c, out, ctype, stream());
+  copy(c, out, ctype, stream());
   if (inputs[0].shape(-1) == 0) {
     return;
   }

mlx/backend/cpu/primitives.cpp

@@ -22,7 +22,7 @@ void reshape(const array& in, array& out) {
   auto [copy_necessary, out_strides] = prepare_reshape(in, out);
   if (copy_necessary) {
     out.set_data(allocator::malloc(out.nbytes()));
-    copy_cpu_inplace(in, out, CopyType::General, out.primitive().stream());
+    copy_inplace(in, out, CopyType::General, out.primitive().stream());
   } else {
     shared_buffer_reshape(in, out_strides, out);
   }
@@ -175,7 +175,7 @@ void AsType::eval_cpu(const std::vector<array>& inputs, array& out) {
   assert(inputs.size() == 1);
   auto& in = inputs[0];
   CopyType ctype = in.flags().contiguous ? CopyType::Vector : CopyType::General;
-  copy_cpu(in, out, ctype, stream());
+  copy(in, out, ctype, stream());
 }
 
 void Concatenate::eval_cpu(const std::vector<array>& inputs, array& out) {
@@ -198,7 +198,7 @@ void Concatenate::eval_cpu(const std::vector<array>& inputs, array& out) {
     size_t data_offset = strides[axis_] * sizes[i];
     out_slice.copy_shared_buffer(
         out, strides, flags, out_slice.size(), data_offset);
-    copy_cpu_inplace(inputs[i], out_slice, CopyType::GeneralGeneral, stream());
+    copy_inplace(inputs[i], out_slice, CopyType::GeneralGeneral, stream());
   }
 }
 
@@ -211,7 +211,7 @@ void Contiguous::eval_cpu(const std::vector<array>& inputs, array& out) {
        (allow_col_major_ && in.flags().col_contiguous))) {
     out.copy_shared_buffer(in);
   } else {
-    copy_cpu(in, out, CopyType::General, stream());
+    copy(in, out, CopyType::General, stream());
   }
 }
 
@@ -235,7 +235,7 @@ void Full::eval_cpu(const std::vector<array>& inputs, array& out) {
   } else {
     ctype = CopyType::General;
   }
-  copy_cpu(in, out, ctype, stream());
+  copy(in, out, ctype, stream());
 }
 
 void Pad::eval_cpu(const std::vector<array>& inputs, array& out) {
@@ -251,7 +251,7 @@ void Pad::eval_cpu(const std::vector<array>& inputs, array& out) {
   assert(val.dtype() == in.dtype() && in.dtype() == out.dtype());
 
   // Fill output with val
-  copy_cpu(val, out, CopyType::Scalar, stream());
+  copy(val, out, CopyType::Scalar, stream());
 
   // Find offset for start of input values
   size_t data_offset = 0;
@@ -266,7 +266,7 @@ void Pad::eval_cpu(const std::vector<array>& inputs, array& out) {
       out, out.strides(), out.flags(), out_slice.size(), data_offset);
 
   // Copy input values into the slice
-  copy_cpu_inplace(in, out_slice, CopyType::GeneralGeneral, stream());
+  copy_inplace(in, out_slice, CopyType::GeneralGeneral, stream());
 }
 
 void RandomBits::eval_cpu(const std::vector<array>& inputs, array& out) {
@@ -340,7 +340,7 @@ void DynamicSlice::eval_cpu(const std::vector<array>& inputs, array& out) {
   out.set_data(allocator::malloc(out.nbytes()));
   auto [in_offset, donated] =
       compute_dynamic_offset(inputs[1], in.strides(), axes_, stream());
-  copy_cpu_inplace(
+  copy_inplace(
       /* const array& src = */ in,
       /* array& dst = */ out,
       /* const Shape& data_shape = */ out.shape(),
@@ -372,11 +372,11 @@ void DynamicSliceUpdate::eval_cpu(
   auto ctype = in.flags().contiguous && in.size() == in.data_size()
       ? CopyType::Vector
       : CopyType::General;
-  copy_cpu(in, out, in.data_size() == 1 ? CopyType::Scalar : ctype, stream());
+  copy(in, out, in.data_size() == 1 ? CopyType::Scalar : ctype, stream());
 
   auto [out_offset, donated] =
       compute_dynamic_offset(inputs[2], out.strides(), axes_, stream());
-  copy_cpu_inplace(
+  copy_inplace(
       /* const array& src = */ upd,
       /* array& dst = */ out,
       /* const std::vector<int>& data_shape = */ upd.shape(),
@@ -412,14 +412,14 @@ void SliceUpdate::eval_cpu(const std::vector<array>& inputs, array& out) {
   auto ctype = in.flags().contiguous && in.size() == in.data_size()
       ? CopyType::Vector
       : CopyType::General;
-  copy_cpu(in, out, in.data_size() == 1 ? CopyType::Scalar : ctype, stream());
+  copy(in, out, in.data_size() == 1 ? CopyType::Scalar : ctype, stream());
 
   // Calculate out strides, initial offset and if copy needs to be made
   auto [data_offset, out_strides] =
       prepare_slice(out, start_indices_, strides_);
 
   // Do copy
-  copy_cpu_inplace(
+  copy_inplace(
       /* const array& src = */ upd,
       /* array& dst = */ out,
       /* const std::vector<int>& data_shape = */ upd.shape(),
@@ -456,9 +456,9 @@ void View::eval_cpu(const std::vector<array>& inputs, array& out) {
     if (in.dtype() == bool_) {
       auto in_tmp = array(in.shape(), uint8, nullptr, {});
       in_tmp.copy_shared_buffer(in);
-      copy_cpu_inplace(in_tmp, tmp, CopyType::General, stream());
+      copy_inplace(in_tmp, tmp, CopyType::General, stream());
     } else {
-      copy_cpu_inplace(in, tmp, CopyType::General, stream());
+      copy_inplace(in, tmp, CopyType::General, stream());
     }
 
     auto flags = out.flags();

mlx/backend/cpu/qrf.cpp

@@ -26,7 +26,7 @@ void qrf_impl(const array& a, array& q, array& r, Stream stream) {
   strides[in.ndim() - 2] = 1;
   strides[in.ndim() - 1] = M;
   in.set_data(allocator::malloc(in.nbytes()), in.nbytes(), strides, flags);
-  copy_cpu_inplace(a, in, CopyType::GeneralGeneral, stream);
+  copy_inplace(a, in, CopyType::GeneralGeneral, stream);
   auto& encoder = cpu::get_command_encoder(stream);
   q.set_data(allocator::malloc(q.nbytes()));
   r.set_data(allocator::malloc(r.nbytes()));

mlx/backend/cpu/quantized.cpp

@@ -529,7 +529,7 @@ void QuantizedMatmul::eval_cpu(const std::vector<array>& inputs, array& out) {
       return arr;
     } else {
       temps.push_back(array(arr.shape(), arr.dtype(), nullptr, {}));
-      copy_cpu(arr, temps.back(), CopyType::General, s);
+      copy(arr, temps.back(), CopyType::General, s);
       return temps.back();
     }
   };
@@ -579,7 +579,7 @@ void GatherQMM::eval_cpu(const std::vector<array>& inputs, array& out) {
       return arr;
     } else {
       temps.push_back(array(arr.shape(), arr.dtype(), nullptr, {}));
-      copy_cpu(arr, temps.back(), CopyType::General, s);
+      copy(arr, temps.back(), CopyType::General, s);
       return temps.back();
     }
   };
@@ -712,7 +712,9 @@ void fast::AffineQuantize::eval_cpu(
     if (arr.flags().row_contiguous) {
       return std::make_pair(arr, false);
     } else {
-      return std::make_pair(contiguous_copy_cpu(arr, s), true);
+      array arr_copy(arr.shape(), arr.dtype(), nullptr, {});
+      copy(arr, arr_copy, CopyType::General, s);
+      return std::make_pair(arr_copy, true);
     }
   };
 

mlx/backend/cpu/reduce.cpp

@@ -491,27 +491,19 @@ void Reduce::eval_cpu(const std::vector<array>& inputs, array& out) {
         switch (in.dtype()) {
           case bool_:
           case uint8:
-            reduce_dispatch_sum_prod<uint8_t>(in, out, reduce_type_, axes_);
-            break;
-          case uint16:
-            reduce_dispatch_sum_prod<uint16_t>(in, out, reduce_type_, axes_);
-            break;
-          case uint32:
-            reduce_dispatch_sum_prod<uint32_t>(in, out, reduce_type_, axes_);
-            break;
-          case uint64:
-            reduce_dispatch_sum_prod<uint64_t>(in, out, reduce_type_, axes_);
-            break;
           case int8:
             reduce_dispatch_sum_prod<int8_t>(in, out, reduce_type_, axes_);
             break;
           case int16:
+          case uint16:
             reduce_dispatch_sum_prod<int16_t>(in, out, reduce_type_, axes_);
             break;
           case int32:
+          case uint32:
             reduce_dispatch_sum_prod<int32_t>(in, out, reduce_type_, axes_);
             break;
           case int64:
+          case uint64:
             reduce_dispatch_sum_prod<int64_t>(in, out, reduce_type_, axes_);
             break;
           case float16:

mlx/backend/cpu/scan.cpp

@@ -250,8 +250,10 @@ void Scan::eval_cpu(const std::vector<array>& inputs, array& out) {
   // Ensure contiguity
   auto in = inputs[0];
   if (!in.flags().row_contiguous) {
-    in = contiguous_copy_cpu(in, stream());
-    encoder.add_temporary(in);
+    array arr_copy(in.shape(), in.dtype(), nullptr, {});
+    copy(in, arr_copy, CopyType::General, stream());
+    in = arr_copy;
+    encoder.add_temporary(arr_copy);
   }
   out.set_data(allocator::malloc(out.nbytes()));
 

mlx/backend/cpu/simd/accelerate_simd.h

@@ -234,7 +234,6 @@ Simd<T, N> remainder(Simd<T, N> a, Simd<T, N> b) {
 
 template <typename MaskT, typename T1, typename T2, int N>
 Simd<T1, N> select(Simd<MaskT, N> mask, Simd<T1, N> x, Simd<T2, N> y) {
-  static_assert(std::is_same_v<MaskT, bool>);
   if constexpr (sizeof(T1) == 1) {
     return asd::bitselect(y.value, x.value, asd::convert<char>(mask.value));
   } else if constexpr (sizeof(T1) == 2) {
@@ -252,13 +251,9 @@ Simd<T, N> pow(Simd<T, N> base, Simd<T, N> exp) {
     return asd::pow(base.value, exp.value);
   } else {
     Simd<T, N> res = 1;
-    // Raising an integer to a negative power is undefined
-    if (any(exp < 0)) {
-      return 0;
-    }
-    while (any(exp > 0)) {
-      res = select((exp & 1) != 0, res * base, res);
-      base = select(exp > 0, base * base, base);
+    while (any(exp)) {
+      res = select(exp & 1, res * base, res);
+      base = select(exp, base * base, base);
       exp = exp >> 1;
     }
     return res;

mlx/backend/cpu/softmax.cpp

@@ -131,7 +131,8 @@ void Softmax::eval_cpu(const std::vector<array>& inputs, array& out) {
       }
       return x;
     } else {
-      array x_copy = contiguous_copy_cpu(x, s);
+      array x_copy(x.shape(), x.dtype(), nullptr, {});
+      copy(x, x_copy, CopyType::General, s);
       out.copy_shared_buffer(x_copy);
       return x_copy;
     }

mlx/backend/cpu/sort.cpp

@@ -8,7 +8,7 @@
 #include "mlx/backend/common/utils.h"
 #include "mlx/backend/cpu/copy.h"
 #include "mlx/backend/cpu/encoder.h"
-#include "mlx/dtype_utils.h"
+
 #include "mlx/primitives.h"
 
 namespace mlx::core {
@@ -333,23 +333,44 @@ void Sort::eval_cpu(const std::vector<array>& inputs, array& out) {
   assert(inputs.size() == 1);
   auto& in = inputs[0];
 
-  int axis = axis_;
-  if (axis < 0) {
-    axis += in.ndim();
-  }
-
   // Copy input to output
-  CopyType ctype = (in.flags().contiguous && in.strides()[axis] != 0)
-      ? CopyType::Vector
-      : CopyType::General;
-  copy_cpu(in, out, ctype, stream());
+  CopyType ctype = in.flags().contiguous ? CopyType::Vector : CopyType::General;
+  copy(in, out, ctype, stream());
 
   auto& encoder = cpu::get_command_encoder(stream());
   encoder.set_output_array(out);
-  encoder.dispatch([out = array::unsafe_weak_copy(out), axis]() mutable {
-    dispatch_all_types(out.dtype(), [&](auto type_tag) {
-      sort<MLX_GET_TYPE(type_tag)>(out, axis);
-    });
+  encoder.dispatch(
+      [out = array::unsafe_weak_copy(out), axis_ = axis_]() mutable {
+        switch (out.dtype()) {
+          case bool_:
+            return sort<bool>(out, axis_);
+          case uint8:
+            return sort<uint8_t>(out, axis_);
+          case uint16:
+            return sort<uint16_t>(out, axis_);
+          case uint32:
+            return sort<uint32_t>(out, axis_);
+          case uint64:
+            return sort<uint64_t>(out, axis_);
+          case int8:
+            return sort<int8_t>(out, axis_);
+          case int16:
+            return sort<int16_t>(out, axis_);
+          case int32:
+            return sort<int32_t>(out, axis_);
+          case int64:
+            return sort<int64_t>(out, axis_);
+          case float32:
+            return sort<float>(out, axis_);
+          case float64:
+            return sort<double>(out, axis_);
+          case float16:
+            return sort<float16_t>(out, axis_);
+          case bfloat16:
+            return sort<bfloat16_t>(out, axis_);
+          case complex64:
+            return sort<complex64_t>(out, axis_);
+        }
       });
 }
 
@@ -405,10 +426,8 @@ void Partition::eval_cpu(const std::vector<array>& inputs, array& out) {
   auto& in = inputs[0];
 
   // Copy input to output
-  CopyType ctype = (in.flags().contiguous && in.strides()[axis_] != 0)
-      ? CopyType::Vector
-      : CopyType::General;
-  copy_cpu(in, out, ctype, stream());
+  CopyType ctype = in.flags().contiguous ? CopyType::Vector : CopyType::General;
+  copy(in, out, ctype, stream());
 
   auto& encoder = cpu::get_command_encoder(stream());
   encoder.set_output_array(out);

mlx/backend/cpu/svd.cpp

@@ -31,7 +31,7 @@ void svd_impl(
 
   // lapack clobbers the input, so we have to make a copy.
   array in(a.shape(), a.dtype(), nullptr, {});
-  copy_cpu(
+  copy(
       a,
       in,
       a.flags().row_contiguous ? CopyType::Vector : CopyType::General,
@@ -81,7 +81,9 @@ void svd_impl(
     // Vᵀ of shape N x N. (M x M in lapack).
     const int ldvt = M;
 
-    auto jobz = (u_ptr) ? "A" : "N";
+    auto job_u = (u_ptr) ? "V" : "N";
+    auto job_vt = (u_ptr) ? "V" : "N";
+    static constexpr auto range = "A";
 
     // Will contain the number of singular values after the call has returned.
     int ns = 0;
@@ -89,20 +91,30 @@ void svd_impl(
 
     // Will contain the indices of eigenvectors that failed to converge (not
     // used here but required by lapack).
-    auto iwork = array::Data{allocator::malloc(sizeof(int) * 8 * K)};
+    auto iwork = array::Data{allocator::malloc(sizeof(int) * 12 * K)};
 
     static const int lwork_query = -1;
 
+    static const int ignored_int = 0;
+    static const T ignored_float = 0;
+
     int info;
 
     // Compute workspace size.
-    gesdd<T>(
-        /* jobz = */ jobz,
+    gesvdx<T>(
+        /* jobu = */ job_u,
+        /* jobvt = */ job_vt,
+        /* range = */ range,
         // M and N are swapped since lapack expects column-major.
         /* m = */ &N,
         /* n = */ &M,
         /* a = */ nullptr,
         /* lda = */ &lda,
+        /* vl = */ &ignored_float,
+        /* vu = */ &ignored_float,
+        /* il = */ &ignored_int,
+        /* iu = */ &ignored_int,
+        /* ns = */ &ns,
         /* s = */ nullptr,
         /* u = */ nullptr,
         /* ldu = */ &ldu,
@@ -124,13 +136,20 @@ void svd_impl(
 
     // Loop over matrices.
     for (int i = 0; i < num_matrices; i++) {
-      gesdd<T>(
-          /* jobz = */ jobz,
+      gesvdx<T>(
+          /* jobu = */ job_u,
+          /* jobvt = */ job_vt,
+          /* range = */ range,
           // M and N are swapped since lapack expects column-major.
           /* m = */ &N,
           /* n = */ &M,
           /* a = */ in_ptr + M * N * i,
           /* lda = */ &lda,
+          /* vl = */ &ignored_float,
+          /* vu = */ &ignored_float,
+          /* il = */ &ignored_int,
+          /* iu = */ &ignored_int,
+          /* ns = */ &ns,
           /* s = */ s_ptr + K * i,
           // According to the identity above, lapack will write Vᵀᵀ as U.
           /* u = */ vt_ptr ? vt_ptr + N * N * i : nullptr,
@@ -148,6 +167,13 @@ void svd_impl(
         ss << "svd_impl: sgesvdx_ failed with code " << info;
         throw std::runtime_error(ss.str());
       }
+
+      if (ns != K) {
+        std::stringstream ss;
+        ss << "svd_impl: expected " << K << " singular values, but " << ns
+           << " were computed.";
+        throw std::runtime_error(ss.str());
+      }
     }
   });
   encoder.add_temporary(in);

mlx/backend/cuda/CMakeLists.txt

@@ -6,8 +6,8 @@
 target_sources(
   mlx
   PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/allocator.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/arange.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/arg_reduce.cu
+          ${CMAKE_CURRENT_SOURCE_DIR}/binary.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/binary_two.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/compiled.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/copy.cu
@@ -15,26 +15,18 @@ target_sources(
           ${CMAKE_CURRENT_SOURCE_DIR}/copy/copy_general.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/copy/copy_general_dynamic.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/copy/copy_general_input.cu
-          ${CMAKE_CURRENT_SOURCE_DIR}/conv.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/conv/gemm_conv.cu
-          ${CMAKE_CURRENT_SOURCE_DIR}/conv/gemm_grouped_conv.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/cuda.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/cudnn_utils.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/custom_kernel.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/device.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/distributed.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/eval.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/event.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/fence.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/gemms/gemv.cu
-          ${CMAKE_CURRENT_SOURCE_DIR}/gemms/cublas_gemm.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/jit_module.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/indexing.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/kernel_utils.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/matmul.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/layer_norm.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/logsumexp.cu
-          ${CMAKE_CURRENT_SOURCE_DIR}/primitives.cpp
+          ${CMAKE_CURRENT_SOURCE_DIR}/primitives.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/random.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/reduce.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/reduce/all_reduce.cu
@@ -43,28 +35,15 @@ target_sources(
           ${CMAKE_CURRENT_SOURCE_DIR}/reduce/row_reduce.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/rms_norm.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/rope.cu
-          ${CMAKE_CURRENT_SOURCE_DIR}/scaled_dot_product_attention.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/scan.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/slicing.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/softmax.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/sort.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/ternary.cu
+          ${CMAKE_CURRENT_SOURCE_DIR}/unary.cu
           ${CMAKE_CURRENT_SOURCE_DIR}/utils.cpp
-          ${CMAKE_CURRENT_SOURCE_DIR}/quantized/affine_quantize.cu
-          ${CMAKE_CURRENT_SOURCE_DIR}/quantized/quantized.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/worker.cpp)
 
-add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/binary)
-add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/unary)
-
-if(CMAKE_CUDA_COMPILER_VERSION VERSION_GREATER_EQUAL 12.9.0)
-  target_sources(
-    mlx PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/gemms/cublas_gemm_batched_12_9.cu)
-else()
-  target_sources(
-    mlx PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/gemms/cublas_gemm_batched_12_0.cpp)
-endif()
-
 target_compile_definitions(mlx PRIVATE MLX_USE_CUDA)
 
 # Embed kernel sources in binary for JIT compilation.
@@ -107,18 +86,11 @@ endif()
 target_compile_options(
   mlx PRIVATE "$<$<COMPILE_LANGUAGE:CUDA>:--Wno-deprecated-gpu-targets>")
 
-# Use stronger binaries compression. This feature was introduced in CUDA 12.8
-# and requires drivers released after CUDA 12.4.
-if(CMAKE_CUDA_COMPILER_VERSION VERSION_GREATER_EQUAL 12.8.0)
-  target_compile_options(
-    mlx PRIVATE "$<$<COMPILE_LANGUAGE:CUDA>:--compress-mode=size>")
-endif()
-
-# Compute capability >= 7.0 is required for synchronization between CPU/GPU with
-# managed memory.
-if(NOT DEFINED MLX_CUDA_ARCHITECTURES)
-  set(MLX_CUDA_ARCHITECTURES "native")
-endif()
+# Compute capability 7 is required for synchronization between CPU/GPU with
+# managed memory. TODO: Add more architectures for potential performance gain.
+set(MLX_CUDA_ARCHITECTURES
+    "70;80"
+    CACHE STRING "CUDA architectures")
 message(STATUS "CUDA architectures: ${MLX_CUDA_ARCHITECTURES}")
 set_target_properties(mlx PROPERTIES CUDA_ARCHITECTURES
                                      "${MLX_CUDA_ARCHITECTURES}")
@@ -150,27 +122,6 @@ target_link_libraries(mlx PRIVATE CUDA::cublasLt)
 # Use NVRTC and driver APIs.
 target_link_libraries(mlx PRIVATE CUDA::nvrtc CUDA::cuda_driver)
 
-# Use the frontend APIs of cuDNN.
-FetchContent_Declare(
-  cudnn
-  GIT_REPOSITORY https://github.com/NVIDIA/cudnn-frontend.git
-  GIT_TAG v1.14.0
-  GIT_SHALLOW TRUE
-  EXCLUDE_FROM_ALL)
-set(CUDNN_FRONTEND_SKIP_JSON_LIB ON)
-set(CUDNN_FRONTEND_BUILD_SAMPLES OFF)
-set(CUDNN_FRONTEND_BUILD_TESTS OFF)
-set(CUDNN_FRONTEND_BUILD_PYTHON_BINDINGS OFF)
-FetchContent_MakeAvailable(cudnn)
-target_link_libraries(mlx PRIVATE cudnn_frontend)
-# Link with the actual cuDNN libraries.
-include(${cudnn_frontend_SOURCE_DIR}/cmake/cuDNN.cmake)
-target_link_libraries(mlx PRIVATE CUDNN::cudnn_all)
-
 # Suppress nvcc warnings on MLX headers.
 target_compile_options(mlx PRIVATE $<$<COMPILE_LANGUAGE:CUDA>:-Xcudafe
                                    --diag_suppress=997>)
-
-# Install CCCL headers for JIT.
-install(DIRECTORY ${cccl_SOURCE_DIR}/include/cuda
-        DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/cccl)

mlx/backend/cuda/allocator.cpp

@@ -2,6 +2,7 @@
 
 #include "mlx/backend/cuda/allocator.h"
 #include "mlx/backend/cuda/utils.h"
+#include "mlx/backend/cuda/worker.h"
 #include "mlx/utils.h"
 
 #include <cuda_runtime.h>
@@ -16,73 +17,14 @@ namespace cu {
 
 constexpr int page_size = 16384;
 
-// Any allocations smaller than this will try to use the small pool
-constexpr int small_block_size = 8;
-
-// The small pool size in bytes. This should be a multiple of the host page
-// size and small_block_size.
-constexpr int small_pool_size = 4 * page_size;
-
-SmallSizePool::SmallSizePool() {
-  auto num_blocks = small_pool_size / small_block_size;
-  buffer_ = new Block[num_blocks];
-
-  next_free_ = buffer_;
-
-  CHECK_CUDA_ERROR(cudaMallocManaged(&data_, small_pool_size));
-#if CUDART_VERSION >= 13000
-  cudaMemLocation loc;
-  loc.type = cudaMemLocationTypeDevice;
-  loc.id = 0;
-#else
-  int loc = 0;
-#endif // CUDART_VERSION >= 13000
-  CHECK_CUDA_ERROR(
-      cudaMemAdvise(data_, small_pool_size, cudaMemAdviseSetReadMostly, loc));
-
-  auto curr = next_free_;
-  for (size_t i = 1; i < num_blocks; ++i) {
-    curr->next = buffer_ + i;
-    curr = curr->next;
-  }
-  curr->next = nullptr;
-}
-
-SmallSizePool::~SmallSizePool() {
-  CHECK_CUDA_ERROR(cudaFree(data_));
-  delete[] buffer_;
-}
-
-CudaBuffer* SmallSizePool::malloc() {
-  if (next_free_ == nullptr) {
-    return nullptr;
-  }
-  Block* b = next_free_;
-  uint64_t i = next_free_ - buffer_;
-  next_free_ = next_free_->next;
-  b->buf.data = static_cast<char*>(data_) + i * small_block_size;
-  b->buf.size = small_block_size;
-  return &b->buf;
-}
-
-void SmallSizePool::free(CudaBuffer* buf) {
-  auto b = reinterpret_cast<Block*>(buf);
-  b->next = next_free_;
-  next_free_ = b;
-}
-
-bool SmallSizePool::in_pool(CudaBuffer* buf) {
-  constexpr int num_blocks = (small_pool_size / small_block_size);
-  auto b = reinterpret_cast<Block*>(buf);
-  int64_t block_num = b - buffer_;
-  return block_num >= 0 && block_num < num_blocks;
-}
-
 CudaAllocator::CudaAllocator()
     : buffer_cache_(
           page_size,
           [](CudaBuffer* buf) { return buf->size; },
-          [this](CudaBuffer* buf) { cuda_free(buf); }) {
+          [this](CudaBuffer* buf) {
+            cuda_free(buf->data);
+            delete buf;
+          }) {
   // TODO: Set memory limit for multi-device.
   size_t free, total;
   CHECK_CUDA_ERROR(cudaMemGetInfo(&free, &total));
@@ -94,9 +36,7 @@ Buffer CudaAllocator::malloc(size_t size) {
   // Find available buffer from cache.
   auto orig_size = size;
   std::unique_lock lock(mutex_);
-  if (size <= small_block_size) {
-    size = 8;
-  } else if (size < page_size) {
+  if (size < page_size) {
     size = next_power_of_2(size);
   } else {
     size = page_size * ((size + page_size - 1) / page_size);
@@ -104,26 +44,20 @@ Buffer CudaAllocator::malloc(size_t size) {
 
   CudaBuffer* buf = buffer_cache_.reuse_from_cache(size);
   if (!buf) {
-    // If we have a lot of memory pressure try to reclaim memory from the cache.
-    int64_t mem_to_free =
-        get_active_memory() + get_cache_memory() + size - memory_limit_;
-    if (mem_to_free > 0) {
-      buffer_cache_.release_cached_buffers(mem_to_free);
+    // If we have a lot of memory pressure or are over the maximum cache size,
+    // try to reclaim memory from the cache.
+    size_t mem_required = get_active_memory() + get_cache_memory() + size;
+    if (mem_required >= memory_limit_) {
+      buffer_cache_.release_cached_buffers(mem_required - memory_limit_);
     }
 
-    // Try the scalar pool first
-    if (size <= small_block_size) {
-      buf = scalar_pool_.malloc();
-    }
     lock.unlock();
-    if (!buf) {
     buf = new CudaBuffer{nullptr, size};
     cudaError_t err = cudaMallocManaged(&buf->data, size);
     if (err != cudaSuccess && err != cudaErrorMemoryAllocation) {
       throw std::runtime_error(fmt::format(
           "cudaMallocManaged failed: {}.", cudaGetErrorString(err)));
     }
-    }
     lock.lock();
   }
   active_memory_ += size;
@@ -133,6 +67,7 @@ Buffer CudaAllocator::malloc(size_t size) {
   if (get_cache_memory() > max_pool_size_) {
     buffer_cache_.release_cached_buffers(get_cache_memory() - max_pool_size_);
   }
+
   return Buffer{buf};
 }
 
@@ -147,7 +82,9 @@ void CudaAllocator::free(Buffer buffer) {
   if (get_cache_memory() < max_pool_size_) {
     buffer_cache_.recycle_to_cache(buf);
   } else {
-    cuda_free(buf);
+    lock.unlock();
+    cuda_free(buf->data);
+    delete buf;
   }
 }
 
@@ -159,14 +96,27 @@ size_t CudaAllocator::size(Buffer buffer) const {
   return buf->size;
 }
 
-// This must be called with mutex_ aquired
-void CudaAllocator::cuda_free(CudaBuffer* buf) {
-  if (scalar_pool_.in_pool(buf)) {
-    scalar_pool_.free(buf);
-  } else {
-    cudaFree(buf->data);
-    delete buf;
+void CudaAllocator::register_this_thread() {
+  std::lock_guard lock(worker_mutex_);
+  allowed_threads_.insert(std::this_thread::get_id());
+}
+
+void CudaAllocator::cuda_free(void* buf) {
+  // If cuda_free() is called from a unregistered thread, reschedule the call to
+  // worker.
+  {
+    std::lock_guard lock(worker_mutex_);
+    if (allowed_threads_.count(std::this_thread::get_id()) == 0) {
+      if (!worker_) {
+        worker_.reset(new Worker);
       }
+      worker_->add_task([this, buf]() { this->cuda_free(buf); });
+      worker_->end_batch();
+      worker_->commit();
+      return;
+    }
+  }
+  cudaFree(buf);
 }
 
 size_t CudaAllocator::get_active_memory() const {

mlx/backend/cuda/allocator.h

@@ -7,10 +7,13 @@
 
 #include <mutex>
 #include <set>
+#include <thread>
 #include <utility>
 
 namespace mlx::core::cu {
 
+class Worker;
+
 using allocator::Buffer;
 
 // Stores cuda-managed unified memory.
@@ -19,35 +22,21 @@ struct CudaBuffer {
   size_t size;
 };
 
-class SmallSizePool {
- private:
-  union Block {
-    Block* next;
-    CudaBuffer buf;
-  };
-
-  Block* buffer_{nullptr};
-  void* data_{nullptr};
-  Block* next_free_{nullptr};
-
- public:
-  SmallSizePool();
-  ~SmallSizePool();
-
-  SmallSizePool(const SmallSizePool&) = delete;
-  SmallSizePool& operator=(const SmallSizePool&) = delete;
-
-  CudaBuffer* malloc();
-  void free(CudaBuffer* buf);
-  bool in_pool(CudaBuffer* buf);
-};
-
 class CudaAllocator : public allocator::Allocator {
  public:
   Buffer malloc(size_t size) override;
   void free(Buffer buffer) override;
   size_t size(Buffer buffer) const override;
 
+  // Register current thread as safe to free buffers.
+  // In cuda freeing a buffer implicitly synchronizes stream, and for threads
+  // that may be waited by gpu stream (for example cpu stream threads), freeing
+  // buffers there would result in dead lock.
+  void register_this_thread();
+
+  // Call cudaFree in the safe thread.
+  void cuda_free(void* buf);
+
   size_t get_active_memory() const;
   size_t get_peak_memory() const;
   void reset_peak_memory();
@@ -58,18 +47,19 @@ class CudaAllocator : public allocator::Allocator {
   void clear_cache();
 
  private:
-  void cuda_free(CudaBuffer* buf);
-
   CudaAllocator();
   friend CudaAllocator& allocator();
 
+  std::mutex worker_mutex_;
+  std::unique_ptr<Worker> worker_;
+  std::set<std::thread::id> allowed_threads_;
+
   std::mutex mutex_;
   size_t memory_limit_;
   size_t max_pool_size_;
   BufferCache<CudaBuffer> buffer_cache_;
   size_t active_memory_{0};
   size_t peak_memory_{0};
-  SmallSizePool scalar_pool_;
 };
 
 CudaAllocator& allocator();

mlx/backend/cuda/arange.cu

@@ -1,69 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/device.h"
-#include "mlx/backend/cuda/device/fp16_math.cuh"
-#include "mlx/backend/cuda/kernel_utils.cuh"
-#include "mlx/dtype_utils.h"
-#include "mlx/primitives.h"
-
-#include <cooperative_groups.h>
-#include <nvtx3/nvtx3.hpp>
-
-namespace mlx::core {
-
-namespace cu {
-
-namespace cg = cooperative_groups;
-
-template <typename T, typename IdxT, int N_WRITES>
-__global__ void arange(T* out, IdxT size, T start, T step) {
-  IdxT index = cg::this_grid().thread_rank();
-
-  if ((index + 1) * N_WRITES > size) {
-    for (IdxT i = index * N_WRITES; i < size; ++i) {
-      out[i] = start + i * step;
-    }
-  } else {
-    AlignedVector<T, N_WRITES> out_vec;
-#pragma unroll
-    for (int i = 0; i < N_WRITES; ++i) {
-      out_vec[i] = start + (index * N_WRITES + i) * step;
-    }
-
-    store_vector<N_WRITES>(out, index, out_vec);
-  }
-}
-
-} // namespace cu
-
-void Arange::eval_gpu(const std::vector<array>& inputs, array& out) {
-  nvtx3::scoped_range r("Arange::eval_gpu");
-  if (out.size() == 0) {
-    return;
-  }
-  out.set_data(allocator::malloc(out.nbytes()));
-
-  auto& encoder = cu::get_command_encoder(stream());
-  encoder.set_output_array(out);
-
-  dispatch_int_float_types(out.dtype(), "Arange", [&](auto type_tag) {
-    using CTYPE = MLX_GET_TYPE(type_tag);
-    using OutType = cuda_type_t<CTYPE>;
-    constexpr int N_WRITES = 16 / sizeof(OutType);
-    dispatch_bool(out.data_size() > INT32_MAX, [&](auto large) {
-      using IdxT = std::conditional_t<large(), int64_t, int32_t>;
-      auto [num_blocks, block_dims] = get_launch_args(out, large(), N_WRITES);
-      encoder.add_kernel_node(
-          cu::arange<OutType, IdxT, N_WRITES>,
-          num_blocks,
-          block_dims,
-          0,
-          out.data<OutType>(),
-          out.data_size(),
-          static_cast<CTYPE>(start_),
-          static_cast<CTYPE>(start_ + step_) - static_cast<CTYPE>(start_));
-    });
-  });
-}
-
-} // namespace mlx::core

mlx/backend/cuda/arg_reduce.cu

@@ -1,8 +1,8 @@
 // Copyright © 2025 Apple Inc.
-
 #include "mlx/backend/common/utils.h"
 #include "mlx/backend/cuda/device.h"
 #include "mlx/backend/cuda/device/fp16_math.cuh"
+#include "mlx/backend/cuda/iterators/strided_iterator.cuh"
 #include "mlx/backend/cuda/kernel_utils.cuh"
 #include "mlx/dtype_utils.h"
 #include "mlx/primitives.h"
@@ -44,11 +44,8 @@ struct ArgMin {
   }
 
   template <int N>
-  __device__ IndexValPair<T> reduce_many(
-      IndexValPair<T> best,
-      const AlignedVector<T, N>& vals,
-      uint32_t offset) {
-#pragma unroll
+  __device__ IndexValPair<T>
+  reduce_many(IndexValPair<T> best, T (&vals)[N], uint32_t offset) {
     for (int i = 0; i < N; i++) {
       if (vals[i] < best.val) {
         best.val = vals[i];
@@ -77,11 +74,8 @@ struct ArgMax {
   }
 
   template <int N>
-  __device__ IndexValPair<T> reduce_many(
-      IndexValPair<T> best,
-      const AlignedVector<T, N>& vals,
-      uint32_t offset) {
-#pragma unroll
+  __device__ IndexValPair<T>
+  reduce_many(IndexValPair<T> best, T (&vals)[N], uint32_t offset) {
     for (int i = 0; i < N; i++) {
       if (vals[i] > best.val) {
         best.val = vals[i];
@@ -112,15 +106,16 @@ __global__ void arg_reduce_general(
 
   int64_t in_idx = elem_to_loc(index, shape.data(), in_strides.data(), ndim);
   int64_t out_idx = elem_to_loc(index, shape.data(), out_strides.data(), ndim);
-  in += in_idx;
 
   Op op;
   T init = op.init();
   IndexValPair<T> best{0, init};
 
   for (int r = 0; r < cuda::ceil_div(axis_size, BLOCK_DIM * N_READS); ++r) {
+    T vals[N_READS];
     auto tid = r * BLOCK_DIM + block.thread_index().x;
-    auto vals = load_vector<N_READS>(in, tid, axis_size, axis_stride, init);
+    cub::LoadDirectBlocked(
+        tid, strided_iterator(in + in_idx, axis_stride), vals, axis_size, init);
     best = op.reduce_many(best, vals, tid * N_READS);
   }
 
@@ -171,7 +166,6 @@ void ArgReduce::eval_gpu(const std::vector<array>& inputs, array& out) {
           kernel,
           num_blocks,
           block_dim(),
-          0,
           in.data<T>(),
           out.data<uint32_t>(),
           out.size(),

mlx/backend/cuda/binary.cu

@@ -3,6 +3,7 @@
 #include "mlx/backend/common/binary.h"
 #include "mlx/backend/cuda/device.h"
 #include "mlx/backend/cuda/device/binary_ops.cuh"
+#include "mlx/backend/cuda/device/cucomplex_math.cuh"
 #include "mlx/backend/cuda/kernel_utils.cuh"
 #include "mlx/dtype_utils.h"
 #include "mlx/primitives.h"
@@ -28,7 +29,7 @@ __global__ void binary_ss(const In* a, const In* b, Out* out, IdxT size) {
     AlignedVector<Out, N_READS> out_vec;
 #pragma unroll
     for (int i = 0; i < N_READS; ++i) {
-      out_vec[i] = Op{}(a[0], b[0]);
+      out_vec.val[i] = Op{}(a[0], b[0]);
     }
 
     store_vector<N_READS>(out, index, out_vec);
@@ -49,7 +50,7 @@ __global__ void binary_sv(const In* a, const In* b, Out* out, IdxT size) {
     AlignedVector<Out, N_READS> out_vec;
 #pragma unroll
     for (int i = 0; i < N_READS; ++i) {
-      out_vec[i] = Op{}(a[0], b_vec[i]);
+      out_vec.val[i] = Op{}(a[0], b_vec.val[i]);
     }
 
     store_vector<N_READS>(out, index, out_vec);
@@ -70,7 +71,7 @@ __global__ void binary_vs(const In* a, const In* b, Out* out, IdxT size) {
     AlignedVector<Out, N_READS> out_vec;
 #pragma unroll
     for (int i = 0; i < N_READS; ++i) {
-      out_vec[i] = Op{}(a_vec[i], b[0]);
+      out_vec.val[i] = Op{}(a_vec.val[i], b[0]);
     }
 
     store_vector<N_READS>(out, index, out_vec);
@@ -92,96 +93,46 @@ __global__ void binary_vv(const In* a, const In* b, Out* out, IdxT size) {
     AlignedVector<Out, N_READS> out_vec;
 #pragma unroll
     for (int i = 0; i < N_READS; ++i) {
-      out_vec[i] = Op{}(a_vec[i], b_vec[i]);
+      out_vec.val[i] = Op{}(a_vec.val[i], b_vec.val[i]);
     }
 
     store_vector<N_READS>(out, index, out_vec);
   }
 }
 
-template <
-    typename Op,
-    typename In,
-    typename Out,
-    typename IdxT,
-    int NDIM,
-    int N_READS>
+template <typename Op, typename In, typename Out, typename IdxT, int NDIM>
 __global__ void binary_g_nd(
     const In* a,
     const In* b,
     Out* out,
-    IdxT size_rest,
+    IdxT size,
     const __grid_constant__ cuda::std::array<int32_t, NDIM> shape,
     const __grid_constant__ cuda::std::array<int64_t, NDIM> a_strides,
     const __grid_constant__ cuda::std::array<int64_t, NDIM> b_strides) {
-  auto block = cg::this_thread_block();
-  auto grid = cg::this_grid();
-  IdxT index_rest =
-      grid.block_index().y * block.dim_threads().y + block.thread_index().y;
-  if (index_rest >= size_rest) {
-    return;
-  }
-
-  auto shape_x = shape[NDIM - 1];
-  auto a_stride_x = a_strides[NDIM - 1];
-  auto b_stride_x = b_strides[NDIM - 1];
-  IdxT index_x =
-      grid.block_index().x * block.dim_threads().x + block.thread_index().x;
+  IdxT index = cg::this_grid().thread_rank();
+  if (index < size) {
     auto [a_idx, b_idx] = elem_to_loc_nd<NDIM>(
-      index_rest * shape_x, shape.data(), a_strides.data(), b_strides.data());
-  auto a_vec =
-      load_vector<N_READS>(a + a_idx, index_x, shape_x, a_stride_x, In(0));
-  auto b_vec =
-      load_vector<N_READS>(b + b_idx, index_x, shape_x, b_stride_x, In(0));
-
-  AlignedVector<Out, N_READS> out_vec;
-#pragma unroll
-  for (int i = 0; i < N_READS; ++i) {
-    out_vec[i] = Op{}(a_vec[i], b_vec[i]);
+        index, shape.data(), a_strides.data(), b_strides.data());
+    out[index] = Op{}(a[a_idx], b[b_idx]);
   }
-  store_vector(out + shape_x * index_rest, index_x, out_vec, shape_x);
 }
 
-template <typename Op, typename In, typename Out, typename IdxT, int N_READS>
+template <typename Op, typename In, typename Out, typename IdxT>
 __global__ void binary_g(
     const In* a,
     const In* b,
     Out* out,
-    IdxT size_rest,
+    IdxT size,
     const __grid_constant__ Shape shape,
     const __grid_constant__ Strides a_strides,
     const __grid_constant__ Strides b_strides,
     int ndim) {
-  auto block = cg::this_thread_block();
-  auto grid = cg::this_grid();
-  IdxT index_rest =
-      grid.block_index().y * block.dim_threads().y + block.thread_index().y;
-  if (index_rest >= size_rest) {
-    return;
+  IdxT index = cg::this_grid().thread_rank();
+  if (index < size) {
+    auto [a_idx, b_idx] = elem_to_loc_4d(
+        index, shape.data(), a_strides.data(), b_strides.data(), ndim);
+    out[index] = Op{}(a[a_idx], b[b_idx]);
   }
-
-  auto shape_x = shape[ndim - 1];
-  auto a_stride_x = a_strides[ndim - 1];
-  auto b_stride_x = b_strides[ndim - 1];
-  IdxT index_x =
-      grid.block_index().x * block.dim_threads().x + block.thread_index().x;
-  auto [a_idx, b_idx] = elem_to_loc(
-      index_rest * shape_x,
-      shape.data(),
-      a_strides.data(),
-      b_strides.data(),
-      ndim);
-  auto a_vec =
-      load_vector<N_READS>(a + a_idx, index_x, shape_x, a_stride_x, In(0));
-  auto b_vec =
-      load_vector<N_READS>(b + b_idx, index_x, shape_x, b_stride_x, In(0));
-
-  AlignedVector<Out, N_READS> out_vec;
-#pragma unroll
-  for (int i = 0; i < N_READS; ++i) {
-    out_vec[i] = Op{}(a_vec[i], b_vec[i]);
-  }
-  store_vector(out + shape_x * index_rest, index_x, out_vec, shape_x);
 }
 
 template <typename Op, typename In, typename Out>
@@ -226,7 +177,7 @@ template <typename Op>
 void binary_op_gpu_inplace(
     const std::vector<array>& inputs,
     array& out,
-    const char* op,
+    std::string_view op,
     const Stream& s) {
   assert(inputs.size() > 1);
   const auto& a = inputs[0];
@@ -259,61 +210,36 @@ void binary_op_gpu_inplace(
                 auto& a_strides = strides[0];
                 auto& b_strides = strides[1];
                 int ndim = shape.size();
-                int work_per_thread = 1;
-                auto dim0 = ndim > 0 ? shape.back() : 1;
-                auto rest = out.size() / dim0;
-                if (dim0 >= 4) {
-                  work_per_thread = 4;
-                }
-                dim0 = (dim0 + work_per_thread - 1) / work_per_thread;
-                auto block_dims = get_block_dims(dim0, rest, 1);
-                uint32_t num_blocks_x = cuda::ceil_div(dim0, block_dims.x);
-                uint32_t num_blocks_y = cuda::ceil_div(rest, block_dims.y);
                 if (ndim <= 3) {
                   dispatch_1_2_3(ndim, [&](auto dims_constant) {
-                    auto kernel = cu::binary_g_nd<
-                        Op,
-                        InType,
-                        OutType,
-                        IdxT,
-                        dims_constant(),
-                        1>;
-                    if (work_per_thread == 4) {
-                      kernel = cu::binary_g_nd<
-                          Op,
-                          InType,
-                          OutType,
-                          IdxT,
-                          dims_constant(),
-                          4>;
-                    }
+                    auto kernel = cu::
+                        binary_g_nd<Op, InType, OutType, IdxT, dims_constant()>;
+                    auto [num_blocks, block_dims] =
+                        get_launch_args(kernel, out, large());
                     encoder.add_kernel_node(
                         kernel,
-                        {num_blocks_x, num_blocks_y},
+                        num_blocks,
                         block_dims,
-                        0,
                         a.data<InType>(),
                         b.data<InType>(),
                         out.data<OutType>(),
-                        rest,
+                        out.size(),
                         const_param<dims_constant()>(shape),
                         const_param<dims_constant()>(a_strides),
                         const_param<dims_constant()>(b_strides));
                   });
                 } else {
-                  auto kernel = cu::binary_g<Op, InType, OutType, IdxT, 1>;
-                  if (work_per_thread == 4) {
-                    kernel = cu::binary_g<Op, InType, OutType, IdxT, 4>;
-                  }
+                  auto kernel = cu::binary_g<Op, InType, OutType, IdxT>;
+                  auto [num_blocks, block_dims] =
+                      get_launch_args(kernel, out, large());
                   encoder.add_kernel_node(
                       kernel,
-                      {num_blocks_x, num_blocks_y},
+                      num_blocks,
                       block_dims,
-                      0,
                       a.data<InType>(),
                       b.data<InType>(),
                       out.data<OutType>(),
-                      rest,
+                      out.size(),
                       const_param(shape),
                       const_param(a_strides),
                       const_param(b_strides),
@@ -323,7 +249,8 @@ void binary_op_gpu_inplace(
         } else {
           dispatch_bool(out.data_size() > UINT32_MAX, [&](auto large) {
             using IdxT = std::conditional_t<large(), int64_t, uint32_t>;
-            constexpr int N_READS = 16 / sizeof(InType);
+            // TODO: Choose optimized value based on type size.
+            constexpr int N_READS = 4;
             auto kernel = cu::binary_ss<Op, InType, OutType, IdxT, N_READS>;
             if (bopt == BinaryOpType::ScalarVector) {
               kernel = cu::binary_sv<Op, InType, OutType, IdxT, N_READS>;
@@ -333,12 +260,16 @@ void binary_op_gpu_inplace(
               kernel = cu::binary_vv<Op, InType, OutType, IdxT, N_READS>;
             }
             auto [num_blocks, block_dims] = get_launch_args(
-                out.data_size(), out.shape(), out.strides(), large(), N_READS);
+                kernel,
+                out.data_size(),
+                out.shape(),
+                out.strides(),
+                large(),
+                N_READS);
             encoder.add_kernel_node(
                 kernel,
                 num_blocks,
                 block_dims,
-                0,
                 a.data<InType>(),
                 b.data<InType>(),
                 out.data<OutType>(),
@@ -360,7 +291,7 @@ template <typename Op>
 void binary_op_gpu(
     const std::vector<array>& inputs,
     array& out,
-    const char* op,
+    std::string_view op,
     const Stream& s) {
   auto& a = inputs[0];
   auto& b = inputs[1];
@@ -373,7 +304,59 @@ void binary_op_gpu(
   void func::eval_gpu(const std::vector<array>& inputs, array& out) {    \
     nvtx3::scoped_range r(#func "::eval_gpu");                           \
     auto& s = out.primitive().stream();                                  \
-    binary_op_gpu<cu::func>(inputs, out, name(), s);                  \
+    binary_op_gpu<cu::func>(inputs, out, get_primitive_string(this), s); \
   }
 
+BINARY_GPU(Add)
+BINARY_GPU(ArcTan2)
+BINARY_GPU(Divide)
+BINARY_GPU(Remainder)
+BINARY_GPU(Greater)
+BINARY_GPU(GreaterEqual)
+BINARY_GPU(Less)
+BINARY_GPU(LessEqual)
+BINARY_GPU(LogicalAnd)
+BINARY_GPU(LogicalOr)
+BINARY_GPU(LogAddExp)
+BINARY_GPU(Maximum)
+BINARY_GPU(Minimum)
+BINARY_GPU(Multiply)
+BINARY_GPU(NotEqual)
+BINARY_GPU(Power)
+BINARY_GPU(Subtract)
+
+void Equal::eval_gpu(const std::vector<array>& inputs, array& out) {
+  nvtx3::scoped_range r("Equal::eval_gpu");
+  auto& s = out.primitive().stream();
+  auto op = get_primitive_string(this);
+  if (equal_nan_) {
+    binary_op_gpu<cu::NaNEqual>(inputs, out, op, s);
+  } else {
+    binary_op_gpu<cu::Equal>(inputs, out, op, s);
+  }
+}
+
+void BitwiseBinary::eval_gpu(const std::vector<array>& inputs, array& out) {
+  nvtx3::scoped_range r("BitwiseBinary::eval_gpu");
+  auto& s = out.primitive().stream();
+  auto op = get_primitive_string(this);
+  switch (op_) {
+    case BitwiseBinary::And:
+      binary_op_gpu<cu::BitwiseAnd>(inputs, out, op, s);
+      break;
+    case BitwiseBinary::Or:
+      binary_op_gpu<cu::BitwiseOr>(inputs, out, op, s);
+      break;
+    case BitwiseBinary::Xor:
+      binary_op_gpu<cu::BitwiseXor>(inputs, out, op, s);
+      break;
+    case BitwiseBinary::LeftShift:
+      binary_op_gpu<cu::LeftShift>(inputs, out, op, s);
+      break;
+    case BitwiseBinary::RightShift:
+      binary_op_gpu<cu::RightShift>(inputs, out, op, s);
+      break;
+  }
+}
+
 } // namespace mlx::core

mlx/backend/cuda/binary/CMakeLists.txt

@@ -1,21 +0,0 @@
-target_sources(
-  mlx
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/add.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/arctan2.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/bitwise_binary.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/divide.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/equal.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/greater.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/greater_equal.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/less.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/less_equal.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/logical_and.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/logical_or.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/log_add_exp.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/minimum.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/maximum.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/multiply.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/power.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/remainder.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/not_equal.cu
-  PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/subtract.cu)

mlx/backend/cuda/binary/add.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(Add)
-} // namespace mlx::core

mlx/backend/cuda/binary/arctan2.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(ArcTan2)
-} // namespace mlx::core

mlx/backend/cuda/binary/bitwise_binary.cu

@@ -1,27 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-void BitwiseBinary::eval_gpu(const std::vector<array>& inputs, array& out) {
-  nvtx3::scoped_range r("BitwiseBinary::eval_gpu");
-  auto& s = out.primitive().stream();
-  switch (op_) {
-    case BitwiseBinary::And:
-      binary_op_gpu<cu::BitwiseAnd>(inputs, out, name(), s);
-      break;
-    case BitwiseBinary::Or:
-      binary_op_gpu<cu::BitwiseOr>(inputs, out, name(), s);
-      break;
-    case BitwiseBinary::Xor:
-      binary_op_gpu<cu::BitwiseXor>(inputs, out, name(), s);
-      break;
-    case BitwiseBinary::LeftShift:
-      binary_op_gpu<cu::LeftShift>(inputs, out, name(), s);
-      break;
-    case BitwiseBinary::RightShift:
-      binary_op_gpu<cu::RightShift>(inputs, out, name(), s);
-      break;
-  }
-}
-} // namespace mlx::core

mlx/backend/cuda/binary/divide.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(Divide)
-} // namespace mlx::core

mlx/backend/cuda/binary/equal.cu

@@ -1,15 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-void Equal::eval_gpu(const std::vector<array>& inputs, array& out) {
-  nvtx3::scoped_range r("Equal::eval_gpu");
-  auto& s = out.primitive().stream();
-  if (equal_nan_) {
-    binary_op_gpu<cu::NaNEqual>(inputs, out, name(), s);
-  } else {
-    binary_op_gpu<cu::Equal>(inputs, out, name(), s);
-  }
-}
-} // namespace mlx::core

mlx/backend/cuda/binary/greater.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(Greater)
-} // namespace mlx::core

mlx/backend/cuda/binary/greater_equal.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(GreaterEqual)
-} // namespace mlx::core

mlx/backend/cuda/binary/less.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(Less)
-} // namespace mlx::core

mlx/backend/cuda/binary/less_equal.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(LessEqual)
-} // namespace mlx::core

mlx/backend/cuda/binary/log_add_exp.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(LogAddExp)
-} // namespace mlx::core

mlx/backend/cuda/binary/logical_and.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(LogicalAnd)
-} // namespace mlx::core

mlx/backend/cuda/binary/logical_or.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(LogicalOr)
-} // namespace mlx::core

mlx/backend/cuda/binary/maximum.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(Maximum)
-} // namespace mlx::core

mlx/backend/cuda/binary/minimum.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(Minimum)
-} // namespace mlx::core

mlx/backend/cuda/binary/multiply.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(Multiply)
-} // namespace mlx::core

mlx/backend/cuda/binary/not_equal.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(NotEqual)
-} // namespace mlx::core

mlx/backend/cuda/binary/power.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(Power)
-} // namespace mlx::core

mlx/backend/cuda/binary/remainder.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(Remainder)
-} // namespace mlx::core

mlx/backend/cuda/binary/subtract.cu

@@ -1,7 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/binary/binary.cuh"
-
-namespace mlx::core {
-BINARY_GPU(Subtract)
-} // namespace mlx::core

mlx/backend/cuda/binary_two.cu

@@ -3,6 +3,7 @@
 #include "mlx/backend/common/binary.h"
 #include "mlx/backend/cuda/device.h"
 #include "mlx/backend/cuda/device/binary_ops.cuh"
+#include "mlx/backend/cuda/device/cucomplex_math.cuh"
 #include "mlx/backend/cuda/kernel_utils.cuh"
 #include "mlx/dtype_utils.h"
 #include "mlx/primitives.h"
@@ -33,8 +34,8 @@ binary_two_ss(const In* a, const In* b, Out* out_a, Out* out_b, IdxT size) {
 #pragma unroll
     for (int i = 0; i < N_READS; ++i) {
       auto out = Op{}(a[0], b[0]);
-      out_a_vec[i] = out[0];
-      out_b_vec[i] = out[1];
+      out_a_vec.val[i] = out[0];
+      out_b_vec.val[i] = out[1];
     }
 
     store_vector<N_READS>(out_a, index, out_a_vec);
@@ -60,9 +61,9 @@ binary_two_sv(const In* a, const In* b, Out* out_a, Out* out_b, IdxT size) {
     AlignedVector<Out, N_READS> out_b_vec;
 #pragma unroll
     for (int i = 0; i < N_READS; ++i) {
-      auto out = Op{}(a[0], b_vec[i]);
-      out_a_vec[i] = out[0];
-      out_b_vec[i] = out[1];
+      auto out = Op{}(a[0], b_vec.val[i]);
+      out_a_vec.val[i] = out[0];
+      out_b_vec.val[i] = out[1];
     }
 
     store_vector<N_READS>(out_a, index, out_a_vec);
@@ -88,9 +89,9 @@ binary_two_vs(const In* a, const In* b, Out* out_a, Out* out_b, IdxT size) {
     AlignedVector<Out, N_READS> out_b_vec;
 #pragma unroll
     for (int i = 0; i < N_READS; ++i) {
-      auto out = Op{}(a_vec[i], b[0]);
-      out_a_vec[i] = out[0];
-      out_b_vec[i] = out[1];
+      auto out = Op{}(a_vec.val[i], b[0]);
+      out_a_vec.val[i] = out[0];
+      out_b_vec.val[i] = out[1];
     }
 
     store_vector<N_READS>(out_a, index, out_a_vec);
@@ -117,9 +118,9 @@ binary_two_vv(const In* a, const In* b, Out* out_a, Out* out_b, IdxT size) {
     AlignedVector<Out, N_READS> out_b_vec;
 #pragma unroll
     for (int i = 0; i < N_READS; ++i) {
-      auto out = Op{}(a_vec[i], b_vec[i]);
-      out_a_vec[i] = out[0];
-      out_b_vec[i] = out[1];
+      auto out = Op{}(a_vec.val[i], b_vec.val[i]);
+      out_a_vec.val[i] = out[0];
+      out_b_vec.val[i] = out[1];
     }
 
     store_vector<N_READS>(out_a, index, out_a_vec);
@@ -127,99 +128,45 @@ binary_two_vv(const In* a, const In* b, Out* out_a, Out* out_b, IdxT size) {
   }
 }
 
-template <
-    typename Op,
-    typename In,
-    typename Out,
-    typename IdxT,
-    int NDIM,
-    int N_READS>
+template <typename Op, typename In, typename Out, typename IdxT, int NDIM>
 __global__ void binary_two_g_nd(
     const In* a,
     const In* b,
     Out* out_a,
     Out* out_b,
-    IdxT size_rest,
+    IdxT size,
     const __grid_constant__ cuda::std::array<int32_t, NDIM> shape,
     const __grid_constant__ cuda::std::array<int64_t, NDIM> a_strides,
     const __grid_constant__ cuda::std::array<int64_t, NDIM> b_strides) {
-  auto block = cg::this_thread_block();
-  auto grid = cg::this_grid();
-  IdxT index_rest =
-      grid.block_index().y * block.dim_threads().y + block.thread_index().y;
-  if (index_rest >= size_rest) {
-    return;
-  }
-
-  auto shape_x = shape[NDIM - 1];
-  auto a_stride_x = a_strides[NDIM - 1];
-  auto b_stride_x = b_strides[NDIM - 1];
-  IdxT index_x =
-      grid.block_index().x * block.dim_threads().x + block.thread_index().x;
+  IdxT index = cg::this_grid().thread_rank();
+  if (index < size) {
     auto [a_idx, b_idx] = elem_to_loc_nd<NDIM>(
-      index_rest * shape_x, shape.data(), a_strides.data(), b_strides.data());
-  auto a_vec =
-      load_vector<N_READS>(a + a_idx, index_x, shape_x, a_stride_x, In(0));
-  auto b_vec =
-      load_vector<N_READS>(b + b_idx, index_x, shape_x, b_stride_x, In(0));
-
-  AlignedVector<Out, N_READS> out_vec_a;
-  AlignedVector<Out, N_READS> out_vec_b;
-#pragma unroll
-  for (int i = 0; i < N_READS; ++i) {
-    auto out = Op{}(a_vec[i], b_vec[i]);
-    out_vec_a[i] = out[0];
-    out_vec_b[i] = out[1];
+        index, shape.data(), a_strides.data(), b_strides.data());
+    auto out = Op{}(a[a_idx], b[b_idx]);
+    out_a[index] = out[0];
+    out_b[index] = out[1];
   }
-  store_vector(out_a + shape_x * index_rest, index_x, out_vec_a, shape_x);
-  store_vector(out_b + shape_x * index_rest, index_x, out_vec_b, shape_x);
 }
 
-template <typename Op, typename In, typename Out, typename IdxT, int N_READS>
+template <typename Op, typename In, typename Out, typename IdxT>
 __global__ void binary_two_g(
     const In* a,
     const In* b,
     Out* out_a,
     Out* out_b,
-    IdxT size_rest,
+    IdxT size,
     const __grid_constant__ Shape shape,
     const __grid_constant__ Strides a_strides,
     const __grid_constant__ Strides b_strides,
     int ndim) {
-  auto block = cg::this_thread_block();
-  auto grid = cg::this_grid();
-  IdxT index_rest =
-      grid.block_index().y * block.dim_threads().y + block.thread_index().y;
-  if (index_rest >= size_rest) {
-    return;
+  IdxT index = cg::this_grid().thread_rank();
+  if (index < size) {
+    auto [a_idx, b_idx] = elem_to_loc_4d(
+        index, shape.data(), a_strides.data(), b_strides.data(), ndim);
+    auto out = Op{}(a[a_idx], b[b_idx]);
+    out_a[index] = out[0];
+    out_b[index] = out[1];
   }
-
-  auto shape_x = shape[ndim - 1];
-  auto a_stride_x = a_strides[ndim - 1];
-  auto b_stride_x = b_strides[ndim - 1];
-  IdxT index_x =
-      grid.block_index().x * block.dim_threads().x + block.thread_index().x;
-  auto [a_idx, b_idx] = elem_to_loc(
-      index_rest * shape_x,
-      shape.data(),
-      a_strides.data(),
-      b_strides.data(),
-      ndim);
-  auto a_vec =
-      load_vector<N_READS>(a + a_idx, index_x, shape_x, a_stride_x, In(0));
-  auto b_vec =
-      load_vector<N_READS>(b + b_idx, index_x, shape_x, b_stride_x, In(0));
-
-  AlignedVector<Out, N_READS> out_vec_a;
-  AlignedVector<Out, N_READS> out_vec_b;
-#pragma unroll
-  for (int i = 0; i < N_READS; ++i) {
-    auto out = Op{}(a_vec[i], b_vec[i]);
-    out_vec_a[i] = out[0];
-    out_vec_b[i] = out[1];
-  }
-  store_vector(out_a + shape_x * index_rest, index_x, out_vec_a, shape_x);
-  store_vector(out_b + shape_x * index_rest, index_x, out_vec_b, shape_x);
 }
 
 template <typename Op, typename In, typename Out>
@@ -237,7 +184,7 @@ template <typename Op>
 void binary_two_op_gpu_inplace(
     const std::vector<array>& inputs,
     std::vector<array>& outputs,
-    const char* op,
+    std::string_view op,
     const Stream& s) {
   assert(inputs.size() > 1);
   const auto& a = inputs[0];
@@ -279,17 +226,6 @@ void binary_two_op_gpu_inplace(
                 auto& a_strides = strides[0];
                 auto& b_strides = strides[1];
                 int ndim = shape.size();
-                int work_per_thread = 1;
-                auto dim0 = ndim > 0 ? shape.back() : 1;
-                auto rest = out_a.size() / dim0;
-                if (dim0 >= 4) {
-                  work_per_thread = 4;
-                }
-                dim0 = (dim0 + work_per_thread - 1) / work_per_thread;
-                auto block_dims = get_block_dims(dim0, rest, 1);
-                uint32_t num_blocks_x = cuda::ceil_div(dim0, block_dims.x);
-                uint32_t num_blocks_y = cuda::ceil_div(rest, block_dims.y);
-
                 if (ndim <= 3) {
                   dispatch_1_2_3(ndim, [&](auto dims_constant) {
                     auto kernel = cu::binary_two_g_nd<
@@ -297,46 +233,35 @@ void binary_two_op_gpu_inplace(
                         InType,
                         OutType,
                         IdxT,
-                        dims_constant(),
-                        1>;
-                    if (work_per_thread == 4) {
-                      kernel = cu::binary_two_g_nd<
-                          Op,
-                          InType,
-                          OutType,
-                          IdxT,
-                          dims_constant(),
-                          4>;
-                    }
+                        dims_constant()>;
+                    auto [num_blocks, block_dims] =
+                        get_launch_args(kernel, out_a, large());
                     encoder.add_kernel_node(
                         kernel,
-                        {num_blocks_x, num_blocks_y},
+                        num_blocks,
                         block_dims,
-                        0,
                         a.data<InType>(),
                         b.data<InType>(),
                         out_a.data<OutType>(),
                         out_b.data<OutType>(),
-                        rest,
+                        out_a.size(),
                         const_param<dims_constant()>(shape),
                         const_param<dims_constant()>(a_strides),
                         const_param<dims_constant()>(b_strides));
                   });
                 } else {
-                  auto kernel = cu::binary_two_g<Op, InType, OutType, IdxT, 1>;
-                  if (work_per_thread == 4) {
-                    kernel = cu::binary_two_g<Op, InType, OutType, IdxT, 4>;
-                  }
+                  auto kernel = cu::binary_two_g<Op, InType, OutType, IdxT>;
+                  auto [num_blocks, block_dims] =
+                      get_launch_args(kernel, out_a, large());
                   encoder.add_kernel_node(
                       kernel,
-                      {num_blocks_x, num_blocks_y},
+                      num_blocks,
                       block_dims,
-                      0,
                       a.data<InType>(),
                       b.data<InType>(),
                       out_a.data<OutType>(),
                       out_b.data<OutType>(),
-                      rest,
+                      out_a.size(),
                       const_param(shape),
                       const_param(a_strides),
                       const_param(b_strides),
@@ -346,7 +271,8 @@ void binary_two_op_gpu_inplace(
         } else {
           dispatch_bool(out_a.data_size() > UINT32_MAX, [&](auto large) {
             using IdxT = std::conditional_t<large(), int64_t, uint32_t>;
-            constexpr int N_READS = 16 / sizeof(InType);
+            // TODO: Choose optimized value based on type size.
+            constexpr int N_READS = 4;
             auto kernel = cu::binary_two_ss<Op, InType, OutType, IdxT, N_READS>;
             if (bopt == BinaryOpType::ScalarVector) {
               kernel = cu::binary_two_sv<Op, InType, OutType, IdxT, N_READS>;
@@ -356,6 +282,7 @@ void binary_two_op_gpu_inplace(
               kernel = cu::binary_two_vv<Op, InType, OutType, IdxT, N_READS>;
             }
             auto [num_blocks, block_dims] = get_launch_args(
+                kernel,
                 out_a.data_size(),
                 out_a.shape(),
                 out_a.strides(),
@@ -365,7 +292,6 @@ void binary_two_op_gpu_inplace(
                 kernel,
                 num_blocks,
                 block_dims,
-                0,
                 a.data<InType>(),
                 b.data<InType>(),
                 out_a.data<OutType>(),
@@ -388,7 +314,7 @@ template <typename Op>
 void binary_two_op_gpu(
     const std::vector<array>& inputs,
     std::vector<array>& outputs,
-    const char* op,
+    std::string_view op,
     const Stream& s) {
   auto& a = inputs[0];
   auto& b = inputs[1];
@@ -403,7 +329,7 @@ void DivMod::eval_gpu(
     std::vector<array>& outputs) {
   nvtx3::scoped_range r("DivMod::eval_gpu");
   auto& s = outputs[0].primitive().stream();
-  binary_two_op_gpu<cu::DivMod>(inputs, outputs, name(), s);
+  binary_two_op_gpu<cu::DivMod>(inputs, outputs, get_primitive_string(this), s);
 }
 
 } // namespace mlx::core

mlx/backend/cuda/compiled.cpp

@@ -53,10 +53,9 @@ struct FusedKernelBuilder {
 
     // Build function signature.
     if (contiguous) {
-      os += "template <typename IdxT = uint32_t, int work_per_thread = 1>\n";
+      os += "template <typename IdxT = uint32_t>\n";
     } else {
-      os +=
-          "template <int NDIM, typename IdxT = uint32_t, int work_per_thread = 1>\n";
+      os += "template <int NDIM, typename IdxT = uint32_t>\n";
     }
     os += fmt::format("__global__ void {}(\n", kernel_name + name);
     for (size_t i = 0; i < params.size(); ++i) {
@@ -68,77 +67,12 @@ struct FusedKernelBuilder {
     }
     os += ") {\n";
 
-    // Index. For non contiguous kernels we create a separate index
-    // variable per variable otherwise everyone uses `index`.
+    // Index.
     os +=
-        "  IdxT index = cg::this_grid().thread_rank() * work_per_thread;\n"
+        "  IdxT index = cg::this_grid().thread_rank();\n"
         "  if (index >= size) {\n"
         "    return;\n"
         "  }\n";
-    if (!contiguous) {
-      for (size_t i = 0; i < inputs.size(); ++i) {
-        const auto& x = inputs[i];
-        const std::string& xname = namer.get_name(x);
-        if (is_scalar(x) || is_constant(i)) {
-          continue;
-        }
-        os += "  IdxT " + xname + "_idx = 0;\n";
-      }
-      os += "  {\n";
-      os += "    IdxT loc = index;\n";
-      os +=
-          "    #pragma unroll\n"
-          "    for (int i = NDIM - 1; i >= 0; i--) {\n";
-      for (size_t i = 0; i < inputs.size(); ++i) {
-        const auto& x = inputs[i];
-        const std::string& xname = namer.get_name(x);
-        if (is_scalar(x) || is_constant(i)) {
-          continue;
-        }
-        os += "      " + xname + "_idx += (loc \% shape[i]) * IdxT(" + xname +
-            "_strides[i]);\n";
-      }
-      os +=
-          "      loc /= shape[i];\n"
-          "    }\n"
-          "  }\n";
-    }
-
-    // Vectorized read loop
-    if (contiguous) {
-      for (size_t i = 0; i < inputs.size(); ++i) {
-        const auto& x = inputs[i];
-        if (is_scalar(x) || is_constant(i)) {
-          continue;
-        }
-        const std::string& xname = namer.get_name(x);
-        std::string type = dtype_to_cuda_type(x.dtype());
-        os += fmt::format(
-            "  auto vec_{0} = load_vector<work_per_thread, {1}>({0} + index, 0, size - index, 0);\n",
-            xname,
-            type);
-      }
-    }
-
-    // Create some space for the outputs
-    for (const auto& x : outputs) {
-      const std::string& xname = namer.get_name(x);
-      std::string type = dtype_to_cuda_type(x.dtype());
-      os += fmt::format(
-          "  AlignedVector<{}, work_per_thread> vec_{};\n", type, xname);
-    }
-
-    // Work loop
-    if (!contiguous) {
-      os +=
-          "\n"
-          "  for (int i = 0; i < work_per_thread && index < size; i++) {\n";
-    } else {
-      os +=
-          "\n"
-          "  #pragma unroll\n"
-          "  for (int i = 0; i < work_per_thread; i++) {\n";
-    }
 
     // Read inputs.
     for (size_t i = 0; i < inputs.size(); ++i) {
@@ -153,9 +87,12 @@ struct FusedKernelBuilder {
       } else if (is_scalar(x)) {
         value = fmt::format("{}[0]", xname);
       } else if (contiguous) {
-        value = fmt::format("vec_{}[i]", xname);
+        value = fmt::format("{}[index]", xname);
       } else {
-        value = fmt::format("{}[{}_idx]", xname, xname);
+        std::string index = fmt::format(
+            "elem_to_loc_nd<NDIM>(index, shape.data(), {}_strides.data())",
+            xname);
+        value = fmt::format("{}[{}]", xname, index);
       }
       os += fmt::format("  {} tmp_{} = {};\n", type, xname, value);
     }
@@ -169,7 +106,9 @@ struct FusedKernelBuilder {
         value = fmt::format(
             "static_cast<{}>(tmp_{})", type, namer.get_name(x.inputs()[0]));
       } else {
-        value = x.primitive().name();
+        std::ostringstream ss;
+        x.primitive().print(ss);
+        value = ss.str();
         value += "{}(";
         for (size_t i = 0; i < x.inputs().size() - 1; ++i) {
           value += fmt::format("tmp_{}, ", namer.get_name(x.inputs()[i]));
@@ -181,28 +120,7 @@ struct FusedKernelBuilder {
 
     // Write output.
     for (const auto& x : outputs) {
-      os += fmt::format("    vec_{0}[i] = tmp_{0};\n", namer.get_name(x));
-    }
-
-    // End of work loop
-    if (!contiguous) {
-      os += "\n";
-      for (size_t i = 0; i < inputs.size(); ++i) {
-        const auto& x = inputs[i];
-        const std::string& xname = namer.get_name(x);
-        if (is_scalar(x) || is_constant(i)) {
-          continue;
-        }
-        os += fmt::format("    {0}_idx += {0}_strides[NDIM - 1];\n", xname);
-      }
-    }
-    os += "  }\n";
-
-    // Store the output to global memory
-    for (const auto& x : outputs) {
-      os += fmt::format(
-          "  store_vector({0} + index, 0, vec_{0}, size - index);\n",
-          namer.get_name(x));
+      os += fmt::format("  {0}[index] = tmp_{0};\n", namer.get_name(x));
     }
 
     os += "}\n";
@@ -228,15 +146,6 @@ void Compiled::eval_gpu(
   nvtx3::scoped_range r("Compiled::eval_gpu");
   auto& s = stream();
 
-  // Determine the work per thread for the vectorized reads/writes. We take it
-  // as 16 over the max itemsize for the outputs. Another heuristic could be
-  // over the max itemsize of all arrays.
-  int max_size = 1;
-  for (const auto& x : outputs) {
-    max_size = (max_size > x.itemsize()) ? max_size : x.itemsize();
-  }
-  int work_per_thread = 16 / max_size;
-
   cu::JitModule& mod = cu::get_jit_module(s.device, lib_name(), [&]() {
     // Build source code.
     cu::FusedKernelBuilder builder{
@@ -249,26 +158,17 @@ void Compiled::eval_gpu(
     builder.build("_strided", false);
     builder.os += "\n} // namespace mlx::core::cu\n";
     // Build kernel names.
-    std::vector<std::string> kernel_names;
-    kernel_names.push_back(fmt::format(
-        "mlx::core::cu::{}_contiguous<uint32_t, {}>",
-        lib_name(),
-        work_per_thread));
-    kernel_names.push_back(fmt::format(
-        "mlx::core::cu::{}_contiguous<int64_t, {}>",
-        lib_name(),
-        work_per_thread));
-    for (auto wpt : std::array<int, 2>{1, work_per_thread}) {
+    std::vector<std::string> kernel_names = {
+        fmt::format("mlx::core::cu::{}_contiguous<uint32_t>", lib_name()),
+        fmt::format("mlx::core::cu::{}_contiguous<int64_t>", lib_name()),
+    };
     for (int i = 1; i <= MAX_NDIM; ++i) {
       kernel_names.push_back(fmt::format(
-            "mlx::core::cu::{}_strided<{}, uint32_t, {}>", lib_name(), i, wpt));
-        kernel_names.push_back(fmt::format(
-            "mlx::core::cu::{}_strided<{}, int64_t, {}>", lib_name(), i, wpt));
+          "mlx::core::cu::{}_strided<{}, uint32_t>", lib_name(), i));
+      kernel_names.push_back(
+          fmt::format("mlx::core::cu::{}_strided<{}, int64_t>", lib_name(), i));
     }
-    }
-
-    return std::make_tuple(
-        false, std::move(builder.os), std::move(kernel_names));
+    return std::make_pair(std::move(builder.os), std::move(kernel_names));
   });
 
   // Collapse contiguous dims to route to a faster kernel if possible. Also
@@ -309,20 +209,13 @@ void Compiled::eval_gpu(
     args.append<uint32_t>(outputs[0].data_size());
   }
 
-  // Choose work per thread
-  if (!contiguous && shape.back() % work_per_thread != 0) {
-    work_per_thread = 1;
-  }
-
   // Launch kernel.
   const char* index_type = large ? "int64_t" : "uint32_t";
   std::string kernel_name = fmt::format("mlx::core::cu::{}", lib_name());
   if (contiguous) {
-    kernel_name +=
-        fmt::format("_contiguous<{}, {}>", index_type, work_per_thread);
+    kernel_name += fmt::format("_contiguous<{}>", index_type);
   } else {
-    kernel_name += fmt::format(
-        "_strided<{}, {}, {}>", shape.size(), index_type, work_per_thread);
+    kernel_name += fmt::format("_strided<{}, {}>", shape.size(), index_type);
   }
   auto& encoder = cu::get_command_encoder(s);
   for (const auto& in : inputs) {
@@ -333,9 +226,8 @@ void Compiled::eval_gpu(
   }
 
   auto kernel = mod.get_kernel(kernel_name);
-  auto [num_blocks, block_dims] =
-      get_launch_args(outputs[0], large, work_per_thread);
-  encoder.add_kernel_node(kernel, num_blocks, block_dims, 0, args.args());
+  auto [num_blocks, block_dims] = get_launch_args(kernel, outputs[0], large);
+  encoder.add_kernel_node(kernel, num_blocks, block_dims, args.args());
 }
 
 } // namespace mlx::core

mlx/backend/cuda/conv.cpp

@@ -1,418 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/conv/conv.h"
-#include "mlx/backend/cuda/cudnn_utils.h"
-#include "mlx/backend/cuda/device.h"
-#include "mlx/backend/cuda/lru_cache.h"
-#include "mlx/backend/gpu/copy.h"
-#include "mlx/primitives.h"
-
-#include <nvtx3/nvtx3.hpp>
-
-#include <cassert>
-
-namespace mlx::core {
-
-namespace {
-
-// Alias for better readability.
-#define CONV_FORWARD CUDNN_BACKEND_OPERATION_CONVOLUTION_FORWARD_DESCRIPTOR
-#define CONV_BACKWARD_INPUT \
-  CUDNN_BACKEND_OPERATION_CONVOLUTION_BACKWARD_DATA_DESCRIPTOR
-#define CONV_BACKWARD_WEIGHT \
-  CUDNN_BACKEND_OPERATION_CONVOLUTION_BACKWARD_FILTER_DESCRIPTOR
-
-// Custom placeholder representing fallback kernel.
-#define CONV_FALLBACK static_cast<cudnnBackendDescriptorType_t>(-1)
-
-struct ConvCacheKey {
-  int device_id;
-  cudnnDataType_t cudnn_dtype;
-  std::array<int, MAX_NDIM> input_shape;
-  std::array<int, MAX_NDIM> weight_shape;
-  std::array<int, MAX_NDIM> stride;
-  std::array<int, MAX_NDIM> padding_lo;
-  std::array<int, MAX_NDIM> padding_hi;
-  std::array<int, MAX_NDIM> dilation;
-  int groups;
-  bool flip;
-  uint8_t input_alignment;
-  uint8_t weight_alignment;
-  uint8_t output_alignment;
-};
-
-auto& conv_cache() {
-  static LRUBytesKeyCache<
-      ConvCacheKey,
-      std::pair<
-          cudnnBackendDescriptorType_t,
-          std::optional<cudnn_frontend::ExecutionPlan>>>
-      cache(/* capacity */ 128);
-  return cache;
-}
-
-auto get_conv_op_settings(
-    cudnnBackendDescriptorType_t backend_type,
-    array& x,
-    array& w,
-    array& y,
-    const std::vector<int>& kernel_strides,
-    const std::vector<int>& padding_lo_,
-    const std::vector<int>& padding_hi_,
-    const std::vector<int>& kernel_dilation,
-    const std::vector<int>& input_dilation) {
-  auto padding_lo = convert_vector<int64_t>(padding_lo_);
-  auto padding_hi = convert_vector<int64_t>(padding_hi_);
-
-  if (backend_type == CONV_BACKWARD_INPUT) {
-    for (int i = 0; i < padding_lo.size(); ++i) {
-      int wt_size = 1 + kernel_dilation[i] * (w.shape(1 + i) - 1);
-      padding_lo[i] = wt_size - padding_lo[i] - 1;
-      int in_size = 1 + kernel_strides[i] * (x.shape(1 + i) - 1);
-      int out_size = 1 + input_dilation[i] * (y.shape(1 + i) - 1);
-      padding_hi[i] = out_size - in_size + padding_hi[i];
-    }
-    return std::make_tuple(
-        convert_vector<int64_t>(input_dilation),
-        std::move(padding_lo),
-        std::move(padding_hi),
-        convert_vector<int64_t>(kernel_dilation));
-
-  } else if (backend_type == CONV_BACKWARD_WEIGHT) {
-    padding_hi = padding_lo;
-    return std::make_tuple(
-        convert_vector<int64_t>(kernel_dilation),
-        std::move(padding_lo),
-        std::move(padding_hi),
-        convert_vector<int64_t>(kernel_strides));
-
-  } else {
-    return std::make_tuple(
-        convert_vector<int64_t>(kernel_strides),
-        std::move(padding_lo),
-        std::move(padding_hi),
-        convert_vector<int64_t>(kernel_dilation));
-  }
-}
-
-std::optional<cudnn_frontend::OperationGraph> build_conv_op_graph(
-    cu::CommandEncoder& encoder,
-    cudnnBackendDescriptorType_t backend_type,
-    Dtype dtype,
-    array& x,
-    array& w,
-    array& y,
-    const SmallVector<int64_t>& stride,
-    const SmallVector<int64_t>& padding_lo,
-    const SmallVector<int64_t>& padding_hi,
-    const SmallVector<int64_t>& dilation) {
-  try {
-    auto compute_dtype = (dtype == float16 || dtype == bfloat16)
-        ? CUDNN_DATA_FLOAT
-        : dtype_to_cudnn_type(dtype);
-    auto conv_desc = cudnn_frontend::ConvDescBuilder()
-                         .setDataType(compute_dtype)
-                         .setMathMode(CUDNN_CROSS_CORRELATION)
-                         .setNDims(stride.size())
-                         .setStrides(stride.size(), stride.data())
-                         .setPrePadding(padding_lo.size(), padding_lo.data())
-                         .setPostPadding(padding_hi.size(), padding_hi.data())
-                         .setDilation(dilation.size(), dilation.data())
-                         .build();
-
-    auto op = cudnn_frontend::OperationBuilder(backend_type)
-                  .setxDesc(build_cudnn_tensor_nchw('x', x))
-                  .setwDesc(build_cudnn_tensor_nchw('w', w))
-                  .setyDesc(build_cudnn_tensor_nchw('y', y))
-                  .setcDesc(conv_desc)
-                  .build();
-
-    std::array<cudnn_frontend::Operation const*, 1> ops = {&op};
-    return cudnn_frontend::OperationGraphBuilder()
-        .setHandle(encoder.device().cudnn_handle())
-        .setOperationGraph(ops.size(), ops.data())
-        .build();
-  } catch (cudnn_frontend::cudnnException& error) {
-    if (error.getCudnnStatus() != CUDNN_STATUS_BAD_PARAM) {
-      throw;
-    }
-    return std::nullopt;
-  }
-}
-
-// Transpose from (C_out, H, W, C_in / groups) to (C_in, H, W, C_out / groups).
-array group_transpose(
-    const array& x,
-    int groups,
-    int group_dim,
-    int axis1,
-    int axis2,
-    Stream s) {
-  if (groups == 1) {
-    return swapaxes_in_eval(x, axis1, axis2);
-  }
-  int ndim = x.ndim();
-  if (group_dim < 0) {
-    group_dim += ndim;
-  }
-  if (axis1 < 0) {
-    axis1 += ndim;
-  }
-  if (axis2 < 0) {
-    axis2 += ndim;
-  }
-  if (group_dim <= axis1) {
-    axis1 += 1;
-  }
-  if (group_dim <= axis2) {
-    axis2 += 1;
-  }
-  auto shape = x.shape();
-  shape.insert(shape.begin() + group_dim, groups);
-  shape[group_dim + 1] = shape[group_dim + 1] / groups;
-  array x_trans = reshape_in_eval(x, std::move(shape), s);
-  x_trans = swapaxes_in_eval(x_trans, axis1, axis2);
-  x_trans = flatten_in_eval(x_trans, group_dim, group_dim + 1, s);
-  return x_trans;
-}
-
-// Do necessary transposes and copies to prepare the inputs and outputs for
-// building the cuDNN conv op. It is safe to be called multiple times in one
-// eval_gpu, with cost of possible redundant copies.
-std::tuple<array, array, array> prepare_args(
-    cu::CommandEncoder& encoder,
-    cudnnBackendDescriptorType_t backend_type,
-    array in,
-    array wt,
-    array out,
-    int groups,
-    Stream s) {
-  // Transpose the args depending on the backend type.
-  // TODO: Handle groups.
-  if (backend_type == CONV_BACKWARD_INPUT) {
-    wt = group_transpose(wt, groups, 0, 0, -1, s);
-  } else if (backend_type == CONV_BACKWARD_WEIGHT) {
-    in = group_transpose(in, groups, -1, 0, -1, s);
-    wt = swapaxes_in_eval(wt, 0, -1);
-    // Create a contiguous array that shares the data with |out|, but with dim
-    // C_in and C_out swapped.
-    Shape shape(out.shape());
-    std::swap(shape.front(), shape.back());
-    Strides strides(shape.size(), 1);
-    for (int i = shape.size() - 2; i >= 0; --i) {
-      strides[i] = shape[i + 1] * strides[i + 1];
-    }
-    array intermediate(std::move(shape), out.dtype(), nullptr, {});
-    intermediate.copy_shared_buffer(
-        out, std::move(strides), {true, true, false}, out.data_size());
-    out = intermediate;
-  }
-
-  // cuDNN requires contiguous input.
-  if (!in.flags().row_contiguous) {
-    in = contiguous_copy_gpu(in, s);
-    encoder.add_temporary(in);
-  }
-  if (!wt.flags().row_contiguous) {
-    wt = contiguous_copy_gpu(wt, s);
-    encoder.add_temporary(wt);
-  }
-
-  return {std::move(in), std::move(wt), std::move(out)};
-}
-
-// Get the x/w/y args from the in/wt/out args depending on backend type.
-inline std::tuple<array&, array&, array&> dispatch_args(
-    cudnnBackendDescriptorType_t backend_type,
-    array& in,
-    array& wt,
-    array& out) {
-  switch (backend_type) {
-    case CONV_BACKWARD_INPUT:
-      return {out, wt, in};
-    case CONV_BACKWARD_WEIGHT:
-      return {in, out, wt};
-    default:
-      return {in, wt, out};
-  }
-}
-
-// Register inputs and outputs before actually running conv op. Can only be
-// called once per eval_gpu.
-void register_args(
-    cu::CommandEncoder& encoder,
-    cudnnBackendDescriptorType_t backend_type,
-    array& in,
-    array& wt,
-    array& intermediate_out,
-    array& final_out) {
-  encoder.set_input_array(in);
-  encoder.set_input_array(wt);
-  encoder.set_output_array(final_out);
-
-  if (backend_type == CONV_BACKWARD_WEIGHT) {
-    // Turn |out| into a strided array, which will have C_in and C_out swapped
-    // in vjp and the final |grad_weight| will then be contiguous.
-    Strides strides = intermediate_out.strides();
-    std::swap(strides.front(), strides.back());
-    final_out.copy_shared_buffer(
-        intermediate_out,
-        std::move(strides),
-        {false, false, false},
-        intermediate_out.data_size());
-  }
-}
-
-} // namespace
-
-void Convolution::eval_gpu(const std::vector<array>& inputs, array& out_) {
-  nvtx3::scoped_range r("Convolution::eval_gpu");
-  if (out_.size() == 0) {
-    return;
-  }
-
-  assert(inputs.size() == 2);
-  array in = inputs[0];
-  array wt = inputs[1];
-  array out = out_;
-  out.set_data(allocator::malloc(out.nbytes()));
-  Dtype dtype = out.dtype();
-
-  auto& s = stream();
-  auto& encoder = cu::get_command_encoder(s);
-
-  // Search cache.
-  ConvCacheKey cache_key{
-      encoder.device().cuda_device(),
-      dtype_to_cudnn_type(dtype),
-      vector_key(in.shape()),
-      vector_key(wt.shape()),
-      vector_key(kernel_strides_),
-      vector_key(padding_lo_),
-      vector_key(padding_hi_),
-      vector_key(kernel_dilation_),
-      groups_,
-      flip_,
-      get_alignment(in),
-      get_alignment(wt),
-      get_alignment(out)};
-  if (auto it = conv_cache().find(cache_key); it != conv_cache().end()) {
-    auto& [backend_type, plan] = it->second;
-    if (plan) {
-      // Run cached plan.
-      std::tie(in, wt, out) =
-          prepare_args(encoder, backend_type, in, wt, out, groups_, s);
-      register_args(encoder, backend_type, in, wt, out, out_);
-      auto [x, w, y] = dispatch_args(backend_type, in, wt, out);
-      if (!encode_cudnn_plan(encoder, *plan, {'x', 'w', 'y'}, x, w, y)) {
-        throw std::runtime_error("[conv] Cached plan failed to execute.");
-      }
-    } else {
-      // Run fallback kernel.
-      gemm_conv(
-          encoder,
-          in,
-          wt,
-          out,
-          kernel_strides_,
-          padding_lo_,
-          kernel_dilation_,
-          input_dilation_,
-          groups_,
-          flip_,
-          s);
-    }
-    return;
-  }
-
-  // There is no reliable way to deduce the proper cuDNN backend for the
-  // convolution, so we make a best guess and then try.
-  SmallVector<cudnnBackendDescriptorType_t, 2> try_backends;
-  if (flip_) {
-    // When weight is flipped, we assume it is backward input convolution.
-    try_backends.push_back(CONV_BACKWARD_INPUT);
-  } else {
-    // Otherwise it could be backward weight convolution or forward convolution,
-    // mathematically there is no difference so we have to use heuristics.
-    // Empirically backward convolutions have large kernel dimensions, and
-    // usually have |in| and |wt| transposed.
-    if (!in.flags().row_contiguous && !wt.flags().row_contiguous &&
-        wt.shape(2) > out.shape(2)) {
-      try_backends = {CONV_BACKWARD_WEIGHT, CONV_FORWARD};
-    } else {
-      try_backends = {CONV_FORWARD, CONV_BACKWARD_WEIGHT};
-    }
-  }
-
-  // Try to build op graph.
-  cudnnBackendDescriptorType_t backend_type;
-  std::optional<cudnn_frontend::OperationGraph> op_graph;
-  for (auto try_backend : try_backends) {
-    auto [in_copy, wt_copy, out_copy] =
-        prepare_args(encoder, try_backend, in, wt, out, groups_, s);
-    auto [x, w, y] = dispatch_args(try_backend, in_copy, wt_copy, out_copy);
-    auto [stride, padding_lo, padding_hi, dilation] = get_conv_op_settings(
-        try_backend,
-        x,
-        w,
-        y,
-        kernel_strides_,
-        padding_lo_,
-        padding_hi_,
-        kernel_dilation_,
-        input_dilation_);
-    op_graph = build_conv_op_graph(
-        encoder,
-        try_backend,
-        dtype,
-        x,
-        w,
-        y,
-        stride,
-        padding_lo,
-        padding_hi,
-        dilation);
-    if (op_graph) {
-      backend_type = try_backend;
-      in = std::move(in_copy);
-      wt = std::move(wt_copy);
-      out = std::move(out_copy);
-      break;
-    }
-  }
-
-  if (op_graph) {
-    // Setup inputs and outputs.
-    register_args(encoder, backend_type, in, wt, out, out_);
-
-    // Find a plan for the graph and execute it.
-    auto plan = find_cudnn_plan_from_op_graph(
-        encoder.device().cudnn_handle(), backend_type, dtype, *op_graph);
-    if (!plan) {
-      throw std::runtime_error("[conv] Unable to find an execution plan.");
-    }
-    auto [x, w, y] = dispatch_args(backend_type, in, wt, out);
-    if (encode_cudnn_plan(encoder, *plan, {'x', 'w', 'y'}, x, w, y)) {
-      conv_cache().emplace(
-          cache_key, std::make_pair(backend_type, std::move(*plan)));
-      return;
-    }
-  }
-
-  // Use fallback kernel for settings not supported by cuDNN.
-  gemm_conv(
-      encoder,
-      in,
-      wt,
-      out,
-      kernel_strides_,
-      padding_lo_,
-      kernel_dilation_,
-      input_dilation_,
-      groups_,
-      flip_,
-      s);
-  conv_cache().emplace(cache_key, std::make_pair(CONV_FALLBACK, std::nullopt));
-}
-
-} // namespace mlx::core

mlx/backend/cuda/conv/conv.h

@@ -1,126 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#pragma once
-
-#include "mlx/backend/cuda/device.h"
-#include "mlx/backend/gpu/copy.h"
-
-namespace mlx::core {
-
-template <int NDIM>
-struct ConvParams {
-  int N; // Batch size
-  int C; // In channels
-  int O; // Out channels
-  int strides[NDIM];
-  int padding[NDIM];
-  int kernel_dilation[NDIM];
-  int input_dilation[NDIM];
-  int groups;
-  bool flip;
-  int in_spatial_dims[NDIM];
-  int wt_spatial_dims[NDIM];
-  int out_spatial_dims[NDIM];
-  int64_t in_strides[NDIM + 2];
-
-  ConvParams(
-      const array& in,
-      const array& wt,
-      const array& out,
-      const std::vector<int>& strides,
-      const std::vector<int>& padding,
-      const std::vector<int>& kernel_dilation,
-      const std::vector<int>& input_dilation,
-      int groups,
-      bool flip)
-      : N(in.shape(0)),
-        C(in.shape(-1)),
-        O(wt.shape(0)),
-        groups(groups),
-        flip(flip) {
-    std::copy_n(strides.begin(), NDIM, this->strides);
-    std::copy_n(padding.begin(), NDIM, this->padding);
-    std::copy_n(kernel_dilation.begin(), NDIM, this->kernel_dilation);
-    std::copy_n(input_dilation.begin(), NDIM, this->input_dilation);
-    std::copy_n(in.shape().begin() + 1, NDIM, this->in_spatial_dims);
-    std::copy_n(wt.shape().begin() + 1, NDIM, this->wt_spatial_dims);
-    std::copy_n(out.shape().begin() + 1, NDIM, this->out_spatial_dims);
-    std::copy_n(in.strides().begin(), NDIM + 2, this->in_strides);
-  }
-};
-
-void gemm_grouped_conv(
-    cu::CommandEncoder& encoder,
-    const array& in,
-    const array& wt,
-    array& out,
-    const std::vector<int>& strides,
-    const std::vector<int>& padding,
-    const std::vector<int>& kernel_dilation,
-    const std::vector<int>& input_dilation,
-    int groups,
-    bool flip,
-    Stream s);
-
-void gemm_conv(
-    cu::CommandEncoder& encoder,
-    const array& in,
-    const array& wt,
-    array& out,
-    const std::vector<int>& strides,
-    const std::vector<int>& padding,
-    const std::vector<int>& kernel_dilation,
-    const std::vector<int>& input_dilation,
-    bool flip,
-    Stream s);
-
-inline void gemm_conv(
-    cu::CommandEncoder& encoder,
-    array in,
-    array wt,
-    array& out,
-    const std::vector<int>& strides,
-    const std::vector<int>& padding,
-    const std::vector<int>& kernel_dilation,
-    const std::vector<int>& input_dilation,
-    int groups,
-    bool flip,
-    Stream s) {
-  if (!in.flags().row_contiguous) {
-    in = contiguous_copy_gpu(in, s);
-    encoder.add_temporary(in);
-  }
-  if (!wt.flags().row_contiguous) {
-    wt = contiguous_copy_gpu(wt, s);
-    encoder.add_temporary(wt);
-  }
-
-  if (groups == 1) {
-    gemm_conv(
-        encoder,
-        in,
-        wt,
-        out,
-        strides,
-        padding,
-        kernel_dilation,
-        input_dilation,
-        flip,
-        s);
-  } else {
-    gemm_grouped_conv(
-        encoder,
-        in,
-        wt,
-        out,
-        strides,
-        padding,
-        kernel_dilation,
-        input_dilation,
-        groups,
-        flip,
-        s);
-  }
-}
-
-} // namespace mlx::core

mlx/backend/cuda/conv/gemm_conv.cu

@@ -1,217 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/conv/conv.h"
-#include "mlx/backend/cuda/gemms/cublas_gemm.h"
-#include "mlx/backend/cuda/kernel_utils.cuh"
-#include "mlx/dtype_utils.h"
-
-#include <cooperative_groups.h>
-
-namespace mlx::core {
-
-namespace cu {
-
-namespace cg = cooperative_groups;
-
-template <typename T, int NDIM>
-__global__ void naive_unfold_nd(
-    const T* in,
-    T* out,
-    int filter_size,
-    int out_pixels,
-    const __grid_constant__ ConvParams<NDIM> params) {
-  auto block = cg::this_thread_block();
-  auto tid = block.group_index();
-  auto lid = block.thread_index();
-
-  int index_batch = tid.z / out_pixels; // [0, N)
-  int index_out_spatial = tid.z % out_pixels; // [0, H_out * W_out)
-  int index_wt_spatial =
-      tid.x * block.dim_threads().x + lid.x; // [0, H_wt * W_wt)
-
-  if (index_wt_spatial >= filter_size / params.C) {
-    return;
-  }
-
-  in += tid.y; // [0, C)
-  out += tid.z * filter_size + index_wt_spatial * params.C + tid.y;
-
-  bool valid = index_batch < params.N;
-
-  // Get the coordinates in input.
-  int index_in[NDIM] = {};
-#pragma unroll
-  for (int i = NDIM - 1; i >= 0; --i) {
-    int index_out = index_out_spatial % params.out_spatial_dims[i];
-    int index_wt = index_wt_spatial % params.wt_spatial_dims[i];
-
-    if (params.flip) {
-      index_wt = params.wt_spatial_dims[i] - index_wt - 1;
-    }
-
-    int index = index_out * params.strides[i] - params.padding[i] +
-        index_wt * params.kernel_dilation[i];
-    int index_max =
-        1 + params.input_dilation[i] * (params.in_spatial_dims[i] - 1);
-
-    valid &= (index >= 0) && (index < index_max) &&
-        (index % params.input_dilation[i] == 0);
-
-    index_in[i] = index / params.input_dilation[i];
-
-    index_out_spatial /= params.out_spatial_dims[i];
-    index_wt_spatial /= params.wt_spatial_dims[i];
-  }
-
-  if (valid) {
-    int in_offset = index_batch * params.in_strides[0];
-#pragma unroll
-    for (int i = 0; i < NDIM; ++i) {
-      in_offset += index_in[i] * params.in_strides[i + 1];
-    }
-    *out = in[in_offset];
-  } else {
-    *out = T{0};
-  }
-}
-
-} // namespace cu
-
-template <int NDIM>
-array unfold_inputs_nd(
-    cu::CommandEncoder& encoder,
-    const array& in,
-    int mat_M,
-    int mat_K,
-    int mat_N,
-    ConvParams<NDIM>& params) {
-  array unfolded({mat_M, mat_K}, in.dtype(), nullptr, {});
-  unfolded.set_data(allocator::malloc(unfolded.nbytes()));
-  encoder.add_temporary(unfolded);
-
-  int filter_size = params.C;
-#pragma unroll
-  for (int i = 0; i < NDIM; ++i) {
-    filter_size *= params.wt_spatial_dims[i];
-  }
-
-  int out_pixels = 1;
-#pragma unroll
-  for (int i = 0; i < NDIM; ++i) {
-    out_pixels *= params.out_spatial_dims[i];
-  }
-
-  int wt_spatial_size = mat_K / params.C;
-  dim3 block_dims;
-  block_dims.x = std::min(std::max(wt_spatial_size, 32), 1024);
-  dim3 num_blocks;
-  num_blocks.x = cuda::ceil_div(wt_spatial_size, block_dims.x);
-  num_blocks.y = params.C;
-  num_blocks.z = mat_M;
-
-  encoder.set_input_array(in);
-  encoder.set_output_array(unfolded);
-  dispatch_float_types(in.dtype(), "unfold", [&](auto type_tag) {
-    using DataType = cuda_type_t<MLX_GET_TYPE(type_tag)>;
-    encoder.add_kernel_node(
-        cu::naive_unfold_nd<DataType, NDIM>,
-        num_blocks,
-        block_dims,
-        0,
-        in.data<DataType>(),
-        unfolded.data<DataType>(),
-        filter_size,
-        out_pixels,
-        params);
-  });
-
-  return unfolded;
-}
-
-template <int NDIM>
-void gemm_conv_nd(
-    cu::CommandEncoder& encoder,
-    const array& in,
-    const array& wt,
-    array& out,
-    ConvParams<NDIM>& params,
-    Stream s) {
-  // Get gemm shapes.
-  int mat_M = out.size() / params.O; // N * H_out * W_out
-  int mat_K = wt.size() / params.O; // C * H_wt * W_wt
-  int mat_N = params.O; // O
-
-  // Unfold input to (N * H_out * W_out, C * H_wt * W_wt) for gemm.
-  array in_unfolded =
-      unfold_inputs_nd<NDIM>(encoder, in, mat_M, mat_K, mat_N, params);
-
-  // Reshape weight to (C * H_wt * W_wt, O) for gemm.
-  array wt_reshaped({mat_K, mat_N}, wt.dtype(), nullptr, {});
-  wt_reshaped.copy_shared_buffer(
-      wt,
-      {1, mat_K},
-      {false, false, /* col_contiguous */ true},
-      wt.data_size());
-
-  // Single batch.
-  Shape batch_shape{1};
-  Strides a_batch_strides{0};
-  Strides b_batch_strides{0};
-
-  // Run matmul.
-  CublasGemm gemm(
-      encoder.device(),
-      in.dtype(),
-      false, // a_transposed
-      mat_M, // a_rows
-      mat_K, // a_cols
-      mat_K, // lda
-      true, // b_transposed
-      mat_K, // b_rows
-      mat_N, // b_cols
-      mat_K, // ldb
-      batch_shape.back(),
-      a_batch_strides.back(),
-      b_batch_strides.back());
-  gemm.run(
-      encoder,
-      out,
-      in_unfolded,
-      wt_reshaped,
-      batch_shape,
-      a_batch_strides,
-      b_batch_strides);
-}
-
-void gemm_conv(
-    cu::CommandEncoder& encoder,
-    const array& in,
-    const array& wt,
-    array& out,
-    const std::vector<int>& strides,
-    const std::vector<int>& padding,
-    const std::vector<int>& kernel_dilation,
-    const std::vector<int>& input_dilation,
-    bool flip,
-    Stream s) {
-  int conv_ndim = in.ndim() - 2;
-  if (conv_ndim < 1 || conv_ndim > 3) {
-    throw std::runtime_error(
-        fmt::format("[conv] Unsupported gemm_conv for {}D conv.", conv_ndim));
-  }
-  dispatch_1_2_3(conv_ndim, [&](auto ndim_constant) {
-    ConvParams<ndim_constant()> params(
-        in,
-        wt,
-        out,
-        strides,
-        padding,
-        kernel_dilation,
-        input_dilation,
-        1, // groups
-        flip);
-    gemm_conv_nd<ndim_constant()>(encoder, in, wt, out, params, s);
-  });
-}
-
-} // namespace mlx::core

mlx/backend/cuda/conv/gemm_grouped_conv.cu

@@ -1,231 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/conv/conv.h"
-#include "mlx/backend/cuda/gemms/cublas_gemm.h"
-#include "mlx/backend/cuda/kernel_utils.cuh"
-#include "mlx/dtype_utils.h"
-
-#include <cooperative_groups.h>
-
-namespace mlx::core {
-
-namespace cu {
-
-namespace cg = cooperative_groups;
-
-template <typename T, int NDIM>
-__global__ void naive_grouped_unfold_transpose_nd(
-    const T* in,
-    T* out,
-    int filter_size,
-    int out_pixels,
-    const __grid_constant__ ConvParams<NDIM> params) {
-  auto block = cg::this_thread_block();
-  auto tid = block.group_index();
-  auto lid = block.thread_index();
-
-  int index_batch = tid.z / out_pixels; // [0, N)
-  int index_out_spatial = tid.z % out_pixels; // [0, H_out * W_out)
-  int index_wt_spatial =
-      tid.x * block.dim_threads().x + lid.x; // [0, H_wt * W_wt)
-
-  if (index_wt_spatial >= filter_size / params.C) {
-    return;
-  }
-
-  in += tid.y; // [0, C)
-  out += tid.z * filter_size + tid.y * (filter_size / params.C);
-
-  bool valid = index_batch < params.N;
-
-  // Get the coordinates in input.
-  int index_in[NDIM] = {};
-  int wt_stride = 1;
-#pragma unroll
-  for (int i = NDIM - 1; i >= 0; --i) {
-    int index_out = index_out_spatial % params.out_spatial_dims[i];
-    int index_wt = index_wt_spatial % params.wt_spatial_dims[i];
-    out += index_wt * wt_stride;
-
-    if (params.flip) {
-      index_wt = params.wt_spatial_dims[i] - index_wt - 1;
-    }
-
-    int index = index_out * params.strides[i] - params.padding[i] +
-        index_wt * params.kernel_dilation[i];
-    int index_max =
-        1 + params.input_dilation[i] * (params.in_spatial_dims[i] - 1);
-
-    valid &= (index >= 0) && (index < index_max) &&
-        (index % params.input_dilation[i] == 0);
-
-    index_in[i] = index / params.input_dilation[i];
-
-    index_out_spatial /= params.out_spatial_dims[i];
-    index_wt_spatial /= params.wt_spatial_dims[i];
-    wt_stride *= params.wt_spatial_dims[i];
-  }
-
-  if (valid) {
-    int in_offset = index_batch * params.in_strides[0];
-#pragma unroll
-    for (int i = 0; i < NDIM; ++i) {
-      in_offset += index_in[i] * params.in_strides[i + 1];
-    }
-    *out = in[in_offset];
-  } else {
-    *out = T{0};
-  }
-}
-
-} // namespace cu
-
-template <int NDIM>
-array grouped_unfold_transpose_inputs_nd(
-    cu::CommandEncoder& encoder,
-    const array& in,
-    int mat_M,
-    int mat_K,
-    int mat_N,
-    ConvParams<NDIM>& params) {
-  array unfolded({mat_M, mat_K * params.groups}, in.dtype(), nullptr, {});
-  unfolded.set_data(allocator::malloc(unfolded.nbytes()));
-  encoder.add_temporary(unfolded);
-
-  int filter_size = params.C;
-#pragma unroll
-  for (int i = 0; i < NDIM; ++i) {
-    filter_size *= params.wt_spatial_dims[i];
-  }
-
-  int out_pixels = 1;
-#pragma unroll
-  for (int i = 0; i < NDIM; ++i) {
-    out_pixels *= params.out_spatial_dims[i];
-  }
-
-  int wt_spatial_size = (mat_K * params.groups) / params.C;
-  dim3 block_dims;
-  block_dims.x = std::min(std::max(wt_spatial_size, 32), 1024);
-  dim3 num_blocks;
-  num_blocks.x = cuda::ceil_div(wt_spatial_size, block_dims.x);
-  num_blocks.y = params.C;
-  num_blocks.z = mat_M;
-
-  encoder.set_input_array(in);
-  encoder.set_output_array(unfolded);
-  dispatch_float_types(in.dtype(), "unfold", [&](auto type_tag) {
-    using DataType = cuda_type_t<MLX_GET_TYPE(type_tag)>;
-    encoder.add_kernel_node(
-        cu::naive_grouped_unfold_transpose_nd<DataType, NDIM>,
-        num_blocks,
-        block_dims,
-        0,
-        in.data<DataType>(),
-        unfolded.data<DataType>(),
-        filter_size,
-        out_pixels,
-        params);
-  });
-
-  return unfolded;
-}
-
-template <int NDIM>
-void gemm_grouped_conv_nd(
-    cu::CommandEncoder& encoder,
-    const array& in,
-    const array& wt,
-    array& out,
-    ConvParams<NDIM>& params,
-    Stream s) {
-  // Get gemm shapes.
-  int C_per_group = params.C / params.groups;
-  int O_per_group = params.O / params.groups;
-  int mat_M = out.size() / params.O; // N * H_out * W_out
-  int mat_K = wt.size() / params.O; // C_per_group * H_wt * W_wt
-  int mat_N = O_per_group; // O_per_group
-
-  // Unfold input to (N * H_out * W_out, C * H_wt * W_wt) for gemm.
-  array in_unfolded = grouped_unfold_transpose_inputs_nd<NDIM>(
-      encoder, in, mat_M, mat_K, mat_N, params);
-
-  // Reshape weight to (O, C_per_group, H_wt * W_wt) for gemm.
-  int wt_spatial_size = (wt.size() / wt.shape(0)) / wt.shape(-1);
-  array wt_view(
-      {params.O, C_per_group, wt_spatial_size}, wt.dtype(), nullptr, {});
-  wt_view.copy_shared_buffer(
-      wt, {wt.strides(0), 1, C_per_group}, wt.flags(), wt.size());
-  array wt_reshaped = contiguous_copy_gpu(wt_view, s);
-
-  // Batch with size of groups.
-  Shape batch_shape{params.groups};
-  Strides a_batch_strides{mat_K};
-  Strides b_batch_strides{mat_N * mat_K};
-
-  // Run matmul.
-  CublasGemm gemm(
-      encoder.device(),
-      in.dtype(),
-      false, // a_transposed
-      mat_M, // a_rows
-      mat_K, // a_cols
-      mat_K * params.groups, // lda
-      true, // b_transposed
-      mat_K, // b_rows
-      mat_N, // b_cols
-      mat_K, // ldb
-      batch_shape.back(),
-      a_batch_strides.back(),
-      b_batch_strides.back());
-  gemm.set_out(
-      out.dtype(),
-      false, // out_transposed
-      mat_M, // out_rows
-      mat_N, // out_cols
-      mat_N * params.groups, // out_ld
-      params.groups, // batch_count
-      mat_N); // batch_stride
-  gemm.run(
-      encoder,
-      out,
-      in_unfolded,
-      wt_reshaped,
-      batch_shape,
-      a_batch_strides,
-      b_batch_strides);
-}
-
-void gemm_grouped_conv(
-    cu::CommandEncoder& encoder,
-    const array& in,
-    const array& wt,
-    array& out,
-    const std::vector<int>& strides,
-    const std::vector<int>& padding,
-    const std::vector<int>& kernel_dilation,
-    const std::vector<int>& input_dilation,
-    int groups,
-    bool flip,
-    Stream s) {
-  int conv_ndim = in.ndim() - 2;
-  if (conv_ndim < 1 || conv_ndim > 3) {
-    throw std::runtime_error(
-        fmt::format("[conv] Unsupported gemm_conv for {}D conv.", conv_ndim));
-  }
-  dispatch_1_2_3(conv_ndim, [&](auto ndim_constant) {
-    ConvParams<ndim_constant()> params(
-        in,
-        wt,
-        out,
-        strides,
-        padding,
-        kernel_dilation,
-        input_dilation,
-        groups,
-        flip);
-    gemm_grouped_conv_nd<ndim_constant()>(encoder, in, wt, out, params, s);
-  });
-}
-
-} // namespace mlx::core

mlx/backend/cuda/copy/copy_contiguous.cu

@@ -22,7 +22,7 @@ __global__ void copy_s(const In* in, Out* out, IdxT size) {
     AlignedVector<Out, N_READS> out_vec;
 #pragma unroll
     for (int i = 0; i < N_READS; ++i) {
-      out_vec[i] = cast_to<Out>(in[0]);
+      out_vec.val[i] = cast_to<Out>(in[0]);
     }
 
     store_vector<N_READS>(out, index, out_vec);
@@ -43,7 +43,7 @@ __global__ void copy_v(const In* in, Out* out, IdxT size) {
     AlignedVector<Out, N_READS> out_vec;
 #pragma unroll
     for (int i = 0; i < N_READS; ++i) {
-      out_vec[i] = cast_to<Out>(in_vec[i]);
+      out_vec.val[i] = cast_to<Out>(in_vec.val[i]);
     }
 
     store_vector<N_READS>(out, index, out_vec);
@@ -65,18 +65,23 @@ void copy_contiguous(
         using InType = cuda_type_t<MLX_GET_TYPE(in_type_tag)>;
         using OutType = cuda_type_t<MLX_GET_TYPE(out_type_tag)>;
         using IdxT = std::conditional_t<large(), int64_t, uint32_t>;
-        constexpr int N_READS = 16 / sizeof(InType);
+        // TODO: Choose optimized value based on type size.
+        constexpr int N_READS = 4;
         auto kernel = cu::copy_s<InType, OutType, IdxT, N_READS>;
         if (ctype == CopyType::Vector) {
           kernel = cu::copy_v<InType, OutType, IdxT, N_READS>;
         }
         auto [num_blocks, block_dims] = get_launch_args(
-            out.data_size(), out.shape(), out.strides(), large(), N_READS);
+            kernel,
+            out.data_size(),
+            out.shape(),
+            out.strides(),
+            large(),
+            N_READS);
         encoder.add_kernel_node(
             kernel,
             num_blocks,
             block_dims,
-            0,
             in.data<InType>() + in_offset,
             out.data<OutType>() + out_offset,
             out.data_size());

mlx/backend/cuda/copy/copy_general.cu

@@ -10,80 +10,37 @@ namespace cu {
 
 namespace cg = cooperative_groups;
 
-template <typename In, typename Out, typename IdxT, int NDIM, int N_READS>
+template <typename In, typename Out, typename IdxT, int NDIM>
 __global__ void copy_gg_nd(
     const In* in,
     Out* out,
-    IdxT size_rest,
+    IdxT size,
     const __grid_constant__ cuda::std::array<int32_t, NDIM> shape,
     const __grid_constant__ cuda::std::array<int64_t, NDIM> strides_in,
     const __grid_constant__ cuda::std::array<int64_t, NDIM> strides_out) {
-  auto block = cg::this_thread_block();
-  auto grid = cg::this_grid();
-  IdxT index_rest =
-      grid.block_index().y * block.dim_threads().y + block.thread_index().y;
-  if (index_rest >= size_rest) {
-    return;
-  }
-
-  auto shape_x = shape[NDIM - 1];
-  auto in_stride_x = strides_in[NDIM - 1];
-  auto out_stride_x = strides_out[NDIM - 1];
-  IdxT index_x =
-      grid.block_index().x * block.dim_threads().x + block.thread_index().x;
+  IdxT index = cg::this_grid().thread_rank();
+  if (index < size) {
     auto [idx_in, idx_out] = elem_to_loc_nd<NDIM>(
-      index_rest * shape_x,
-      shape.data(),
-      strides_in.data(),
-      strides_out.data());
-
-  auto in_vec =
-      load_vector<N_READS>(in + idx_in, index_x, shape_x, in_stride_x, In(0));
-  AlignedVector<Out, N_READS> out_vec;
-#pragma unroll
-  for (int i = 0; i < N_READS; ++i) {
-    out_vec[i] = CastOp<In, Out>{}(in_vec[i]);
+        index, shape.data(), strides_in.data(), strides_out.data());
+    out[idx_out] = CastOp<In, Out>{}(in[idx_in]);
   }
-  store_vector(out + idx_out, index_x, out_vec, shape_x, out_stride_x);
 }
 
-template <typename In, typename Out, typename IdxT, int N_READS>
+template <typename In, typename Out, typename IdxT>
 __global__ void copy_gg(
     const In* in,
     Out* out,
-    IdxT size_rest,
+    IdxT size,
     const __grid_constant__ Shape shape,
     const __grid_constant__ Strides strides_in,
     const __grid_constant__ Strides strides_out,
     int ndim) {
-  auto block = cg::this_thread_block();
-  auto grid = cg::this_grid();
-  IdxT index_rest =
-      grid.block_index().y * block.dim_threads().y + block.thread_index().y;
-  if (index_rest >= size_rest) {
-    return;
+  IdxT index = cg::this_grid().thread_rank();
+  if (index < size) {
+    auto [idx_in, idx_out] = elem_to_loc_4d(
+        index, shape.data(), strides_in.data(), strides_out.data(), ndim);
+    out[idx_out] = CastOp<In, Out>{}(in[idx_in]);
   }
-
-  auto shape_x = shape[ndim - 1];
-  auto in_stride_x = strides_in[ndim - 1];
-  auto out_stride_x = strides_out[ndim - 1];
-  IdxT index_x =
-      grid.block_index().x * block.dim_threads().x + block.thread_index().x;
-  auto [idx_in, idx_out] = elem_to_loc(
-      index_rest * shape_x,
-      shape.data(),
-      strides_in.data(),
-      strides_out.data(),
-      ndim);
-
-  auto in_vec =
-      load_vector<N_READS>(in + idx_in, index_x, shape_x, in_stride_x, In(0));
-  AlignedVector<Out, N_READS> out_vec;
-#pragma unroll
-  for (int i = 0; i < N_READS; ++i) {
-    out_vec[i] = CastOp<In, Out>{}(in_vec[i]);
-  }
-  store_vector(out + idx_out, index_x, out_vec, shape_x, out_stride_x);
 }
 
 } // namespace cu
@@ -112,52 +69,34 @@ void copy_general(
             size_t data_size = 1;
             for (auto& s : shape)
               data_size *= s;
-
-            int work_per_thread = 1;
-            auto dim0 = ndim > 0 ? shape.back() : 1;
-            auto rest = data_size / dim0;
-            if (dim0 >= 4) {
-              work_per_thread = 4;
-            }
-
-            dim0 = (dim0 + work_per_thread - 1) / work_per_thread;
-            auto block_dims = get_block_dims(dim0, rest, 1);
-            uint32_t num_blocks_x = cuda::ceil_div(dim0, block_dims.x);
-            uint32_t num_blocks_y = cuda::ceil_div(rest, block_dims.y);
-
             if (ndim <= 3) {
               dispatch_1_2_3(ndim, [&](auto ndim_constant) {
                 auto kernel =
-                    cu::copy_gg_nd<InType, OutType, IdxT, ndim_constant(), 1>;
-                if (work_per_thread == 4) {
-                  kernel =
-                      cu::copy_gg_nd<InType, OutType, IdxT, ndim_constant(), 4>;
-                }
+                    cu::copy_gg_nd<InType, OutType, IdxT, ndim_constant()>;
+                auto [num_blocks, block_dims] = get_launch_args(
+                    kernel, data_size, shape, out.strides(), large());
                 encoder.add_kernel_node(
                     kernel,
-                    {num_blocks_x, num_blocks_y},
+                    num_blocks,
                     block_dims,
-                    0,
                     in_ptr,
                     out_ptr,
-                    rest,
+                    data_size,
                     const_param<ndim_constant()>(shape),
                     const_param<ndim_constant()>(strides_in),
                     const_param<ndim_constant()>(strides_out));
               });
             } else { // ndim >= 4
-              auto kernel = cu::copy_gg<InType, OutType, IdxT, 1>;
-              if (work_per_thread == 4) {
-                kernel = cu::copy_gg<InType, OutType, IdxT, 4>;
-              }
+              auto kernel = cu::copy_gg<InType, OutType, IdxT>;
+              auto [num_blocks, block_dims] = get_launch_args(
+                  kernel, data_size, shape, out.strides(), large());
               encoder.add_kernel_node(
                   kernel,
-                  {num_blocks_x, num_blocks_y},
+                  num_blocks,
                   block_dims,
-                  0,
                   in_ptr,
                   out_ptr,
-                  rest,
+                  data_size,
                   const_param(shape),
                   const_param(strides_in),
                   const_param(strides_out),

mlx/backend/cuda/copy/copy_general_dynamic.cu

@@ -41,7 +41,7 @@ __global__ void copy_gg_dynamic(
     const int64_t* offset_out) {
   IdxT index = cg::this_grid().thread_rank();
   if (index < size) {
-    auto [idx_in, idx_out] = elem_to_loc(
+    auto [idx_in, idx_out] = elem_to_loc_4d(
         index, shape.data(), strides_in.data(), strides_out.data(), ndim);
     out[idx_out + *offset_out] = CastOp<In, Out>{}(in[idx_in + *offset_in]);
   }
@@ -74,16 +74,14 @@ void copy_general_dynamic(
             int ndim = shape.size();
             if (ndim <= 3) {
               dispatch_1_2_3(ndim, [&](auto dims_constant) {
-                auto [num_blocks, block_dims] = get_launch_args(out, large());
+                auto kernel = cu::
+                    copy_gg_dynamic_nd<InType, OutType, IdxT, dims_constant()>;
+                auto [num_blocks, block_dims] =
+                    get_launch_args(kernel, out, large());
                 encoder.add_kernel_node(
-                    cu::copy_gg_dynamic_nd<
-                        InType,
-                        OutType,
-                        IdxT,
-                        dims_constant()>,
+                    kernel,
                     num_blocks,
                     block_dims,
-                    0,
                     in_ptr,
                     out_ptr,
                     out.size(),
@@ -94,12 +92,13 @@ void copy_general_dynamic(
                     dynamic_offset_out.data<int64_t>());
               });
             } else { // ndim >= 4
-              auto [num_blocks, block_dims] = get_launch_args(out, large());
+              auto kernel = cu::copy_gg_dynamic<InType, OutType, IdxT>;
+              auto [num_blocks, block_dims] =
+                  get_launch_args(kernel, out, large());
               encoder.add_kernel_node(
-                  cu::copy_gg_dynamic<InType, OutType, IdxT>,
+                  kernel,
                   num_blocks,
                   block_dims,
-                  0,
                   in_ptr,
                   out_ptr,
                   out.size(),

mlx/backend/cuda/copy/copy_general_input.cu

@@ -10,67 +10,33 @@ namespace cu {
 
 namespace cg = cooperative_groups;
 
-template <typename In, typename Out, typename IdxT, int NDIM, int N_READS>
+template <typename In, typename Out, typename IdxT, int NDIM>
 __global__ void copy_g_nd(
     const In* in,
     Out* out,
-    IdxT size_rest,
+    IdxT size,
     const __grid_constant__ cuda::std::array<int32_t, NDIM> shape,
-    const __grid_constant__ cuda::std::array<int64_t, NDIM> strides) {
-  auto block = cg::this_thread_block();
-  auto grid = cg::this_grid();
-  IdxT index_rest =
-      grid.block_index().y * block.dim_threads().y + block.thread_index().y;
-  if (index_rest >= size_rest) {
-    return;
+    const __grid_constant__ cuda::std::array<int64_t, NDIM> strides_in) {
+  IdxT index = cg::this_grid().thread_rank();
+  if (index < size) {
+    IdxT idx_in = elem_to_loc_nd<NDIM>(index, shape.data(), strides_in.data());
+    out[index] = CastOp<In, Out>{}(in[idx_in]);
   }
-
-  auto shape_x = shape[NDIM - 1];
-  auto stride_x = strides[NDIM - 1];
-  IdxT index_x =
-      grid.block_index().x * block.dim_threads().x + block.thread_index().x;
-  auto idx =
-      elem_to_loc_nd<NDIM>(index_rest * shape_x, shape.data(), strides.data());
-  auto in_vec =
-      load_vector<N_READS>(in + idx, index_x, shape_x, stride_x, In(0));
-  AlignedVector<Out, N_READS> out_vec;
-#pragma unroll
-  for (int i = 0; i < N_READS; ++i) {
-    out_vec[i] = CastOp<In, Out>{}(in_vec[i]);
-  }
-  store_vector(out + shape_x * index_rest, index_x, out_vec, shape_x);
 }
 
-template <typename In, typename Out, typename IdxT, int N_READS>
+template <typename In, typename Out, typename IdxT>
 __global__ void copy_g(
     const In* in,
     Out* out,
-    IdxT size_rest,
+    IdxT size,
     const __grid_constant__ Shape shape,
-    const __grid_constant__ Strides strides,
+    const __grid_constant__ Strides strides_in,
     int ndim) {
-  auto block = cg::this_thread_block();
-  auto grid = cg::this_grid();
-  IdxT index_rest =
-      grid.block_index().y * block.dim_threads().y + block.thread_index().y;
-  if (index_rest >= size_rest) {
-    return;
+  IdxT index = cg::this_grid().thread_rank();
+  if (index < size) {
+    IdxT idx_in = elem_to_loc_4d(index, shape.data(), strides_in.data(), ndim);
+    out[index] = CastOp<In, Out>{}(in[idx_in]);
   }
-
-  auto shape_x = shape[ndim - 1];
-  auto stride_x = strides[ndim - 1];
-  IdxT index_x =
-      grid.block_index().x * block.dim_threads().x + block.thread_index().x;
-  auto idx =
-      elem_to_loc(index_rest * shape_x, shape.data(), strides.data(), ndim);
-  auto in_vec =
-      load_vector<N_READS>(in + idx, index_x, shape_x, stride_x, In(0));
-  AlignedVector<Out, N_READS> out_vec;
-#pragma unroll
-  for (int i = 0; i < N_READS; ++i) {
-    out_vec[i] = CastOp<In, Out>{}(in_vec[i]);
-  }
-  store_vector(out + shape_x * index_rest, index_x, out_vec, shape_x);
 }
 
 } // namespace cu
@@ -95,49 +61,33 @@ void copy_general_input(
             const InType* in_ptr = in.data<InType>() + offset_in;
             OutType* out_ptr = out.data<OutType>() + offset_out;
             int ndim = shape.size();
-            int work_per_thread = 1;
-            auto dim0 = ndim > 0 ? shape.back() : 1;
-            auto rest = out.size() / dim0;
-            if (dim0 >= 4) {
-              work_per_thread = 4;
-            }
-            dim0 = (dim0 + work_per_thread - 1) / work_per_thread;
-            auto block_dims = get_block_dims(dim0, rest, 1);
-            uint32_t num_blocks_x = cuda::ceil_div(dim0, block_dims.x);
-            uint32_t num_blocks_y = cuda::ceil_div(rest, block_dims.y);
-
             if (ndim <= 3) {
               dispatch_1_2_3(ndim, [&](auto dims_constant) {
                 auto kernel =
-                    cu::copy_g_nd<InType, OutType, IdxT, dims_constant(), 1>;
-                if (work_per_thread == 4) {
-                  kernel =
-                      cu::copy_g_nd<InType, OutType, IdxT, dims_constant(), 4>;
-                }
+                    cu::copy_g_nd<InType, OutType, IdxT, dims_constant()>;
+                auto [num_blocks, block_dims] =
+                    get_launch_args(kernel, out, large());
                 encoder.add_kernel_node(
                     kernel,
-                    {num_blocks_x, num_blocks_y},
+                    num_blocks,
                     block_dims,
-                    0,
                     in_ptr,
                     out_ptr,
-                    rest,
+                    out.size(),
                     const_param<dims_constant()>(shape),
                     const_param<dims_constant()>(strides_in));
               });
             } else { // ndim >= 4
-              auto kernel = cu::copy_g<InType, OutType, IdxT, 1>;
-              if (work_per_thread == 4) {
-                kernel = cu::copy_g<InType, OutType, IdxT, 4>;
-              }
+              auto kernel = cu::copy_g<InType, OutType, IdxT>;
+              auto [num_blocks, block_dims] =
+                  get_launch_args(kernel, out, large());
               encoder.add_kernel_node(
                   kernel,
-                  {num_blocks_x, num_blocks_y},
+                  num_blocks,
                   block_dims,
-                  0,
                   in_ptr,
                   out_ptr,
-                  rest,
+                  out.size(),
                   const_param(shape),
                   const_param(strides_in),
                   ndim);

mlx/backend/cuda/cudnn_utils.cpp

@@ -1,272 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include "mlx/backend/cuda/cudnn_utils.h"
-#include "mlx/backend/cuda/device.h"
-
-namespace mlx::core {
-
-namespace {
-
-// Create a cudnn tensor descriptor.
-template <typename Vec>
-inline cudnn_frontend::Tensor build_cudnn_tensor(
-    int64_t id,
-    const array& x,
-    const Vec& shape,
-    const Vec& strides) {
-  return cudnn_frontend::TensorBuilder()
-      .setDim(shape.size(), shape.data())
-      .setStrides(strides.size(), strides.data())
-      .setId(id)
-      .setAlignment(get_alignment(x))
-      .setDataType(dtype_to_cudnn_type(x.dtype()))
-      .build();
-}
-
-// In MLX a singleton dim (shape[dim] == 1) can have any stride, but in cuDNN
-// whether a tensor is contiguous is determined with:
-// shape[dim] == shape[dim + 1] * strides[dim + 1]
-// So a contiguous array with singleton dims in MLX may be mistakenly treated
-// as strided in cuDNN, and we work around it by normalizing the strides.
-Strides normalized_strides(const array& x) {
-  if (!x.flags().row_contiguous || x.ndim() < 2) {
-    return x.strides();
-  }
-  Strides strides = x.strides();
-  for (int i = x.ndim() - 2; i >= 0; --i) {
-    if (x.shape(i) == 1) {
-      strides[i] = x.shape(i + 1) * strides[i + 1];
-    }
-  }
-  return strides;
-}
-
-// Return the shape and strides after transposing from NHWC to NCHW.
-auto nhwc_to_nchw(SmallVector<int64_t> shape, SmallVector<int64_t> strides) {
-  assert(shape.size() >= 3);
-  shape.insert(shape.begin() + 1, shape.back());
-  shape.erase(shape.end() - 1);
-  strides.insert(strides.begin() + 1, strides.back());
-  strides.erase(strides.end() - 1);
-  return std::make_tuple(std::move(shape), std::move(strides));
-}
-
-inline auto nhwc_to_nchw(const array& x) {
-  return nhwc_to_nchw(
-      convert_vector<int64_t>(x.shape()), normalized_strides(x));
-}
-
-// Return available engines for a |op_graph|.
-cudnn_frontend::EngineConfigList get_cudnn_engine_configs(
-    cudnnBackendDescriptorType_t backend_type,
-    Dtype dtype,
-    cudnn_frontend::OperationGraph& op_graph,
-    bool use_fallback = true) {
-  SmallVector<cudnn_frontend::GeneratorSource, 2> sources;
-  sources.push_back([](auto& op_graph) {
-    auto heuristics = cudnn_frontend::EngineHeuristicsBuilder()
-                          .setOperationGraph(op_graph)
-                          .setHeurMode(CUDNN_HEUR_MODE_A)
-                          .build();
-    return heuristics.getEngineConfig(heuristics.getEngineConfigCount());
-  });
-  if (use_fallback) {
-    sources.push_back([&backend_type](auto& op_graph) {
-      auto fallback = cudnn_frontend::EngineFallbackListBuilder()
-                          .setOperationGraph(op_graph)
-                          .setOperation(backend_type)
-                          .build();
-      return fallback.getFallbackList();
-    });
-  }
-
-  auto configs =
-      cudnn_frontend::EngineConfigGenerator(sources.size(), sources.data())
-          .generate_engine_config(op_graph);
-
-  cudnn_frontend::EngineConfigList filtered_configs;
-  cudnn_frontend::filter(configs, filtered_configs, [dtype](auto c) {
-    if (cudnn_frontend::hasNumericalNote<
-            CUDNN_NUMERICAL_NOTE_DOWN_CONVERT_INPUTS>(c)) {
-      return true;
-    }
-    if (cudnn_frontend::hasNumericalNote<CUDNN_NUMERICAL_NOTE_TENSOR_CORE>(c) &&
-        dtype == float32 && !env::enable_tf32()) {
-      return true;
-    }
-    return false;
-  });
-  return filtered_configs;
-}
-
-// Take |engine_configs| and |op_graph| and find a working execution plans
-// from them.
-std::optional<cudnn_frontend::ExecutionPlan>
-find_cudnn_plan_from_engine_configs(
-    cudnnHandle_t handle,
-    const cudnn_frontend::EngineConfigList& engine_configs,
-    const cudnn_frontend::OperationGraph& op_graph) {
-  auto op_graph_tag = op_graph.getTag();
-  for (const auto& config : engine_configs) {
-    try {
-      return cudnn_frontend::ExecutionPlanBuilder()
-          .setHandle(handle)
-          .setEngineConfig(config, op_graph_tag)
-          .build();
-    } catch (cudnn_frontend::cudnnException& error) {
-      if (error.getCudnnStatus() != CUDNN_STATUS_NOT_SUPPORTED) {
-        throw;
-      }
-    }
-  }
-  return std::nullopt;
-}
-
-// Prepare workspace and args to execute plan.
-template <typename F>
-bool prepare_cudnn_plan(
-    cu::CommandEncoder& encoder,
-    cudnn_frontend::ExecutionPlan& plan,
-    int num_args,
-    const int64_t* uids,
-    void** data_ptrs,
-    F&& execute) {
-  int workspace_size = plan.getWorkspaceSize();
-  array workspace(
-      workspace_size > 0 ? allocator::malloc(workspace_size)
-                         : allocator::Buffer(nullptr),
-      {workspace_size},
-      uint8);
-
-  auto args = cudnn_frontend::VariantPackBuilder()
-                  .setWorkspacePointer(workspace.data<void>())
-                  .setDataPointers(num_args, data_ptrs)
-                  .setUids(num_args, uids)
-                  .build();
-
-  auto handle = encoder.device().cudnn_handle();
-  cudnnSetStream(handle, encoder.stream());
-
-  if (!execute(handle, plan.get_raw_desc(), args.get_raw_desc())) {
-    return false;
-  }
-
-  encoder.add_temporary(workspace);
-  return true;
-}
-
-} // namespace
-
-cudnn_frontend::Tensor build_cudnn_tensor(int64_t id, const array& x) {
-  auto shape = convert_vector<int64_t>(x.shape());
-  return build_cudnn_tensor(id, x, shape, normalized_strides(x));
-}
-
-cudnn_frontend::Tensor build_cudnn_tensor_nchw(int64_t id, const array& x) {
-  auto [shape, strides] = nhwc_to_nchw(x);
-  return build_cudnn_tensor(id, x, shape, strides);
-}
-
-cudnn_frontend::Tensor build_cudnn_tensor_4d_nchw(int64_t id, const array& x) {
-  if (x.ndim() == 0) {
-    SmallVector<int64_t, 4> scalar_dims = {1, 1, 1, 1};
-    return build_cudnn_tensor(id, x, scalar_dims, scalar_dims);
-  }
-  if (x.ndim() == 1) {
-    int64_t s = x.shape(0);
-    SmallVector<int64_t, 4> shape = {1, x.shape(0), 1, 1};
-    SmallVector<int64_t, 4> strides = {s, 1, s, s};
-    return build_cudnn_tensor(id, x, shape, strides);
-  }
-  if (x.ndim() == 2) {
-    int64_t s =
-        x.flags().row_contiguous ? x.shape(1) * x.strides(1) : x.strides(0);
-    SmallVector<int64_t, 4> shape = {x.shape(0), x.shape(1), 1, 1};
-    SmallVector<int64_t, 4> strides = {s, x.strides(1), s, s};
-    return build_cudnn_tensor(id, x, shape, strides);
-  }
-  if (x.ndim() == 3 || x.ndim() == 4) {
-    return build_cudnn_tensor_nchw(id, x);
-  }
-  throw std::runtime_error(
-      fmt::format("Unsupported array with {} dims.", x.ndim()));
-}
-
-cudnn_frontend::Tensor build_cudnn_scalar_4d(int64_t id, Dtype dtype) {
-  SmallVector<int64_t, 4> scalar_dims = {1, 1, 1, 1};
-  return cudnn_frontend::TensorBuilder()
-      .setDim(scalar_dims.size(), scalar_dims.data())
-      .setStrides(scalar_dims.size(), scalar_dims.data())
-      .setId(id)
-      .setAlignment(16)
-      .setDataType(dtype_to_cudnn_type(dtype))
-      .setByValue(true)
-      .build();
-}
-
-std::optional<cudnn_frontend::ExecutionPlan> find_cudnn_plan_from_op_graph(
-    cudnnHandle_t handle,
-    cudnnBackendDescriptorType_t backend_type,
-    Dtype dtype,
-    cudnn_frontend::OperationGraph& op_graph) {
-  auto engine_configs = get_cudnn_engine_configs(backend_type, dtype, op_graph);
-  return find_cudnn_plan_from_engine_configs(handle, engine_configs, op_graph);
-}
-
-bool encode_cudnn_plan_with_capturing(
-    cu::CommandEncoder& encoder,
-    cudnn_frontend::ExecutionPlan& plan,
-    int num_args,
-    const int64_t* uids,
-    void** data_ptrs) {
-  return prepare_cudnn_plan(
-      encoder,
-      plan,
-      num_args,
-      uids,
-      data_ptrs,
-      [&](auto handle, auto plan, auto args) {
-        auto capture = encoder.capture_context();
-        if (cudnnBackendExecute(handle, plan, args) != CUDNN_STATUS_SUCCESS) {
-          // Discard the captured graph when failed.
-          capture.discard = true;
-          return false;
-        }
-        return true;
-      });
-}
-
-#if CUDNN_VERSION >= 90500
-bool encode_cudnn_plan_with_graph_api(
-    cu::CommandEncoder& encoder,
-    cudnn_frontend::ExecutionPlan& plan,
-    CudaGraph& graph,
-    int num_args,
-    const int64_t* uids,
-    void** data_ptrs) {
-  return prepare_cudnn_plan(
-      encoder,
-      plan,
-      num_args,
-      uids,
-      data_ptrs,
-      [&](auto handle, auto plan, auto args) {
-        if (!graph) {
-          graph = CudaGraph(encoder.device());
-          if (cudnnBackendPopulateCudaGraph(handle, plan, args, graph) !=
-              CUDNN_STATUS_SUCCESS) {
-            return false;
-          }
-        } else {
-          if (cudnnBackendUpdateCudaGraph(handle, plan, args, graph) !=
-              CUDNN_STATUS_SUCCESS) {
-            return false;
-          }
-        }
-        encoder.add_graph_node(graph);
-        return true;
-      });
-}
-#endif
-
-} // namespace mlx::core

mlx/backend/cuda/cudnn_utils.h

@@ -1,164 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#pragma once
-
-#include "mlx/array.h"
-#include "mlx/backend/cuda/device/config.h"
-#include "mlx/backend/cuda/utils.h"
-#include "mlx/dtype_utils.h"
-
-#include <cudnn_frontend.h>
-#include <cudnn_frontend_find_plan.h>
-#include <fmt/format.h>
-
-#include <algorithm>
-#include <array>
-
-namespace mlx::core {
-
-namespace cu {
-class CommandEncoder;
-}
-
-// Return pointer alignment of |x|'s data.
-inline uint8_t get_alignment(const array& x) {
-  uint8_t alignment = 1;
-  uintptr_t address = reinterpret_cast<uintptr_t>(x.data<void>());
-  for (; alignment < 32; alignment *= 2) {
-    if (address % (alignment * 2)) {
-      return alignment;
-    }
-  }
-  return alignment;
-}
-
-// Convert the type of elements in |vec| to |T|.
-template <typename T, typename Vec>
-inline SmallVector<T> convert_vector(const Vec& vec) {
-  return SmallVector<T>(vec.begin(), vec.end());
-}
-
-// Return an array that can be used as map key for |vec| with size <= MAX_NDIM.
-//
-// There are 2 differences from the const_param util from kernel_utils.cuh:
-// 1. The rest of array is filled with 0.
-// 2. This util can be used in .cpp files.
-template <typename T, template <typename U> class Vec>
-inline std::array<T, MAX_NDIM> vector_key(const Vec<T>& vec) {
-  if (vec.size() > MAX_NDIM) {
-    throw std::runtime_error(
-        fmt::format("ndim can not be larger than {}.", MAX_NDIM));
-  }
-  std::array<T, MAX_NDIM> result = {};
-  std::copy_n(vec.begin(), vec.size(), result.begin());
-  return result;
-}
-
-// Helpers used by get_data_ptrs to get pointers.
-inline void* get_data_ptr(const array& arr) {
-  return const_cast<void*>(arr.data<void>());
-}
-
-template <typename T, typename = std::enable_if_t<std::is_scalar_v<T>>>
-inline void* get_data_ptr(T& scalar) {
-  return &scalar;
-}
-
-// Return an array filled with data pointers of args.
-template <typename... Args>
-inline std::array<void*, sizeof...(Args)> get_data_ptrs(Args&... args) {
-  return {get_data_ptr(args)...};
-}
-
-// Map dtype to cudnn data type.
-inline cudnnDataType_t dtype_to_cudnn_type(Dtype dtype) {
-  switch (dtype) {
-    case int8:
-      return CUDNN_DATA_INT8;
-    case int32:
-      return CUDNN_DATA_INT32;
-    case uint8:
-      return CUDNN_DATA_UINT8;
-    case float16:
-      return CUDNN_DATA_HALF;
-    case bfloat16:
-      return CUDNN_DATA_BFLOAT16;
-    case float32:
-      return CUDNN_DATA_FLOAT;
-    case float64:
-      return CUDNN_DATA_DOUBLE;
-    default:
-      throw std::runtime_error(fmt::format(
-          "Unsupported dtype in Convolution: {}.", dtype_to_string(dtype)));
-  }
-}
-
-// Create a tensor descriptor from |x|.
-cudnn_frontend::Tensor build_cudnn_tensor(int64_t id, const array& x);
-
-// Create a tensor descriptor from |x|, and transpose from NHWC to NCHW.
-cudnn_frontend::Tensor build_cudnn_tensor_nchw(int64_t id, const array& x);
-
-// Create a tensor descriptor from |x|, make sure it is 4D, and transpose it
-// from NHWC to NCHW.
-cudnn_frontend::Tensor build_cudnn_tensor_4d_nchw(int64_t id, const array& x);
-
-// Create a 4D scalar tensor descriptor, which is passed by value.
-cudnn_frontend::Tensor build_cudnn_scalar_4d(int64_t id, Dtype dtype);
-
-// Find a working plan for |op_graph|.
-std::optional<cudnn_frontend::ExecutionPlan> find_cudnn_plan_from_op_graph(
-    cudnnHandle_t handle,
-    cudnnBackendDescriptorType_t backend_type,
-    Dtype dtype,
-    cudnn_frontend::OperationGraph& op_graph);
-
-// Encode the plan to command buffer by capturing.
-bool encode_cudnn_plan_with_capturing(
-    cu::CommandEncoder& encoder,
-    cudnn_frontend::ExecutionPlan& plan,
-    int num_args,
-    const int64_t* uids,
-    void** data_ptrs);
-
-#if CUDNN_VERSION >= 90500
-// Encode the plan to command buffer by using native graph api of cudnn. If the
-// |graph| is empty it will be populated, otherwise it will be updated.
-bool encode_cudnn_plan_with_graph_api(
-    cu::CommandEncoder& encoder,
-    cudnn_frontend::ExecutionPlan& plan,
-    CudaGraph& graph,
-    int num_args,
-    const int64_t* uids,
-    void** data_ptrs);
-#endif
-
-// Helpers to make calls like encode_cudnn_plan(..., {'x', 'y', 'z'}, x, y, z).
-template <typename... Args>
-bool encode_cudnn_plan(
-    cu::CommandEncoder& encoder,
-    cudnn_frontend::ExecutionPlan& plan,
-    std::initializer_list<int64_t> uids,
-    Args&... args) {
-  assert(uids.size() == sizeof...(args));
-  auto data_ptrs = get_data_ptrs(args...);
-  return encode_cudnn_plan_with_capturing(
-      encoder, plan, uids.size(), uids.begin(), data_ptrs.data());
-}
-
-#if CUDNN_VERSION >= 90500
-template <typename... Args>
-bool encode_cudnn_plan(
-    cu::CommandEncoder& encoder,
-    cudnn_frontend::ExecutionPlan& plan,
-    CudaGraph& graph,
-    std::initializer_list<int64_t> uids,
-    Args&... args) {
-  assert(uids.size() == sizeof...(args));
-  auto data_ptrs = get_data_ptrs(args...);
-  return encode_cudnn_plan_with_graph_api(
-      encoder, plan, graph, uids.size(), uids.begin(), data_ptrs.data());
-}
-#endif
-
-} // namespace mlx::core

mlx/backend/cuda/custom_kernel.cpp

@@ -1,379 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#include <iostream>
-
-#include "mlx/backend/common/compiled.h"
-#include "mlx/backend/cuda/jit_module.h"
-#include "mlx/backend/cuda/utils.h"
-#include "mlx/backend/gpu/copy.h"
-#include "mlx/fast.h"
-#include "mlx/fast_primitives.h"
-
-#include <fmt/format.h>
-#include <nvtx3/nvtx3.hpp>
-
-namespace mlx::core::fast {
-
-namespace {
-
-constexpr const char* default_header = R"(
-#include "mlx/backend/cuda/device/utils.cuh"
-
-#include <cooperative_groups.h>
-
-#define inf cuda::std::numeric_limits<float>::infinity()
-
-)";
-
-std::string template_arguments_hash(
-    const std::vector<std::pair<std::string, TemplateArg>>& template_args) {
-  if (template_args.empty()) {
-    return "";
-  }
-
-  std::string hash;
-  hash.reserve(512);
-
-  for (const auto& [name, arg] : template_args) {
-    if (std::holds_alternative<int>(arg)) {
-      hash += fmt::format("_{}", std::get<int>(arg));
-    } else if (std::holds_alternative<bool>(arg)) {
-      hash += (std::get<bool>(arg)) ? "_t" : "_f";
-    } else if (std::holds_alternative<Dtype>(arg)) {
-      hash += "_";
-      hash += get_type_string(std::get<Dtype>(arg));
-    }
-  }
-
-  return hash;
-}
-
-std::string build_kernel(
-    const std::string& func_name,
-    const std::string& header,
-    const std::string& source,
-    const std::vector<std::string>& input_names,
-    const std::vector<array>& inputs,
-    const std::vector<std::string>& output_names,
-    const std::vector<Dtype>& output_dtypes,
-    const std::vector<std::pair<std::string, TemplateArg>>& template_args,
-    const std::vector<CustomKernelShapeInfo>& shape_infos) {
-  std::string kernel_source;
-  kernel_source.reserve(header.size() + source.size() + 8192);
-  kernel_source += default_header;
-  kernel_source += header;
-  kernel_source +=
-      "namespace mlx::core::cu {\n\n"
-      "namespace cg = cooperative_groups;\n\n";
-
-  kernel_source += "__global__ void ";
-  kernel_source += func_name;
-  kernel_source += "(\n";
-
-  // Add inputs
-  for (int i = 0; i < inputs.size(); ++i) {
-    const auto& name = input_names[i];
-    const auto& arr = inputs[i];
-    kernel_source += "    const ";
-    kernel_source += dtype_to_cuda_type(arr.dtype());
-    kernel_source += "* ";
-    kernel_source += name;
-    kernel_source += ",\n";
-    // Add input shape, strides and ndim if present in the source
-    if (arr.ndim() > 0) {
-      if (shape_infos[i].shape) {
-        kernel_source += "    const __grid_constant__ Shape ";
-        kernel_source += name;
-        kernel_source += "_shape,\n";
-      }
-      if (shape_infos[i].strides) {
-        kernel_source += "    const __grid_constant__ Strides ";
-        kernel_source += name;
-        kernel_source += "_strides,\n";
-      }
-      if (shape_infos[i].ndim) {
-        kernel_source += "    const __grid_constant__ int ";
-        kernel_source += name;
-        kernel_source += "_ndim,\n";
-      }
-    }
-  }
-
-  // Add outputs
-  for (int i = 0; i < output_names.size(); ++i) {
-    const auto& name = output_names[i];
-    const auto& dtype = output_dtypes[i];
-    kernel_source += "    ";
-    kernel_source += dtype_to_cuda_type(dtype);
-    kernel_source += "* ";
-    kernel_source += name;
-    if (i < output_names.size() - 1) {
-      kernel_source += ",\n";
-    } else {
-      kernel_source += ") {\n";
-    }
-  }
-
-  // Set compile time constants
-  if (!template_args.empty()) {
-    for (const auto& [name, arg] : template_args) {
-      if (std::holds_alternative<int>(arg)) {
-        kernel_source +=
-            fmt::format("  constexpr int {} = {};\n", name, std::get<int>(arg));
-      } else if (std::holds_alternative<bool>(arg)) {
-        kernel_source += fmt::format(
-            "  constexpr bool {} = {};\n", name, std::get<bool>(arg));
-      } else {
-        kernel_source += fmt::format(
-            "  using {} = {};\n",
-            name,
-            dtype_to_cuda_type(std::get<Dtype>(arg)));
-      }
-    }
-    kernel_source += "\n";
-  }
-
-  kernel_source += source;
-  kernel_source += "\n}\n\n} // namespace mlx::core::cu\n";
-
-  return kernel_source;
-}
-
-} // namespace
-
-CustomKernelFunction cuda_kernel(
-    const std::string& name,
-    const std::vector<std::string>& input_names,
-    const std::vector<std::string>& output_names,
-    const std::string& source,
-    const std::string& header,
-    bool ensure_row_contiguous,
-    int shared_memory) {
-  if (output_names.empty()) {
-    throw std::invalid_argument(
-        "[custom_kernel] Must specify at least one output.");
-  }
-
-  std::vector<CustomKernelShapeInfo> shape_infos;
-  for (auto& n : input_names) {
-    CustomKernelShapeInfo shape_info;
-    shape_info.shape = source.find(n + "_shape") != std::string::npos;
-    shape_info.strides = source.find(n + "_strides") != std::string::npos;
-    shape_info.ndim = source.find(n + "_ndim") != std::string::npos;
-    shape_infos.push_back(shape_info);
-  }
-
-  return [=, shape_infos = std::move(shape_infos)](
-             const std::vector<array>& inputs,
-             const std::vector<Shape>& output_shapes,
-             const std::vector<Dtype>& output_dtypes,
-             std::tuple<int, int, int> grid,
-             std::tuple<int, int, int> threadgroup,
-             const std::vector<std::pair<std::string, TemplateArg>>&
-                 template_args = {},
-             std::optional<float> init_value = std::nullopt,
-             bool verbose = false,
-             StreamOrDevice s_ = {}) {
-    if (inputs.size() != input_names.size()) {
-      std::ostringstream msg;
-      msg << "[custom_kernel] Expected `inputs` to have size "
-          << input_names.size() << " but got size " << inputs.size() << "."
-          << std::endl;
-      throw std::invalid_argument(msg.str());
-    }
-    if (output_shapes.size() != output_names.size()) {
-      std::ostringstream msg;
-      msg << "[custom_kernel] Expected `output_shapes` to have size "
-          << output_names.size() << " but got size " << output_shapes.size()
-          << "." << std::endl;
-      throw std::invalid_argument(msg.str());
-    }
-    if (output_dtypes.size() != output_names.size()) {
-      std::ostringstream msg;
-      msg << "[custom_kernel] Expected `output_dtypes` to have size "
-          << output_names.size() << " but got size " << output_dtypes.size()
-          << "." << std::endl;
-      throw std::invalid_argument(msg.str());
-    }
-
-    auto s = to_stream(s_);
-    if (s.device != Device::gpu) {
-      throw std::invalid_argument("[custom_kernel] Only supports the GPU.");
-    }
-
-    std::string kernel_name =
-        "custom_kernel_" + name + template_arguments_hash(template_args);
-    std::string kernel_source = build_kernel(
-        kernel_name,
-        header,
-        source,
-        input_names,
-        inputs,
-        output_names,
-        output_dtypes,
-        template_args,
-        shape_infos);
-
-    if (verbose) {
-      std::cout << "Generated source code for `" << kernel_name
-                << "`:" << std::endl
-                << "```" << std::endl
-                << kernel_source << std::endl
-                << "```" << std::endl;
-    }
-
-    return array::make_arrays(
-        std::move(output_shapes),
-        std::move(output_dtypes),
-        std::make_shared<CustomKernel>(
-            s,
-            std::move(kernel_name),
-            std::move(kernel_source),
-            grid,
-            threadgroup,
-            shape_infos,
-            ensure_row_contiguous,
-            init_value,
-            std::vector<ScalarArg>{},
-            false,
-            shared_memory),
-        std::move(inputs));
-  };
-}
-
-std::vector<array> precompiled_cuda_kernel(
-    const std::string& name,
-    const std::string& compiled_source,
-    const std::vector<array>& inputs,
-    const std::vector<Shape>& output_shapes,
-    const std::vector<Dtype>& output_dtypes,
-    const std::vector<ScalarArg>& scalars,
-    std::tuple<int, int, int> grid,
-    std::tuple<int, int, int> threadgroup,
-    int shared_memory,
-    std::optional<float> init_value,
-    bool ensure_row_contiguous,
-    StreamOrDevice s) {
-  std::vector<CustomKernelShapeInfo> shape_infos(
-      inputs.size(), CustomKernelShapeInfo{false, false, false});
-  return array::make_arrays(
-      output_shapes,
-      output_dtypes,
-      std::make_shared<CustomKernel>(
-          to_stream(s),
-          name,
-          compiled_source,
-          grid,
-          threadgroup,
-          shape_infos,
-          ensure_row_contiguous,
-          init_value,
-          scalars,
-          true,
-          shared_memory),
-      inputs);
-}
-
-void CustomKernel::eval_gpu(
-    const std::vector<array>& inputs,
-    std::vector<array>& outputs) {
-  nvtx3::scoped_range r("CustomKernel::eval_gpu");
-  auto& s = stream();
-
-  std::vector<array> copies;
-
-  // Allocate and initialize the output arrays
-  for (auto& out : outputs) {
-    if (init_value_) {
-      copies.emplace_back(init_value_.value(), out.dtype());
-      fill_gpu(copies.back(), out, s);
-    } else {
-      out.set_data(allocator::malloc(out.nbytes()));
-    }
-  }
-
-  // Create the input arrays and copy if needed
-  auto check_input = [&copies, &s, this](const array& x) -> const array {
-    bool no_copy = x.flags().row_contiguous;
-    if (!ensure_row_contiguous_ || no_copy) {
-      return x;
-    } else {
-      copies.push_back(array(x.shape(), x.dtype(), nullptr, {}));
-      copy_gpu(x, copies.back(), CopyType::General, s);
-      return copies.back();
-    }
-  };
-  std::vector<array> checked_inputs;
-  for (const array& in : inputs) {
-    checked_inputs.push_back(check_input(in));
-  }
-
-  // Compile the custom kernel
-  std::string kernel_name =
-      (is_precompiled_) ? name_ : "mlx::core::cu::" + name_;
-  cu::JitModule& mod = cu::get_jit_module(
-      s.device,
-      name_,
-      [&]() {
-        return std::make_tuple(
-            is_precompiled_, source_, std::vector{kernel_name});
-      },
-      false);
-
-  // Make the arguments
-  cu::KernelArgs args;
-  for (int i = 0; i < checked_inputs.size(); i++) {
-    const array& in = checked_inputs[i];
-    auto& shape_info = shape_infos_[i];
-    args.append(in);
-    if (shape_info.shape) {
-      args.append_ndim(in.shape());
-    }
-    if (shape_info.strides) {
-      args.append_ndim(in.strides());
-    }
-    if (shape_info.ndim) {
-      args.append<int32_t>(in.ndim());
-    }
-  }
-  for (auto& out : outputs) {
-    args.append(out);
-  }
-  for (auto& s : scalar_arguments_) {
-    if (std::holds_alternative<bool>(s)) {
-      args.append(std::get<bool>(s));
-    } else if (std::holds_alternative<int>(s)) {
-      args.append(std::get<int>(s));
-    } else if (std::holds_alternative<float>(s)) {
-      args.append(std::get<float>(s));
-    }
-  }
-
-  // Make the grid
-  const auto [tx, ty, tz] = threadgroup_;
-  const auto [gx, gy, gz] = grid_;
-  dim3 block(std::min(tx, gx), std::min(ty, gy), std::min(tz, gz));
-  dim3 grid((gx + tx - 1) / tx, (gy + ty - 1) / ty, (gz + tz - 1) / tz);
-
-  // Call the kernel
-  auto& encoder = cu::get_command_encoder(s);
-  for (const auto& in : checked_inputs) {
-    encoder.set_input_array(in);
-  }
-  for (const auto& out : outputs) {
-    encoder.set_output_array(out);
-  }
-  for (const auto& t : copies) {
-    encoder.add_temporary(t);
-  }
-  auto kernel =
-      mod.get_kernel(kernel_name, [smem = shared_memory_](CUfunction kernel) {
-        if (smem > 0 && smem > 48000) {
-          cuFuncSetAttribute(
-              kernel, CU_FUNC_ATTRIBUTE_MAX_DYNAMIC_SHARED_SIZE_BYTES, smem);
-        }
-      });
-  encoder.add_kernel_node(kernel, grid, block, shared_memory_, args.args());
-}
-
-} // namespace mlx::core::fast

mlx/backend/cuda/device.cpp

@@ -1,7 +1,6 @@
 // Copyright © 2025 Apple Inc.
 
 #include "mlx/backend/cuda/device.h"
-#include "mlx/backend/cuda/jit_module.h"
 #include "mlx/backend/cuda/worker.h"
 #include "mlx/utils.h"
 
@@ -10,23 +9,12 @@
 #include <future>
 #include <unordered_set>
 
-namespace mlx::core::cu {
-
-namespace {
+namespace mlx::core {
 
 // Can be tuned with MLX_MAX_OPS_PER_BUFFER
 // This should be less than 255
 constexpr int default_max_nodes_per_graph = 20;
 
-#define CHECK_CUDNN_ERROR(cmd) check_cudnn_error(#cmd, (cmd))
-
-void check_cudnn_error(const char* name, cudnnStatus_t err) {
-  if (err != CUDNN_STATUS_SUCCESS) {
-    throw std::runtime_error(
-        fmt::format("{} failed: {}.", name, cudnnGetErrorString(err)));
-  }
-}
-
 int cuda_graph_cache_size() {
   static int cache_size = []() {
     return env::get_var("MLX_CUDA_GRAPH_CACHE_SIZE", 100);
@@ -34,7 +22,7 @@ int cuda_graph_cache_size() {
   return cache_size;
 }
 
-} // namespace
+namespace cu {
 
 Device::Device(int device) : device_(device) {
   CHECK_CUDA_ERROR(cudaDeviceGetAttribute(
@@ -52,18 +40,11 @@ Device::Device(int device) : device_(device) {
   }
   // The cublasLt handle is used by matmul.
   make_current();
-  CHECK_CUBLAS_ERROR(cublasLtCreate(&lt_));
-  // The cudnn handle is used by Convolution.
-  CHECK_CUDNN_ERROR(cudnnCreate(&cudnn_));
-
-  // Initialize the jit module cache here ensures it is not
-  // unloaded before any evaluation is done
-  get_jit_module_cache();
+  cublasLtCreate(&lt_);
 }
 
 Device::~Device() {
-  CHECK_CUDNN_ERROR(cudnnDestroy(cudnn_));
-  CHECK_CUBLAS_ERROR(cublasLtDestroy(lt_));
+  cublasLtDestroy(lt_);
 }
 
 void Device::make_current() {
@@ -85,17 +66,30 @@ CommandEncoder& Device::get_command_encoder(Stream s) {
 }
 
 CommandEncoder::CaptureContext::CaptureContext(CommandEncoder& enc) : enc(enc) {
-  enc.device().make_current();
+  CHECK_CUDA_ERROR(cudaGraphCreate(&graph, 0));
   CHECK_CUDA_ERROR(
       cudaStreamBeginCapture(enc.stream(), cudaStreamCaptureModeGlobal));
 }
 
 CommandEncoder::CaptureContext::~CaptureContext() {
-  graph.end_capture(enc.stream());
-  if (discard) {
-    return;
+  CHECK_CUDA_ERROR(cudaStreamEndCapture(enc.stream(), &graph));
+  size_t num_nodes;
+  CHECK_CUDA_ERROR(cudaGraphGetNodes(graph, NULL, &num_nodes));
+  if (num_nodes == 1) {
+    cudaGraphNode_t captured_node;
+    CHECK_CUDA_ERROR(cudaGraphGetNodes(graph, &captured_node, &num_nodes));
+    CUDA_KERNEL_NODE_PARAMS params;
+    CHECK_CUDA_ERROR(cuGraphKernelNodeGetParams(captured_node, &params));
+    cudaGraphNode_t node;
+    CHECK_CUDA_ERROR(cuGraphAddKernelNode(&node, enc.graph_, NULL, 0, &params));
+    enc.insert_graph_dependencies(GraphNode{node, 'K'});
+  } else {
+    cudaGraphNode_t node;
+    CHECK_CUDA_ERROR(
+        cudaGraphAddChildGraphNode(&node, enc.graph_, NULL, 0, graph));
+    enc.insert_graph_dependencies(GraphNode{node, 'G'});
   }
-  enc.add_graph_node(graph);
+  CHECK_CUDA_ERROR(cudaGraphDestroy(graph));
 }
 
 CommandEncoder::ConcurrentContext::ConcurrentContext(CommandEncoder& enc)
@@ -182,11 +176,20 @@ void CommandEncoder::insert_graph_dependencies(std::vector<GraphNode> nodes) {
   }
 }
 
-CommandEncoder::CommandEncoder(Device& d)
-    : device_(d),
-      stream_(d),
-      graph_(d),
-      graph_cache_(cuda_graph_cache_size()) {}
+CommandEncoder::CommandEncoder(Device& d) : device_(d), stream_(d) {
+  CHECK_CUDA_ERROR(cudaGraphCreate(&graph_, 0));
+}
+
+void clear_graphs(std::unordered_map<std::string, cudaGraphExec_t>& graphs) {
+  for (auto& [_, graph_exec] : graphs) {
+    CHECK_CUDA_ERROR(cudaGraphExecDestroy(graph_exec));
+  }
+  graphs.clear();
+}
+
+CommandEncoder::~CommandEncoder() {
+  clear_graphs(graph_cache_);
+}
 
 void CommandEncoder::add_completed_handler(std::function<void()> task) {
   worker_.add_task(std::move(task));
@@ -213,22 +216,22 @@ void CommandEncoder::add_kernel_node(
     void* func,
     dim3 grid_dim,
     dim3 block_dim,
-    uint32_t smem_bytes,
     void** params) {
   cudaKernelNodeParams kernel_params = {0};
   kernel_params.func = func;
   kernel_params.gridDim = grid_dim;
   kernel_params.blockDim = block_dim;
   kernel_params.kernelParams = params;
-  kernel_params.sharedMemBytes = smem_bytes;
-  add_kernel_node(kernel_params);
+  cudaGraphNode_t node;
+  CHECK_CUDA_ERROR(
+      cudaGraphAddKernelNode(&node, graph_, NULL, 0, &kernel_params));
+  insert_graph_dependencies(GraphNode{node, 'K'});
 }
 
 void CommandEncoder::add_kernel_node(
     CUfunction func,
     dim3 grid_dim,
     dim3 block_dim,
-    uint32_t smem_bytes,
     void** params) {
   CUDA_KERNEL_NODE_PARAMS kernel_params = {0};
   kernel_params.func = func;
@@ -239,43 +242,20 @@ void CommandEncoder::add_kernel_node(
   kernel_params.blockDimY = block_dim.y;
   kernel_params.blockDimZ = block_dim.z;
   kernel_params.kernelParams = params;
-  kernel_params.sharedMemBytes = smem_bytes;
-  add_kernel_node(kernel_params);
-}
-
-void CommandEncoder::add_kernel_node(const cudaKernelNodeParams& params) {
-  cudaGraphNode_t node;
-  CHECK_CUDA_ERROR(cudaGraphAddKernelNode(&node, graph_, NULL, 0, &params));
-  insert_graph_dependencies(GraphNode{node, 'K'});
-}
-
-void CommandEncoder::add_kernel_node(const CUDA_KERNEL_NODE_PARAMS& params) {
   CUgraphNode node;
-  CHECK_CUDA_ERROR(cuGraphAddKernelNode(&node, graph_, NULL, 0, &params));
+  CHECK_CUDA_ERROR(
+      cuGraphAddKernelNode(&node, graph_, NULL, 0, &kernel_params));
   insert_graph_dependencies(GraphNode{node, 'K'});
 }
 
-void CommandEncoder::add_graph_node(cudaGraph_t child) {
-  cudaGraphNode_t node;
-  CHECK_CUDA_ERROR(cudaGraphAddChildGraphNode(&node, graph_, NULL, 0, child));
-  insert_graph_dependencies(GraphNode{node, 'G'});
-}
-
 void CommandEncoder::commit() {
-  nvtx3::scoped_range r("CommandEncoder::commit");
   if (!temporaries_.empty()) {
     add_completed_handler([temporaries = std::move(temporaries_)]() {});
   }
   if (node_count_ > 0) {
     if (!from_nodes_.empty()) {
       CHECK_CUDA_ERROR(cudaGraphAddDependencies(
-          graph_,
-          from_nodes_.data(),
-          to_nodes_.data(),
-#if CUDART_VERSION >= 13000
-          nullptr, // edgeData
-#endif // CUDART_VERSION >= 13000
-          from_nodes_.size()));
+          graph_, from_nodes_.data(), to_nodes_.data(), from_nodes_.size()));
     }
 
     graph_key_ += ".";
@@ -285,7 +265,7 @@ void CommandEncoder::commit() {
     graph_key_ += ".";
     graph_key_ += std::to_string(empty_node_count_);
 
-    CudaGraphExec& graph_exec = graph_cache_[graph_key_];
+    cudaGraphExec_t& graph_exec = graph_cache_[graph_key_];
 
     if (graph_exec != nullptr) {
       cudaGraphExecUpdateResult update_result;
@@ -299,27 +279,35 @@ void CommandEncoder::commit() {
 #endif // CUDART_VERSION >= 12000
       if (update_result != cudaGraphExecUpdateSuccess) {
         cudaGetLastError(); // reset error
-        graph_exec.reset();
+        CHECK_CUDA_ERROR(cudaGraphExecDestroy(graph_exec));
+        graph_exec = nullptr;
       }
     }
     if (graph_exec == nullptr) {
-      graph_exec.instantiate(graph_);
+      CHECK_CUDA_ERROR(
+          cudaGraphInstantiate(&graph_exec, graph_, NULL, NULL, 0));
     }
     device_.make_current();
     CHECK_CUDA_ERROR(cudaGraphLaunch(graph_exec, stream_));
 
+    // TODO smarter cache policy
+    if (graph_cache_.size() > cuda_graph_cache_size()) {
+      clear_graphs(graph_cache_);
+    }
+
     // Reset state
     node_count_ = 0;
     graph_node_count_ = 0;
-    empty_node_count_ = 0;
     from_nodes_.clear();
     to_nodes_.clear();
     graph_key_.clear();
     node_map_.clear();
-    graph_ = CudaGraph(device_);
+    CHECK_CUDA_ERROR(cudaGraphDestroy(graph_));
+    CHECK_CUDA_ERROR(cudaGraphCreate(&graph_, 0));
   }
 
   // Put completion handlers in a batch.
+  worker_.end_batch();
   worker_.commit(stream_);
 }
 
@@ -328,6 +316,7 @@ void CommandEncoder::synchronize() {
   auto p = std::make_shared<std::promise<void>>();
   std::future<void> f = p->get_future();
   add_completed_handler([p = std::move(p)]() { p->set_value(); });
+  worker_.end_batch();
   commit();
   f.wait();
 }
@@ -345,4 +334,6 @@ CommandEncoder& get_command_encoder(Stream s) {
   return device(s.device).get_command_encoder(s);
 }
 
-} // namespace mlx::core::cu
+} // namespace cu
+
+} // namespace mlx::core

mlx/backend/cuda/device.h

@@ -3,13 +3,11 @@
 #pragma once
 
 #include "mlx/array.h"
-#include "mlx/backend/cuda/lru_cache.h"
 #include "mlx/backend/cuda/worker.h"
 #include "mlx/stream.h"
 
 #include <cublasLt.h>
 #include <cuda.h>
-#include <cudnn.h>
 #include <thrust/execution_policy.h>
 
 #include <unordered_map>
@@ -21,9 +19,8 @@ class CommandEncoder {
   struct CaptureContext {
     CaptureContext(CommandEncoder& enc);
     ~CaptureContext();
-    CudaGraph graph;
+    cudaGraph_t graph;
     CommandEncoder& enc;
-    bool discard{false};
   };
   struct ConcurrentContext {
     ConcurrentContext(CommandEncoder& enc);
@@ -32,6 +29,7 @@ class CommandEncoder {
   };
 
   explicit CommandEncoder(Device& d);
+  ~CommandEncoder();
 
   CommandEncoder(const CommandEncoder&) = delete;
   CommandEncoder& operator=(const CommandEncoder&) = delete;
@@ -47,39 +45,25 @@ class CommandEncoder {
   void set_output_array(const array& arr);
 
   template <typename F, typename... Params>
-  void add_kernel_node(
-      F* func,
-      dim3 grid_dim,
-      dim3 block_dim,
-      uint32_t smem_bytes,
-      Params&&... params) {
+  void
+  add_kernel_node(F* func, dim3 grid_dim, dim3 block_dim, Params&&... params) {
     constexpr size_t num = sizeof...(Params);
     void* ptrs[num];
     size_t i = 0;
     ([&](auto&& p) { ptrs[i++] = static_cast<void*>(&p); }(
          std::forward<Params>(params)),
      ...);
-    add_kernel_node((void*)func, grid_dim, block_dim, smem_bytes, ptrs);
+    add_kernel_node((void*)func, grid_dim, block_dim, ptrs);
   }
 
   void add_kernel_node(
       CUfunction func,
       dim3 grid_dim,
       dim3 block_dim,
-      uint32_t smem_bytes,
       void** params);
 
-  void add_kernel_node(
-      void* func,
-      dim3 grid_dim,
-      dim3 block_dim,
-      uint32_t smem_bytes,
-      void** params);
-
-  // Low-level graph helpers.
-  void add_kernel_node(const cudaKernelNodeParams& params);
-  void add_kernel_node(const CUDA_KERNEL_NODE_PARAMS& params);
-  void add_graph_node(cudaGraph_t child);
+  void
+  add_kernel_node(void* func, dim3 grid_dim, dim3 block_dim, void** params);
 
   void add_temporary(const array& arr) {
     temporaries_.push_back(arr.data_shared_ptr());
@@ -89,10 +73,6 @@ class CommandEncoder {
   void maybe_commit();
   void commit();
 
-  Device& device() {
-    return device_;
-  }
-
   CudaStream& stream() {
     return stream_;
   }
@@ -115,7 +95,7 @@ class CommandEncoder {
 
   Device& device_;
   CudaStream stream_;
-  CudaGraph graph_;
+  cudaGraph_t graph_;
   Worker worker_;
   char node_count_{0};
   char graph_node_count_{0};
@@ -126,7 +106,7 @@ class CommandEncoder {
   std::string graph_key_;
   std::vector<GraphNode> concurrent_nodes_;
   std::vector<std::shared_ptr<array::Data>> temporaries_;
-  LRUCache<std::string, CudaGraphExec> graph_cache_;
+  std::unordered_map<std::string, cudaGraphExec_t> graph_cache_;
   std::vector<std::uintptr_t> active_deps_;
   std::vector<std::uintptr_t> active_outputs_;
   std::unordered_map<std::uintptr_t, GraphNode> node_map_;
@@ -157,16 +137,12 @@ class Device {
   cublasLtHandle_t lt_handle() const {
     return lt_;
   }
-  cudnnHandle_t cudnn_handle() const {
-    return cudnn_;
-  }
 
  private:
   int device_;
   int compute_capability_major_;
   int compute_capability_minor_;
   cublasLtHandle_t lt_;
-  cudnnHandle_t cudnn_;
   std::unordered_map<int, CommandEncoder> encoders_;
 };
 

mlx/backend/cuda/device/arange.cuh

@@ -0,0 +1,15 @@
+// Copyright © 2025 Apple Inc.
+
+namespace mlx::core::cu {
+
+template <typename T>
+struct Arange {
+  const T start;
+  const T step;
+
+  __device__ T operator()(uint32_t i) const {
+    return start + i * step;
+  }
+};
+
+} // namespace mlx::core::cu

mlx/backend/cuda/device/atomic_ops.cuh

@@ -2,7 +2,7 @@
 
 #pragma once
 
-#include "mlx/backend/cuda/device/complex.cuh"
+#include "mlx/backend/cuda/device/cucomplex_math.cuh"
 #include "mlx/backend/cuda/device/fp16_math.cuh"
 
 #include <cuda/atomic>
@@ -48,16 +48,25 @@ inline __device__ void atomic_add(__half* out, __half val) {
   atomicAdd(out, val);
 }
 
-inline __device__ void atomic_add(complex64_t* out, complex64_t val) {
-  atomic_add_general(out, val);
-}
-
-inline __device__ void atomic_add(__nv_bfloat16* out, __nv_bfloat16 val) {
-#if __CUDA_ARCH__ < 800
+inline __device__ void atomic_add(cuComplex* out, cuComplex val) {
+#if __CUDA_ARCH__ < 900
   atomic_add_general(out, val);
 #else
   atomicAdd(out, val);
 #endif
 }
 
+inline __device__ void atomic_add(__nv_bfloat16* out, __nv_bfloat16 val) {
+#if __CUDA_ARCH__ < 800
+#if CCCL_VERSION >= 2008000
+  atomic_add_general(out, val);
+#else
+  bool cccl_version_too_old_for_bfloat16_atomic_add = false;
+  assert(cccl_version_too_old_for_bfloat16_atomic_add);
+#endif
+#else
+  atomicAdd(out, val);
+#endif
+}
+
 } // namespace mlx::core::cu

mlx/backend/cuda/device/binary_ops.cuh

@@ -44,7 +44,7 @@ struct Remainder {
       } else {
         return x % y;
       }
-    } else if constexpr (is_complex_v<T>) {
+    } else if constexpr (cuda::std::is_same_v<T, cuComplex>) {
       return x % y;
     } else {
       T r = fmod(x, y);
@@ -66,12 +66,14 @@ struct Equal {
 struct NaNEqual {
   template <typename T>
   __device__ bool operator()(T x, T y) {
-    if constexpr (is_complex_v<T>) {
+    if constexpr (std::is_same_v<T, cuComplex>) {
       return x == y ||
-          (isnan(x.real()) && isnan(y.real()) && isnan(x.imag()) &&
-           isnan(y.imag())) ||
-          (x.real() == y.real() && isnan(x.imag()) && isnan(y.imag())) ||
-          (isnan(x.real()) && isnan(y.real()) && x.imag() == y.imag());
+          (isnan(cuCrealf(x)) && isnan(cuCrealf(y)) && isnan(cuCimagf(x)) &&
+           isnan(cuCimagf(y))) ||
+          (cuCrealf(x) == cuCrealf(y) && isnan(cuCimagf(x)) &&
+           isnan(cuCimagf(y))) ||
+          (isnan(cuCrealf(x)) && isnan(cuCrealf(y)) &&
+           cuCimagf(x) == cuCimagf(y));
     } else {
       return x == y || (isnan(x) && isnan(y));
     }
@@ -109,17 +111,17 @@ struct LessEqual {
 struct LogAddExp {
   template <typename T>
   __device__ T operator()(T x, T y) {
-    if constexpr (is_complex_v<T>) {
-      if (isnan(x.real()) || isnan(x.imag()) || isnan(y.real()) ||
-          isnan(y.imag())) {
+    if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      if (isnan(cuCrealf(x)) || isnan(cuCimagf(x)) || isnan(cuCrealf(y)) ||
+          isnan(cuCimagf(y))) {
         return {
             cuda::std::numeric_limits<float>::quiet_NaN(),
             cuda::std::numeric_limits<float>::quiet_NaN()};
       }
-      auto max = x.real() > y.real() ? x : y;
-      auto min = x.real() < y.real() ? x : y;
-      auto min_real = min.real();
-      auto max_real = max.real();
+      auto max = cuCrealf(x) > cuCrealf(y) ? x : y;
+      auto min = cuCrealf(x) < cuCrealf(y) ? x : y;
+      auto min_real = cuCrealf(min);
+      auto max_real = cuCrealf(max);
       if (!isfinite(min_real) && (min_real == max_real)) {
         if (min_real < 0) {
           return min;
@@ -148,8 +150,8 @@ struct Maximum {
   __device__ T operator()(T x, T y) {
     if constexpr (cuda::std::is_integral_v<T>) {
       return max(x, y);
-    } else if constexpr (is_complex_v<T>) {
-      if (isnan(x.real()) || isnan(x.imag())) {
+    } else if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      if (isnan(cuCrealf(x)) || isnan(cuCimagf(x))) {
         return x;
       }
       return x > y ? x : y;
@@ -167,8 +169,8 @@ struct Minimum {
   __device__ T operator()(T x, T y) {
     if constexpr (cuda::std::is_integral_v<T>) {
       return min(x, y);
-    } else if constexpr (is_complex_v<T>) {
-      if (isnan(x.real()) || isnan(x.imag())) {
+    } else if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      if (isnan(cuCrealf(x)) || isnan(cuCimagf(x))) {
         return x;
       }
       return x < y ? x : y;
@@ -191,8 +193,8 @@ struct Multiply {
 struct NotEqual {
   template <typename T>
   __device__ bool operator()(T x, T y) {
-    if constexpr (is_complex_v<T>) {
-      return x.real() != y.real() || x.imag() != y.imag();
+    if constexpr (std::is_same_v<T, cuComplex>) {
+      return cuCrealf(x) != cuCrealf(y) || cuCimagf(x) != cuCimagf(y);
     } else {
       return x != y;
     }
@@ -204,12 +206,6 @@ struct Power {
   __device__ T operator()(T base, T exp) {
     if constexpr (cuda::std::is_integral_v<T>) {
       T res = 1;
-      // Raising an integer to a negative power is undefined
-      if constexpr (cuda::std::is_signed_v<T>) {
-        if (exp < 0) {
-          return 0;
-        }
-      }
       while (exp) {
         if (exp & 1) {
           res *= base;
@@ -218,8 +214,19 @@ struct Power {
         base *= base;
       }
       return res;
-    } else if constexpr (is_complex_v<T>) {
-      return pow(base, exp);
+    } else if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      if (base.y == 0 && base.x == 0) {
+        if (isnan(exp.x) || isnan(exp.y)) {
+          auto nan = cuda::std::numeric_limits<float>::quiet_NaN();
+          return make_cuFloatComplex(nan, nan);
+        }
+        return make_cuFloatComplex(0.0, 0.0);
+      }
+      auto x_theta = atan2f(base.y, base.x);
+      auto x_ln_r = 0.5 * logf(base.x * base.x + base.y * base.y);
+      auto mag = expf(exp.x * x_ln_r - exp.y * x_theta);
+      auto phase = exp.y * x_ln_r + exp.x * x_theta;
+      return make_cuFloatComplex(mag * cosf(phase), mag * sinf(phase));
     } else {
       return powf(base, exp);
     }

mlx/backend/cuda/device/cast_op.cuh

@@ -2,10 +2,10 @@
 
 #pragma once
 
-#include "mlx/backend/cuda/device/complex.cuh"
-
+#include <cuComplex.h>
 #include <cuda_bf16.h>
 #include <cuda_fp16.h>
+#include <thrust/iterator/transform_iterator.h>
 
 namespace mlx::core::cu {
 
@@ -20,43 +20,50 @@ struct CastOp {
 };
 
 // Castings between complex and boolean.
-template <typename T>
-struct CastOp<complex_t<T>, bool> {
+// TODO: Should make a custom complex type.
+template <>
+struct CastOp<cuComplex, bool> {
   static constexpr bool is_castable = true;
 
-  __device__ bool operator()(complex_t<T> x) {
-    return x.real() != 0 && x.imag() != 0;
+  __device__ bool operator()(cuComplex x) {
+    return x.x != 0 && x.y != 0;
   }
 };
 
-template <typename T>
-struct CastOp<bool, complex_t<T>> {
+template <>
+struct CastOp<bool, cuComplex> {
   static constexpr bool is_castable = true;
 
-  __device__ complex_t<T> operator()(bool x) {
-    return x ? complex_t<T>{1, 1} : complex_t<T>{0, 0};
+  __device__ cuComplex operator()(bool x) {
+    return x ? make_cuFloatComplex(1, 1) : make_cuFloatComplex(0, 0);
   }
 };
 
 // Converting a complex number to real number discards the imaginary part.
-template <typename T, typename DstT>
-struct CastOp<complex_t<T>, DstT, cuda::std::enable_if_t<!is_complex_v<DstT>>> {
-  static constexpr bool is_castable = cuda::std::is_convertible_v<T, DstT>;
+template <typename DstT>
+struct CastOp<
+    cuComplex,
+    DstT,
+    cuda::std::enable_if_t<!cuda::std::is_same_v<cuComplex, DstT>>> {
+  static constexpr bool is_castable = cuda::std::is_convertible_v<float, DstT>;
 
-  __device__ DstT operator()(complex_t<T> x) {
-    static_assert(!is_complex_v<DstT>);
-    return static_cast<DstT>(x.real());
+  __device__ DstT operator()(cuComplex x) {
+    static_assert(!cuda::std::is_same_v<cuComplex, DstT>);
+    return static_cast<DstT>(cuCrealf(x));
   }
 };
 
 // Allow converting a real number to complex number.
-template <typename SrcT, typename T>
-struct CastOp<SrcT, complex_t<T>, cuda::std::enable_if_t<!is_complex_v<SrcT>>> {
-  static constexpr bool is_castable = cuda::std::is_convertible_v<SrcT, T>;
+template <typename SrcT>
+struct CastOp<
+    SrcT,
+    cuComplex,
+    cuda::std::enable_if_t<!cuda::std::is_same_v<SrcT, cuComplex>>> {
+  static constexpr bool is_castable = cuda::std::is_convertible_v<SrcT, float>;
 
-  __device__ complex_t<T> operator()(SrcT x) {
-    static_assert(!is_complex_v<SrcT>);
-    return complex_t<T>{static_cast<T>(x), 0};
+  __device__ cuComplex operator()(SrcT x) {
+    static_assert(!cuda::std::is_same_v<SrcT, cuComplex>);
+    return cuComplex{static_cast<float>(x), 0};
   }
 };
 
@@ -81,7 +88,8 @@ struct CastOp<
     SrcT,
     DstT,
     cuda::std::enable_if_t<
-        !cuda::std::is_convertible_v<SrcT, DstT> && !is_complex_v<SrcT> &&
+        !cuda::std::is_convertible_v<SrcT, DstT> &&
+        !cuda::std::is_same_v<SrcT, cuComplex> &&
         (cuda::std::is_same_v<DstT, __half> ||
          cuda::std::is_same_v<DstT, __nv_bfloat16>)>> {
   static constexpr bool is_castable = true;
@@ -96,7 +104,8 @@ struct CastOp<
     SrcT,
     DstT,
     cuda::std::enable_if_t<
-        !cuda::std::is_convertible_v<SrcT, DstT> && !is_complex_v<SrcT> &&
+        !cuda::std::is_convertible_v<SrcT, DstT> &&
+        !cuda::std::is_same_v<DstT, cuComplex> &&
         !cuda::std::is_same_v<DstT, __half> &&
         !cuda::std::is_same_v<DstT, __nv_bfloat16> &&
         (cuda::std::is_same_v<SrcT, __half> ||
@@ -115,4 +124,15 @@ inline __host__ __device__ auto cast_to(SrcT x) {
   return CastOp<SrcT, DstT>{}(x);
 }
 
+// Return an iterator that cast the value to DstT using CastOp.
+template <typename DstT, typename Iterator>
+inline __host__ __device__ auto make_cast_iterator(Iterator it) {
+  using SrcT = typename cuda::std::iterator_traits<Iterator>::value_type;
+  if constexpr (std::is_same_v<SrcT, DstT>) {
+    return it;
+  } else {
+    return thrust::make_transform_iterator(it, CastOp<SrcT, DstT>{});
+  }
+}
+
 } // namespace mlx::core::cu

mlx/backend/cuda/device/cexpf.cuh

@@ -0,0 +1,138 @@
+// Copyright © 2025 Apple Inc.
+// Copyright © 2008-2013 NVIDIA Corporation
+// Copyright © 2013 Filipe RNC Maia
+//
+// 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.
+//
+// Forked from
+// https://github.com/NVIDIA/cccl/blob/main/thrust/thrust/detail/complex/cexpf.h
+
+// TODO: We should use thrust::exp but the thrust header in old CUDA versions
+// can not be used in JIT.
+
+#pragma once
+
+#include <cuComplex.h>
+#include <cuda/std/cstdint>
+
+namespace mlx::core::cu::detail {
+
+using ieee_float_shape_type = union {
+  float value;
+  uint32_t word;
+};
+
+inline __device__ void get_float_word(uint32_t& i, float d) {
+  ieee_float_shape_type gf_u;
+  gf_u.value = (d);
+  (i) = gf_u.word;
+}
+
+inline __device__ void get_float_word(int32_t& i, float d) {
+  ieee_float_shape_type gf_u;
+  gf_u.value = (d);
+  (i) = gf_u.word;
+}
+
+inline __device__ void set_float_word(float& d, uint32_t i) {
+  ieee_float_shape_type sf_u;
+  sf_u.word = (i);
+  (d) = sf_u.value;
+}
+
+inline __device__ float frexp_expf(float x, int* expt) {
+  const uint32_t k = 235;
+  const float kln2 = 162.88958740F;
+
+  float exp_x;
+  uint32_t hx;
+
+  exp_x = expf(x - kln2);
+  get_float_word(hx, exp_x);
+  *expt = (hx >> 23) - (0x7f + 127) + k;
+  set_float_word(exp_x, (hx & 0x7fffff) | ((0x7f + 127) << 23));
+  return exp_x;
+}
+
+inline __device__ cuComplex ldexp_cexpf(cuComplex z, int expt) {
+  float x, y, exp_x, scale1, scale2;
+  int ex_expt, half_expt;
+
+  x = cuCrealf(z);
+  y = cuCimagf(z);
+  exp_x = frexp_expf(x, &ex_expt);
+  expt += ex_expt;
+
+  half_expt = expt / 2;
+  set_float_word(scale1, (0x7f + half_expt) << 23);
+  half_expt = expt - half_expt;
+  set_float_word(scale2, (0x7f + half_expt) << 23);
+
+  return cuComplex{
+      cosf(y) * exp_x * scale1 * scale2, sinf(y) * exp_x * scale1 * scale2};
+}
+
+inline __device__ cuComplex cexpf(const cuComplex& z) {
+  float x, y, exp_x;
+  uint32_t hx, hy;
+
+  const uint32_t exp_ovfl = 0x42b17218, cexp_ovfl = 0x43400074;
+
+  x = cuCrealf(z);
+  y = cuCimagf(z);
+
+  get_float_word(hy, y);
+  hy &= 0x7fffffff;
+
+  /* cexp(x + I 0) = exp(x) + I 0 */
+  if (hy == 0) {
+    return cuComplex{expf(x), y};
+  }
+  get_float_word(hx, x);
+  /* cexp(0 + I y) = cos(y) + I sin(y) */
+  if ((hx & 0x7fffffff) == 0) {
+    return cuComplex{cosf(y), sinf(y)};
+  }
+  if (hy >= 0x7f800000) {
+    if ((hx & 0x7fffffff) != 0x7f800000) {
+      /* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
+      return cuComplex{y - y, y - y};
+    } else if (hx & 0x80000000) {
+      /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
+      return cuComplex{0.0, 0.0};
+    } else {
+      /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
+      return cuComplex{x, y - y};
+    }
+  }
+
+  if (hx >= exp_ovfl && hx <= cexp_ovfl) {
+    /*
+     * x is between 88.7 and 192, so we must scale to avoid
+     * overflow in expf(x).
+     */
+    return ldexp_cexpf(z, 0);
+  } else {
+    /*
+     * Cases covered here:
+     *  -  x < exp_ovfl and exp(x) won't overflow (common case)
+     *  -  x > cexp_ovfl, so exp(x) * s overflows for all s > 0
+     *  -  x = +-Inf (generated by exp())
+     *  -  x = NaN (spurious inexact exception from y)
+     */
+    exp_x = expf(x);
+    return cuComplex{exp_x * cosf(y), exp_x * sinf(y)};
+  }
+}
+
+} // namespace mlx::core::cu::detail

mlx/backend/cuda/device/complex.cuh

@@ -1,60 +0,0 @@
-// Copyright © 2025 Apple Inc.
-
-#pragma once
-
-// Make multiplication and division faster.
-#define LIBCUDACXX_ENABLE_SIMPLIFIED_COMPLEX_OPERATIONS
-
-#include <cuda/std/complex>
-#include <cuda/std/type_traits>
-
-namespace mlx::core::cu {
-
-// TODO: Consider using a faster implementation as cuda::std::complex has to
-// conform to C++ standard.
-template <typename T>
-using complex_t = cuda::std::complex<T>;
-
-using complex64_t = complex_t<float>;
-using complex128_t = complex_t<double>;
-
-template <typename T>
-struct is_complex : cuda::std::false_type {};
-
-template <typename T>
-struct is_complex<cuda::std::complex<T>> : cuda::std::true_type {};
-
-template <typename T>
-inline constexpr bool is_complex_v = is_complex<T>::value;
-
-// cuda::std::complex is missing some operators.
-template <typename T>
-inline __host__ __device__ complex_t<T> operator%(
-    complex_t<T> a,
-    complex_t<T> b) {
-  T r = a.real() - floor(a.real() / b.real()) * b.real();
-  T i = a.imag() - floor(a.imag() / b.imag()) * b.imag();
-  return complex_t<T>{r, i};
-}
-
-template <typename T>
-inline __host__ __device__ bool operator>(complex_t<T> a, complex_t<T> b) {
-  return (a.real() > b.real()) || (a.real() == b.real() && a.imag() > b.imag());
-}
-
-template <typename T>
-inline __host__ __device__ bool operator<(complex_t<T> a, complex_t<T> b) {
-  return operator>(b, a);
-}
-
-template <typename T>
-inline __host__ __device__ bool operator<=(complex_t<T> a, complex_t<T> b) {
-  return !(a > b);
-}
-
-template <typename T>
-inline __host__ __device__ bool operator>=(complex_t<T> a, complex_t<T> b) {
-  return !(a < b);
-}
-
-} // namespace mlx::core::cu

mlx/backend/cuda/device/cucomplex_math.cuh

@@ -0,0 +1,240 @@
+// Copyright © 2025 Apple Inc.
+// Copyright © 2017-2024 The Simons Foundation, Inc.
+//
+// FINUFFT is 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.
+//
+// Forked from
+// https://github.com/flatironinstitute/finufft/blob/main/include/cufinufft/contrib/helper_math.h
+
+#pragma once
+
+#include <cuComplex.h>
+
+// This header provides some helper functions for cuComplex types.
+// It mainly wraps existing CUDA implementations to provide operator overloads
+// e.g. cuAdd, cuSub, cuMul, cuDiv, cuCreal, cuCimag, cuCabs, cuCarg, cuConj are
+// all provided by CUDA
+
+__forceinline__ __host__ __device__ cuDoubleComplex
+operator+(const cuDoubleComplex& a, const cuDoubleComplex& b) {
+  return cuCadd(a, b);
+}
+
+__forceinline__ __host__ __device__ cuDoubleComplex
+operator-(const cuDoubleComplex& a, const cuDoubleComplex& b) {
+  return cuCsub(a, b);
+}
+
+__forceinline__ __host__ __device__ cuDoubleComplex
+operator*(const cuDoubleComplex& a, const cuDoubleComplex& b) {
+  return cuCmul(a, b);
+}
+
+__forceinline__ __host__ __device__ cuDoubleComplex
+operator/(const cuDoubleComplex& a, const cuDoubleComplex& b) {
+  return cuCdiv(a, b);
+}
+
+__forceinline__ __host__ __device__ cuDoubleComplex
+operator%(const cuDoubleComplex& a, const cuDoubleComplex& b) {
+  double r = cuCreal(a) - (floorf(cuCreal(a) / cuCreal(b)) * cuCreal(b));
+  double i = cuCimag(a) - (floorf(cuCimag(a) / cuCimag(b)) * cuCimag(b));
+  return make_cuDoubleComplex(r, i);
+}
+
+__forceinline__ __host__ __device__ bool operator==(
+    const cuDoubleComplex& a,
+    const cuDoubleComplex& b) {
+  return cuCreal(a) == cuCreal(b) && cuCimag(a) == cuCimag(b);
+}
+
+__forceinline__ __host__ __device__ bool operator!=(
+    const cuDoubleComplex& a,
+    const cuDoubleComplex& b) {
+  return !(a == b);
+}
+
+__forceinline__ __host__ __device__ bool operator>(
+    const cuDoubleComplex& a,
+    const cuDoubleComplex& b) {
+  double mag_a = sqrt(cuCreal(a) * cuCreal(a) + cuCimag(a) * cuCimag(a));
+  double mag_b = sqrt(cuCreal(b) * cuCreal(b) + cuCimag(b) * cuCimag(b));
+  return mag_a > mag_b;
+}
+
+__forceinline__ __host__ __device__ bool operator>=(
+    const cuDoubleComplex& a,
+    const cuDoubleComplex& b) {
+  return a > b || a == b;
+}
+
+__forceinline__ __host__ __device__ bool operator<(
+    const cuDoubleComplex& a,
+    const cuDoubleComplex& b) {
+  return b > a;
+}
+
+__forceinline__ __host__ __device__ bool operator<=(
+    const cuDoubleComplex& a,
+    const cuDoubleComplex& b) {
+  return b > a || a == b;
+}
+
+__forceinline__ __host__ __device__ cuDoubleComplex
+operator+(const cuDoubleComplex& a, double b) {
+  return make_cuDoubleComplex(cuCreal(a) + b, cuCimag(a));
+}
+
+__forceinline__ __host__ __device__ cuDoubleComplex
+operator+(double a, const cuDoubleComplex& b) {
+  return make_cuDoubleComplex(a + cuCreal(b), cuCimag(b));
+}
+
+__forceinline__ __host__ __device__ cuDoubleComplex
+operator-(const cuDoubleComplex& a, double b) {
+  return make_cuDoubleComplex(cuCreal(a) - b, cuCimag(a));
+}
+
+__forceinline__ __host__ __device__ cuDoubleComplex
+operator-(double a, const cuDoubleComplex& b) {
+  return make_cuDoubleComplex(a - cuCreal(b), -cuCimag(b));
+}
+
+__forceinline__ __host__ __device__ cuDoubleComplex
+operator*(const cuDoubleComplex& a, double b) {
+  return make_cuDoubleComplex(cuCreal(a) * b, cuCimag(a) * b);
+}
+
+__forceinline__ __host__ __device__ cuDoubleComplex
+operator*(double a, const cuDoubleComplex& b) {
+  return make_cuDoubleComplex(a * cuCreal(b), a * cuCimag(b));
+}
+
+__forceinline__ __host__ __device__ cuDoubleComplex
+operator/(const cuDoubleComplex& a, double b) {
+  return make_cuDoubleComplex(cuCreal(a) / b, cuCimag(a) / b);
+}
+
+__forceinline__ __host__ __device__ cuDoubleComplex
+operator/(double a, const cuDoubleComplex& b) {
+  double denom = cuCreal(b) * cuCreal(b) + cuCimag(b) * cuCimag(b);
+  return make_cuDoubleComplex(
+      (a * cuCreal(b)) / denom, (-a * cuCimag(b)) / denom);
+}
+
+__forceinline__ __host__ __device__ cuFloatComplex
+operator+(const cuFloatComplex& a, const cuFloatComplex& b) {
+  return cuCaddf(a, b);
+}
+
+__forceinline__ __host__ __device__ cuFloatComplex
+operator-(const cuFloatComplex& a, const cuFloatComplex& b) {
+  return cuCsubf(a, b);
+}
+
+__forceinline__ __host__ __device__ cuFloatComplex
+operator*(const cuFloatComplex& a, const cuFloatComplex& b) {
+  return cuCmulf(a, b);
+}
+
+__forceinline__ __host__ __device__ cuFloatComplex
+operator/(const cuFloatComplex& a, const cuFloatComplex& b) {
+  return cuCdivf(a, b);
+}
+
+__forceinline__ __host__ __device__ cuFloatComplex
+operator%(const cuFloatComplex& a, const cuFloatComplex& b) {
+  float r = cuCrealf(a) - (floorf(cuCrealf(a) / cuCrealf(b)) * cuCrealf(b));
+  float i = cuCimagf(a) - (floorf(cuCimagf(a) / cuCimagf(b)) * cuCimagf(b));
+  return make_cuFloatComplex(r, i);
+}
+
+__forceinline__ __host__ __device__ bool operator==(
+    const cuFloatComplex& a,
+    const cuFloatComplex& b) {
+  return cuCrealf(a) == cuCrealf(b) && cuCimagf(a) == cuCimagf(b);
+}
+
+__forceinline__ __host__ __device__ bool operator!=(
+    const cuFloatComplex& a,
+    const cuFloatComplex& b) {
+  return !(a == b);
+}
+
+__forceinline__ __host__ __device__ bool operator>(
+    const cuFloatComplex& a,
+    const cuFloatComplex& b) {
+  float mag_a = sqrt(cuCrealf(a) * cuCrealf(a) + cuCimagf(a) * cuCimagf(a));
+  float mag_b = sqrt(cuCrealf(b) * cuCrealf(b) + cuCimagf(b) * cuCimagf(b));
+  return mag_a > mag_b;
+}
+
+__forceinline__ __host__ __device__ bool operator>=(
+    const cuFloatComplex& a,
+    const cuFloatComplex& b) {
+  return a > b || a == b;
+}
+
+__forceinline__ __host__ __device__ bool operator<(
+    const cuFloatComplex& a,
+    const cuFloatComplex& b) {
+  return b > a;
+}
+
+__forceinline__ __host__ __device__ bool operator<=(
+    const cuFloatComplex& a,
+    const cuFloatComplex& b) {
+  return b > a || a == b;
+}
+
+__forceinline__ __host__ __device__ cuFloatComplex
+operator+(const cuFloatComplex& a, float b) {
+  return make_cuFloatComplex(cuCrealf(a) + b, cuCimagf(a));
+}
+
+__forceinline__ __host__ __device__ cuFloatComplex
+operator+(float a, const cuFloatComplex& b) {
+  return make_cuFloatComplex(a + cuCrealf(b), cuCimagf(b));
+}
+
+__forceinline__ __host__ __device__ cuFloatComplex
+operator-(const cuFloatComplex& a, float b) {
+  return make_cuFloatComplex(cuCrealf(a) - b, cuCimagf(a));
+}
+
+__forceinline__ __host__ __device__ cuFloatComplex
+operator-(float a, const cuFloatComplex& b) {
+  return make_cuFloatComplex(a - cuCrealf(b), -cuCimagf(b));
+}
+
+__forceinline__ __host__ __device__ cuFloatComplex
+operator*(const cuFloatComplex& a, float b) {
+  return make_cuFloatComplex(cuCrealf(a) * b, cuCimagf(a) * b);
+}
+
+__forceinline__ __host__ __device__ cuFloatComplex
+operator*(float a, const cuFloatComplex& b) {
+  return make_cuFloatComplex(a * cuCrealf(b), a * cuCimagf(b));
+}
+
+__forceinline__ __host__ __device__ cuFloatComplex
+operator/(const cuFloatComplex& a, float b) {
+  return make_cuFloatComplex(cuCrealf(a) / b, cuCimagf(a) / b);
+}
+
+__forceinline__ __host__ __device__ cuFloatComplex
+operator/(float a, const cuFloatComplex& b) {
+  float denom = cuCrealf(b) * cuCrealf(b) + cuCimagf(b) * cuCimagf(b);
+  return make_cuFloatComplex(
+      (a * cuCrealf(b)) / denom, (-a * cuCimagf(b)) / denom);
+}

mlx/backend/cuda/device/unary_ops.cuh

@@ -2,6 +2,8 @@
 
 #pragma once
 
+#include "mlx/backend/cuda/device/cexpf.cuh"
+#include "mlx/backend/cuda/device/cucomplex_math.cuh"
 #include "mlx/backend/cuda/device/fp16_math.cuh"
 #include "mlx/backend/cuda/device/utils.cuh"
 
@@ -14,6 +16,8 @@ struct Abs {
   __device__ T operator()(T x) {
     if constexpr (cuda::std::is_unsigned_v<T>) {
       return x;
+    } else if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      return {sqrt(cuCrealf(x) * cuCrealf(x) + cuCimagf(x) * cuCimagf(x)), 0};
     } else {
       return abs(x);
     }
@@ -25,6 +29,8 @@ struct ArcCos {
   __device__ T operator()(T x) {
     return acos(x);
   }
+
+  __device__ cuComplex operator()(cuComplex x);
 };
 
 struct ArcCosh {
@@ -39,6 +45,8 @@ struct ArcSin {
   __device__ T operator()(T x) {
     return asin(x);
   }
+
+  __device__ cuComplex operator()(cuComplex x);
 };
 
 struct ArcSinh {
@@ -53,6 +61,8 @@ struct ArcTan {
   __device__ T operator()(T x) {
     return atan(x);
   }
+
+  __device__ cuComplex operator()(cuComplex x);
 };
 
 struct ArcTanh {
@@ -74,8 +84,6 @@ struct Ceil {
   __device__ T operator()(T x) {
     if constexpr (cuda::std::is_integral_v<T>) {
       return x;
-    } else if constexpr (is_complex_v<T>) {
-      return T{ceil(x.real()), ceil(x.imag())};
     } else {
       return ceil(x);
     }
@@ -83,24 +91,35 @@ struct Ceil {
 };
 
 struct Conjugate {
-  template <typename T>
-  __device__ complex_t<T> operator()(complex_t<T> x) {
-    return conj(x);
+  __device__ cuComplex operator()(cuComplex x) {
+    return {cuCrealf(x), -cuCimagf(x)};
   }
 };
 
 struct Cos {
   template <typename T>
   __device__ T operator()(T x) {
+    if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      return {
+          cos(cuCrealf(x)) * cosh(cuCimagf(x)),
+          -sin(cuCrealf(x)) * sinh(cuCimagf(x))};
+    } else {
       return cos(x);
     }
+  }
 };
 
 struct Cosh {
   template <typename T>
   __device__ T operator()(T x) {
+    if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      return {
+          cosh(cuCrealf(x)) * cos(cuCimagf(x)),
+          sinh(cuCrealf(x)) * sin(cuCimagf(x))};
+    } else {
       return cosh(x);
     }
+  }
 };
 
 struct Erf {
@@ -132,8 +151,12 @@ struct ErfInv {
 struct Exp {
   template <typename T>
   __device__ T operator()(T x) {
+    if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      return detail::cexpf(x);
+    } else {
       return exp(x);
     }
+  }
 };
 
 struct Expm1 {
@@ -154,8 +177,6 @@ struct Floor {
   __device__ T operator()(T x) {
     if constexpr (cuda::std::is_integral_v<T>) {
       return x;
-    } else if constexpr (is_complex_v<T>) {
-      return T{floor(x.real()), floor(x.imag())};
     } else {
       return floor(x);
     }
@@ -163,25 +184,30 @@ struct Floor {
 };
 
 struct Imag {
-  template <typename T>
-  __device__ auto operator()(complex_t<T> x) {
-    return x.imag();
+  __device__ float operator()(cuComplex x) {
+    return cuCimagf(x);
   }
 };
 
 struct Log {
   template <typename T>
   __device__ T operator()(T x) {
+    if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      auto r = log(cuCrealf(Abs{}(x)));
+      auto i = atan2f(cuCimagf(x), cuCrealf(x));
+      return {r, i};
+    } else {
       return log(x);
     }
+  }
 };
 
 struct Log2 {
   template <typename T>
   __device__ T operator()(T x) {
-    if constexpr (is_complex_v<T>) {
+    if constexpr (cuda::std::is_same_v<T, cuComplex>) {
       auto y = Log{}(x);
-      return {y.real() / CUDART_LN2_F, y.imag() / CUDART_LN2_F};
+      return {cuCrealf(y) / CUDART_LN2_F, cuCimagf(y) / CUDART_LN2_F};
     } else {
       return log2(x);
     }
@@ -191,17 +217,23 @@ struct Log2 {
 struct Log10 {
   template <typename T>
   __device__ T operator()(T x) {
+    if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      auto y = Log{}(x);
+      return {cuCrealf(y) / CUDART_LNT_F, cuCimagf(y) / CUDART_LNT_F};
+      return y;
+    } else {
       return log10(x);
     }
+  }
 };
 
 struct Log1p {
   template <typename T>
   __device__ T operator()(T z) {
-    if constexpr (is_complex_v<T>) {
-      float x = z.real();
-      float y = z.imag();
-      float zabs = Abs{}(z).real();
+    if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      float x = cuCrealf(z);
+      float y = cuCimagf(z);
+      float zabs = cuCrealf(Abs{}(z));
       float theta = atan2f(y, x + 1);
       if (zabs < 0.5f) {
         float r = x * (2 + x) + y * y;
@@ -228,8 +260,8 @@ struct LogicalNot {
 struct Negative {
   template <typename T>
   __device__ T operator()(T x) {
-    if constexpr (is_complex_v<T>) {
-      return T{0, 0} - x;
+    if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      return 0 - x;
     } else {
       return -x;
     }
@@ -237,17 +269,16 @@ struct Negative {
 };
 
 struct Real {
-  template <typename T>
-  __device__ auto operator()(complex_t<T> x) {
-    return x.real();
+  __device__ float operator()(cuComplex x) {
+    return cuCrealf(x);
   }
 };
 
 struct Round {
   template <typename T>
   __device__ T operator()(T x) {
-    if constexpr (is_complex_v<T>) {
-      return {rint(x.real()), rint(x.imag())};
+    if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      return {rint(cuCrealf(x)), rint(cuCimagf(x))};
     } else {
       return rint(x);
     }
@@ -267,8 +298,8 @@ struct Sign {
   __device__ T operator()(T x) {
     if constexpr (cuda::std::is_unsigned_v<T>) {
       return x != 0;
-    } else if constexpr (is_complex_v<T>) {
-      if (x.real() == 0 && x.imag() == 0) {
+    } else if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      if (cuCrealf(x) == 0 && cuCimagf(x) == 0) {
         return x;
       } else {
         return x / Abs()(x);
@@ -284,15 +315,27 @@ struct Sign {
 struct Sin {
   template <typename T>
   __device__ T operator()(T x) {
+    if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      return {
+          sin(cuCrealf(x)) * cosh(cuCimagf(x)),
+          cos(cuCrealf(x)) * sinh(cuCimagf(x))};
+    } else {
       return sin(x);
     }
+  }
 };
 
 struct Sinh {
   template <typename T>
   __device__ T operator()(T x) {
+    if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      return {
+          sinh(cuCrealf(x)) * cos(cuCimagf(x)),
+          cosh(cuCrealf(x)) * sin(cuCimagf(x))};
+    } else {
       return sinh(x);
     }
+  }
 };
 
 struct Square {
@@ -307,31 +350,77 @@ struct Sqrt {
   __device__ T operator()(T x) {
     return sqrt(x);
   }
+
+  __device__ cuComplex operator()(cuComplex x) {
+    auto xr = cuCrealf(x);
+    auto xi = cuCimagf(x);
+    if (xr == 0.0f && xi == 0.0f) {
+      return {0.0f, 0.0f};
+    }
+    auto r = cuCrealf(Abs{}(x));
+    auto a = sqrt((r + xr) / 2.0f);
+    auto b_abs = sqrt((r - xr) / 2.0f);
+    auto b = copysign(b_abs, xi);
+    return {a, b};
+  }
 };
 
 struct Rsqrt {
   template <typename T>
   __device__ T operator()(T x) {
-    if constexpr (is_complex_v<T>) {
-      return 1.0f / Sqrt{}(x);
-    } else {
     return rsqrt(x);
   }
+  __device__ cuComplex operator()(cuComplex x) {
+    return 1.0f / Sqrt{}(x);
   }
 };
 
 struct Tan {
   template <typename T>
   __device__ T operator()(T x) {
+    if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      float tan_a = tan(cuCrealf(x));
+      float tanh_b = tanh(cuCimagf(x));
+      float t1 = tan_a * tanh_b;
+      float denom = 1. + t1 * t1;
+      return {(tan_a - tanh_b * t1) / denom, (tanh_b + tan_a * t1) / denom};
+    } else {
       return tan(x);
     }
+  }
 };
 
 struct Tanh {
   template <typename T>
   __device__ T operator()(T x) {
+    if constexpr (cuda::std::is_same_v<T, cuComplex>) {
+      float tanh_a = tanh(cuCrealf(x));
+      float tan_b = tan(cuCimagf(x));
+      float t1 = tanh_a * tan_b;
+      float denom = 1. + t1 * t1;
+      return {(tanh_a + tan_b * t1) / denom, (tan_b - tanh_a * t1) / denom};
+    } else {
       return tanh(x);
     }
+  }
+};
+
+inline __device__ cuComplex ArcCos::operator()(cuComplex x) {
+  auto i = cuComplex{0.0, 1.0};
+  auto y = Log{}(x + i * Sqrt{}(1.0 - x * x));
+  return {cuCimagf(y), -cuCrealf(y)};
+};
+
+inline __device__ cuComplex ArcSin::operator()(cuComplex x) {
+  auto i = cuComplex{0.0f, 1.0f};
+  auto y = Log{}(i * x + Sqrt{}(1.0f - x * x));
+  return {cuCimagf(y), -cuCrealf(y)};
+};
+
+inline __device__ cuComplex ArcTan::operator()(cuComplex x) {
+  auto i = cuComplex{0.0f, 1.0f};
+  auto ix = i * x;
+  return (1.0f / cuComplex{0.0f, 2.0f}) * Log{}((1.0f + ix) / (1.0f - ix));
 };
 
 } // namespace mlx::core::cu

mlx/backend/cuda/device/utils.cuh

@@ -8,9 +8,9 @@
 
 #pragma once
 
-#include "mlx/backend/cuda/device/complex.cuh"
 #include "mlx/backend/cuda/device/config.h"
 
+#include <cuComplex.h>
 #include <cuda_bf16.h>
 #include <cuda_fp16.h>
 #include <cuda/std/array>
@@ -32,135 +32,21 @@ using Strides = cuda::std::array<int64_t, MAX_NDIM>;
 template <typename T, int N>
 struct alignas(sizeof(T) * N) AlignedVector {
   T val[N];
-
-  __device__ T& operator[](int i) {
-    return val[i];
-  }
-
-  __device__ T operator[](int i) const {
-    return val[i];
-  }
 };
 
-template <int N, typename T>
-inline __host__ __device__ bool is_aligned(T* x) {
-  return (reinterpret_cast<uintptr_t>(x) % (N * sizeof(T))) == 0;
-}
-
-template <int N, typename T>
-inline __device__ AlignedVector<T, N> unsafe_load_vector(
-    const T* ptr,
-    uint32_t offset) {
-  auto* from = reinterpret_cast<const AlignedVector<T, N>*>(ptr);
-  return from[offset];
-}
-
 template <int N, typename T>
 inline __device__ AlignedVector<T, N> load_vector(
     const T* ptr,
     uint32_t offset) {
-  if (is_aligned<N>(ptr)) {
   auto* from = reinterpret_cast<const AlignedVector<T, N>*>(ptr);
   return from[offset];
-  } else {
-    AlignedVector<T, N> v;
-#pragma unroll
-    for (int i = 0; i < N; ++i) {
-      v[i] = ptr[offset * N + i];
-    }
-    return v;
-  }
-}
-
-template <int N, typename T, typename SizeT>
-inline __device__ AlignedVector<T, N>
-load_vector(const T* ptr, uint32_t offset, SizeT size, T fallback) {
-  if (is_aligned<N>(ptr) && (offset + 1) * N <= size) {
-    auto* from = reinterpret_cast<const AlignedVector<T, N>*>(ptr);
-    return from[offset];
-  } else {
-    AlignedVector<T, N> v;
-#pragma unroll
-    for (int i = 0; i < N; ++i) {
-      v[i] = (N * offset + i) < size ? ptr[offset * N + i] : fallback;
-    }
-    return v;
-  }
-}
-
-template <int N, typename T, typename SizeT>
-inline __device__ AlignedVector<T, N> load_vector(
-    const T* ptr,
-    uint32_t offset,
-    SizeT size,
-    int64_t stride,
-    T fallback) {
-  if (is_aligned<N>(ptr) && stride == 1 && (offset + 1) * N <= size) {
-    auto* from = reinterpret_cast<const AlignedVector<T, N>*>(ptr);
-    return from[offset];
-  } else {
-    AlignedVector<T, N> v;
-#pragma unroll
-    for (int i = 0; i < N; ++i) {
-      v[i] =
-          (N * offset + i) < size ? ptr[stride * (offset * N + i)] : fallback;
-    }
-    return v;
-  }
-}
-
-template <int N, typename T>
-inline __device__ void
-unsafe_store_vector(T* ptr, uint32_t offset, const AlignedVector<T, N>& vec) {
-  auto* to = reinterpret_cast<AlignedVector<T, N>*>(ptr);
-  to[offset] = vec;
 }
 
 template <int N, typename T>
 inline __device__ void
 store_vector(T* ptr, uint32_t offset, const AlignedVector<T, N>& vec) {
-  if (is_aligned<N>(ptr)) {
   auto* to = reinterpret_cast<AlignedVector<T, N>*>(ptr);
   to[offset] = vec;
-  } else {
-#pragma unroll
-    for (int i = 0; i < N; ++i) {
-      ptr[offset * N + i] = vec[i];
-    }
-  }
-}
-
-template <int N, typename T, typename SizeT>
-inline __device__ void store_vector(
-    T* ptr,
-    uint32_t offset,
-    const AlignedVector<T, N>& vec,
-    SizeT size) {
-  if (is_aligned<N>(ptr) && (offset + 1) * N <= size) {
-    auto* to = reinterpret_cast<AlignedVector<T, N>*>(ptr);
-    to[offset] = vec;
-  } else {
-    for (int i = 0; (offset * N + i) < size && i < N; ++i) {
-      ptr[offset * N + i] = vec[i];
-    }
-  }
-}
-
-template <int N, typename T, typename SizeT>
-inline __device__ void store_vector(
-    T* ptr,
-    uint32_t offset,
-    const AlignedVector<T, N>& vec,
-    SizeT size,
-    int64_t stride) {
-  if (is_aligned<N>(ptr) && (offset + 1) * N <= size && stride == 1) {
-    auto* to = reinterpret_cast<AlignedVector<T, N>*>(ptr);
-    to[offset] = vec;
-  } else {
-    for (int i = 0; (offset * N + i) < size && i < N; ++i) {
-      ptr[stride * (offset * N + i)] = vec[i];
-    }
-  }
 }
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -241,13 +127,13 @@ struct Limits<bool> {
   }
 };
 
-template <typename T>
-struct Limits<complex_t<T>> {
-  static constexpr __host__ __device__ complex_t<T> max() {
-    return {Limits<T>::max(), Limits<T>::max()};
+template <>
+struct Limits<cuComplex> {
+  static constexpr __host__ __device__ cuComplex max() {
+    return {Limits<float>::max(), Limits<float>::max()};
   }
-  static constexpr __host__ __device__ complex_t<T> min() {
-    return {Limits<T>::min(), Limits<T>::min()};
+  static constexpr __host__ __device__ cuComplex min() {
+    return {Limits<float>::min(), Limits<float>::min()};
   }
 };
 
@@ -318,8 +204,20 @@ inline __host__ __device__ cuda::std::tuple<IdxT, IdxT, IdxT> elem_to_loc_nd(
   return cuda::std::make_tuple(a_loc, b_loc, c_loc);
 }
 
+// Optimized version when ndim is larger than 4.
 template <typename IdxT = int64_t>
-inline __host__ __device__ cuda::std::tuple<IdxT, IdxT> elem_to_loc(
+inline __host__ __device__ IdxT
+elem_to_loc_4d(IdxT elem, const int* shape, const int64_t* strides, int ndim) {
+  IdxT loc = 0;
+  for (int i = ndim - 1; i >= 0; --i) {
+    loc += (elem % shape[i]) * IdxT(strides[i]);
+    elem /= shape[i];
+  }
+  return loc;
+}
+
+template <typename IdxT = int64_t>
+inline __host__ __device__ cuda::std::tuple<IdxT, IdxT> elem_to_loc_4d(
     IdxT elem,
     const int* shape,
     const int64_t* a_strides,
@@ -337,7 +235,7 @@ inline __host__ __device__ cuda::std::tuple<IdxT, IdxT> elem_to_loc(
 }
 
 template <typename IdxT = int64_t>
-inline __host__ __device__ cuda::std::tuple<IdxT, IdxT, IdxT> elem_to_loc(
+inline __host__ __device__ cuda::std::tuple<IdxT, IdxT, IdxT> elem_to_loc_4d(
     IdxT elem,
     const int* shape,
     const int64_t* a_strides,