diff --git a/.githash b/.githash
index 6cdfeea5..2c5d87fb 100644
--- a/.githash
+++ b/.githash
@@ -1 +1 @@
-dec047f1bd1c8287513c6c437f946982e516ccd4
+e4cc4b6a778aa828e5647803f594222053797f54
diff --git a/Project.toml b/Project.toml
index 535b8060..fb5c3074 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,8 +4,12 @@ version = "0.1.0"
 
 [deps]
 CNPreferences = "3e078157-ea10-49d5-bf32-908f777cd46f"
+CUDACore = "bd0ed864-bdfe-4181-a5ed-ce625a5fdea2"
+CUDATools = "9ec180c6-1c07-47c7-9e6e-ebefa4d1f6d0"
 CxxWrap = "1f15a43c-97ca-5a2a-ae31-89f07a497df4"
+ExpressionExplorer = "21656369-7473-754a-2065-74616d696c43"
 JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
+KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
 Legate = "1238f2cf-6593-4d60-9aca-2f5364e49909"
 LegatePreferences = "8028f36a-2b64-49e9-aa04-2d0933fd2ed9"
 Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
@@ -20,17 +24,14 @@ cunumeric_jl_wrapper_jll = "49048992-29d2-5fd1-994f-9cecf112d624"
 cupynumeric_jll = "2862d674-414d-5b0b-a494-b21f8deca547"
 libcxxwrap_julia_jll = "3eaa8342-bff7-56a5-9981-c04077f7cee7"
 
-[weakdeps]
-CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
-
-[extensions]
-CUDAExt = "CUDA"
-
 [compat]
 CNPreferences = "0.1.2"
-CUDA = "5.9"
+CUDACore = "6.0.0"
+CUDATools = "6.0.0"
 CxxWrap = "0.17"
+ExpressionExplorer = "1.1.4"
 JuliaFormatter = "2.3.0"
+KernelAbstractions = "0.9.41"
 Legate = "0.1.2"
 LegatePreferences = "0.1.6"
 MacroTools = "0.5.16"
diff --git a/benchmark/gray_scott_fusion.jl b/benchmark/gray_scott_fusion.jl
new file mode 100644
index 00000000..3ce91151
--- /dev/null
+++ b/benchmark/gray_scott_fusion.jl
@@ -0,0 +1,159 @@
+using cuNumeric
+using Random
+using StatsBase
+
+cuNumeric.versioninfo()
+
+Random.seed!(1234)
+
+struct Params{T}
+    dx::T
+    dt::T
+    c_u::T
+    c_v::T
+    f::T
+    k::T
+
+    function Params(dx=1.0f0, c_u=1.0f0, c_v=0.3f0, f=0.03f0, k=0.06f0)
+        new{Float32}(dx, dx/5, c_u, c_v, f, k)
+    end
+end
+
+function bc!(u_new, v_new, u, v)
+    u_new[:, 1] = u[:, end - 1]
+    u_new[:, end] = u[:, 2]
+    u_new[1, :] = u[end - 1, :]
+    u_new[end, :] = u[2, :]
+    v_new[:, 1] = v[:, end - 1]
+    v_new[:, end] = v[:, 2]
+    v_new[1, :] = v[end - 1, :]
+    v_new[end, :] = v[2, :]
+end
+
+function step!(u, v, u_new, v_new, args::Params)
+
+    dx_sq = args.dx^2
+
+    # AVOID SOME ISSUES WITH ^ operator
+    v_sq = v[2:(end - 1), 2:(end - 1)] .* v[2:(end - 1), 2:(end - 1)]
+
+    # calculate F_u and F_v functions
+    F_u = (
+        (-u[2:(end - 1), 2:(end - 1)] .* (v_sq)) .+
+        args.f*(1.0f0 .- u[2:(end - 1), 2:(end - 1)])
+    )
+    F_v = (
+        (u[2:(end - 1), 2:(end - 1)] .* (v_sq)) .-
+        (args.f+args.k) .* v[2:(end - 1), 2:(end - 1)]
+    )
+    # 2-D Laplacian of f using array slicing, excluding boundaries
+    # For an N x N array f, f_lap is the Nend x Nend array in the "middle"
+    u_lap = (
+        (u[3:end, 2:(end - 1)] - 2*u[2:(end - 1), 2:(end - 1)] + u[1:(end - 2), 2:(end - 1)]) ./
+        dx_sq
+        +
+        (u[2:(end - 1), 3:end] - 2*u[2:(end - 1), 2:(end - 1)] + u[2:(end - 1), 1:(end - 2)]) ./
+        dx_sq
+    )
+    v_lap = (
+        (v[3:end, 2:(end - 1)] - 2*v[2:(end - 1), 2:(end - 1)] + v[1:(end - 2), 2:(end - 1)]) ./
+        dx_sq
+        +
+        (v[2:(end - 1), 3:end] - 2*v[2:(end - 1), 2:(end - 1)] + v[2:(end - 1), 1:(end - 2)]) ./
+        dx_sq
+    )
+
+    # Forward-Euler time step for all points except the boundaries
+    u_new[2:(end - 1), 2:(end - 1)] =
+        ((args.c_u * u_lap) + F_u) * args.dt + u[2:(end - 1), 2:(end - 1)]
+    v_new[2:(end - 1), 2:(end - 1)] =
+        ((args.c_v * v_lap) + F_v) * args.dt + v[2:(end - 1), 2:(end - 1)]
+
+    # Apply periodic boundary conditions
+    bc!(u_new, v_new, u, v)
+end
+
+function initialize(N)
+    u = cuNumeric.ones(N, N)
+    v = cuNumeric.zeros(N, N)
+    u_new = cuNumeric.zeros(N, N)
+    v_new = cuNumeric.zeros(N, N)
+
+    rand_frac = round(Int, 0.15 * N)
+
+    # Make julia arrays first for reproducability, I don't think
+    # we have a way to set cuNumeric's random seed.
+    u_tmp = rand(Float32, rand_frac, rand_frac)
+    v_tmp = rand(Float32, rand_frac, rand_frac)
+
+    # Copy over values to cuNumeric arrays
+    allowscalar() do
+        for i in 1:rand_frac, j in 1:rand_frac
+            u[i,j] = u_tmp[i,j]
+            v[i,j] = v_tmp[i,j]
+        end
+    end
+
+    GC.gc() # remove any intermediates
+
+    return u, v, u_new, v_new
+end
+
+# Re-iplement loop without plotting, and add timing
+function gray_scott(N::Int, n_steps::Int, n_warmup::Int)
+
+    args = Params()
+    u, v, u_new, v_new = initialize(N)
+    @info "Initialized Arrays"
+
+    for n in 1:n_warmup
+        @info "Warmup Step $n / $n_warmup"
+        step!(u, v, u_new, v_new, args)
+        u, u_new = u_new, u
+        v, v_new = v_new, v
+    end
+
+    step_times_us = Vector{Float64}(undef, n_steps)
+    for n in 1:n_steps
+        @info "Step $n / $n_steps"
+        start_time = get_time_microseconds()
+        step!(u, v, u_new, v_new, args)
+        u, u_new = u_new, u
+        v, v_new = v_new, v
+        step_times_us[n] = get_time_microseconds() - start_time
+
+        if n%10 == 0
+            GC.gc()
+        end
+    end
+
+    return u, v, step_times_us
+end
+
+gpus = parse(Int, ARGS[1])
+N = parse(Int, ARGS[2])
+n_steps = parse(Int, ARGS[3])
+n_warmup = parse(Int, ARGS[4])
+
+println(
+    "[cuNumeric]  Gray-Scott benchmark on $(N)x$(N) grid for $(n_steps) iterations, $(n_warmup) warmups"
+)
+
+u, v, step_times_us = gray_scott(N, n_steps, n_warmup)
+
+step_times_ms = step_times_us ./ 1000.0
+mean_time_ms = mean(step_times_ms) 
+std_time_ms = std(step_times_ms)
+
+total_time_s = sum(step_times_ms) / 1000.0
+iter_per_second = n_steps / total_time_s
+
+println("Mean Step Time: $(mean_time_ms) ms")
+println("StdDev Step Time: $(std_time_ms) ms")
+println("Min Step Time: $(minimum(step_times_ms)) ms")
+println("Max Step Time: $(maximum(step_times_ms)) ms")
+println("Throughput: $(iter_per_second) ")
+
+# open("./gray_scott.csv", "a") do io
+#     @printf(io, "%s,%d,%d,%d,%.6f,%.6f\n", "cunumeric", gpus, N, M, mean_time_ms, gflops)
+# end
\ No newline at end of file
diff --git a/deps/build.jl b/deps/build.jl
index d53e576f..7899822c 100644
--- a/deps/build.jl
+++ b/deps/build.jl
@@ -20,6 +20,11 @@
 using Preferences
 using Legate
 using CNPreferences
+using CUDACore: CUDACore
+
+# Maybe needed as build deps
+using cupynumeric_jll: cupynumeric_jll
+using OpenBLAS32_jll: OpenBLAS32_jll
 
 include("version.jl")
 
@@ -62,17 +67,22 @@ end
 function build_jlcxxwrap(repo_root, cupynumeric_root)
     build_libcxxwrap = joinpath(repo_root, "scripts/install_cxxwrap.sh")
     version_path = joinpath(DEPOT_PATH[1], "dev/libcxxwrap_julia_jll/override/LEGATE_INSTALL.txt")
-    if isfile(version_path)
-        version = VersionNumber(strip(read(version_path, String)))
-        @info "libcxxwrap: Found cuNumeric $version"
-        if is_supported_version(version)
-            @info "libcxxwrap: Found supported version built with cuNumeric.jl: $version"
-            return nothing
+    jlcxxwrap_libpath = joinpath(DEPOT_PATH[1], "dev/libcxxwrap_julia_jll/override/lib/libjlcxxwrap_julia.so")
+    if isfile(jlcxxwrap_libpath)
+        if isfile(version_path)
+            version = VersionNumber(strip(read(version_path, String)))
+            @info "libcxxwrap: Found cuNumeric $version"
+            if is_supported_version(version)
+                @info "libcxxwrap: Found supported version built with cuNumeric.jl: $version"
+                return nothing
+            else
+                @info "libcxxwrap: Unsupported version found: $version. Rebuilding..."
+            end
         else
-            @info "libcxxwrap: Unsupported version found: $version. Rebuilding..."
+            @info "libcxxwrap: No version file found. Starting build..."
         end
     else
-        @info "libcxxwrap: No version file found. Starting build..."
+        @info "libcxxwrap: No jlcxxwrap_julia library found. Starting build..."
     end
 
     @info "libcxxwrap: Running build script: $build_libcxxwrap"
@@ -111,13 +121,6 @@ function build_cpp_wrapper(
     run_sh(`bash $bld_command`, "cpp_wrapper")
 end
 
-function _find_jll_artifact_dir(jll)
-    eval(:(using $(jll)))
-    jll_mod = getfield(Main, jll)
-    root = jll_mod.artifact_dir
-    return root
-end
-
 function _start_build()
     pkg_root = up_dir(@__DIR__)
     deps_dir = joinpath(@__DIR__)
@@ -131,6 +134,35 @@ function _start_build()
     return pkg_root
 end
 
+function _cuda_include_dir_from_toolkit(cuda_toolkit_root)
+    cuda_include_dir = joinpath(cuda_toolkit_root, "include")
+
+    isfile(joinpath(cuda_include_dir, "cuda.h")) ||
+        error("Could not find cuda.h at $(joinpath(cuda_include_dir, "cuda.h"))")
+
+    return cuda_include_dir
+end
+
+function _cuda_driver_lib_from_toolkit(cuda_toolkit_root)
+    candidates = [
+        joinpath(cuda_toolkit_root, "lib64", "stubs", "libcuda.so"),
+        joinpath(cuda_toolkit_root, "lib", "stubs", "libcuda.so"),
+        joinpath(cuda_toolkit_root, "lib64", "libcuda.so"),
+        joinpath(cuda_toolkit_root, "lib", "libcuda.so"),
+    ]
+
+    for path in candidates
+        isfile(path) && return path
+    end
+
+    error("""
+    Could not find libcuda.so under CUDA_TOOLKIT_ROOT=$cuda_toolkit_root.
+
+    Tried:
+    $(join(candidates, "\n"))
+    """)
+end
+
 """
     build CxxWrap and cunumeric_jl_wrapper
 """
@@ -147,7 +179,7 @@ function build_deps(pkg_root, cupynumeric_root, blas_root)
     build_jlcxxwrap(pkg_root, cupynumeric_root)
     build_cpp_wrapper(
         pkg_root, cupynumeric_root, up_dir(legate_lib), blas_root,
-        install_lib,
+        install_lib
     ) # $pkg_root/lib/cunumeric_jl_wrapper
 end
 
@@ -161,9 +193,17 @@ function build(::CNPreferences.Conda)
     pkg_root = _start_build()
 
     cupynumeric_root = load_preference(CNPreferences, "cunumeric_conda_env", nothing)
+    cuda_toolkit_root = load_preference(CNPreferences, "CUDA_TOOLKIT_ROOT", nothing)
     if isnothing(cupynumeric_root)
         error("This shouldn't happen. cunumeric_conda_env = nothing?")
     end
+    if isnothing(cuda_toolkit_root)
+        error(
+            "CUDA_TOOLKIT_ROOT must be set by CNPreferences to point to the CUDA linked in your cupynumeric build."
+        )
+    end
+
+    #!TODO SET LocalPreferences.toml to use local CUDA libraries
 
     is_cupynumeric_installed(cupynumeric_root; throw_errors=true)
     build_deps(pkg_root, cupynumeric_root, cupynumeric_root) # blas is same root as cupynumeric
@@ -172,22 +212,40 @@ end
 function build(::CNPreferences.Developer)
     pkg_root = _start_build()
 
-    # can be nothing so this errors if not set
     cupynumeric_root = load_preference(CNPreferences, "cunumeric_path", nothing)
-    blas_lib = load_preference(CNPreferences, "BLAS_LIB", nothing)
+    blas_lib_dir = load_preference(CNPreferences, "BLAS_LIB_DIR", nothing)
+    cuda_toolkit_root = load_preference(CNPreferences, "CUDA_TOOLKIT_ROOT", nothing)
+
+    cuda_header_path = nothing
+    cuda_driver_path = nothing
     if isnothing(cupynumeric_root)
-        # we are using cupynumeric_jll
-        cupynumeric_root = _find_jll_artifact_dir(:cupynumeric_jll)
+        # We are using cupynumeric_jll.
+        cupynumeric_root = cupynumeric_jll.artifact_dir
+        # cuda_header_path = dirname(cupynumeric_jll.cuda_header)
+        # cuda_driver_path = joinpath(CUDACore.CUDA_Driver_jll.artifact_dir, "lib", "libcuda.so")
     else
-        # this means we have a custom path set
+        # User provided a custom cupynumeric install.
         is_cupynumeric_installed(cupynumeric_root; throw_errors=true)
+
+        isnothing(cuda_toolkit_root) &&
+            error(
+                "CUDA_TOOLKIT_ROOT must be set by CNPreferences when not using cupynumeric_jll in developer mode"
+            )
+
+        # Local full CUDA toolkit path.
+        # cuda_header_path = _cuda_include_dir_from_toolkit(cuda_toolkit_root)
+        # cuda_driver_path = _cuda_driver_lib_from_toolkit(cuda_toolkit_root)
     end
 
-    if isnothing(blas_lib)
-        blas_lib = _find_jll_artifact_dir(:OpenBLAS32_jll)
+    if isnothing(blas_lib_dir)
+        blas_lib_dir = joinpath(OpenBLAS32_jll.artifact_dir, "lib")
     end
 
-    build_deps(pkg_root, cupynumeric_root, up_dir(blas_lib))
+    build_deps(
+        pkg_root,
+        cupynumeric_root,
+        blas_lib_dir
+    )
 end
 
 const mode_str = load_preference(CNPreferences, "cunumeric_mode", CNPreferences.MODE_JLL)
diff --git a/ext/CUDAExt/CUDAExt.jl b/ext/CUDAExt/CUDAExt.jl
deleted file mode 100644
index ed9c07c6..00000000
--- a/ext/CUDAExt/CUDAExt.jl
+++ /dev/null
@@ -1,26 +0,0 @@
-module CUDAExt
-
-using Random
-using CUDA
-using Legate: Legate
-using CxxWrap: CxxWrap
-using cuNumeric: cuNumeric
-import cuNumeric:
-    @cuda_task, @launch, NDArray, assert_experimental
-
-const KERNEL_OFFSET = sizeof(CUDA.KernelState)
-
-include("cuda.jl")
-
-function __init__()
-    if CUDA.functional()
-        # in cuda.jl to notify /wrapper/src/cuda.cpp about CUDA.jl kernel state size
-        cuNumeric.register_kernel_state_size(UInt64(KERNEL_OFFSET))
-        # in /wrapper/src/cuda.cpp
-        cuNumeric.register_tasks();
-    else
-        @warn "CUDA.jl is not functional; skipping CUDA kernel registration."
-    end
-end
-
-end # module CUDAExt
diff --git a/lib/cunumeric_jl_wrapper/CMakeLists.txt b/lib/cunumeric_jl_wrapper/CMakeLists.txt
index 2da543ab..87848a77 100644
--- a/lib/cunumeric_jl_wrapper/CMakeLists.txt
+++ b/lib/cunumeric_jl_wrapper/CMakeLists.txt
@@ -95,4 +95,4 @@ if(HAVE_CUDA)
     target_include_directories(${C_INTERFACE_LIB} PRIVATE ${CUDAToolkit_INCLUDE_DIRS})
 endif()
 
-install(TARGETS ${C_INTERFACE_LIB} DESTINATION lib)
+install(TARGETS ${C_INTERFACE_LIB} DESTINATION lib)
\ No newline at end of file
diff --git a/lib/cunumeric_jl_wrapper/include/cuda_macros.h b/lib/cunumeric_jl_wrapper/include/cuda_macros.h
new file mode 100644
index 00000000..de9a9f2f
--- /dev/null
+++ b/lib/cunumeric_jl_wrapper/include/cuda_macros.h
@@ -0,0 +1,74 @@
+#pragma once
+
+#define ERROR_CHECK(x)                                                 \
+  {                                                                    \
+    cudaError_t status = x;                                            \
+    if (status != cudaSuccess) {                                       \
+      fprintf(stderr, "CUDA Error at %s:%d: %s\n", __FILE__, __LINE__, \
+              cudaGetErrorString(status));                             \
+      if (stream_) cudaStreamDestroy(stream_);                         \
+      exit(-1);                                                        \
+    }                                                                  \
+  }
+
+#define DRIVER_ERROR_CHECK(x)                                                 \
+  {                                                                           \
+    CUresult status = x;                                                      \
+    if (status != CUDA_SUCCESS) {                                             \
+      const char *err_str = nullptr;                                          \
+      cuGetErrorString(status, &err_str);                                     \
+      fprintf(stderr, "CUDA Driver Error at %s:%d: %s\n", __FILE__, __LINE__, \
+              err_str);                                                       \
+      if (stream_) cudaStreamDestroy(stream_);                                \
+      exit(-1);                                                               \
+    }                                                                         \
+  }
+
+#define TEST_PRINT_DEBUG(dev_ptr, N, T, format, stream, message)            \
+  {                                                                         \
+    std::vector<T> host_arr(N);                                             \
+    ERROR_CHECK(cudaMemcpy(host_arr.data(),                                 \
+                           reinterpret_cast<const T *>(dev_ptr),            \
+                           sizeof(T) * N, cudaMemcpyDeviceToHost));         \
+    ERROR_CHECK(cudaStreamSynchronize(stream));                             \
+    fprintf(stderr, "[TEST_PRINT] %s: " format "\n", message, host_arr[0]); \
+  }
+
+#ifdef CUDA_DEBUG
+#define CUDA_DEBUG_PRINT(x) \
+  do {                      \
+    x;                      \
+  } while (0)
+#else
+#define CUDA_DEBUG_PRINT(x) \
+  do {                      \
+  } while (0)
+#endif
+
+#define CUDA_DEVICE_ARRAY_ARG(MODE, ACCESSOR_CALL)                             \
+  template <                                                                   \
+      typename T, int D,                                                       \
+      typename std::enable_if<(D >= 1 && D <= REALM_MAX_DIM), int>::type = 0>  \
+  void cuda_device_array_arg_##MODE(char *&p,                                  \
+                                    const legate::PhysicalArray &rf) {         \
+    auto shp = rf.shape<D>();                                                  \
+    auto acc = rf.data().ACCESSOR_CALL<T, D>();                                \
+    CUDA_DEBUG_PRINT(std::cerr << "[RunPTXTask] " #MODE " accessor shape: "    \
+                               << shp.lo << " - " << shp.hi << ", dim: " << D  \
+                               << std::endl;                                   \
+                     std::cerr << "[RunPTXTask] " #MODE " accessor strides: "  \
+                               << acc.accessor.strides << std::endl;);         \
+    void *dev_ptr = const_cast<void *>(/*.lo to ensure multiple GPU support*/  \
+                                       static_cast<const void *>(              \
+                                           acc.ptr(Realm::Point<D>(shp.lo)))); \
+    auto extents = shp.hi - shp.lo + legate::Point<D>::ONES();                 \
+    CuDeviceArray<D> desc;                                                     \
+    desc.ptr = dev_ptr;                                                        \
+    desc.maxsize = shp.volume() * sizeof(T);                                   \
+    for (size_t i = 0; i < D; ++i) {                                           \
+      desc.dims[i] = extents[i];                                               \
+    }                                                                          \
+    desc.length = shp.volume();                                                \
+    memcpy(p, &desc, sizeof(CuDeviceArray<D>));                                \
+    p += sizeof(CuDeviceArray<D>);                                             \
+  }
diff --git a/lib/cunumeric_jl_wrapper/include/ufi.h b/lib/cunumeric_jl_wrapper/include/ufi.h
index f5558c16..6a2f04f1 100644
--- a/lib/cunumeric_jl_wrapper/include/ufi.h
+++ b/lib/cunumeric_jl_wrapper/include/ufi.h
@@ -28,6 +28,7 @@ namespace ufi {
 enum TaskIDs {
   LOAD_PTX_TASK = 143432,
   RUN_PTX_TASK = 143433,
+  RUN_PTX_BROADCAST_TASK = 143434,
 };
 
 class LoadPTXTask : public legate::LegateTask<LoadPTXTask> {
@@ -46,6 +47,14 @@ class RunPTXTask : public legate::LegateTask<RunPTXTask> {
   static void gpu_variant(legate::TaskContext context);
 };
 
+class RunPTXBroadcastTask : public legate::LegateTask<RunPTXBroadcastTask> {
+ public:
+  static inline const auto TASK_CONFIG =
+      legate::TaskConfig{legate::LocalTaskID{ufi::RUN_PTX_BROADCAST_TASK}};
+
+  static void gpu_variant(legate::TaskContext context);
+};
+
 }  // namespace ufi
 void wrap_cuda_methods(jlcxx::Module& mod);
 
diff --git a/lib/cunumeric_jl_wrapper/src/cuda.cpp b/lib/cunumeric_jl_wrapper/src/cuda.cpp
index d5e779d0..3582b8a6 100644
--- a/lib/cunumeric_jl_wrapper/src/cuda.cpp
+++ b/lib/cunumeric_jl_wrapper/src/cuda.cpp
@@ -29,7 +29,7 @@
 #include "types.h"
 #include "ufi.h"
 
-#define CUDA_DEBUG 0
+// #define CUDA_DEBUG
 
 #define BLOCK_START 1
 #define THREAD_START 4
@@ -175,37 +175,40 @@ struct ufiFunctor {
   }
 };
 
-// https://github.com/nv-legate/legate.pandas/blob/branch-22.01/src/udf/eval_udf_gpu.cc
-/*static*/ void RunPTXTask::gpu_variant(legate::TaskContext context) {
-  cudaStream_t stream_ = context.get_task_stream();
-  std::string kernel_name = context.scalar(0).value<std::string>();  // 0
-
-  std::uint32_t bx =
-      context.scalar(BLOCK_START + 0).value<std::uint32_t>();  // 1
-  std::uint32_t by =
-      context.scalar(BLOCK_START + 1).value<std::uint32_t>();  // 2
-  std::uint32_t bz =
-      context.scalar(BLOCK_START + 2).value<std::uint32_t>();  // 3
-
-  std::uint32_t tx =
-      context.scalar(THREAD_START + 0).value<std::uint32_t>();  // 4
-  std::uint32_t ty =
-      context.scalar(THREAD_START + 1).value<std::uint32_t>();  // 5
-  std::uint32_t tz =
-      context.scalar(THREAD_START + 2).value<std::uint32_t>();  // 6
+struct PTXLaunchParams {
+  cudaStream_t stream;
+  CUstream custream;
+  CUfunction func;
+  std::string kernel_name;
+  std::uint32_t bx, by, bz;
+  std::uint32_t tx, ty, tz;
+};
+
+// Reads common scalars (kernel_name, blocks, threads) and looks up the
+// compiled CUfunction. Shared by RunPTXTask and RunPTXBroadcastTask.
+static PTXLaunchParams read_launch_params(legate::TaskContext &context) {
+  PTXLaunchParams p;
+  p.stream = context.get_task_stream();
+  p.kernel_name = context.scalar(0).value<std::string>();
+
+  p.bx = context.scalar(BLOCK_START + 0).value<std::uint32_t>();
+  p.by = context.scalar(BLOCK_START + 1).value<std::uint32_t>();
+  p.bz = context.scalar(BLOCK_START + 2).value<std::uint32_t>();
+
+  p.tx = context.scalar(THREAD_START + 0).value<std::uint32_t>();
+  p.ty = context.scalar(THREAD_START + 1).value<std::uint32_t>();
+  p.tz = context.scalar(THREAD_START + 2).value<std::uint32_t>();
 
   CUcontext ctx;
-  cuStreamGetCtx(stream_, &ctx);
+  cuStreamGetCtx(p.stream, &ctx);
 
-  FunctionKey key = {ctx, kernel_name};
+  FunctionKey key = {ctx, p.kernel_name};
   assert(cufunction_ptr.has_value());
   FunctionMap &fmap = cufunction_ptr.get();
-
   auto it = fmap.find(key);
 
 #ifdef CUDA_DEBUG
   if (it == fmap.end()) {
-    // for DEBUG output
     std::cerr << "[RunPTXTask] Could not find key: " << key_to_string(key)
               << std::endl;
     for (const auto &[k, v] : fmap) {
@@ -216,17 +219,51 @@ struct ufiFunctor {
 #endif
 
   assert(it != fmap.end());
-  CUfunction func = it->second;
+  p.func = it->second;
+  p.custream = reinterpret_cast<CUstream>(p.stream);
+  return p;
+}
+
+// Launch the kernel with the filled arg_buffer.
+static void launch_kernel(const PTXLaunchParams &lp,
+                          std::vector<char> &arg_buffer,
+                          std::size_t buffer_size) {
+  cudaStream_t stream_ = lp.stream;  // alias for DRIVER_ERROR_CHECK macro
+  void *config[] = {
+      CU_LAUNCH_PARAM_BUFFER_POINTER,
+      static_cast<void *>(arg_buffer.data()),
+      CU_LAUNCH_PARAM_BUFFER_SIZE,
+      &buffer_size,
+      CU_LAUNCH_PARAM_END,
+  };
+
+#ifdef CUDA_DEBUG
+  std::cerr << "[RunPTXTask] Launching kernel " << lp.kernel_name
+            << " with blocks (" << lp.bx << "," << lp.by << "," << lp.bz
+            << ") and threads (" << lp.tx << "," << lp.ty << "," << lp.tz << ")"
+            << std::endl;
+#endif
+
+  DRIVER_ERROR_CHECK(cuLaunchKernel(lp.func, lp.bx, lp.by, lp.bz, lp.tx, lp.ty,
+                                    lp.tz, 0, lp.custream, nullptr, config));
+}
+
+// Helper: align pointer to 8-byte boundary.
+static inline void align8(char *&ptr) {
+  std::uintptr_t addr = reinterpret_cast<std::uintptr_t>(ptr);
+  ptr = reinterpret_cast<char *>((addr + 7) & ~std::uintptr_t(7));
+}
+
+// RunPTXTask: user-defined @cuda_task kernels
+// Arg buffer: [kernel_state | inputs... | outputs... | scalars...]
+// https://github.com/nv-legate/legate.pandas/blob/branch-22.01/src/udf/eval_udf_gpu.cc
+/*static*/ void RunPTXTask::gpu_variant(legate::TaskContext context) {
+  auto lp = read_launch_params(context);
 
   const std::size_t num_inputs = context.num_inputs();
   const std::size_t num_outputs = context.num_outputs();
   const std::size_t num_scalars = context.num_scalars();
-  const std::size_t num_reductions =
-      context.num_reductions();  // unused for now
 
-  // compute total size: all device arrays + all scalars
-  // skip scalar 0-2 (kernel_name, threads, blocks)
-  // we allocate extra to decrease looping and dynamic dispatching on dim
   std::size_t max_buffer_size =
       padded_bytes_kernel_state +
       (num_inputs + num_outputs) * sizeof(CuDeviceArray<REALM_MAX_DIM>);
@@ -238,26 +275,13 @@ struct ufiFunctor {
 
   for (std::size_t i = 0; i < num_inputs; ++i) {
     auto ps = context.input(i);
-    auto code = ps.type().code();
-    auto dim = ps.dim();
-#ifdef CUDA_DEBUG
-    std::cerr << "[RunPTXTask] Input " << i << " type: " << code
-              << ", dim: " << dim << std::endl;
-#endif
-    // dispatch on dim and code with ufiFunctor operator()
-    legate::double_dispatch(dim, code, ufiFunctor{}, ufi::AccessMode::READ, p,
-                            ps);
+    legate::double_dispatch(ps.dim(), ps.type().code(), ufiFunctor{},
+                            ufi::AccessMode::READ, p, ps);
   }
   for (std::size_t i = 0; i < num_outputs; ++i) {
     auto ps = context.output(i);
-    auto code = ps.type().code();
-    auto dim = ps.dim();
-#ifdef CUDA_DEBUG
-    std::cerr << "[RunPTXTask] Output " << i << " type: " << code
-              << ", dim: " << dim << std::endl;
-#endif
-    legate::double_dispatch(dim, code, ufiFunctor{}, ufi::AccessMode::WRITE, p,
-                            ps);
+    legate::double_dispatch(ps.dim(), ps.type().code(), ufiFunctor{},
+                            ufi::AccessMode::WRITE, p, ps);
   }
   for (std::size_t i = ARG_OFFSET; i < num_scalars; ++i) {
     const auto &scalar = context.scalar(i);
@@ -265,29 +289,78 @@ struct ufiFunctor {
     p += scalar.size();
   }
 
-  std::size_t buffer_size = p - arg_buffer.data();  // calc used buffer
+  launch_kernel(lp, arg_buffer, p - arg_buffer.data());
+}
 
-  void *config[] = {
-      CU_LAUNCH_PARAM_BUFFER_POINTER,
-      static_cast<void *>(arg_buffer.data()),
-      CU_LAUNCH_PARAM_BUFFER_SIZE,
-      &buffer_size,
-      CU_LAUNCH_PARAM_END,
-  };
+// RunPTXBroadcastTask: broadcast fusion kernels
+// Arg buffer: [kernel_state | ctx | arg_map-driven args...]
+//
+// Scalars after ARG_OFFSET:
+//   [7]      = ctx (CompilerMetadata, raw bytes)
+//   [8]      = num_kernel_args (Int32)
+//   [9..8+N] = arg_map entries (Int32 each)
+//   [9+N..]  = actual scalar values
+//
+// arg_map encoding:
+//   val >= 0, val < num_outputs  → output[val] (write CuDeviceArray)
+//   val >= num_outputs           → input[val - num_outputs] (read
+//   CuDeviceArray) val < 0                      → scalar at index -(val + 1) in
+//   trailing scalars
+/*static*/ void RunPTXBroadcastTask::gpu_variant(legate::TaskContext context) {
+  auto lp = read_launch_params(context);
 
-  CUstream custream_ = reinterpret_cast<CUstream>(stream_);
+  const std::size_t num_inputs = context.num_inputs();
+  const std::size_t num_outputs = context.num_outputs();
+  const std::size_t num_scalars = context.num_scalars();
+  // Read num_kernel_args first so we can size the buffer precisely
+  std::int32_t num_kernel_args =
+      context.scalar(ARG_OFFSET + 1).value<std::int32_t>();
+  std::size_t map_start = ARG_OFFSET + 2;
+  std::size_t scalar_values_start = map_start + num_kernel_args;
 
-#ifdef CUDA_DEBUG
-  std::cerr << "[RunPTXTask] Launching kernel " << kernel_name
-            << " with blocks (" << bx << "," << by << "," << bz
-            << ") and threads (" << tx << "," << ty << "," << tz << ")"
-            << " on CUcontext " << context_to_string(ctx) << std::endl;
-#endif
-  // Launch the kernel
-  DRIVER_ERROR_CHECK(cuLaunchKernel(func, bx, by, bz, tx, ty, tz, 0, custream_,
-                                    nullptr, config));
+  std::size_t max_buffer_size =
+      padded_bytes_kernel_state +
+      context.scalar(ARG_OFFSET).size() +  // ctx (CompilerMetadata)
+      num_kernel_args *
+          (sizeof(CuDeviceArray<REALM_MAX_DIM>) +
+           8);  // worst case: all args are CuDeviceArrays + alignment
+  for (std::size_t i = scalar_values_start; i < num_scalars; ++i) {
+    max_buffer_size += context.scalar(i).size();
+  }
 
-  // DRIVER_ERROR_CHECK(cuStreamSynchronize(stream_));
+  std::vector<char> arg_buffer(max_buffer_size);
+  char *p = arg_buffer.data() + padded_bytes_kernel_state;
+
+  // 1. Write ctx (CompilerMetadata)
+  if (num_scalars > ARG_OFFSET) {
+    const auto &ctx_scalar = context.scalar(ARG_OFFSET);
+    memcpy(p, ctx_scalar.ptr(), ctx_scalar.size());
+    p += ctx_scalar.size();
+  }
+
+  // 2. Read arg_map and reconstruct arg buffer
+  for (std::int32_t i = 0; i < num_kernel_args; ++i) {
+    std::int32_t val = context.scalar(map_start + i).value<std::int32_t>();
+
+    if (val >= 0 && val < static_cast<std::int32_t>(num_outputs)) {
+      align8(p);
+      auto ps = context.output(val);
+      legate::double_dispatch(ps.dim(), ps.type().code(), ufiFunctor{},
+                              ufi::AccessMode::WRITE, p, ps);
+    } else if (val >= static_cast<std::int32_t>(num_outputs)) {
+      align8(p);
+      auto ps = context.input(val - num_outputs);
+      legate::double_dispatch(ps.dim(), ps.type().code(), ufiFunctor{},
+                              ufi::AccessMode::READ, p, ps);
+    } else {
+      std::size_t scalar_idx = static_cast<std::size_t>(-(val + 1));
+      const auto &scalar = context.scalar(scalar_values_start + scalar_idx);
+      memcpy(p, scalar.ptr(), scalar.size());
+      p += scalar.size();
+    }
+  }
+
+  launch_kernel(lp, arg_buffer, p - arg_buffer.data());
 }
 
 // https://github.com/nv-legate/legate.pandas/blob/branch-22.01/src/udf/load_ptx.cc
@@ -389,6 +462,13 @@ inline void add_xyz_scalars(legate::AutoTask &task,
   task.add_scalar_arg(legate::Scalar(xyz[2]));
 }
 
+inline void add_scalar_from_ptr(legate::AutoTask &task, void *ptr,
+                                size_t size) {
+  uint8_t *byte_ptr = static_cast<uint8_t *>(ptr);
+  std::vector<uint8_t> vec(byte_ptr, byte_ptr + size);
+  task.add_scalar_arg(legate::Scalar(vec));
+}
+
 void gpu_sync() {
   cudaStream_t stream_ = nullptr;
   ERROR_CHECK(cudaDeviceSynchronize());
@@ -412,9 +492,12 @@ void register_kernel_state_size(uint64_t st_size) {
 
 void wrap_cuda_methods(jlcxx::Module &mod) {
   mod.method("add_xyz_scalars", &add_xyz_scalars);
+  mod.method("add_scalar_from_ptr", &add_scalar_from_ptr);
   mod.method("register_kernel_state_size", &register_kernel_state_size);
   mod.method("gpu_sync", &gpu_sync);
   mod.method("extract_kernel_name", &extract_kernel_name);
   mod.set_const("LOAD_PTX", legate::LocalTaskID{ufi::TaskIDs::LOAD_PTX_TASK});
   mod.set_const("RUN_PTX", legate::LocalTaskID{ufi::TaskIDs::RUN_PTX_TASK});
+  mod.set_const("RUN_PTX_BROADCAST",
+                legate::LocalTaskID{ufi::TaskIDs::RUN_PTX_BROADCAST_TASK});
 }
diff --git a/lib/cunumeric_jl_wrapper/src/wrapper.cpp b/lib/cunumeric_jl_wrapper/src/wrapper.cpp
index 29334333..1c5ddb96 100644
--- a/lib/cunumeric_jl_wrapper/src/wrapper.cpp
+++ b/lib/cunumeric_jl_wrapper/src/wrapper.cpp
@@ -61,6 +61,7 @@ void register_tasks() {
   auto library = get_lib();
   ufi::LoadPTXTask::register_variants(library);
   ufi::RunPTXTask::register_variants(library);
+  ufi::RunPTXBroadcastTask::register_variants(library);
 }
 #endif
 
diff --git a/scripts/build_cpp_wrapper.sh b/scripts/build_cpp_wrapper.sh
index aae6f432..6bd7e73d 100755
--- a/scripts/build_cpp_wrapper.sh
+++ b/scripts/build_cpp_wrapper.sh
@@ -60,4 +60,4 @@ else
     echo "Skipping configure (already done in $BUILD_DIR)"
 fi
 
-cmake --build "$BUILD_DIR" --parallel "$NTHREADS" --verbose
+cmake --build "$BUILD_DIR" --parallel "$NTHREADS" --verbose
\ No newline at end of file
diff --git a/scripts/install_cxxwrap.sh b/scripts/install_cxxwrap.sh
index ba430d50..8335858c 100755
--- a/scripts/install_cxxwrap.sh
+++ b/scripts/install_cxxwrap.sh
@@ -69,7 +69,7 @@ JULIA_CXXWRAP=$JULIA_CXXWRAP_DEV/override
 cd $CUNUMERIC_ROOT_DIR
 [ -f Manifest.toml ] && rm Manifest.toml
 rm -rf $JULIA_CXXWRAP_DEV
-julia -e 'using Pkg; Pkg.activate("."); Pkg.add(url="https://github.com/JuliaLegate/Legate.jl")'
+julia -e 'using Pkg; Pkg.activate("."); Pkg.add("Legate")'
 julia -e 'using Pkg; Pkg.activate("."); Pkg.precompile(["CxxWrap"])'
 
 # https://github.com/JuliaInterop/libcxxwrap-julia/tree/v0.13.3?tab=readme-ov-file#preparing-the-install-location
diff --git a/src/cuNumeric.jl b/src/cuNumeric.jl
index 3a9280e6..23389a38 100644
--- a/src/cuNumeric.jl
+++ b/src/cuNumeric.jl
@@ -27,6 +27,12 @@ using Legate
 using Libdl
 using CxxWrap
 
+using CUDATools: CUDATools
+using CUDACore: CUDACore
+import CUDACore: CuArray
+import KernelAbstractions: @kernel, @index
+import KernelAbstractions as KA
+
 using cupynumeric_jll
 using cunumeric_jl_wrapper_jll
 
@@ -138,21 +144,23 @@ include("warnings.jl")
 # NDArray internal
 include("ndarray/detail/ndarray.jl")
 
-# NDArray interface
+# Utilities
+include("cuda/cuda_util.jl")
+include("utilities/version.jl")
+include("util.jl")
+
+const FUSE_BROADCAST_EXPRS = load_preference(CNPreferences, "FUSE_BROADCAST_EXPRS", true)
+
+# Functionality
 include("ndarray/promotion.jl")
+include("cuda/cuda_ptx_task.jl")
+include("ndarray/broadcast_fusion.jl")
 include("ndarray/broadcast.jl")
 include("ndarray/ndarray.jl")
 include("ndarray/unary.jl")
 include("ndarray/binary.jl")
-
-# scoping macro
 include("scoping.jl")
 
-# Utilities
-include("utilities/version.jl")
-include("utilities/cuda_stubs.jl")
-include("util.jl")
-
 # From https://github.com/JuliaGraphics/QML.jl/blob/dca239404135d85fe5d4afe34ed3dc5f61736c63/src/QML.jl#L147
 mutable struct ArgcArgv
     argv
@@ -230,11 +238,14 @@ function __init__()
     _is_precompiling() && return nothing
 
     # Cannot set LEGATE_CONFIG on CI machines used
-    # to register packages. So we will just skip starting 
+    # to register packages. So we will just skip starting
     # legate/cunumeric when using registry CI machines.
     get(ENV, "JULIA_REGISTRYCI_AUTOMERGE", false) == "true" && return nothing
 
     ensure_runtime!()
+
+    # Requries runtime to be started
+    _setup_cuda_tasking()
 end
 
 end #module cuNumeric
diff --git a/ext/CUDAExt/cuda.jl b/src/cuda/cuda_ptx_task.jl
similarity index 54%
rename from ext/CUDAExt/cuda.jl
rename to src/cuda/cuda_ptx_task.jl
index 05d7e9f1..3c2594b8 100644
--- a/ext/CUDAExt/cuda.jl
+++ b/src/cuda/cuda_ptx_task.jl
@@ -1,19 +1,12 @@
-function ndarray_cuda_type(A::NDArray{T,N}) where {T,N}
-    if N == 1
-        CuDeviceVector{T,1}
-    elseif N == 2
-        CuDeviceMatrix{T,1}
-    else
-        CuDeviceArray{T,N,1}
-    end
-end
+export @cuda_task, @launch, CUDATask
 
-function ndarray_cuda_type(arg::T) where {T}
-    Base.isbits(arg) || throw(ArgumentError("Unsupported argument type: $(typeof(arg))"))
-    typeof(arg)
+struct CUDATask
+    func::String
+    argtypes::NTuple{N,Type} where {N} #! THIS IS TYPE UNSTABLE
 end
 
-cuNumeric.map_ndarray_cuda_types(args...) = tuple(map(ndarray_cuda_type, args)...)
+#! JUST PASS TYPES HERE INSTEAD OF CALLING typeof()
+map_ndarray_cuda_types(args...) = tuple(map(ndarray_cuda_type, typeof.(args))...)
 
 function to_stdvec(::Type{T}, vec) where {T}
     stdvec = CxxWrap.StdVector{T}()
@@ -79,8 +72,9 @@ function nda_to_logical_array(arr::NDArray{T,N}) where {T,N}
     return Legate.LogicalArray{T,N}(st_handle[], size(arr))
 end
 
-function Launch(kernel::cuNumeric.CUDATask, inputs::Tuple{Vararg{NDArray}},
-    outputs::Tuple{Vararg{NDArray}}, scalars::Tuple{Vararg{Any}}; blocks, threads)
+function Launch(kernel::CUDATask, inputs::Tuple{Vararg{NDArray}},
+    outputs::Tuple{Vararg{NDArray}}, scalars::Tuple{Vararg{Any}};
+    blocks, threads, taskid=cuNumeric.RUN_PTX, ctx=nothing)
 
     # we find the largest input/output.
     ndarrays = vcat(inputs..., outputs...)
@@ -91,7 +85,6 @@ function Launch(kernel::cuNumeric.CUDATask, inputs::Tuple{Vararg{NDArray}},
 
     rt = Legate.get_runtime()
     lib = cuNumeric.get_lib()
-    taskid = cuNumeric.RUN_PTX
     task = Legate.create_auto_task(rt, lib, taskid)
 
     input_vars = Vector{Legate.Variable}()
@@ -110,14 +103,22 @@ function Launch(kernel::cuNumeric.CUDATask, inputs::Tuple{Vararg{NDArray}},
         push!(output_vars, p)
     end
 
-    # next 3 lines are reserved scalars in the RUN_PTX task
+    # Reserved scalars: kernel_name (0), blocks (1,2,3), threads (4,5,6)
     Legate.add_scalar(task, Legate.string_to_scalar(kernel.func)) # 0
-    cuNumeric.add_xyz_scalars(task, to_stdvec(UInt32, blocks))  # bx,by,bz 1,2,3
-    cuNumeric.add_xyz_scalars(task, to_stdvec(UInt32, threads)) # tx,ty,tz 4,5,6
+    cuNumeric.add_xyz_scalars(task, to_stdvec(UInt32, blocks))    # 1,2,3
+    cuNumeric.add_xyz_scalars(task, to_stdvec(UInt32, threads))   # 4,5,6
+
+    # CompilerMetadata ctx for broadcast tasks (scalar 7, raw bytes)
+    if !isnothing(ctx)
+        ref = Ref(ctx)
+        GC.@preserve ref begin
+            cuNumeric.add_scalar_from_ptr(task, Base.unsafe_convert(Ptr{Cvoid}, ref), sizeof(ctx))
+        end
+    end
 
-    # any user defined scalars in the launch macro
+    # User-defined scalars
     for s in scalars
-        Legate.add_scalar(task, Legate.Scalar(s)) # 7+ -> ARG_OFFSET
+        Legate.add_scalar(task, Legate.Scalar(s))
     end
 
     # all inputs are aligned with all outputs
@@ -125,17 +126,20 @@ function Launch(kernel::cuNumeric.CUDATask, inputs::Tuple{Vararg{NDArray}},
     Legate.submit_auto_task(rt, task)
 end
 
-function cuNumeric.launch(kernel::cuNumeric.CUDATask, inputs, outputs, scalars; blocks, threads)
+function launch(kernel::CUDATask, inputs, outputs, scalars;
+    blocks, threads, taskid=cuNumeric.RUN_PTX, ctx=nothing)
     Launch(kernel,
         isa(inputs, Tuple) ? inputs : (inputs,),
         isa(outputs, Tuple) ? outputs : (outputs,),
         isa(scalars, Tuple) ? scalars : (scalars,);
         blocks=isa(blocks, Tuple) ? blocks : (blocks,),
         threads=isa(threads, Tuple) ? threads : (threads,),
+        taskid=taskid,
+        ctx=ctx,
     )
 end
 
-function cuNumeric.ptx_task(ptx::String, kernel_name)
+function ptx_task(ptx::String, kernel_name)
     rt = Legate.get_runtime()
     lib = cuNumeric.get_lib() # grab lib of legate app
     # this taskid is directly tied to cpp code in our setup
@@ -147,6 +151,42 @@ function cuNumeric.ptx_task(ptx::String, kernel_name)
     Legate.submit_auto_task(rt, task)
 end
 
+"""
+    @cuda_task(f(args...))
+
+Compile a Julia GPU kernel to PTX, register it with the Legate runtime,
+and return a `CUDATask` object for later launch.
+
+# Arguments
+- `f` — The name of the Julia CUDA.jl GPU kernel function to compile.
+- `args...` — Example arguments to the kernel, used to determine the
+  argument type signature when generating PTX.
+
+# Description
+This macro automates the process of:
+1. Inferring the CUDA argument types for the given `args` using
+   `map_ndarray_cuda_types`.
+2. Using `CUDA.code_ptx` to compile the specified GPU kernel
+   (`f`) into raw PTX text for the inferred types.
+3. Extracting the kernel's function symbol name from the PTX using
+   `extract_kernel_name`.
+4. Registering the compiled PTX and kernel name with the Legate runtime
+   via `ptx_task`, making it available for GPU execution.
+5. Returning a `CUDATask` struct that stores the kernel name and type signature,
+   which can be used to configure and launch the kernel later.
+
+# Notes
+- The `args...` are not executed; they are used solely for type inference.
+- This macro is intended for use with the Legate runtime and
+  assumes a CUDA context is available.
+- Make sure your kernel code is GPU-compatible and does not rely on
+  unsupported Julia features.
+
+# Example
+```julia
+mytask = @cuda_task my_kernel(A, B, C)
+```
+"""
 macro cuda_task(call_expr)
     cuNumeric.assert_experimental()
 
@@ -155,21 +195,58 @@ macro cuda_task(call_expr)
 
     esc(quote
         local _buf = IOBuffer()
-        local _types = cuNumeric.map_ndarray_cuda_types($(fargs...))
+        local _types = map_ndarray_cuda_types($(fargs...))
         # generate ptx using CUDA.jl
-        CUDA.code_ptx(_buf, $fname, _types; raw=false, kernel=true)
+        CUDATools.code_ptx(_buf, $fname, _types; raw=false, kernel=true)
 
         local _ptx = String(take!(_buf))
-        local _func_name = cuNumeric.extract_kernel_name(_ptx)
+        local _func_name = extract_kernel_name(_ptx)
 
         # issue ptx_task within legate runtime to register cufunction ptr with cucontext
-        cuNumeric.ptx_task(_ptx, _func_name)
+        ptx_task(_ptx, _func_name)
 
         # create a cuNumeric.CUDAtask that stores some info for a launch config
-        cuNumeric.CUDATask(_func_name, _types)
+        CUDATask(_func_name, _types)
     end)
 end
 
+"""
+    @launch(; task, blocks=(1,), threads=(256,), inputs=(), outputs=(), scalars=())
+
+Launch a GPU kernel (previously registered via [`@cuda_task`](@ref))  through the Legate runtime.
+
+# Keywords
+- `task` — A `CUDATask` object, typically returned by [`@cuda_task`](@ref).
+- `blocks`  — Tuple or single element specifying the CUDA grid dimensions. Defaults to `(1,)`.
+- `threads` — Tuple or single element specifying the CUDA block dimensions. Defaults to `(256,)`.
+- `inputs`  — Tuple or single element of input NDArray objects.
+- `outputs` — Tuple or single element of output NDArray objects.
+- `scalars` — Tuple or single element of scalar values.
+
+# Description
+The `@launch` macro validates the provided keywords, ensuring only
+the allowed set (`:task`, `:blocks`, `:threads`, `:inputs`, `:outputs`, `:scalars`)
+are present. It then expands to a call to `cuNumeric.launch`,
+passing the given arguments to the Legate runtime for execution.
+
+This macro is meant to provide a concise, declarative syntax for
+launching GPU kernels, separating kernel compilation (via `@cuda_task`)
+from execution configuration.
+
+# Notes
+- `task` **must** be a kernel registered with the runtime, usually from `@cuda_task`.
+- All keyword arguments must be specified as assignments, e.g. `blocks=(2,2)` not positional arguments.
+- Defaults are chosen for single-block, 256-thread 1D launches.
+- The macro escapes its body so that the values of inputs/outputs/scalars are captured
+  from the surrounding scope at macro expansion time.
+
+# Example
+```julia
+mytask = @cuda_task my_kernel(A, B, C)
+
+@launch task=mytask blocks=(8,8) threads=(32,32) inputs=(A, B) outputs=(C)
+```
+"""
 macro launch(args...)
     cuNumeric.assert_experimental()
 
@@ -203,7 +280,8 @@ macro launch(args...)
     esc(
         quote
             cuNumeric.launch(
-                $task, $inputs, $outputs, $scalars; blocks=($blocks), threads=($threads)
+                $task, $inputs, $outputs, $scalars;
+                blocks=($blocks), threads=($threads),
             )
         end,
     )
diff --git a/src/cuda/cuda_util.jl b/src/cuda/cuda_util.jl
new file mode 100644
index 00000000..70760bcb
--- /dev/null
+++ b/src/cuda/cuda_util.jl
@@ -0,0 +1,45 @@
+const KERNEL_OFFSET = sizeof(CUDACore.KernelState)
+
+function _setup_cuda_tasking()
+    if CUDACore.functional()
+        # in cuda.jl to notify /wrapper/src/cuda.cpp about CUDA.jl kernel state size
+        register_kernel_state_size(UInt64(KERNEL_OFFSET))
+        # in /wrapper/src/cuda.cpp
+        register_tasks()
+    else
+        @warn "CUDA.jl is not functional; skipping CUDA kernel registration."
+    end
+end
+
+# Other memeory types here: https://github.com/JuliaGPU/CUDA.jl/blob/345c1600ebd561135148bb04ee2657f521a40e25/CUDACore/src/device/pointer.jl#L7
+function ndarray_cuda_type(::Type{<:NDArray{T,N}}) where {T,N}
+    CUDACore.CuDeviceArray{T,N,CUDACore.AS.Global}
+end
+
+function ndarray_cuda_type(::Type{T}) where {T}
+    Base.isbitstype(T) || throw(ArgumentError("Unsupported argument type: $(T)"))
+    return T
+end
+
+"""
+    map_cuda_type(::Type{T})::Type
+
+Recursively rewrite cuNumeric broadcast-related types so they can be treated as CUDA-friendly
+types for code generation (e.g. mapping `NDArray{...}` to `CuDeviceArray{...}` inside
+`Base.Broadcast.Broadcasted{...}` type parameters).
+"""
+map_cuda_type(::Type{T}) where {T} = T
+
+map_cuda_type(::Type{<:NDArray{T,N}}) where {T,N} = CUDACore.CuDeviceArray{T,N,CUDACore.AS.Global}
+
+function map_cuda_type(::Type{T}) where {T<:Tuple}
+    return Tuple{map_cuda_type.(T.parameters)...}
+end
+
+function map_cuda_type(::Type{Base.Broadcast.Broadcasted{S,Ax,F,Args}}) where {S,Ax,F,Args}
+    return Base.Broadcast.Broadcasted{map_cuda_type(S),Ax,F,map_cuda_type(Args)}
+end
+
+function map_cuda_type(::Type{Base.Broadcast.Extruded{X,K,D}}) where {X,K,D}
+    return Base.Broadcast.Extruded{map_cuda_type(X),K,D}
+end
diff --git a/src/ndarray/broadcast.jl b/src/ndarray/broadcast.jl
index fcf0a3e5..1f0eceeb 100644
--- a/src/ndarray/broadcast.jl
+++ b/src/ndarray/broadcast.jl
@@ -4,6 +4,11 @@ struct NDArrayStyle{N} <: AbstractArrayStyle{N} end
 Base.BroadcastStyle(::Type{<:NDArray{<:Any,N}}) where {N} = NDArrayStyle{N}()
 Base.BroadcastStyle(::NDArrayStyle{N}, ::NDArrayStyle{M}) where {N,M} = NDArrayStyle{max(N, M)}()
 
+# Some other functions in cuda_util.jl
+function map_cuda_type(::Type{cuNumeric.NDArrayStyle{N}}) where {N}
+    CUDACore.CuArrayStyle{N,CUDACore.DeviceMemory}
+end # Also can be HostMemory or UnifiedMemory
+
 function _nd_forbid_mix()
     throw(
         ArgumentError(
@@ -38,6 +43,7 @@ function Base.similar(bc::Broadcasted{NDArrayStyle{N}}, ::Type{ElType}) where {N
 end
 
 function __broadcast(f::Function, _, args...)
+    #! WITH FUSION I THINK WE CAN SUPPORT THIS BY JUST CALLING MAP or MAP!
     error(
         """
         Tried to broadcast $(f). cuNumeric.jl does not support broadcasting user-defined functions yet. Please re-define \
@@ -52,18 +58,27 @@ end
 bcast_depth(bc::Base.Broadcast.Broadcasted) = maximum(bcast_depth, bc.args, init=0) + 1;
 bcast_depth(::Any) = 0
 
-function Base.Broadcast.materialize(bc::Broadcasted{<:NDArrayStyle})
+struct BrokenBroadcast{T} end
+Base.convert(::Type{BrokenBroadcast{T}}, x) where {T} = BrokenBroadcast{T}()
+Base.convert(::Type{BrokenBroadcast{T}}, x::BrokenBroadcast{T}) where {T} = x
+Base.eltype(::Type{BrokenBroadcast{T}}) where {T} = T
+
+function Broadcast.copy(bc::Broadcasted{<:NDArrayStyle{0}})
     ElType = Broadcast.combine_eltypes(bc.f, bc.args)
-    if ElType == Union{} || !Base.allocatedinline(ElType)
-        error("Cannot broadcast $(bc.f) over NDArrays with eltypes: $(eltype.(bc.args))")
+    if ElType == Union{}
+        ElType = Nothing
     end
+    dest = copyto!(similar(bc, ElType), bc)
+    #! CHECK THIS DOESNT CAUSE ISSUES DUE TO BLOCKING NATURE
+    return @allowscalar dest[CartesianIndex()]
+end
 
-    #* This be the place to inject kernel fusion via CUDA.jl
-    #* Use the function in Base.Broadcast.flatten(bc).
-    #* How can we check all the funcs in this expr
-    #* are supported by CUDA?
-
-    return unravel_broadcast_tree(bc)
+@inline function Broadcast.copy(bc::Broadcasted{<:NDArrayStyle})
+    ElType = Broadcast.combine_eltypes(bc.f, bc.args)
+    if ElType == Union{} || !Base.allocatedinline(ElType)
+        ElType = BrokenBroadcast{ElType}
+    end
+    copyto!(similar(bc, ElType), bc)
 end
 
 # Recursion base cases
@@ -84,6 +99,7 @@ function __materialize(bc::Broadcasted{<:NDArrayStyle})
     unravel_broadcast_tree(bc)
 end
 
+# Un-fused implementation of broadcast tree
 function unravel_broadcast_tree(bc::Broadcasted)
 
     # Recursively materialize/unravel any nested broadcasts
@@ -109,10 +125,32 @@ function unravel_broadcast_tree(bc::Broadcasted)
     return __broadcast(bc.f, out, in_args...)
 end
 
-# Support .=
-function Base.copyto!(dest::NDArray{T,N}, bc::Broadcasted{<:NDArrayStyle{N}}) where {T,N}
-    # Moves result from broadcast (src) to dest. src array is no longer valid
-    #! THIS ENABLES FOOT GUN IF USER SPECIFIES INTEGER ARRAY AT OUTPUT
-    nda_move(dest, checked_promote_arr(Base.Broadcast.materialize(bc), T))
-    return dest
+@inline function _copyto!(dest::NDArray, bc::Broadcasted)
+    axes(dest) == axes(bc) || Broadcast.throwdm(axes(dest), axes(bc))
+    isempty(dest) && return dest
+    if eltype(dest) <: BrokenBroadcast
+        throw(
+            ArgumentError(
+                "Broadcast operation resulting in $(eltype(eltype(dest))) is not NDArray compatible"
+            ),
+        )
+    end
+
+    #! IF THIS IS KNOWN AT COMPILE TIME WE CAN GENERATE
+    #! THIS CODE WITHOUT THE IF STATEMENT
+    if FUSE_BROADCAST_EXPRS
+        #! DO I NEED TO DO TYPE PROMOTION CHECKS BEFORE RETURNING?
+        #! WE MIGHT NEED TO CHECK IF ON GPU OR CPU AND FALLBACK
+        return fuse_broadcast_tree!(dest, bc)
+    else
+        temp_result = unravel_broadcast_tree(bc)
+        nda_move(dest, checked_promote_arr(temp_result, eltype(dest)))
+        return dest
+    end
 end
+
+# Support .=
+@inline Base.copyto!(dest::NDArray, bc::Broadcasted{Nothing}) = _copyto!(dest, bc)
+@inline Base.copyto!(dest::NDArray, bc::Broadcasted{<:NDArrayStyle}) = _copyto!(dest, bc)
+
+#! TODO ADD MAP FUSED IMPLEMENTATIONS
diff --git a/src/ndarray/broadcast_fusion.jl b/src/ndarray/broadcast_fusion.jl
new file mode 100644
index 00000000..6961d6a9
--- /dev/null
+++ b/src/ndarray/broadcast_fusion.jl
@@ -0,0 +1,372 @@
+
+struct RuntimeBroadcastArg{J} end
+struct StaticBroadcastArg{J} end
+
+Base.@propagate_inbounds @inline function _gpu_broadcast_getindex(x, I)
+    return @inbounds Base.Broadcast._broadcast_getindex(x, I)
+end
+
+Base.@propagate_inbounds @inline function _gpu_broadcast_getindex(x::Number, I)
+    return x
+end
+
+Base.@propagate_inbounds @inline function _materialize_broadcast_arg(
+    ::RuntimeBroadcastArg{J},
+    runtime_args,
+    static_args,
+    I,
+) where {J}
+    arg = getfield(runtime_args, J)
+    return @inbounds _gpu_broadcast_getindex(arg, I)
+end
+
+Base.@propagate_inbounds @inline function _materialize_broadcast_arg(
+    ::StaticBroadcastArg{J},
+    runtime_args,
+    static_args,
+    I,
+) where {J}
+    return getfield(static_args, J)
+end
+
+Base.@propagate_inbounds @inline function _materialize_broadcast_args(
+    ::Tuple{},
+    runtime_args,
+    static_args,
+    I,
+)
+    return ()
+end
+
+Base.@propagate_inbounds @inline function _materialize_broadcast_args(
+    plan::Tuple,
+    runtime_args,
+    static_args,
+    I,
+)
+    head = getfield(plan, 1)
+    tail = Base.tail(plan)
+
+    return (
+        @inbounds(_materialize_broadcast_arg(head, runtime_args, static_args, I)),
+        @inbounds(_materialize_broadcast_args(tail, runtime_args, static_args, I))...,
+    )
+end
+
+_is_runtime_broadcast_arg(x::NDArray) = true
+_is_runtime_broadcast_arg(x::Base.Broadcast.Extruded) = true
+_is_runtime_broadcast_arg(x::Number) = true
+_is_runtime_broadcast_arg(x) = false
+
+function _push_runtime_arg!(runtime_args, arg_plan, x)
+    push!(runtime_args, x)
+    push!(arg_plan, RuntimeBroadcastArg{length(runtime_args)}())
+    return nothing
+end
+
+function _push_static_arg!(static_args, arg_plan, x)
+    isbits(x) || throw(
+        ArgumentError(
+            "Broadcast fusion cannot statically capture non-isbits broadcast leaf " *
+            "$(repr(x)) of type $(typeof(x))",
+        ),
+    )
+
+    push!(static_args, x)
+    push!(arg_plan, StaticBroadcastArg{length(static_args)}())
+    return nothing
+end
+
+function split_broadcast_args_for_kernel(args::Tuple)
+    runtime_args = Any[]
+    static_args = Any[]
+    arg_plan = Any[]
+
+    for arg in args
+        if arg isa Base.RefValue
+            value = arg[]
+
+            # Keep ordinary numeric broadcast scalars dynamic so scalar values
+            # do not cause recompilation.
+            if value isa Number
+                _push_runtime_arg!(runtime_args, arg_plan, value)
+            else
+                # Function singletons, Val{N}(), etc.
+                _push_static_arg!(static_args, arg_plan, value)
+            end
+
+        elseif _is_runtime_broadcast_arg(arg)
+            _push_runtime_arg!(runtime_args, arg_plan, arg)
+
+        elseif isbits(arg)
+            # Conservative fallback for singleton/isbits scalar broadcast leaves.
+            _push_static_arg!(static_args, arg_plan, arg)
+
+        else
+            throw(
+                ArgumentError(
+                    "Broadcast fusion does not know how to lower argument " *
+                    "$(repr(arg)) of type $(typeof(arg))",
+                ),
+            )
+        end
+    end
+
+    return tuple(runtime_args...), tuple(static_args...), tuple(arg_plan...)
+end
+
+##############
+
+# The indexing in these kernels is based off this internal function from julia/broadcast.jl
+
+# Base.@propagate_inbounds function _broadcast_getindex(bc::Broadcasted{<:Any,<:Any,<:Any,<:Any}, I)
+#     args = _getindex(bc.args, I)
+#     return _broadcast_getindex_evalf(bc.f, args...)
+# end
+
+function make_linear_kernel(dest, bc::Base.Broadcast.Broadcasted, arg_plan, static_args)
+    f = bc.f
+
+    @kernel function broadcast_kernel_linear_splat(dest, runtime_args...)
+        I = @index(Global, Linear)
+        @inbounds args_modified = _materialize_broadcast_args(
+            arg_plan, runtime_args, static_args, I
+        )
+        @inbounds dest[I] = Base.Broadcast._broadcast_getindex_evalf(f, args_modified...)
+    end
+
+    return broadcast_kernel_linear_splat
+end
+
+function make_cartesian_kernel(dest, bc::Base.Broadcast.Broadcasted, arg_plan, static_args)
+    f = bc.f
+
+    @kernel function broadcast_kernel_cartesian_splat(dest, runtime_args...)
+        I = @index(Global, Cartesian)
+        @inbounds args_modified = _materialize_broadcast_args(
+            arg_plan, runtime_args, static_args, I
+        )
+        @inbounds dest[I] = Base.Broadcast._broadcast_getindex_evalf(f, args_modified...)
+    end
+
+    return broadcast_kernel_cartesian_splat
+end
+
+struct FusedBroadcastMetadata
+    ctx::Any # KA.CompilerMetadata
+    threads::Int
+    blocks::Int
+    cuda_task::CUDATask
+end
+
+const _BCAST_PTX_CACHE = Dict{Tuple{Any,DataType,DataType,Any},FusedBroadcastMetadata}()
+const _BCAST_PTX_CACHE_LOCK = ReentrantLock()
+
+cudevice_array_offset(::Type{T}) where {T<:CUDACore.CuDeviceArray} = 0
+cudevice_array_offset(::Type{T}) where {T<:Base.Broadcast.Extruded} = Int(fieldoffset(T, 1))
+
+stores_cudevicearray(::Type{T}) where {T<:CUDACore.CuDeviceArray} = true
+stores_cudevicearray(::Type{T}) where {T<:Base.Broadcast.Extruded} = true
+stores_cudevicearray(::Type{T}) where {T<:Number} = false
+stores_cudevicearray(::Type{T}) where {T<:Bool} = false
+
+function stores_cudevicearray(::Type{T}) where {T}
+    throw(error("Broadcast fusion. Don't know what to do with type: $T"))
+end
+
+function find_cudevicearray_offsets_and_indices(::Type{BC_ARGS}) where {BC_ARGS}
+    offsets = Vector{Int}()
+    indices = Vector{Int}()
+    for (i, T) in enumerate(fieldtypes(BC_ARGS))
+        if stores_cudevicearray(T)
+            offset = fieldoffset(BC_ARGS, i) + cudevice_array_offset(T)
+            push!(offsets, offset)
+            push!(indices, i)
+        end
+    end
+    return tuple(offsets...), tuple(indices...)
+end
+
+function find_scalar_offsets_and_indices(::Type{BC_ARGS}) where {BC_ARGS}
+    offsets = Vector{Int}()
+    indices = Vector{Int}()
+    for (i, T) in enumerate(fieldtypes(BC_ARGS))
+        if (T <: Number)
+            offset = fieldoffset(BC_ARGS, i)
+            push!(offsets, offset)
+            push!(indices, i)
+        end
+    end
+    return tuple(offsets...), tuple(indices...)
+end
+
+get_ndarray(x::T) where {T<:NDArray} = x
+get_ndarray(x::T) where {T<:Base.Broadcast.Extruded} = x.x
+get_ndarray(x::T) where {T} = throw(error("Broadcast fusion. Don't know what to do with type: $T"))
+
+function _threads_from_occupancy(
+    obj::KA.Kernel{CUDACore.CUDAKernels.CUDABackend},
+    ::Type{DEST_T},
+    ARG_TYPES...;
+    ndrange,
+) where {DEST_T}
+    backend = KA.backend(obj)
+
+    ndrange, workgroupsize, iterspace, dynamic = KA.launch_config(obj, ndrange, nothing)
+    ctx = KA.mkcontext(obj, ndrange, iterspace)
+
+    # If the kernel is statically sized we can tell the compiler about that
+    maxthreads =
+        if KA.workgroupsize(obj) <: KA.StaticSize
+            prod(KA.get(KA.workgroupsize(obj)))
+        else
+            nothing
+        end
+
+    # Determine threads via occupancy
+    tt = Base.to_tuple_type((typeof(ctx), DEST_T, ARG_TYPES...))
+    host_kernel = CUDACore.cufunction(
+        obj.f,
+        tt;
+        kernel=true,
+        maxthreads=maxthreads,
+        always_inline=backend.always_inline,
+    )
+    config = CUDACore.launch_configuration(host_kernel.fun; max_threads=prod(ndrange))
+    threads = config.threads
+
+    workgroupsize = CUDACore.CUDAKernels.threads_to_workgroupsize(threads, ndrange)
+    iterspace, dynamic = KA.partition(obj, ndrange, workgroupsize)
+    ctx = KA.mkcontext(obj, ndrange, iterspace)
+
+    blocks = length(KA.blocks(iterspace))
+    threads = length(KA.workitems(iterspace))
+
+    return threads, blocks, ctx
+end
+
+"""
+    get_ptx(obj::KA.Kernel{CUDABackend}, ::Type{DEST_T}, ::Type{BC_T};
+                         ndrange) -> (ptx::String, threads::Int, blocks::Int)
+
+Compile a KA CUDA kernel (kernel-body `obj.f(ctx, ...)`) using *types only* for `DEST_T` and `BC_T`,
+choose a workgroup size (threads) using CUDA occupancy when possible, and return the generated PTX.
+"""
+function get_ptx(
+    obj::KA.Kernel{CUDACore.CUDAKernels.CUDABackend},
+    ::Type{DEST_T},
+    arg_types...;
+    ndrange,
+) where {DEST_T,BC_T}
+    # println(Base.isbitstype.(arg_types))
+    threads, blocks, ctx = _threads_from_occupancy(obj, DEST_T, arg_types...; ndrange=ndrange)
+    blocks == 0 && return "", 0, 0, ctx #! MAYBE ERROR HERE?
+
+    # Generate PTX
+    buf = IOBuffer()
+    #!TODO REMOVE THE MANUAL PTX VERSION HERE!
+    CUDATools.code_ptx(buf, obj.f, (typeof(ctx), DEST_T, arg_types...);
+        raw=false, kernel=true, ptx=v"7.8")
+
+    return String(take!(buf)), threads, blocks, ctx
+end
+
+function get_cuda_task(
+    obj::KA.Kernel{CUDACore.CUDAKernels.CUDABackend},
+    dest::D,
+    runtime_args::RT,
+    ndrange,
+) where {D<:NDArray,RT<:Tuple}
+    DEST_T = map_cuda_type(D)
+    ARG_TYPES = map_cuda_type.(typeof.(runtime_args))
+
+    key = (obj, D, RT, ndrange)
+
+    lock(_BCAST_PTX_CACHE_LOCK) do
+        return get!(_BCAST_PTX_CACHE, key) do
+            ptx, threads, blocks, ctx = get_ptx(obj, DEST_T, ARG_TYPES...; ndrange=ndrange)
+
+            orig_name = extract_kernel_name(ptx)
+            unique_name = orig_name * "_" * string(hash(ptx); base=16)
+            ptx = replace(ptx, orig_name => unique_name)
+
+            ptx_task(ptx, unique_name)
+            cuda_task = CUDATask(unique_name, (DEST_T, ARG_TYPES...))
+
+            FusedBroadcastMetadata(ctx, threads, blocks, cuda_task)
+        end
+    end
+end
+
+function fuse_broadcast_tree!(dest::D, bc::B) where {D<:NDArray,B<:Base.Broadcast.Broadcasted}
+    bc = Base.Broadcast.preprocess(dest, bc)
+    bc = Base.Broadcast.instantiate(bc)
+    bc = Base.Broadcast.flatten(bc)
+
+    # Things like exponentiation generate arguments like Base.RefValue
+    # which do not work with our pattern for making CUDA kernels as they are
+    # not is-bits types. We split these out manually into static args and handle
+    # them separately from runtime args (i.e., arrays, scalars)
+    runtime_args, static_args, arg_plan = split_broadcast_args_for_kernel(bc.args)
+
+    broadcast_kernel =
+        if ndims(dest) == 1 ||
+            (isa(IndexStyle(dest), IndexLinear) &&
+            isa(IndexStyle(bc), IndexLinear))
+            make_linear_kernel(dest, bc, arg_plan, static_args)
+        else
+            make_cartesian_kernel(dest, bc, arg_plan, static_args)
+        end
+
+    ndrange = ndims(dest) > 0 ? size(dest) : (1,)
+
+    bck_cuda = broadcast_kernel(CUDACore.CUDAKernels.CUDABackend())
+
+    fkm = get_cuda_task(bck_cuda, dest, runtime_args, ndrange)
+
+    num_outputs = 1
+
+    unique_ndarrays = NDArray[]
+    ndarray_to_input_idx = Dict{UInt,Int}()
+
+    arg_map = Int32[]
+    actual_scalars = Any[]
+
+    # First PTX argument after ctx is dest.
+    push!(arg_map, Int32(0))
+
+    # Now map only runtime args, not original bc.args.
+    for arg in runtime_args
+        if stores_cudevicearray(map_cuda_type(typeof(arg)))
+            nda = get_ndarray(arg)
+            oid = objectid(nda)
+
+            if !haskey(ndarray_to_input_idx, oid)
+                push!(unique_ndarrays, nda)
+                ndarray_to_input_idx[oid] = length(unique_ndarrays) - 1
+            end
+
+            input_idx = ndarray_to_input_idx[oid]
+            push!(arg_map, Int32(num_outputs + input_idx))
+        else
+            push!(arg_map, Int32(-1 - length(actual_scalars)))
+            push!(actual_scalars, arg)
+        end
+    end
+
+    input_ndarrays = tuple(unique_ndarrays...)
+
+    launch(
+        fkm.cuda_task,
+        input_ndarrays,
+        (dest,),
+        (Int32(length(arg_map)), arg_map..., actual_scalars...);
+        blocks=fkm.blocks,
+        threads=fkm.threads,
+        taskid=cuNumeric.RUN_PTX_BROADCAST,
+        ctx=fkm.ctx,
+    )
+
+    #! PROMOTION CHECKS?
+    return dest
+end
diff --git a/src/ndarray/ndarray.jl b/src/ndarray/ndarray.jl
index d58f6572..edfed03f 100644
--- a/src/ndarray/ndarray.jl
+++ b/src/ndarray/ndarray.jl
@@ -214,6 +214,7 @@ size(arr, 2)
 """
 Base.size(arr::NDArray{<:Any,N}) where {N} = cuNumeric.shape(arr)
 Base.size(arr::NDArray, dim::Int) = Base.size(arr)[dim]
+Base.isempty(arr::NDArray) = any(==(0), size(arr))
 
 @doc"""
     Base.firstindex(arr::NDArray, dim::Int)
@@ -240,9 +241,8 @@ Base.view(arr::NDArray, inds...) = arr[inds...] # NDArray slices are views by de
 Base.IndexStyle(::NDArray) = IndexCartesian()
 
 function Base.show(io::IO, arr::NDArray{T,0}) where {T}
-    println(io, "0-dimensional NDArray{$(T),0}")
     allowscalar() do
-        print(io, arr[])
+        print(io, "NDArray{$(T),0}(", repr(arr[]), ")")
     end
 end
 
@@ -253,13 +253,13 @@ function Base.show(io::IO, ::MIME"text/plain", arr::NDArray{T,0}) where {T}
     end
 end
 
-function Base.show(io::IO, arr::NDArray{T,D}) where {T,D}
-    println(io, "NDArray{$(T),$(D)}")
-    Base.print_array(io, Array(arr))
+function Base.show(io::IO, arr::NDArray{T,N}) where {T,N}
+    print(io, "NDArray{$(T),$(N)} with size ", size(arr))
 end
 
-function Base.show(io::IO, ::MIME"text/plain", arr::NDArray{T}) where {T}
-    Base.show(io, arr)
+function Base.show(io::IO, ::MIME"text/plain", arr::NDArray{T,N}) where {T,N}
+    println(io, "NDArray{$(T),$(N)} with size ", size(arr))
+    Base.print_array(io, Array(arr))
 end
 
 function Base.print(arr::NDArray{T}) where {T}
diff --git a/src/utilities/cuda_stubs.jl b/src/utilities/cuda_stubs.jl
deleted file mode 100644
index caead31a..00000000
--- a/src/utilities/cuda_stubs.jl
+++ /dev/null
@@ -1,92 +0,0 @@
-## This file contains stubs for methods implemented in
-## the CUDA package extensions not implemented
-## elsewhere in the package.
-
-export @cuda_task, @launch
-export map_ndarray_cuda_types, launch, CUDATask
-
-function ptx_task end
-function map_ndarray_cuda_types end
-function launch end
-
-struct CUDATask
-    func::String
-    argtypes::NTuple{N,Type} where {N}
-end
-
-"""
-    @cuda_task(f(args...))
-
-Compile a Julia GPU kernel to PTX, register it with the Legate runtime,
-and return a `CUDATask` object for later launch.
-
-# Arguments
-- `f` — The name of the Julia CUDA.jl GPU kernel function to compile.
-- `args...` — Example arguments to the kernel, used to determine the
-  argument type signature when generating PTX.
-
-# Description
-This macro automates the process of:
-1. Inferring the CUDA argument types for the given `args` using
-   `map_ndarray_cuda_types`.
-2. Using `CUDA.code_ptx` to compile the specified GPU kernel
-   (`f`) into raw PTX text for the inferred types.
-3. Extracting the kernel's function symbol name from the PTX using
-   `extract_kernel_name`.
-4. Registering the compiled PTX and kernel name with the Legate runtime
-   via `ptx_task`, making it available for GPU execution.
-5. Returning a `CUDATask` struct that stores the kernel name and type signature,
-   which can be used to configure and launch the kernel later.
-
-# Notes
-- The `args...` are not executed; they are used solely for type inference.
-- This macro is intended for use with the Legate runtime and
-  assumes a CUDA context is available.
-- Make sure your kernel code is GPU-compatible and does not rely on
-  unsupported Julia features.
-
-# Example
-```julia
-mytask = @cuda_task my_kernel(A, B, C)
-```
-"""
-macro cuda_task end
-
-"""
-    @launch(; task, blocks=(1,), threads=(256,), inputs=(), outputs=(), scalars=())
-
-Launch a GPU kernel (previously registered via [`@cuda_task`](@ref))  through the Legate runtime.
-
-# Keywords
-- `task` — A `CUDATask` object, typically returned by [`@cuda_task`](@ref).
-- `blocks`  — Tuple or single element specifying the CUDA grid dimensions. Defaults to `(1,)`.
-- `threads` — Tuple or single element specifying the CUDA block dimensions. Defaults to `(256,)`.
-- `inputs`  — Tuple or single element of input NDArray objects.
-- `outputs` — Tuple or single element of output NDArray objects.
-- `scalars` — Tuple or single element of scalar values.
-
-# Description
-The `@launch` macro validates the provided keywords, ensuring only
-the allowed set (`:task`, `:blocks`, `:threads`, `:inputs`, `:outputs`, `:scalars`)
-are present. It then expands to a call to `cuNumeric.launch`,
-passing the given arguments to the Legate runtime for execution.
-
-This macro is meant to provide a concise, declarative syntax for
-launching GPU kernels, separating kernel compilation (via `@cuda_task`)
-from execution configuration.
-
-# Notes
-- `task` **must** be a kernel registered with the runtime, usually from `@cuda_task`.
-- All keyword arguments must be specified as assignments, e.g. `blocks=(2,2)` not positional arguments.
-- Defaults are chosen for single-block, 256-thread 1D launches.
-- The macro escapes its body so that the values of inputs/outputs/scalars are captured
-  from the surrounding scope at macro expansion time.
-
-# Example
-```julia
-mytask = @cuda_task my_kernel(A, B, C)
-
-@launch task=mytask blocks=(8,8) threads=(32,32) inputs=(A, B) outputs=(C)
-```
-"""
-macro launch end
diff --git a/test/runtests.jl b/test/runtests.jl
index d2396b3f..7aa99150 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,4 +1,4 @@
-#= Copyright 2026 Northwestern University, 
+#= Copyright 2026 Northwestern University,
  *                   Carnegie Mellon University University
  *
  * Licensed under the Apache License, Version 2.0 (the "License");
@@ -57,6 +57,7 @@ include("tests/unary_tests.jl")
 include("tests/binary_tests.jl")
 include("tests/scoping.jl")
 include("tests/scoping-advanced.jl")
+include("tests/broadcast_fusion_tests.jl")
 
 @testset verbose = true "AXPY" begin
     N = 100
@@ -430,11 +431,15 @@ end
 end
 
 if run_gpu_tests
+    @testset verbose = true "Broadcast Fusion" begin
+        test_broadcast_fusion()
+    end
+
     # @testset verbose = true "CUDA Tests" begin
     #     cuda_unaryop(rtol(Float32))
     #     cuda_binaryop(rtol(Float32))
     # end
-    @warn "CUDA tests are turned off inside Pkg.test for now. --check-bounds=yes causes issues."
+    @warn "CUDA @cuda_task tests are turned off inside Pkg.test for now. --check-bounds=yes causes issues."
 else
     @warn "The CUDA tests will not be run as a CUDA-enabled device is not available"
 end
diff --git a/test/tests/broadcast_fusion_tests.jl b/test/tests/broadcast_fusion_tests.jl
new file mode 100644
index 00000000..9c5be8a2
--- /dev/null
+++ b/test/tests/broadcast_fusion_tests.jl
@@ -0,0 +1,185 @@
+#= Copyright 2026 Northwestern University,
+ *                   Carnegie Mellon University University
+ *
+ * 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.
+ *
+ * Author(s): David Krasowska <krasow@u.northwestern.edu>
+ *            Ethan Meitz <emeitz@andrew.cmu.edu>
+=#
+
+#= Broadcast fusion arg buffer tests.
+ * Tests every combination of array and scalar argument positions
+ * to verify the arg_map protocol correctly reconstructs the CUDA
+ * kernel arg buffer.
+ *
+ * Pattern naming convention:
+ *   A = NDArray input
+ *   S = scalar input
+ *   Position = left-to-right in the broadcast expression
+ *   "same" = same array appears multiple times (deduplication test)
+=#
+
+function test_broadcast_fusion(; T=Float32, N=100, atol=1e-5, rtol=1e-5)
+    # Create test arrays with known non-zero values
+    julia_a = rand(T, N)
+    julia_b = rand(T, N)
+    julia_c = rand(T, N)
+
+    a = @allowscalar NDArray(julia_a)
+    b = @allowscalar NDArray(julia_b)
+    c = @allowscalar NDArray(julia_c)
+
+    s1 = T(2.5)
+    s2 = T(1.0)
+    s3 = T(0.5)
+
+    # two different arrays
+    @testset "A + B (two different arrays)" begin
+        expected = julia_a .+ julia_b
+        result = a .+ b
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # same array, deduplication
+    @testset "A + A (same array twice)" begin
+        expected = julia_a .+ julia_a
+        result = a .+ a
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # array then scalar
+    @testset "A + scalar (array first)" begin
+        expected = julia_a .+ s1
+        result = a .+ s1
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # scalar then array
+    @testset "scalar + A (scalar first)" begin
+        expected = s1 .+ julia_a
+        result = s1 .+ a
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # array, two scalars, fused
+    @testset "A * scalar - scalar (fused)" begin
+        expected = julia_a .* s1 .- s2
+        result = a .* s1 .- s2
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # scalar-array-scalar
+    @testset "scalar * A + scalar" begin
+        expected = s1 .* julia_a .+ s2
+        result = s1 .* a .+ s2
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # two arrays then scalar
+    @testset "A + B + scalar" begin
+        expected = julia_a .+ julia_b .+ s1
+        result = a .+ b .+ s1
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # scalar then two arrays
+    @testset "scalar + A + B" begin
+        expected = s1 .+ julia_a .+ julia_b
+        result = s1 .+ a .+ b
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # same array twice + scalar (dedup)
+    @testset "A + A + scalar (dedup + scalar)" begin
+        expected = julia_a .+ julia_a .+ s1
+        result = a .+ a .+ s1
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # three different arrays
+    @testset "A + B + C (three arrays)" begin
+        expected = julia_a .+ julia_b .+ julia_c
+        result = a .+ b .+ c
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # same array three times (triple dedup)
+    @testset "A + A + A (triple dedup)" begin
+        expected = julia_a .+ julia_a .+ julia_a
+        result = a .+ a .+ a
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # two scalars then array
+    @testset "scalar * scalar + A" begin
+        expected = s1 .* s2 .+ julia_a
+        result = s1 .* s2 .+ a
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # multiply, two arrays (different PTX kernel name collision test)
+    @testset "A * B (multiply)" begin
+        expected = julia_a .* julia_b
+        result = a .* b
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # same array squared, dedup
+    @testset "A * A (self multiply, dedup)" begin
+        expected = julia_a .* julia_a
+        result = a .* a
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # subtraction, order matters
+    @testset "A - B (subtraction)" begin
+        expected = julia_a .- julia_b
+        result = a .- b
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # scalar minus array
+    @testset "scalar - A" begin
+        expected = s1 .- julia_a
+        result = s1 .- a
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # three arrays, mixed ops
+    @testset "A * B + C (three arrays, mixed ops)" begin
+        expected = julia_a .* julia_b .+ julia_c
+        result = a .* b .+ c
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # two arrays fused, then scaled
+    @testset "(A + B) * scalar" begin
+        expected = (julia_a .+ julia_b) .* s1
+        result = (a .+ b) .* s1
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # scalar-array pairs
+    @testset "scalar * A + scalar * B" begin
+        expected = s1 .* julia_a .+ s2 .* julia_b
+        result = s1 .* a .+ s2 .* b
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+
+    # same array subtracted, should be all zeros
+    @testset "A - A (same array, expect zeros)" begin
+        expected = julia_a .- julia_a
+        result = a .- a
+        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+    end
+end