diff --git a/.buildkite/developer.pipeline.yml b/.buildkite/developer.pipeline.yml
index 26f4917c4..340cf6644 100644
--- a/.buildkite/developer.pipeline.yml
+++ b/.buildkite/developer.pipeline.yml
@@ -22,9 +22,9 @@ steps:
           LD_LIBRARY_PATH: ""
           LEGATE_AUTO_CONFIG: "0"
           LEGATE_SHOW_CONFIG: "1"
-          GPUTESTS: "1"
           LEGATE_TEST: "1"
-          LEGATE_CONFIG: "--cpus 1 --gpus 1 --utility 1 --fbmem 500 --logging legate=info,level=2 --log-to-file 1"
+          LEGATE_SKIP_RUNTIME: "true"
+          LEGATE_CONFIG: "--cpus 1 --gpus 1 --utility 1 --fbmem 250 --logging legate=info,level=2 --log-to-file 1"
           CUNUMERIC_FUSION: "{{matrix.fusion}}"
         matrix:
           setup:
diff --git a/.buildkite/jll.pipeline.yml b/.buildkite/jll.pipeline.yml
index 7c0a3d312..c76c5f8f3 100644
--- a/.buildkite/jll.pipeline.yml
+++ b/.buildkite/jll.pipeline.yml
@@ -23,7 +23,7 @@ steps:
                 fi
                 julia --project -e 'using Pkg; Pkg.resolve(); Pkg.instantiate()'
           - JuliaCI/julia-test#v1:
-              test_args: "--quickfail"
+              test_args: "--quickfail --jobs=8 --verbose"
           - JuliaCI/julia-coverage#v1:
               dirs:
                 - src
@@ -42,7 +42,8 @@ steps:
           LEGATE_SHOW_CONFIG: "1"
           GPUTESTS: "1"
           LEGATE_TEST: "1"
-          LEGATE_CONFIG: "--cpus 1 --gpus 1 --utility 1 --fbmem 500 --logging legate=info,level=2 --log-to-file 1"
+          LEGATE_SKIP_RUNTIME: "true"
+          LEGATE_CONFIG: "--cpus 1 --gpus 1 --utility 1 --fbmem 250 --logging legate=info,level=2 --log-to-file 1"
         matrix:
           setup:
             julia:
diff --git a/.buildkite/run_developer_ci.sh b/.buildkite/run_developer_ci.sh
index 5fe11c045..537494b28 100755
--- a/.buildkite/run_developer_ci.sh
+++ b/.buildkite/run_developer_ci.sh
@@ -83,5 +83,5 @@ cp LocalPreferences.toml test/LocalPreferences.toml
 
 julia --color=yes --project=. -e '
     using Pkg
-    Pkg.test("cuNumeric"; test_args = ["--quickfail"])
+    Pkg.test("cuNumeric"; test_args = ["--quickfail", "--jobs=8", "--verbose"])
 '
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 9cf4df1de..72d381fbb 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -121,8 +121,8 @@ jobs:
 
       - name: Run tests
         env:
-          GPUTESTS: "0" # parsed by runtests.jl
           LEGATE_SHOW_CONFIG: "1"
           LEGATE_AUTO_CONFIG: "0"
-          LEGATE_CONFIG: "--cpus 1 --utility 1 --sysmem 4000"
-        run: julia --project -e 'using Pkg; Pkg.test(test_args=["--quickfail"])'
+          LEGATE_SKIP_RUNTIME: "true"
+          LEGATE_CONFIG: "--cpus 1 --utility 1 --sysmem 500"
+        run: julia --project -e 'using Pkg; Pkg.test(test_args=["--quickfail", "--jobs=2", "--verbose"])'
diff --git a/.github/workflows/developer.yml b/.github/workflows/developer.yml
index e0494bb22..93cd090ce 100644
--- a/.github/workflows/developer.yml
+++ b/.github/workflows/developer.yml
@@ -125,8 +125,8 @@ jobs:
       # Develop into the workspace root so its members (incl. dev/) resolve to the override.
       - name: Build and test cuNumeric (developer mode)
         env:
-          GPUTESTS: "0"
           LEGATE_CHECKOUT: /tmp/Legate.jl
+          LEGATE_SKIP_RUNTIME: "true"
         run: |
           rm -f Manifest.toml test/Manifest.toml dev/Manifest.toml \
                 LocalPreferences.toml test/LocalPreferences.toml
@@ -157,4 +157,4 @@ jobs:
 
           cp LocalPreferences.toml test/LocalPreferences.toml
 
-          julia --color=yes --project=. -e 'using Pkg; Pkg.test("cuNumeric")'
+          julia --color=yes --project=. -e 'using Pkg; Pkg.test("cuNumeric"; test_args=["--quickfail", "--jobs=2", "--verbose"])'
diff --git a/deps/build.jl b/deps/build.jl
index 68ac071ce..efcc23683 100644
--- a/deps/build.jl
+++ b/deps/build.jl
@@ -42,7 +42,7 @@ function build_cpp_wrapper(
 )
     @info "libcunumeric_jl_wrapper: Building C++ Wrapper Library"
     isdir(install_root) && (rm(install_root; recursive=true); mkdir(install_root))
-    bld_command = `$(joinpath(repo_root, "scripts/build_cpp_wrapper.sh")) $repo_root $cupynumeric_loc $legate_loc $blas_loc $install_root $(Threads.nthreads())`
+    bld_command = `$(joinpath(repo_root, "scripts/build_cpp_wrapper.sh")) $repo_root $cupynumeric_loc $legate_loc $blas_loc $install_root 8`
     return BuildTools.run_build_wrapper_script(
         repo_root, bld_command; cuda_root, cuda_enabled, log_dir=@__DIR__
     )
diff --git a/src/cuNumeric.jl b/src/cuNumeric.jl
index 2a5bcf653..bb27774a6 100644
--- a/src/cuNumeric.jl
+++ b/src/cuNumeric.jl
@@ -260,6 +260,8 @@ function __init__()
     # 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
+    # skip runtime here as well
+    get(ENV, "LEGATE_SKIP_RUNTIME", false) == "true" && return nothing
 
     # Start runtime, but only if not pre-compiling
     ensure_runtime!()
diff --git a/src/utilities/version.jl b/src/utilities/version.jl
index 75d9cbca9..6f05603c6 100644
--- a/src/utilities/version.jl
+++ b/src/utilities/version.jl
@@ -65,6 +65,9 @@ function versioninfo(io::IO=stdout)
     Legate Auto Config: $is_auto_config
     Legate Config:      $legate_config
 
+    Brodcast Fusion:   $(FUSE_BROADCAST_EXPRS)
+    Brodcast Min Ops:  $(FUSE_BROADCAST_MIN_OPS)
+
     Hostname:         $hostname
     Julia Version:    $(VERSION)
     C++ Compiler:     $compiler
diff --git a/test/Project.toml b/test/Project.toml
index 82e5ef337..73f4b047f 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -2,7 +2,10 @@
 CNPreferences = "3e078157-ea10-49d5-bf32-908f777cd46f"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 CUDA_Driver_jll = "4ee394cb-3365-5eb0-8335-949819d2adfc"
+InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/tests/lifetime.jl b/test/analysis/lifetime.jl
similarity index 100%
rename from test/tests/lifetime.jl
rename to test/analysis/lifetime.jl
diff --git a/test/analysis/promotion.jl b/test/analysis/promotion.jl
new file mode 100644
index 000000000..7b1fb39a7
--- /dev/null
+++ b/test/analysis/promotion.jl
@@ -0,0 +1,21 @@
+@testset "Type and Shape Promotion" begin
+    N = 100
+    cunumeric_arr1 = cuNumeric.zeros(Float64, N)
+    cunumeric_arr3 = cuNumeric.zeros(Float32, N)
+    cunumeric_int64 = cuNumeric.zeros(Int64, N)
+    cunumeric_int32 = cuNumeric.zeros(Int32, N)
+    cunumeric_arr5 = cuNumeric.zeros(Float64, N - 1, N - 1)
+    @test_throws "Implicit promotion" cunumeric_arr3 .+ cunumeric_arr1
+    @test_throws "Implicit promotion" map(+, cunumeric_arr3, cunumeric_arr1)
+    @test_throws DimensionMismatch cunumeric_arr1 .+ cunumeric_arr5
+    @test_throws DimensionMismatch cunumeric_arr1 ./ cunumeric_arr5
+
+    allowscalar() do
+        @test safe_compare(
+            cunumeric_arr1, cunumeric_int64 .+ cunumeric_arr1, atol(Float64), rtol(Float64)
+        )
+        r1 = @allowpromotion cunumeric_arr3 .+ cunumeric_arr1
+        r2 = @allowpromotion map(+, cunumeric_arr3, cunumeric_arr1)
+        @test safe_compare(r1, r2, atol(Float64), rtol(Float64))
+    end
+end
diff --git a/test/tests/stability.jl b/test/analysis/type_stability.jl
similarity index 100%
rename from test/tests/stability.jl
rename to test/analysis/type_stability.jl
diff --git a/test/tests/axpy_advanced.jl b/test/array/axpy.jl
similarity index 80%
rename from test/tests/axpy_advanced.jl
rename to test/array/axpy.jl
index 18950f7c1..6386d8814 100644
--- a/test/tests/axpy_advanced.jl
+++ b/test/array/axpy.jl
@@ -35,8 +35,8 @@ function axpy_advanced(T, N)
     α = T(56.6)
 
     # base Julia arrays
-    x_cpu = Base.zeros(T, dims);
-    y_cpu = Base.zeros(T, dims);
+    x_cpu = Base.zeros(T, dims)
+    y_cpu = Base.zeros(T, dims)
 
     # cunumeric arrays
     x = cuNumeric.zeros(T, dims)
@@ -113,3 +113,39 @@ function axpy_advanced(T, N)
         @test safe_compare(result_cpu, result, atol(T), rtol(T))
     end
 end
+
+function axpy_basic(T, N)
+    α = T(56.6)
+    dims = (N, N)
+
+    # Base julia arrays
+    x_cpu = rand(T, dims)
+    y_cpu = rand(T, dims)
+
+    # cunumeric arrays
+    x = cuNumeric.zeros(T, dims)
+    y = cuNumeric.zeros(T, dims)
+
+    # Initialize NDArrays with same random values as julia arrays
+    @allowscalar for i in 1:N
+        for j in 1:N
+            x[i, j] = x_cpu[i, j]
+            y[i, j] = y_cpu[i, j]
+        end
+    end
+
+    result = α .* x .+ y
+    result_cpu = α .* x_cpu .+ y_cpu
+    allowscalar() do
+        @test safe_compare(result, result_cpu, atol(T), rtol(T))
+        @test safe_compare(result_cpu, result, atol(T), rtol(T))
+    end
+end
+
+@testset "AXPY" begin
+    N = 100
+    @testset verbose = true for T in Base.uniontypes(cuNumeric.SUPPORTED_FLOAT_TYPES)
+        @testset "basic" axpy_basic(T, N)
+        @testset "advanced" axpy_advanced(T, N)
+    end
+end
diff --git a/test/array/binary/float.jl b/test/array/binary/float.jl
new file mode 100644
index 000000000..490d92773
--- /dev/null
+++ b/test/array/binary/float.jl
@@ -0,0 +1,27 @@
+#= 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>
+=#
+
+include(joinpath(@__DIR__, "tests.jl"))
+
+run_binary_ops_tests(
+    filter(
+        T -> T <: AbstractFloat || T <: Complex,
+        Base.uniontypes(cuNumeric.SUPPORTED_ARRAY_TYPES),
+    ),
+)
diff --git a/test/array/binary/int.jl b/test/array/binary/int.jl
new file mode 100644
index 000000000..822526e60
--- /dev/null
+++ b/test/array/binary/int.jl
@@ -0,0 +1,28 @@
+#= 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>
+=#
+
+include(joinpath(@__DIR__, "tests.jl"))
+
+run_binary_ops_tests(
+    filter(
+        T -> T === Bool || T <: Signed,
+        Base.uniontypes(cuNumeric.SUPPORTED_ARRAY_TYPES),
+    ),
+)
+run_binary_copyto_tests()
diff --git a/test/array/binary/powers.jl b/test/array/binary/powers.jl
new file mode 100644
index 000000000..e7d96cfb9
--- /dev/null
+++ b/test/array/binary/powers.jl
@@ -0,0 +1,111 @@
+#= 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>
+=#
+
+@testset "Powers" begin
+    N = 9
+
+    get_pwrs(::Type{I}) where {I<:Integer} = I.([-10, -5, -2, -1, 0, 1, 2, 5, 10])
+    get_pwrs(::Type{F}) where {F<:AbstractFloat} = F.([-3.141, -2, -1, 0, 1, 2, 3.2, 4.41, 6.233])
+    get_pwrs(::Type{Bool}) = [true, false, true, false, false, true, false, true, true]
+
+    # TYPES = Base.uniontypes(cuNumeric.SUPPORTED_ARRAY_TYPES)
+    TYPES = Base.uniontypes(cuNumeric.SUPPORTED_FLOAT_TYPES)
+
+    @testset "$(BT) ^ $(PT)" for (BT, PT) in Iterators.product(TYPES, TYPES)
+        base_jl = my_rand(BT, N)
+
+        if BT <: Union{Bool,Int32} && PT == Int32
+            # julia doesnt like Int32 powers
+            pwrs = Float32.(get_pwrs(PT))
+        elseif BT <: Union{Bool,Int32,Int64} && PT <: Union{Int32,Int64}
+            # julia doesnt like Int64 powers on bool
+            pwrs = Float64.(get_pwrs(PT))
+        elseif (PT <: AbstractFloat) && (BT <: AbstractFloat || BT <: Signed)
+            # Things like -387 ^ 3.2 will be Complex and error
+            pwrs = get_pwrs(PT)
+            base_jl = abs.(base_jl)
+        else
+            pwrs = get_pwrs(PT)
+        end
+
+        base_cn = @allowscalar NDArray(base_jl)
+        pwrs_cn = @allowscalar NDArray(pwrs)
+
+        # we deviate a bit form Julia here
+        if (PT <: Union{Int32,Int64}) && (BT <: Union{Bool,Int32,Int64})
+            T_OUT = cuNumeric.__my_promote_type(typeof(^), BT, PT)
+        else
+            T_OUT = Base.promote_op(Base.:(^), BT, PT)
+        end
+
+        TEST_BROKEN = (BT <: Union{Int32,Int64} && PT == Bool)
+
+        allowpromotion(true) do
+            allowscalar() do
+                # Power is array
+                @test safe_compare(
+                    base_jl .^ pwrs, base_cn .^ pwrs_cn, atol(T_OUT), rtol(T_OUT)
+                ) skip=TEST_BROKEN
+
+                # Power is scalar
+                for p in pwrs
+                    @test safe_compare(base_jl .^ p, base_cn .^ p, atol(T_OUT), rtol(T_OUT))
+                end
+            end
+        end
+    end
+
+    @testset verbose = true "Reciprocal" begin
+        @testset for T in TYPES
+            arr_jl = Random.rand(T, N)
+            arr_cn = @allowscalar NDArray(arr_jl)
+
+            # Differ from Julia here
+            T_OUT = cuNumeric.__recip_type(T)
+
+            # Cast julia result to whatever we do
+            res_jl = T_OUT.(arr_jl .^ -1)
+            allowpromotion(true) do
+                res_cn = arr_cn .^ -1
+                res_cn2 = inv.(arr_cn)
+                allowscalar() do
+                    @test safe_compare(res_jl, res_cn, atol(T_OUT), rtol(T_OUT))
+                    @test safe_compare(res_jl, res_cn2, atol(T_OUT), rtol(T_OUT))
+                end
+            end
+        end
+    end
+
+    @testset verbose = true "Square" begin
+        @testset for T in TYPES
+            arr_jl = Random.rand(T, N)
+            arr_cn = @allowscalar NDArray(arr_jl)
+
+            T_OUT = Base.promote_op(Base.:(^), T, Int64)
+            res_jl = arr_jl .^ 2
+            res_cn = arr_cn .^ 2
+
+            allowpromotion(true) do
+                allowscalar() do
+                    @test safe_compare(res_jl, res_cn, atol(T_OUT), rtol(T_OUT))
+                end
+            end
+        end
+    end
+end
diff --git a/test/tests/binary_tests.jl b/test/array/binary/tests.jl
similarity index 52%
rename from test/tests/binary_tests.jl
rename to test/array/binary/tests.jl
index d4c9e64d3..794c2e19a 100644
--- a/test/tests/binary_tests.jl
+++ b/test/array/binary/tests.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");
@@ -60,3 +60,60 @@ function test_binary_function_set(func_dict, T, N)
         test_binary_operation(func, arrs_jl[3:4]..., arrs_cunum[3:4]..., T)
     end
 end
+
+function run_binary_ops_tests(types)
+    @testset "Binary Ops" begin
+        N = 100
+
+        @testset for T in types
+            allowpromotion(true) do
+                test_binary_function_set(cuNumeric.floaty_binary_op_map, T, N)
+                return test_binary_function_set(cuNumeric.binary_op_map, T, N)
+            end
+
+            arr_jl = my_rand(T, N)
+            arr_jl2 = my_rand(T, N)
+            arr_cn = @allowscalar NDArray(arr_jl)
+            arr_cn2 = @allowscalar NDArray(arr_jl2)
+
+            # lcm/gcd require specific handling for integers and avoid overflow
+            if T <: cuNumeric.SUPPORTED_INT_TYPES && T != Bool
+                range_limit = (T == Int8 || T == UInt8) ? 10 : 100
+                arr_jl_small = my_rand(T, N; L=1, R=range_limit)
+                arr_jl2_small = my_rand(T, N; L=1, R=range_limit)
+                arr_cn_small = @allowscalar NDArray(arr_jl_small)
+                arr_cn2_small = @allowscalar NDArray(arr_jl2_small)
+
+                allowscalar() do
+                    @test safe_compare(
+                        lcm.(arr_jl_small, arr_jl2_small), lcm.(arr_cn_small, arr_cn2_small),
+                        atol(T),
+                        rtol(T),
+                    )
+                    @test safe_compare(
+                        gcd.(arr_jl_small, arr_jl2_small), gcd.(arr_cn_small, arr_cn2_small),
+                        atol(T),
+                        rtol(T),
+                    )
+                end
+            end
+
+            allowscalar() do
+                @test unwrap(arr_cn == arr_cn)
+                @test !unwrap(arr_cn == arr_cn2)
+                @test unwrap(arr_cn != arr_cn2)
+                @test !unwrap(arr_cn != arr_cn)
+                @test unwrap(all(arr_cn .== arr_cn))
+            end
+        end
+    end
+end
+
+function run_binary_copyto_tests()
+    @testset "Copy-To" begin
+        a = cuNumeric.zeros(2, 2)
+        b = cuNumeric.ones(2, 2)
+        copyto!(a, b)
+        @test is_same(a, b)
+    end
+end
diff --git a/test/array/binary/uint.jl b/test/array/binary/uint.jl
new file mode 100644
index 000000000..2f1cba577
--- /dev/null
+++ b/test/array/binary/uint.jl
@@ -0,0 +1,27 @@
+#= 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>
+=#
+
+include(joinpath(@__DIR__, "tests.jl"))
+
+run_binary_ops_tests(
+    filter(
+        T -> T <: Unsigned,
+        Base.uniontypes(cuNumeric.SUPPORTED_ARRAY_TYPES),
+    ),
+)
diff --git a/test/tests/elementwise.jl b/test/array/broadcast_basic.jl
similarity index 53%
rename from test/tests/elementwise.jl
rename to test/array/broadcast_basic.jl
index 49f22bb2d..e9eb95c04 100644
--- a/test/tests/elementwise.jl
+++ b/test/array/broadcast_basic.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");
@@ -21,7 +21,29 @@
     - [./ .- .+ .*] operators are tested
 =#
 
-function elementwise()
+function operator(u, v)
+    dx=0.1
+    dt = dx/5
+    c_u=1.0
+    c_v=0.3
+    f=0.03
+    k=0.06
+
+    # calculate F_u and F_v functions
+    F_u = (-u*(v * v)) + f*(1 .- u)
+    F_v = (u*(v * v)) - (f+k)*v
+
+    # 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 - 2*u + u) ./ dx^2 + (u - 2*u + u) ./ dx^2
+    v_lap = (v - 2*v + v) ./ dx^2 + (v - 2*v + v) ./ dx^2
+
+    # Forward-Euler time step for all points except the boundaries
+    u = ((c_u * u_lap) + F_u) * dt + u
+    return v = ((c_v * v_lap) + F_v) * dt + v
+end
+
+@testset "[./ .- .+ .*]" begin
     N = 100
     dims = (N, N)
 
@@ -110,24 +132,62 @@ function elementwise()
     end
 end
 
-function operator(u, v)
-    dx=0.1;
-    dt = dx/5;
-    c_u=1.0;
-    c_v=0.3;
-    f=0.03;
-    k=0.06;
-
-    # calculate F_u and F_v functions
-    F_u = (-u*(v * v)) + f*(1 .- u)
-    F_v = (u*(v * v)) - (f+k)*v
-
-    # 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 - 2*u + u) ./ dx^2 + (u - 2*u + u) ./ dx^2
-    v_lap = (v - 2*v + v) ./ dx^2 + (v - 2*v + v) ./ dx^2
+#TODO LOOP BINARY OPS WITH SCALARS
+@testset verbose = true "Scalars" begin
+    N = 10
+
+    for T in Base.uniontypes(cuNumeric.SUPPORTED_NUMERIC_TYPES)
+        julia_arr = my_rand(T, N)
+        julia_arr_2D = my_rand(T, N, N)
+
+        s = Random.rand(T)
+
+        allowscalar() do
+            cunumeric_arr = NDArray(julia_arr)
+            cunumeric_arr_2D = NDArray(julia_arr_2D)
+            allowpromotion(true) do
+                for cn_arr in (cunumeric_arr, cunumeric_arr_2D)
+                    @test safe_compare(s * julia_arr, s * cunumeric_arr, atol(T), rtol(T))
+                    @test safe_compare(julia_arr * s, cunumeric_arr * s, atol(T), rtol(T))
+                    @test safe_compare(s .* julia_arr, s .* cunumeric_arr, atol(T), rtol(T))
+                    @test safe_compare(julia_arr .* s, cunumeric_arr .* s, atol(T), rtol(T))
+                    @test safe_compare(s .+ julia_arr, s .+ cunumeric_arr, atol(T), rtol(T))
+                    @test safe_compare(julia_arr .+ s, cunumeric_arr .+ s, atol(T), rtol(T))
+                    @test safe_compare(s .- julia_arr, s .- cunumeric_arr, atol(T), rtol(T))
+                    @test safe_compare(julia_arr .- s, cunumeric_arr .- s, atol(T), rtol(T))
+                    @test safe_compare(s ./ julia_arr, s ./ cunumeric_arr, atol(T), rtol(T))
+                    @test safe_compare(
+                        s .* julia_arr .+ s, s .* cunumeric_arr .+ s, atol(T), rtol(T)
+                    )
+                    @test s + s ≈ (NDArray(s) + NDArray(s))[] # atol=atol(T) r_tol=atol(T)
+                    # @test s * s ≈ (NDArray(s) * NDArray(s))[] # atol=atol(T) r_tol=atol(T)
+                end
+            end
+        end
+    end
 
-    # Forward-Euler time step for all points except the boundaries
-    u = ((c_u * u_lap) + F_u) * dt + u
-    v = ((c_v * v_lap) + F_v) * dt + v
+    # Boolean things
+    allowpromotion() do
+        allowscalar() do
+            julia_arr = [-3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
+            cunumeric_arr = NDArray(julia_arr)
+
+            @test safe_compare(true * julia_arr, true * cunumeric_arr, atol(Int), rtol(Int))
+            @test safe_compare(false * julia_arr, false * cunumeric_arr, atol(Int), rtol(Int))
+            @test safe_compare(true .* julia_arr, true .* cunumeric_arr, atol(Int), rtol(Int))
+            @test safe_compare(
+                false .* julia_arr, false .* cunumeric_arr, atol(Int), rtol(Int)
+            )
+
+            julia_arr = [true, false, true, false, false, true, true]
+            cunumeric_arr = NDArray(julia_arr)
+
+            @test safe_compare(4 * julia_arr, 4 * cunumeric_arr, atol(Int), rtol(Int))
+            @test safe_compare(4 .* julia_arr, 4 .* cunumeric_arr, atol(Int), rtol(Int))
+            @test safe_compare(4 .+ julia_arr, 4 .+ cunumeric_arr, atol(Int), rtol(Int))
+            @test safe_compare(
+                julia_arr ./ 3, cunumeric_arr ./ 3, atol(Float64), rtol(Float64)
+            )
+        end
+    end
 end
diff --git a/test/tests/linalg.jl b/test/array/linalg.jl
similarity index 71%
rename from test/tests/linalg.jl
rename to test/array/linalg.jl
index 10703f7c2..836d47249 100644
--- a/test/tests/linalg.jl
+++ b/test/array/linalg.jl
@@ -18,6 +18,81 @@
  *            Nader Rahal <naderrahhal2026@u.northwestern.edu>
 =#
 
+function gemm(N, M, T, max_diff)
+    if T == Bool
+        a = cuNumeric.trues(5, 5)
+        b = cuNumeric.as_type(cuNumeric.trues(5, 5), Float32)
+        c = cuNumeric.as_type(cuNumeric.trues(5, 5), Float64)
+        @test_throws ArgumentError a * a # Bool * Bool not supported
+        @allowpromotion d = a * b
+        @allowpromotion e = a * c
+        @test @allowscalar safe_compare(5 * ones(Float32, 5, 5), d, 0.0, max_diff)
+        @test @allowscalar safe_compare(5 * ones(Float64, 5, 5), e, 0.0, max_diff)
+        return nothing
+    end
+
+    if T <: Integer
+        a = cuNumeric.ones(Int32, 5, 5)
+        a_jl = ones(Int32, 5, 5)
+        b = cuNumeric.ones(Float32, 5, 5)
+        b_jl = ones(Float32, 5, 5)
+        @test_throws ArgumentError a * a
+        @test @allowscalar safe_compare(a_jl * b_jl, a * b, 0.0, max_diff)
+        return nothing
+    end
+
+    dims_to_test = [(N, N), (N, M), (M, N)]
+
+    @testset for dims in dims_to_test
+        # Base julia arrays
+        A_cpu = rand(T, dims[1], dims[2])
+        B_cpu = rand(T, dims[2], dims[1])
+        C_out_cpu = zeros(T, dims[1], dims[1])
+
+        # cunumeric arrays
+        A = cuNumeric.NDArray(A_cpu)
+        B = cuNumeric.NDArray(B_cpu)
+        C_out = cuNumeric.zeros(T, dims[1], dims[1])
+
+        # Julia result
+        C_cpu = A_cpu * B_cpu
+        LinearAlgebra.mul!(C_out_cpu, A_cpu, B_cpu)
+
+        @test C_cpu == C_out_cpu # really just making sure test is written right...
+
+        A = cuNumeric.as_type(A, T)
+        B = cuNumeric.as_type(B, T)
+        C = cuNumeric.zeros(T, N, N)
+
+        C = A * B
+        LinearAlgebra.mul!(C_out, A, B)
+
+        allowscalar() do
+            @test isapprox(C, C_cpu, rtol=max_diff)
+            @test isapprox(C, C_out, rtol=max_diff)
+
+            if T != Float64
+                C_wider = cuNumeric.zeros(Float64, dims[1], dims[1])
+                @test_throws "Implicit promotion" LinearAlgebra.mul!(C_wider, A, B)
+            end
+        end
+
+        # Integer output with FP input
+        if !(T <: Integer)
+            bad = cuNumeric.zeros(Int, dims[1], dims[1])
+            @test_throws ArgumentError mul!(bad, A, B)
+        end
+    end
+end
+
+@testset "GEMM" begin
+    N = 50
+    M = 25
+    @testset verbose = true for T in Base.uniontypes(cuNumeric.SUPPORTED_NUMERIC_TYPES)
+        gemm(N, M, T, rtol(T))
+    end
+end
+
 @testset "transpose" begin
     @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_NUMERIC_TYPES)
         A = my_rand(T, 4, 3)
@@ -27,7 +102,7 @@
         out = cuNumeric.transpose(nda)
 
         allowscalar() do
-            @test cuNumeric.compare(ref, out, atol(T), rtol(T))
+            @test safe_compare(ref, out, atol(T), rtol(T))
         end
     end
 end
@@ -38,7 +113,7 @@ end
         ref = Matrix{T}(I, n, n)
         out = cuNumeric.eye(T, n)
         allowscalar() do
-            @test cuNumeric.compare(ref, out, atol(T), rtol(T))
+            @test safe_compare(ref, out, atol(T), rtol(T))
         end
     end
 end
@@ -81,7 +156,7 @@ end
             out = cuNumeric.diag(nda; k=k)
 
             allowscalar() do
-                @test cuNumeric.compare(ref, out, atol(T), rtol(T))
+                @test safe_compare(ref, out, atol(T), rtol(T))
             end
         end
     end
@@ -95,7 +170,7 @@ end
 #     out = cuNumeric.ravel(nda)
 
 #     allowscalar() do
-#         @test cuNumeric.compare(ref, out, atol(Int32), rtol(Int32))
+#         @test safe_compare(ref, out, atol(Int32), rtol(Int32))
 #     end
 # end
 
@@ -122,7 +197,7 @@ end
         end
         x = cuNumeric.solve(A, b)
         allowscalar() do
-            @test cuNumeric.compare(fill(T(0.25), n, 1), x, atol(T), rtol(T))
+            @test safe_compare(fill(T(0.25), n, 1), x, atol(T), rtol(T))
         end
     end
 end
@@ -135,7 +210,7 @@ end
         x = cuNumeric.solve(A, b)
         ref = reshape(T.(collect(1:n)), n, 1)
         allowscalar() do
-            @test cuNumeric.compare(ref, x, atol(T), rtol(T))
+            @test safe_compare(ref, x, atol(T), rtol(T))
         end
     end
 end
@@ -149,7 +224,7 @@ end
         x = cuNumeric.solve(A, b)
         ref = A_ref \ b_ref
         allowscalar() do
-            @test cuNumeric.compare(ref, x, atol(T), rtol(T))
+            @test safe_compare(ref, x, atol(T), rtol(T))
         end
     end
 end
@@ -162,7 +237,7 @@ end
         @test ndims(x) == 1
         ref = A_ref \ b_ref
         allowscalar() do
-            @test cuNumeric.compare(ref, x, atol(T), rtol(T))
+            @test safe_compare(ref, x, atol(T), rtol(T))
         end
     end
 end
@@ -180,27 +255,26 @@ end
             x = cuNumeric.solve(A, b)
             ref = Float64[1 0; 0 1] \ Float64[1; 1;;]
             allowscalar() do
-                @test cuNumeric.compare(ref, x, atol(Float64), rtol(Float64))
+                @test safe_compare(ref, x, atol(Float64), rtol(Float64))
             end
         end
     end
 end
 
-
 function check_svd_reconstruction(ref_A::AbstractMatrix, u, s, vh, tol_a, tol_r)
-    U  = Array(u)
-    S  = Array(s)
+    U = Array(u)
+    S = Array(s)
     Vh = Array(vh)
     A_rec = U * Diagonal(S) * Vh
     return isapprox(ref_A, A_rec; atol=tol_a, rtol=tol_r)
 end
 
 function check_svd_orthonormality(u, vh, tol_a, tol_r)
-    U  = Array(u)
+    U = Array(u)
     Vh = Array(vh)
     ku = size(U, 2)
     kv = size(Vh, 1)
-    ok_u  = isapprox(U'  * U,  Matrix{eltype(U)}(I, ku, ku);  atol=tol_a, rtol=tol_r)
+    ok_u = isapprox(U' * U, Matrix{eltype(U)}(I, ku, ku); atol=tol_a, rtol=tol_r)
     ok_vh = isapprox(Vh * Vh', Matrix{eltype(Vh)}(I, kv, kv); atol=tol_a, rtol=tol_r)
     return ok_u && ok_vh
 end
@@ -208,7 +282,7 @@ end
 @testset "svd square matrix" begin
     @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_SVD_TYPES)
         A_ref = my_rand(T, 5, 5)
-        nda   = cuNumeric.NDArray(A_ref)
+        nda = cuNumeric.NDArray(A_ref)
         u, s, vh = cuNumeric.svd(nda)
         allowscalar() do
             @test check_svd_reconstruction(A_ref, u, s, vh, atol(T), rtol(T))
@@ -220,7 +294,7 @@ end
 @testset "svd tall matrix (m > n)" begin
     @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_SVD_TYPES)
         A_ref = my_rand(T, 6, 4)
-        nda   = cuNumeric.NDArray(A_ref)
+        nda = cuNumeric.NDArray(A_ref)
         u, s, vh = cuNumeric.svd(nda, false)  # thin SVD for reconstruction test
         allowscalar() do
             @test check_svd_reconstruction(A_ref, u, s, vh, atol(T), rtol(T))
@@ -228,41 +302,41 @@ end
         end
     end
 end
- 
+
 @testset "svd thin output shapes (full_matrices=false)" begin
     @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_SVD_TYPES)
         m, n = 6, 4
-        k     = min(m, n)
+        k = min(m, n)
         A_ref = my_rand(T, m, n)
-        nda   = cuNumeric.NDArray(A_ref)
+        nda = cuNumeric.NDArray(A_ref)
         u, s, vh = cuNumeric.svd(nda, false)
         allowscalar() do
-            @test size(Array(u))  == (m, k)
-            @test size(Array(s))  == (k,)
+            @test size(Array(u)) == (m, k)
+            @test size(Array(s)) == (k,)
             @test size(Array(vh)) == (k, n)
             @test check_svd_reconstruction(A_ref, u, s, vh, atol(T), rtol(T))
         end
     end
 end
- 
+
 @testset "svd full output shapes (full_matrices=true)" begin
     @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_SVD_TYPES)
         m, n = 6, 4
         A_ref = my_rand(T, m, n)
-        nda   = cuNumeric.NDArray(A_ref)
+        nda = cuNumeric.NDArray(A_ref)
         u, s, vh = cuNumeric.svd(nda, true)
         allowscalar() do
-            @test size(Array(u))  == (m, m)
-            @test size(Array(s))  == (min(m, n),)
+            @test size(Array(u)) == (m, m)
+            @test size(Array(s)) == (min(m, n),)
             @test size(Array(vh)) == (n, n)
         end
     end
 end
- 
+
 @testset "svd singular values non-negative and sorted" begin
     @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_SVD_TYPES)
         A_ref = my_rand(T, 5, 5)
-        nda   = cuNumeric.NDArray(A_ref)
+        nda = cuNumeric.NDArray(A_ref)
         _, s, _ = cuNumeric.svd(nda)
         allowscalar() do
             sv = Array(s)
@@ -271,27 +345,27 @@ end
         end
     end
 end
- 
+
 @testset "svd identity matrix" begin
     @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_SVD_TYPES)
-        n     = 4
+        n = 4
         A_ref = Matrix{T}(I, n, n)
-        nda   = cuNumeric.NDArray(A_ref)
+        nda = cuNumeric.NDArray(A_ref)
         _, s, _ = cuNumeric.svd(nda)
         allowscalar() do
-            @test cuNumeric.compare(ones(T, n), s, atol(T), rtol(T))
+            @test safe_compare(ones(T, n), s, atol(T), rtol(T))
         end
     end
 end
- 
+
 @testset "svd rank-1 matrix" begin
     @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_SVD_TYPES)
         # outer product of two vectors: exactly one nonzero singular value
         # 5x4 satisfies the M >= N constraint
-        v1    = T.(collect(1:5))
-        v2    = T.(collect(1:4))
+        v1 = T.(collect(1:5))
+        v2 = T.(collect(1:4))
         A_ref = v1 * v2'
-        nda   = cuNumeric.NDArray(A_ref)
+        nda = cuNumeric.NDArray(A_ref)
         _, s, _ = cuNumeric.svd(nda)
         allowscalar() do
             sv = Array(s)
@@ -316,8 +390,8 @@ end
 @testset "qr reconstruction" begin
     @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_QR_TYPES)
         A_ref = my_rand(T, 6, 4)
-        nda   = cuNumeric.NDArray(A_ref)
-        q, r  = cuNumeric.qr(nda)
+        nda = cuNumeric.NDArray(A_ref)
+        q, r = cuNumeric.qr(nda)
         allowscalar() do
             Q = Array(q)
             R = Array(r)
@@ -332,8 +406,8 @@ end
 @testset "qr square matrix" begin
     @testset verbose=true for T in Base.uniontypes(cuNumeric.SUPPORTED_QR_TYPES)
         A_ref = my_rand(T, 5, 5)
-        nda   = cuNumeric.NDArray(A_ref)
-        q, r  = cuNumeric.qr(nda)
+        nda = cuNumeric.NDArray(A_ref)
+        q, r = cuNumeric.qr(nda)
         allowscalar() do
             Q = Array(q)
             R = Array(r)
@@ -353,8 +427,10 @@ end
             q, r = cuNumeric.qr(A)
             allowscalar() do
                 @test eltype(Array(q)) == Float64
-                @test isapprox(Float64.(vals), Array(q) * Array(r); atol=atol(Float64), rtol=rtol(Float64))
+                @test isapprox(
+                    Float64.(vals), Array(q) * Array(r); atol=atol(Float64), rtol=rtol(Float64)
+                )
             end
         end
     end
-end
\ No newline at end of file
+end
diff --git a/test/tests/slicing.jl b/test/array/slicing.jl
similarity index 81%
rename from test/tests/slicing.jl
rename to test/array/slicing.jl
index 0a8c439c0..2e2a539a4 100644
--- a/test/tests/slicing.jl
+++ b/test/array/slicing.jl
@@ -32,7 +32,7 @@ struct Params{T}
     function Params(
         ::Type{T}, dx=T(0.1), c_u=T(1.0), c_v=T(0.3), f=T(0.03), k=T(0.06)
     ) where {T<:AbstractFloat}
-        new{T}(dx, dx/T(5), c_u, c_v, f, k)
+        return new{T}(dx, dx/T(5), c_u, c_v, f, k)
     end
 end
 
@@ -44,8 +44,8 @@ function slicing(T, N)
     u = cuNumeric.zeros(T, dims)
     v = cuNumeric.zeros(T, dims)
 
-    u_cpu = rand(T, dims);
-    v_cpu = rand(T, dims);
+    u_cpu = rand(T, dims)
+    v_cpu = rand(T, dims)
 
     @allowscalar for i in 1:N
         for j in 1:N
@@ -64,10 +64,10 @@ function slicing(T, N)
     step(u, v, u_new, v_new, args)
 
     allowscalar() do
-        @test cuNumeric.compare(u, u_cpu, atol(T), rtol(T))
-        @test cuNumeric.compare(v, v_cpu, atol(T), rtol(T))
-        @test cuNumeric.compare(u_new, u_new_cpu, atol(T), rtol(T))
-        @test cuNumeric.compare(v_new, v_new_cpu, atol(T), rtol(T))
+        @test safe_compare(u, u_cpu, atol(T), rtol(T))
+        @test safe_compare(v, v_cpu, atol(T), rtol(T))
+        @test safe_compare(u_new, u_new_cpu, atol(T), rtol(T))
+        @test safe_compare(v_new, v_new_cpu, atol(T), rtol(T))
     end
 end
 
@@ -122,5 +122,24 @@ function step(u, v, u_new, v_new, args::Params)
     v_new[:, 1] = v[:, end - 1]
     v_new[:, end] = v[:, 2]
     v_new[1, :] = v[end - 1, :]
-    v_new[end, :] = v[2, :]
+    return v_new[end, :] = v[2, :]
+end
+
+@testset "Array Slices" begin
+    N = 100
+    @testset for T in Base.uniontypes(cuNumeric.SUPPORTED_FLOAT_TYPES)
+        slicing(T, N)
+    end
+end
+
+@testset "Basic Array Accessors" begin
+    for T in Base.uniontypes(cuNumeric.SUPPORTED_ARRAY_TYPES)
+        value = rand(T)
+        arr = cuNumeric.zeros(T, 2, 2)
+
+        allowscalar() do
+            arr[1, 2] = value
+            @test arr[1, 2] == value
+        end
+    end
 end
diff --git a/test/array/unary/float.jl b/test/array/unary/float.jl
new file mode 100644
index 000000000..242b3e70a
--- /dev/null
+++ b/test/array/unary/float.jl
@@ -0,0 +1,27 @@
+#= 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>
+=#
+
+include(joinpath(@__DIR__, "tests.jl"))
+
+run_unary_tests(
+    filter(
+        T -> T <: AbstractFloat || T <: Complex,
+        Base.uniontypes(cuNumeric.SUPPORTED_ARRAY_TYPES),
+    ),
+)
diff --git a/test/array/unary/int.jl b/test/array/unary/int.jl
new file mode 100644
index 000000000..f9f9f3632
--- /dev/null
+++ b/test/array/unary/int.jl
@@ -0,0 +1,28 @@
+#= 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>
+=#
+
+include(joinpath(@__DIR__, "tests.jl"))
+
+run_unary_tests(
+    filter(
+        T -> T === Bool || T <: Signed,
+        Base.uniontypes(cuNumeric.SUPPORTED_ARRAY_TYPES),
+    );
+    include_bool_reductions=true,
+)
diff --git a/test/array/unary/tests.jl b/test/array/unary/tests.jl
new file mode 100644
index 000000000..558a2c877
--- /dev/null
+++ b/test/array/unary/tests.jl
@@ -0,0 +1,239 @@
+#= 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>
+=#
+
+# Map functions to their required domains
+const SPECIAL_DOMAINS = Dict(
+    Base.acosh => :greater_than_one,
+    Base.log => :positive,
+    Base.log10 => :positive,
+    Base.log2 => :positive,
+    Base.log1p => :positive, # technicaly anything > -1
+    Base.sqrt => :positive,
+)
+
+function test_unary_operation(func, julia_arr, cunumeric_arr, T)
+    T_OUT = Base.promote_op(func, T)
+
+    # Pre-allocate output arrays
+    cunumeric_in_place = cuNumeric.zeros(T_OUT, size(julia_arr)...)
+
+    # Compute results using different methods
+    julia_res = func.(julia_arr)
+
+    cunumeric_res = func.(cunumeric_arr)
+    cunumeric_in_place .= func.(cunumeric_arr)
+    cunumeric_res2 = map(func, cunumeric_arr)
+
+    allowscalar() do
+        @test safe_compare(julia_res, cunumeric_in_place, atol(T_OUT), rtol(T_OUT))
+        @test safe_compare(julia_res, cunumeric_res, atol(T_OUT), rtol(T_OUT))
+        @test safe_compare(julia_res, cunumeric_res2, atol(T_OUT), rtol(T_OUT))
+    end
+end
+
+skip_on_integer = (Base.acosh, Base.atanh, Base.atan, Base.acos, Base.asin)
+skip_on_bool = (Base.:(-), skip_on_integer...)
+skip_on_complex = (
+    Base.tanh,
+    Base.deg2rad, Base.rad2deg, Base.sign, Base.cbrt,
+    Base.exp2, Base.expm1, Base.log10, Base.log1p, Base.log2,
+    Base.acos, Base.asin, Base.atan, Base.acosh, Base.asinh, Base.atanh,
+)
+
+function test_unary_function_set(func_dict, T, N)
+    default_generator = (T == Bool) ? :uniform : :unit_interval
+
+    @testset "$func" for func in keys(func_dict)
+        if func in skip_on_complex && (T <: Complex)
+            continue
+        end
+
+        # The are only defined for like 3 integers (-1, 0, 1) so just skip them
+        if func in skip_on_integer && (T <: Integer)
+            continue
+        end
+
+        if func in skip_on_bool && (T == Bool)
+            continue
+        end
+
+        domain_type = get(SPECIAL_DOMAINS, func, default_generator)
+
+        # :uniform is the only generator capable of generating bits
+        skip = (T == Bool && domain_type != :uniform)
+        skip && continue
+
+        julia_arr_1D, julia_arr_2D = make_julia_arrays(T, N, domain_type)
+        cunumeric_arr_1D, cunumeric_arr_2D = make_cunumeric_arrays(
+            [julia_arr_1D], [julia_arr_2D], T, N
+        )
+
+        test_unary_operation(func, julia_arr_1D, cunumeric_arr_1D, T)
+        test_unary_operation(func, julia_arr_2D, cunumeric_arr_2D, T)
+    end
+end
+
+function test_unary_reduction_dims(
+    func, julia_arr::AbstractArray{T,N}, cunumeric_arr::NDArray{T,N}
+) where {T,N}
+    allowpromotion(true) do
+        for d in 1:N
+            julia_res = func(julia_arr; dims=d)
+            cunumeric_res = func(cunumeric_arr; dims=d)
+            n = size(julia_arr, d)
+            # Input magnitude for Higham-style absolute floor (see reduction_atol).
+            scale = maximum(abs, julia_arr)
+            allowscalar() do
+                @test safe_compare(
+                    julia_res, cunumeric_res, reduction_atol(T, n, scale), reduction_rtol(T, n)
+                )
+            end
+        end
+
+        # we are testing a multi axis reduction. This will throw a runtime error.
+        # https://github.com/nv-legate/cupynumeric/blob/main/src/cupynumeric/ndarray.cc#L1132
+        if N >= 2
+            @test_throws Exception func(cunumeric_arr, dims=(1, 2))
+        end
+    end
+end
+
+function run_unary_tests(types; include_bool_reductions::Bool=false)
+    @testset verbose = true "Unary Ops w/o Args" begin
+        N = 100 # keep as perfect square
+
+        @testset for T in types
+            allowpromotion(true) do
+                return test_unary_function_set(cuNumeric.floaty_unary_ops_no_args, T, N)
+            end
+
+            allowpromotion(T == Bool) do
+                return test_unary_function_set(cuNumeric.unary_op_map_no_args, T, N)
+            end
+            # Special cases for unary ops that dont use . syntax
+            @testset "- (Negation)" begin
+                arr = my_rand(T, N)
+                arr_cn = @allowscalar NDArray(arr)
+
+                allowscalar() do
+                    allowpromotion(T == Bool) do
+                        T_OUT = T == Bool ? cuNumeric.DEFAULT_INT : T
+                        @test safe_compare(T_OUT.(-arr), -arr_cn, atol(T), rtol(T))
+                    end
+                end
+            end
+
+            # Special cases for complex-related unary ops
+            @testset "Complex Unary Ops (real, imag, conj)" begin
+                if T <: Complex
+                    arr = my_rand(T, N)
+                    arr_cn = NDArray(arr)
+
+                    allowscalar() do
+                        allowpromotion(true) do
+                            @test safe_compare(real(arr), real(arr_cn), atol(T), rtol(T))
+                            @test safe_compare(imag(arr), imag(arr_cn), atol(T), rtol(T))
+                            @test safe_compare(conj(arr), conj(arr_cn), atol(T), rtol(T))
+
+                            @test safe_compare(real.(arr), real.(arr_cn), atol(T), rtol(T))
+                            @test safe_compare(imag.(arr), imag.(arr_cn), atol(T), rtol(T))
+                            @test safe_compare(conj.(arr), conj.(arr_cn), atol(T), rtol(T))
+                        end
+                    end
+                end
+            end
+        end
+    end
+
+    @testset verbose = true "Unary Reductions" begin
+        N = 100
+
+        @testset for T in types
+            julia_arr = my_rand(T, N)
+            cunumeric_arr = @allowscalar NDArray(julia_arr)
+
+            @testset "$(reduction)" for reduction in keys(cuNumeric.unary_reduction_map)
+                # Skip reductions not supported by the cuNumeric backend for complex types
+                if T <: Complex && (
+                    reduction == Base.maximum ||
+                    reduction == Base.minimum ||
+                    reduction == Base.prod
+                )
+                    continue
+                end
+
+                allowpromotion(true) do
+                    cunumeric_res = reduction(cunumeric_arr)
+                    julia_res = reduction(julia_arr)
+
+                    n = length(julia_arr)
+                    scale = maximum(abs, julia_arr)
+                    allowscalar() do
+                        # assumes 0D result
+                        @test isapprox(
+                            julia_res, cunumeric_res[];
+                            atol=reduction_atol(T, n, scale), rtol=reduction_rtol(T, n),
+                        )
+                    end
+                end
+            end
+        end
+
+        if include_bool_reductions
+            # Test things that only work on Booleans
+            julia_bools = rand(Bool, N)
+            allowscalar() do
+                cunumeric_bools = NDArray(julia_bools)
+                @test any(julia_bools) == any(cunumeric_bools)[]
+                @test all(julia_bools) == all(cunumeric_bools)[]
+            end
+        end
+    end
+
+    @testset verbose=true "Unary Reductions with Dims" begin
+        N = 100
+
+        @testset for T in types
+            julia_arr_1D = my_rand(T, N)
+            julia_arr_2D = my_rand(T, isqrt(N), isqrt(N))
+
+            cunumeric_arr_1D = @allowscalar NDArray(julia_arr_1D)
+            cunumeric_arr_2D = @allowscalar NDArray(julia_arr_2D)
+
+            @testset "$(func)" for (func, _) in cuNumeric.unary_reduction_map
+                # Skip reductions not supported by the cuNumeric backend for complex types
+                if T <: Complex && (
+                    func == Base.maximum ||
+                    func == Base.minimum ||
+                    func == Base.prod
+                )
+                    continue
+                end
+
+                ## TODO Int8 min/max along an axis is broken on GPU
+                if cuNumeric.HAS_CUDA && T == Int8 && (func == Base.minimum || func == Base.maximum)
+                    continue
+                end
+
+                test_unary_reduction_dims(func, julia_arr_1D, cunumeric_arr_1D)
+                test_unary_reduction_dims(func, julia_arr_2D, cunumeric_arr_2D)
+            end
+        end
+    end
+end
diff --git a/test/array/unary/uint.jl b/test/array/unary/uint.jl
new file mode 100644
index 000000000..24afcd2c2
--- /dev/null
+++ b/test/array/unary/uint.jl
@@ -0,0 +1,27 @@
+#= 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>
+=#
+
+include(joinpath(@__DIR__, "tests.jl"))
+
+run_unary_tests(
+    filter(
+        T -> T <: Unsigned,
+        Base.uniontypes(cuNumeric.SUPPORTED_ARRAY_TYPES),
+    ),
+)
diff --git a/test/tests/cuda/fusion_compare.jl b/test/defunct/fusion_compare.jl
similarity index 100%
rename from test/tests/cuda/fusion_compare.jl
rename to test/defunct/fusion_compare.jl
diff --git a/test/tests/cuda/fusion_compare_1d.jl b/test/defunct/fusion_compare_1d.jl
similarity index 100%
rename from test/tests/cuda/fusion_compare_1d.jl
rename to test/defunct/fusion_compare_1d.jl
diff --git a/test/tests/cuda/vecadd.jl b/test/defunct/vecadd.jl
similarity index 100%
rename from test/tests/cuda/vecadd.jl
rename to test/defunct/vecadd.jl
diff --git a/test/tests/broadcast_fusion_tests.jl b/test/gpu_only/broadcast_fusion.jl
similarity index 84%
rename from test/tests/broadcast_fusion_tests.jl
rename to test/gpu_only/broadcast_fusion.jl
index 48c16cd9c..5c92911c8 100644
--- a/test/tests/broadcast_fusion_tests.jl
+++ b/test/gpu_only/broadcast_fusion.jl
@@ -31,7 +31,12 @@
 
 _broadcast_fusion_user_add(x, y) = x + y
 
-function test_broadcast_fusion(; T=Float32, N=100, atol=1e-5, rtol=1e-5)
+@testset "Broadcast Fusion" begin
+    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)
@@ -72,126 +77,126 @@ function test_broadcast_fusion(; T=Float32, N=100, atol=1e-5, rtol=1e-5)
     @testset "A + B (two different arrays)" begin
         expected = julia_a .+ julia_b
         result = a .+ b
-        @allowscalar @test cuNumeric.compare(expected, result, atol, rtol)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_compare(expected, result, atol, rtol)
     end
 
     if cuNumeric.FUSE_BROADCAST_EXPRS
@@ -201,7 +206,7 @@ function test_broadcast_fusion(; T=Float32, N=100, atol=1e-5, rtol=1e-5)
             @analyze_lifetimes begin
                 z .= T(2.0) .* _broadcast_fusion_user_add.(a, b)
             end
-            @allowscalar @test cuNumeric.compare(expected, z, atol, rtol)
+            @allowscalar @test safe_compare(expected, z, atol, rtol)
         end
     end
 
@@ -209,14 +214,14 @@ function test_broadcast_fusion(; T=Float32, N=100, atol=1e-5, rtol=1e-5)
     @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)
+        @allowscalar @test safe_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)
+        @allowscalar @test safe_compare(expected, result, atol, rtol)
     end
 end
 
@@ -224,7 +229,10 @@ end
  * Complements `test_broadcast_fusion` with size extremes, 2D same-shape,
  * fusion gating for shape mismatch, 0-d fallback, and dest/input aliasing.
 =#
-function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
+@testset "Broadcast Fusion Edge Cases" begin
+    T=Float32
+    atol=1e-5
+    rtol=1e-5
     s1 = T(2.5)
     s2 = T(1.0)
 
@@ -234,7 +242,7 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
         a = @allowscalar NDArray(ja)
         b = @allowscalar NDArray(jb)
         result = a .+ b .* s1 .- s2
-        @allowscalar @test cuNumeric.compare(ja .+ jb .* s1 .- s2, result, atol, rtol)
+        @allowscalar @test safe_compare(ja .+ jb .* s1 .- s2, result, atol, rtol)
     end
 
     @testset "very small 1D (N=2)" begin
@@ -243,9 +251,9 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
         a = @allowscalar NDArray(ja)
         b = @allowscalar NDArray(jb)
         result = a .+ b
-        @allowscalar @test cuNumeric.compare(ja .+ jb, result, atol, rtol)
+        @allowscalar @test safe_compare(ja .+ jb, result, atol, rtol)
         result = s1 .* a .- b
-        @allowscalar @test cuNumeric.compare(s1 .* ja .- jb, result, atol, rtol)
+        @allowscalar @test safe_compare(s1 .* ja .- jb, result, atol, rtol)
     end
 
     @testset "empty / zero-size 1D" begin
@@ -274,7 +282,7 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
         a = @allowscalar NDArray(ja)
         b = @allowscalar NDArray(jb)
         result = a .+ b .* s1 .- s2
-        @allowscalar @test cuNumeric.compare(ja .+ jb .* s1 .- s2, result, atol, rtol)
+        @allowscalar @test safe_compare(ja .+ jb .* s1 .- s2, result, atol, rtol)
         # Gate: same-shape leaves should be fusible.
         dest = cuNumeric.zeros(T, N)
         bc = Base.Broadcast.instantiate(Base.broadcasted(+, a, b))
@@ -288,9 +296,9 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
         a = @allowscalar NDArray(ja)
         b = @allowscalar NDArray(jb)
         result = a .+ b .* s1
-        @allowscalar @test cuNumeric.compare(ja .+ jb .* s1, result, atol, rtol)
+        @allowscalar @test safe_compare(ja .+ jb .* s1, result, atol, rtol)
         result = a .+ a .* b .- s2
-        @allowscalar @test cuNumeric.compare(ja .+ ja .* jb .- s2, result, atol, rtol)
+        @allowscalar @test safe_compare(ja .+ ja .* jb .- s2, result, atol, rtol)
         dest = cuNumeric.zeros(T, M, N)
         bc = Base.Broadcast.instantiate(Base.broadcasted(+, a, b))
         @test cuNumeric.can_fuse_linear_broadcast(dest, bc)
@@ -311,7 +319,7 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
             a = @allowscalar NDArray(ja)
             b = @allowscalar NDArray(jb)
             result = a .+ b .* s1
-            @allowscalar @test cuNumeric.compare(ja .+ jb .* s1, result, atol, rtol)
+            @allowscalar @test safe_compare(ja .+ jb .* s1, result, atol, rtol)
         end
     end
 
@@ -322,7 +330,7 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
             a = @allowscalar NDArray(ja)
             b = @allowscalar NDArray(jb)
             result = a .+ b .* s1
-            @allowscalar @test cuNumeric.compare(ja .+ jb .* s1, result, atol, rtol)
+            @allowscalar @test safe_compare(ja .+ jb .* s1, result, atol, rtol)
         end
 
         dims = (17, 19, 23)
@@ -332,11 +340,11 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
         b = @allowscalar NDArray(jb)
 
         result = a .+ b .* s1
-        @allowscalar @test cuNumeric.compare(ja .+ jb .* s1, result, atol, rtol)
+        @allowscalar @test safe_compare(ja .+ jb .* s1, result, atol, rtol)
 
         a_alias = @allowscalar NDArray(copy(ja))
         a_alias .= a_alias .* s1 .+ b
-        @allowscalar @test cuNumeric.compare(ja .* s1 .+ jb, a_alias, atol, rtol)
+        @allowscalar @test safe_compare(ja .* s1 .+ jb, a_alias, atol, rtol)
 
         parent_dims = (19, 23, 29)
         jpa = rand(T, parent_dims...)
@@ -354,7 +362,7 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
 
         sliced_result = av .+ bv .* s1
         expected_sliced = jpa[2:18, 3:21, 4:26] .+ jpb[1:17, 2:20, 3:25] .* s1
-        @allowscalar @test cuNumeric.compare(
+        @allowscalar @test safe_compare(
             expected_sliced, sliced_result, atol, rtol
         )
 
@@ -365,7 +373,7 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
         out_view .= av .+ bv .* s1
         expected_out = zeros(T, parent_dims...)
         expected_out[2:18, 3:21, 4:26] = expected_sliced
-        @allowscalar @test cuNumeric.compare(expected_out, out, atol, rtol)
+        @allowscalar @test safe_compare(expected_out, out, atol, rtol)
 
         a4 = @allowscalar NDArray(rand(T, 2, 3, 4, 5))
         @test_throws ArgumentError cuNumeric.nda_get_slice(
@@ -379,9 +387,9 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
         a = @allowscalar NDArray(ja)
         # Literal power uses RefValue{Val} / static-arg lowering.
         result = a .^ 2
-        @allowscalar @test cuNumeric.compare(ja .^ 2, result, atol, rtol)
+        @allowscalar @test safe_compare(ja .^ 2, result, atol, rtol)
         result = s1 .* a .* s2
-        @allowscalar @test cuNumeric.compare(s1 .* ja .* s2, result, atol, rtol)
+        @allowscalar @test safe_compare(s1 .* ja .* s2, result, atol, rtol)
     end
 
     # Gray-Scott-style slice stencils: strided views must fuse correctly via
@@ -399,7 +407,7 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
         expected_x =
             ja[3:end, 2:(end - 1)] .- two .* ja[2:(end - 1), 2:(end - 1)] .+
             ja[1:(end - 2), 2:(end - 1)]
-        @allowscalar @test cuNumeric.compare(expected_x, result_x, atol, rtol)
+        @allowscalar @test safe_compare(expected_x, result_x, atol, rtol)
 
         dest_x = similar(result_x)
         @test !cuNumeric._is_ndarray_slice(dest_x)
@@ -423,7 +431,7 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
         expected_y =
             ja[2:(end - 1), 3:end] .- two .* ja[2:(end - 1), 2:(end - 1)] .+
             ja[2:(end - 1), 1:(end - 2)]
-        @allowscalar @test cuNumeric.compare(expected_y, result_y, atol, rtol)
+        @allowscalar @test safe_compare(expected_y, result_y, atol, rtol)
 
         # Single-op fallback must update the slice's parent.
         out = @allowscalar NDArray(zeros(T, N, N))
@@ -435,7 +443,7 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
         out_interior .= result_x .+ result_y
         expected_out = zeros(T, N, N)
         expected_out[2:(end - 1), 2:(end - 1)] = expected_x .+ expected_y
-        @allowscalar @test cuNumeric.compare(expected_out, out, atol, rtol)
+        @allowscalar @test safe_compare(expected_out, out, atol, rtol)
 
         # Multi-op fused write into a slice.
         out_fused = @allowscalar NDArray(zeros(T, N, N))
@@ -443,7 +451,7 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
         fused_interior .= result_x .+ result_y .* two
         expected_fused = zeros(T, N, N)
         expected_fused[2:(end - 1), 2:(end - 1)] = expected_x .+ expected_y .* two
-        @allowscalar @test cuNumeric.compare(expected_fused, out_fused, atol, rtol)
+        @allowscalar @test safe_compare(expected_fused, out_fused, atol, rtol)
     end
 
     @testset "fused/unfused scalar promotion parity" begin
@@ -456,8 +464,8 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
         bc_i64 = Base.Broadcast.instantiate(Base.broadcasted(*, 2, a))
         @test cuNumeric.can_fuse_linear_broadcast(dest, bc_i64)
         result = 2 .* a
-        @allowscalar @test cuNumeric.compare(T(2) .* ja, result, atol, rtol)
-        @allowscalar @test cuNumeric.compare(
+        @allowscalar @test safe_compare(T(2) .* ja, result, atol, rtol)
+        @allowscalar @test safe_compare(
             result,
             cuNumeric.unravel_broadcast_tree(bc_i64),
             atol,
@@ -538,11 +546,11 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
         a = @allowscalar NDArray(copy(ja))
         b = @allowscalar NDArray(jb)
         a .+= b
-        @allowscalar @test cuNumeric.compare(ja .+ jb, a, atol, rtol)
+        @allowscalar @test safe_compare(ja .+ jb, a, atol, rtol)
 
         a2 = @allowscalar NDArray(copy(ja))
         a2 .= a2 .* s1 .+ b
-        @allowscalar @test cuNumeric.compare(ja .* s1 .+ jb, a2, atol, rtol)
+        @allowscalar @test safe_compare(ja .* s1 .+ jb, a2, atol, rtol)
 
         M, N2 = 37, 65
         j2a = rand(T, M, N2)
@@ -550,7 +558,7 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
         a2d = @allowscalar NDArray(copy(j2a))
         b2d = @allowscalar NDArray(j2b)
         a2d .= a2d .* s1 .+ b2d
-        @allowscalar @test cuNumeric.compare(j2a .* s1 .+ j2b, a2d, atol, rtol)
+        @allowscalar @test safe_compare(j2a .* s1 .+ j2b, a2d, atol, rtol)
     end
 
     @testset "cross-statement fusion into a slice" begin
@@ -564,7 +572,7 @@ function test_broadcast_fusion_edge_cases(; T=Float32, atol=1e-5, rtol=1e-5)
         end
         expected = zeros(T, N + 2, N + 2)
         expected[2:(end - 1), 2:(end - 1)] = ja .* s1 .+ s2
-        @allowscalar @test cuNumeric.compare(expected, out, atol, rtol)
+        @allowscalar @test safe_compare(expected, out, atol, rtol)
     end
 end
 
@@ -575,7 +583,10 @@ end
  * With `FUSE_BROADCAST_MIN_OPS > 1`, single-op exprs are unfused — tests
  * should set min ops to 1 (LocalPreferences / ENV) to exercise the cache.
 =#
-function test_broadcast_fusion_ptx_cache(; T=Float32, N=64)
+@testset "Broadcast Fusion PTX Cache" begin
+    T=Float32
+    N=64
+
     if !(cuNumeric.FUSE_BROADCAST_EXPRS && cuNumeric.HAS_CUDA)
         @info "Skipping PTX cache tests (need FUSE_BROADCAST_EXPRS && HAS_CUDA)"
         return nothing
diff --git a/test/tests/hdf5.jl b/test/io/hdf5.jl
similarity index 83%
rename from test/tests/hdf5.jl
rename to test/io/hdf5.jl
index 159152b01..a35cfb59b 100644
--- a/test/tests/hdf5.jl
+++ b/test/io/hdf5.jl
@@ -36,13 +36,21 @@ function test_hdf5_roundtrip(::Type{T}, shape::Tuple) where {T}
         output = cuNumeric.h5read(path, dataset; layout=:row)
         @test eltype(output) == T
         @test size(output) == shape
-        @allowscalar @test cuNumeric.compare(expected, output, 0, 0)
+        @allowscalar @test safe_compare(expected, output, 0, 0)
 
         if length(shape) > 1
             col_output = cuNumeric.h5read(path, dataset; layout=:col)
             col_expected = permutedims(expected, reverse(1:length(shape)))
             @test size(col_output) == reverse(shape)
-            @allowscalar @test cuNumeric.compare(col_expected, col_output, 0, 0)
+            @allowscalar @test safe_compare(col_expected, col_output, 0, 0)
+        end
+    end
+end
+
+@testset "HDF5" begin
+    for T in (Float32, Float64, Int32, Int64)
+        @testset "$T $shape" for shape in ((7,), (3, 4), (2, 3, 4))
+            test_hdf5_roundtrip(T, shape)
         end
     end
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index df280a3e8..279554fcd 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,507 +1,58 @@
-#= 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>
-=#
-
-using Test
-using LinearAlgebra
-using Random
-import Random: rand
-
-const VERBOSE = get(ENV, "VERBOSE", "1") != "0"
-const run_gpu_tests = (get(ENV, "GPUTESTS", "1") != "0") && (get(ENV, "NO_CUDA", "OFF") != "ON")
-@info "Run GPU Tests: $(run_gpu_tests)"
-
-if run_gpu_tests
-    using CUDA
-    import CUDA: i32
-    VERBOSE && println(CUDA.versioninfo())
-end
-
-if run_gpu_tests && !CUDA.functional()
-    error(
-        "You asked for CUDA tests, but they are disabled because no functional CUDA device was detected."
-    )
-end
-
 using cuNumeric
-VERBOSE && cuNumeric.versioninfo()
-
-@info "Broadcast fusion: FUSE_BROADCAST_EXPRS=$(cuNumeric.FUSE_BROADCAST_EXPRS) FUSE_BROADCAST_MIN_OPS=$(cuNumeric.FUSE_BROADCAST_MIN_OPS)"
-
-# TODO
-# After loading cuNumeric, we should verify that the Legate config has set a GPU device
-# Right now, if you have a gpu device, but your LEGATE_CONFIG is cpu only,
-# @cuda_task will not be defined and tests will fail confusingly.
-# We should error out more gracefully in this situation.
-
-include("tests/util.jl")
-include("tests/axpy.jl")
-include("tests/axpy_advanced.jl")
-include("tests/elementwise.jl")
-include("tests/slicing.jl")
-include("tests/gemm.jl")
-include("tests/unary_tests.jl")
-include("tests/binary_tests.jl")
-include("tests/scoping.jl")
-include("tests/scoping-advanced.jl")
-include("tests/hdf5.jl")
-include("tests/broadcast_fusion_tests.jl")
-
-@testset "Array Accessors" begin
-    for T in Base.uniontypes(cuNumeric.SUPPORTED_ARRAY_TYPES)
-        value = rand(T)
-        arr = cuNumeric.zeros(T, 2, 2)
-
-        allowscalar() do
-            arr[1, 2] = value
-            @test arr[1, 2] == value
-        end
-    end
-end
-
-@testset verbose = true "AXPY" begin
-    N = 100
-    @testset verbose = true for T in Base.uniontypes(cuNumeric.SUPPORTED_FLOAT_TYPES)
-        @testset "basic" axpy_basic(T, N)
-        @testset "advanced" axpy_advanced(T, N)
-    end
-end
-
-@testset verbose = true "Operators" begin
-    @testset elementwise()
-end
-
-@testset verbose = true "Linear Algebra Tests" begin
-    include("tests/linalg.jl")
-end
-
-@testset verbose = true "GEMM" begin
-    N = 50
-    M = 25
-    @testset verbose = true for T in Base.uniontypes(cuNumeric.SUPPORTED_NUMERIC_TYPES)
-        gemm(N, M, T, rtol(T))
-    end
-end
-
-#* TODO TEST VARIANT OVER DIMS
-@testset verbose = true "Unary Ops w/o Args" begin
-    N = 100 # keep as perfect square
-
-    @testset for T in Base.uniontypes(cuNumeric.SUPPORTED_ARRAY_TYPES)
-        allowpromotion(true) do
-            return test_unary_function_set(cuNumeric.floaty_unary_ops_no_args, T, N)
-        end
-
-        allowpromotion(T == Bool) do
-            return test_unary_function_set(cuNumeric.unary_op_map_no_args, T, N)
-        end
-        # Special cases for unary ops that dont use . syntax
-        @testset "- (Negation)" begin
-            arr = my_rand(T, N)
-            arr_cn = @allowscalar NDArray(arr)
-
-            allowscalar() do
-                allowpromotion(T == Bool) do
-                    T_OUT = T == Bool ? cuNumeric.DEFAULT_INT : T
-                    @test cuNumeric.compare(T_OUT.(-arr), -arr_cn, atol(T), rtol(T))
-                end
-            end
-        end
-
-        # Special cases for complex-related unary ops
-        @testset "Complex Unary Ops (real, imag, conj)" begin
-            if T <: Complex
-                arr = my_rand(T, N)
-                arr_cn = NDArray(arr)
-
-                allowscalar() do
-                    allowpromotion(true) do
-                        @test cuNumeric.compare(real(arr), real(arr_cn), atol(T), rtol(T))
-                        @test cuNumeric.compare(imag(arr), imag(arr_cn), atol(T), rtol(T))
-                        @test cuNumeric.compare(conj(arr), conj(arr_cn), atol(T), rtol(T))
-
-                        @test cuNumeric.compare(real.(arr), real.(arr_cn), atol(T), rtol(T))
-                        @test cuNumeric.compare(imag.(arr), imag.(arr_cn), atol(T), rtol(T))
-                        @test cuNumeric.compare(conj.(arr), conj.(arr_cn), atol(T), rtol(T))
-                    end
-                end
-            end
-        end
-    end
-end
-
-@testset verbose = true "Unary Reductions" begin
-    N = 100
-
-    @testset for T in Base.uniontypes(cuNumeric.SUPPORTED_ARRAY_TYPES)
-        julia_arr = my_rand(T, N)
-        cunumeric_arr = @allowscalar NDArray(julia_arr)
-
-        @testset "$(reduction)" for reduction in keys(cuNumeric.unary_reduction_map)
-            # Skip reductions not supported by the cuNumeric backend for complex types
-            if T <: Complex && (
-                reduction == Base.maximum ||
-                reduction == Base.minimum ||
-                reduction == Base.prod
-            )
-                continue
-            end
-
-            allowpromotion(true) do
-                cunumeric_res = reduction(cunumeric_arr)
-                julia_res = reduction(julia_arr)
-
-                n = length(julia_arr)
-                allowscalar() do
-                    # assumes 0D result
-                    @test isapprox(
-                        julia_res, cunumeric_res[];
-                        atol=reduction_atol(T, n), rtol=reduction_rtol(T, n),
-                    )
-                end
-            end
-        end
-    end
-
-    # Test things that only work on Booleans
-    julia_bools = rand(Bool, N)
-    allowscalar() do
-        cunumeric_bools = NDArray(julia_bools)
-        @test any(julia_bools) == any(cunumeric_bools)[]
-        @test all(julia_bools) == all(cunumeric_bools)[]
-    end
-end
-
-@testset verbose=true "Unary Reductions with Dims" begin
-    N = 100
-
-    @testset for T in Base.uniontypes(cuNumeric.SUPPORTED_ARRAY_TYPES)
-        julia_arr_1D = my_rand(T, N)
-        julia_arr_2D = my_rand(T, isqrt(N), isqrt(N))
-
-        cunumeric_arr_1D = @allowscalar NDArray(julia_arr_1D)
-        cunumeric_arr_2D = @allowscalar NDArray(julia_arr_2D)
-
-        @testset "$(func)" for (func, _) in cuNumeric.unary_reduction_map
-            # Skip reductions not supported by the cuNumeric backend for complex types
-            if T <: Complex && (
-                func == Base.maximum ||
-                func == Base.minimum ||
-                func == Base.prod
-            )
-                continue
-            end
-
-            ## TODO Int8 min/max along an axis is broken on GPU
-            if cuNumeric.HAS_CUDA && T == Int8 && (func == Base.minimum || func == Base.maximum)
-                continue
-            end
-
-            test_unary_reduction_dims(func, julia_arr_1D, cunumeric_arr_1D)
-            test_unary_reduction_dims(func, julia_arr_2D, cunumeric_arr_2D)
-        end
-    end
-end
+using CUDA: CUDA
+using ParallelTestRunner
+using Pkg
+using InteractiveUtils: versioninfo
 
-@testset verbose = true "Binary Ops" begin
-    N = 100
+run_gpu_tests = CUDA.functional()
 
-    @testset for T in Base.uniontypes(cuNumeric.SUPPORTED_ARRAY_TYPES)
-        allowpromotion(true) do
-            test_binary_function_set(cuNumeric.floaty_binary_op_map, T, N)
-            return test_binary_function_set(cuNumeric.binary_op_map, T, N)
-        end
+@info "Julia information:\n" * sprint(io -> versioninfo(io))
+run_gpu_tests && @info "CUDA information:\n" * sprint(io -> CUDA.versioninfo(io))
+@info "cuNumeric information:\n" * sprint(io -> cuNumeric.versioninfo(io))
 
-        arr_jl = my_rand(T, N)
-        arr_jl2 = my_rand(T, N)
-        arr_cn = @allowscalar NDArray(arr_jl)
-        arr_cn2 = @allowscalar NDArray(arr_jl2)
+# Forcibly precompile the current environment in parallel: Pkg sometimes ignores
+# dependencies pointed through via `[sources]`
+Pkg.precompile()
 
-        # lcm/gcd require specific handling for integers and avoid overflow
-        if T <: cuNumeric.SUPPORTED_INT_TYPES && T != Bool
-            range_limit = (T == Int8 || T == UInt8) ? 10 : 100
-            arr_jl_small = my_rand(T, N; L=1, R=range_limit)
-            arr_jl2_small = my_rand(T, N; L=1, R=range_limit)
-            arr_cn_small = @allowscalar NDArray(arr_jl_small)
-            arr_cn2_small = @allowscalar NDArray(arr_jl2_small)
-
-            allowscalar() do
-                @test cuNumeric.compare(
-                    lcm.(arr_jl_small, arr_jl2_small), lcm.(arr_cn_small, arr_cn2_small), atol(T),
-                    rtol(T),
-                )
-                @test cuNumeric.compare(
-                    gcd.(arr_jl_small, arr_jl2_small), gcd.(arr_cn_small, arr_cn2_small), atol(T),
-                    rtol(T),
-                )
-            end
-        end
-
-        allowscalar() do
-            @test unwrap(arr_cn == arr_cn)
-            @test !unwrap(arr_cn == arr_cn2)
-            @test unwrap(arr_cn != arr_cn2)
-            @test !unwrap(arr_cn != arr_cn)
-            @test unwrap(all(arr_cn .== arr_cn))
-        end
-    end
-end
-
-@testset "Type and Shape Promotion" begin
-    N = 100
-    cunumeric_arr1 = cuNumeric.zeros(Float64, N)
-    cunumeric_arr3 = cuNumeric.zeros(Float32, N)
-    cunumeric_int64 = cuNumeric.zeros(Int64, N)
-    cunumeric_int32 = cuNumeric.zeros(Int32, N)
-    cunumeric_arr5 = cuNumeric.zeros(Float64, N - 1, N - 1)
-    @test_throws "Implicit promotion" cunumeric_arr3 .+ cunumeric_arr1
-    @test_throws "Implicit promotion" map(+, cunumeric_arr3, cunumeric_arr1)
-    @test_throws DimensionMismatch cunumeric_arr1 .+ cunumeric_arr5
-    @test_throws DimensionMismatch cunumeric_arr1 ./ cunumeric_arr5
-
-    allowscalar() do
-        @test cuNumeric.compare(
-            cunumeric_arr1, cunumeric_int64 .+ cunumeric_arr1, atol(Float64), rtol(Float64)
-        )
-        r1 = @allowpromotion cunumeric_arr3 .+ cunumeric_arr1
-        r2 = @allowpromotion map(+, cunumeric_arr3, cunumeric_arr1)
-        @test cuNumeric.compare(r1, r2, atol(Float64), rtol(Float64))
-    end
-end
-
-@testset "Copy-To" begin
-    a = cuNumeric.zeros(2, 2)
-    b = cuNumeric.ones(2, 2)
-    copyto!(a, b)
-    @test is_same(a, b)
-end
-
-#TODO LOOP BINARY OPS WITH SCALARS
-@testset verbose = true "Scalars" begin
-    N = 10
-
-    for T in Base.uniontypes(cuNumeric.SUPPORTED_NUMERIC_TYPES)
-        julia_arr = my_rand(T, N)
-        julia_arr_2D = my_rand(T, N, N)
-
-        s = Random.rand(T)
-
-        allowscalar() do
-            cunumeric_arr = NDArray(julia_arr)
-            cunumeric_arr_2D = NDArray(julia_arr_2D)
-            allowpromotion(true) do
-                for cn_arr in (cunumeric_arr, cunumeric_arr_2D)
-                    @test cuNumeric.compare(s * julia_arr, s * cunumeric_arr, atol(T), rtol(T))
-                    @test cuNumeric.compare(julia_arr * s, cunumeric_arr * s, atol(T), rtol(T))
-                    @test cuNumeric.compare(s .* julia_arr, s .* cunumeric_arr, atol(T), rtol(T))
-                    @test cuNumeric.compare(julia_arr .* s, cunumeric_arr .* s, atol(T), rtol(T))
-                    @test cuNumeric.compare(s .+ julia_arr, s .+ cunumeric_arr, atol(T), rtol(T))
-                    @test cuNumeric.compare(julia_arr .+ s, cunumeric_arr .+ s, atol(T), rtol(T))
-                    @test cuNumeric.compare(s .- julia_arr, s .- cunumeric_arr, atol(T), rtol(T))
-                    @test cuNumeric.compare(julia_arr .- s, cunumeric_arr .- s, atol(T), rtol(T))
-                    @test cuNumeric.compare(s ./ julia_arr, s ./ cunumeric_arr, atol(T), rtol(T))
-                    @test cuNumeric.compare(
-                        s .* julia_arr .+ s, s .* cunumeric_arr .+ s, atol(T), rtol(T)
-                    )
-                    @test s + s ≈ (NDArray(s) + NDArray(s))[] # atol=atol(T) r_tol=atol(T)
-                    # @test s * s ≈ (NDArray(s) * NDArray(s))[] # atol=atol(T) r_tol=atol(T)
-                end
-            end
-        end
-    end
-
-    # Boolean things
-    allowpromotion() do
-        allowscalar() do
-            julia_arr = [-3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
-            cunumeric_arr = NDArray(julia_arr)
-
-            @test cuNumeric.compare(true * julia_arr, true * cunumeric_arr, atol(Int), rtol(Int))
-            @test cuNumeric.compare(false * julia_arr, false * cunumeric_arr, atol(Int), rtol(Int))
-            @test cuNumeric.compare(true .* julia_arr, true .* cunumeric_arr, atol(Int), rtol(Int))
-            @test cuNumeric.compare(
-                false .* julia_arr, false .* cunumeric_arr, atol(Int), rtol(Int)
-            )
-
-            julia_arr = [true, false, true, false, false, true, true]
-            cunumeric_arr = NDArray(julia_arr)
-
-            @test cuNumeric.compare(4 * julia_arr, 4 * cunumeric_arr, atol(Int), rtol(Int))
-            @test cuNumeric.compare(4 .* julia_arr, 4 .* cunumeric_arr, atol(Int), rtol(Int))
-            @test cuNumeric.compare(4 .+ julia_arr, 4 .+ cunumeric_arr, atol(Int), rtol(Int))
-            @test cuNumeric.compare(
-                julia_arr ./ 3, cunumeric_arr ./ 3, atol(Float64), rtol(Float64)
-            )
-        end
+cuda_init = if run_gpu_tests
+    quote
+        using CUDA
+        import CUDA: i32
     end
+else
+    :()
 end
 
-@testset verbose = true "Powers" begin
-    N = 9
-
-    get_pwrs(::Type{I}) where {I<:Integer} = I.([-10, -5, -2, -1, 0, 1, 2, 5, 10])
-    get_pwrs(::Type{F}) where {F<:AbstractFloat} = F.([-3.141, -2, -1, 0, 1, 2, 3.2, 4.41, 6.233])
-    get_pwrs(::Type{Bool}) = [true, false, true, false, false, true, false, true, true]
-
-    # TYPES = Base.uniontypes(cuNumeric.SUPPORTED_ARRAY_TYPES)
-    TYPES = Base.uniontypes(cuNumeric.SUPPORTED_FLOAT_TYPES)
-
-    @testset "$(BT) ^ $(PT)" for (BT, PT) in Iterators.product(TYPES, TYPES)
-        base_jl = my_rand(BT, N)
-
-        if BT <: Union{Bool,Int32} && PT == Int32
-            # julia doesnt like Int32 powers
-            pwrs = Float32.(get_pwrs(PT))
-        elseif BT <: Union{Bool,Int32,Int64} && PT <: Union{Int32,Int64}
-            # julia doesnt like Int64 powers on bool
-            pwrs = Float64.(get_pwrs(PT))
-        elseif (PT <: AbstractFloat) && (BT <: AbstractFloat || BT <: Signed)
-            # Things like -387 ^ 3.2 will be Complex and error
-            pwrs = get_pwrs(PT)
-            base_jl = abs.(base_jl)
-        else
-            pwrs = get_pwrs(PT)
-        end
-
-        base_cn = @allowscalar NDArray(base_jl)
-        pwrs_cn = @allowscalar NDArray(pwrs)
+const init_code = quote
+    using LinearAlgebra
+    using Random
+    import Random: rand
 
-        # we deviate a bit form Julia here
-        if (PT <: Union{Int32,Int64}) && (BT <: Union{Bool,Int32,Int64})
-            T_OUT = cuNumeric.__my_promote_type(typeof(^), BT, PT)
-        else
-            T_OUT = Base.promote_op(Base.:(^), BT, PT)
-        end
+    ENV["LEGATE_SKIP_RUNTIME"] = "false"
+    using cuNumeric
 
-        TEST_BROKEN = (BT <: Union{Int32,Int64} && PT == Bool)
+    $cuda_init
 
-        allowpromotion(true) do
-            allowscalar() do
-                # Power is array
-                @test cuNumeric.compare(
-                    base_jl .^ pwrs, base_cn .^ pwrs_cn, atol(T_OUT), rtol(T_OUT)
-                ) skip=TEST_BROKEN
-
-                # Power is scalar
-                for p in pwrs
-                    @test cuNumeric.compare(base_jl .^ p, base_cn .^ p, atol(T_OUT), rtol(T_OUT))
-                end
-            end
-        end
-    end
-
-    @testset verbose = true "Reciprocal" begin
-        @testset for T in TYPES
-            arr_jl = Random.rand(T, N)
-            arr_cn = @allowscalar NDArray(arr_jl)
-
-            # Differ from Julia here
-            T_OUT = cuNumeric.__recip_type(T)
-
-            # Cast julia result to whatever we do
-            res_jl = T_OUT.(arr_jl .^ -1)
-            allowpromotion(true) do
-                res_cn = arr_cn .^ -1
-                res_cn2 = inv.(arr_cn)
-                allowscalar() do
-                    @test cuNumeric.compare(res_jl, res_cn, atol(T_OUT), rtol(T_OUT))
-                    @test cuNumeric.compare(res_jl, res_cn2, atol(T_OUT), rtol(T_OUT))
-                end
-            end
-        end
-    end
-
-    @testset verbose = true "Square" begin
-        @testset for T in TYPES
-            arr_jl = Random.rand(T, N)
-            arr_cn = @allowscalar NDArray(arr_jl)
-
-            T_OUT = Base.promote_op(Base.:(^), T, Int64)
-            res_jl = arr_jl .^ 2
-            res_cn = arr_cn .^ 2
-
-            allowpromotion(true) do
-                allowscalar() do
-                    @test cuNumeric.compare(res_jl, res_cn, atol(T_OUT), rtol(T_OUT))
-                end
-            end
-        end
-    end
+    include("util.jl")
 end
 
-@testset verbose = true "Slicing Tests" begin
-    N = 100
-    @testset for T in Base.uniontypes(cuNumeric.SUPPORTED_FLOAT_TYPES)
-        slicing(T, N)
-    end
-end
+# Find all tests, remove ones that are not relevant to the current configuration
 
-@testset verbose = true "Type Stability" begin
-    include("tests/stability.jl")
-end
+testsuite = find_tests(@__DIR__)
+delete!(testsuite, "util")
+delete!(testsuite, "array/unary/tests")
+delete!(testsuite, "array/binary/tests")
 
-@testset verbose = true "NDArray Lifetime and Zero-Copy" begin
-    include("tests/lifetime.jl")
+if !run_gpu_tests
+    @warn "CUDA GPU not available, skipping GPU-only tests"
+    filter!(test -> !startswith(first(test), "gpu_only/"), testsuite)
 end
 
-@testset verbose = true "Scoping" begin
-    N = 100
-
-    @testset "Rewrite pipeline" test_scoping_rewrite_pipeline()
-
-    @testset verbose = true for T in Base.uniontypes(cuNumeric.SUPPORTED_FLOAT_TYPES)
-        allowscalar() do
-            results = run_all_ops(T, N)
-            for (name, (c_base, c_scoped)) in results
-                @test cuNumeric.compare(c_base, c_scoped, atol(T), rtol(T))
-            end
-
-            u_rand = cuNumeric.rand(T, (15, 15))
-            v_rand = cuNumeric.rand(T, (15, 15))
-
-            u, v = gray_scott_base(T, N, u_rand, v_rand)
-            u_scoped, v_scoped = gray_scott(T, N, u_rand, v_rand)
-
-            @test cuNumeric.compare(u, u_scoped, atol(T) * N, rtol(T) * 10)
-        end
-    end
+if !run_gpu_tests || !cuNumeric.FUSE_BROADCAST_EXPRS
+    @warn "Broadcast fusion is disabled, skipping fusion tests"
+    filter!(test -> !startswith(first(test), "gpu_only/broadcast_fusion"), testsuite)
 end
 
-@testset verbose = true "HDF5" begin
-    for T in (Float32, Float64, Int32, Int64)
-        @testset "$T $shape" for shape in ((7,), (3, 4), (2, 3, 4))
-            test_hdf5_roundtrip(T, shape)
-        end
-    end
-end
-
-if run_gpu_tests
-    @testset verbose = true "Broadcast Fusion" begin
-        test_broadcast_fusion()
-        test_broadcast_fusion_edge_cases()
-        test_broadcast_fusion_ptx_cache()
-    end
+filter!(test -> !startswith(first(test), "defunct/"), testsuite)
 
-    # @testset verbose = true "CUDA Tests" begin
-    #     cuda_unaryop(rtol(Float32))
-    #     cuda_binaryop(rtol(Float32))
-    # end
-    @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
+runtests(cuNumeric, ARGS; testsuite, init_code)
diff --git a/test/tests/axpy.jl b/test/tests/axpy.jl
deleted file mode 100644
index 10135fc4f..000000000
--- a/test/tests/axpy.jl
+++ /dev/null
@@ -1,53 +0,0 @@
-#= 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>
-=#
-
-#= Purpose of test: daxpy
-    -- Focused on double 2 dimenional. Does not test other types or dims.
-    -- NDArray intialization
-    -- NDArray writing and reading scalar indexing
-    --          shows both [i, j] and [(i, j)] working
-    -- NDArray addition and multiplication
-=#
-function axpy_basic(T, N)
-    α = T(56.6)
-    dims = (N, N)
-
-    # Base julia arrays
-    x_cpu = rand(T, dims);
-    y_cpu = rand(T, dims);
-
-    # cunumeric arrays
-    x = cuNumeric.zeros(T, dims)
-    y = cuNumeric.zeros(T, dims)
-
-    # Initialize NDArrays with same random values as julia arrays
-    @allowscalar for i in 1:N
-        for j in 1:N
-            x[i, j] = x_cpu[i, j]
-            y[i, j] = y_cpu[i, j]
-        end
-    end
-
-    result = α .* x .+ y
-    result_cpu = α .* x_cpu .+ y_cpu
-    allowscalar() do
-        @test cuNumeric.compare(result, result_cpu, atol(T), rtol(T))
-        @test cuNumeric.compare(result_cpu, result, atol(T), rtol(T))
-    end
-end
diff --git a/test/tests/gemm.jl b/test/tests/gemm.jl
deleted file mode 100644
index 646eb0e08..000000000
--- a/test/tests/gemm.jl
+++ /dev/null
@@ -1,85 +0,0 @@
-#= 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>
-=#
-
-function gemm(N, M, T, max_diff)
-    if T == Bool
-        a = cuNumeric.trues(5, 5)
-        b = cuNumeric.as_type(cuNumeric.trues(5, 5), Float32)
-        c = cuNumeric.as_type(cuNumeric.trues(5, 5), Float64)
-        @test_throws ArgumentError a * a # Bool * Bool not supported
-        @allowpromotion d = a * b
-        @allowpromotion e = a * c
-        @test @allowscalar cuNumeric.compare(5 * ones(Float32, 5, 5), d, 0.0, max_diff)
-        @test @allowscalar cuNumeric.compare(5 * ones(Float64, 5, 5), e, 0.0, max_diff)
-        return nothing
-    end
-
-    if T <: Integer
-        a = cuNumeric.ones(Int32, 5, 5)
-        a_jl = ones(Int32, 5, 5)
-        b = cuNumeric.ones(Float32, 5, 5)
-        b_jl = ones(Float32, 5, 5)
-        @test_throws ArgumentError a * a
-        @test @allowscalar cuNumeric.compare(a_jl * b_jl, a * b, 0.0, max_diff)
-        return nothing
-    end
-
-    dims_to_test = [(N, N), (N, M), (M, N)]
-
-    @testset for dims in dims_to_test
-        # Base julia arrays
-        A_cpu = rand(T, dims[1], dims[2]);
-        B_cpu = rand(T, dims[2], dims[1]);
-        C_out_cpu = zeros(T, dims[1], dims[1])
-
-        # cunumeric arrays
-        A = cuNumeric.NDArray(A_cpu)
-        B = cuNumeric.NDArray(B_cpu)
-        C_out = cuNumeric.zeros(T, dims[1], dims[1])
-
-        # Julia result
-        C_cpu = A_cpu * B_cpu
-        LinearAlgebra.mul!(C_out_cpu, A_cpu, B_cpu)
-
-        @test C_cpu == C_out_cpu # really just making sure test is written right...
-
-        A = cuNumeric.as_type(A, T)
-        B = cuNumeric.as_type(B, T)
-        C = cuNumeric.zeros(T, N, N)
-
-        C = A * B
-        LinearAlgebra.mul!(C_out, A, B)
-
-        allowscalar() do
-            @test isapprox(C, C_cpu, rtol=max_diff)
-            @test isapprox(C, C_out, rtol=max_diff)
-
-            if T != Float64
-                C_wider = cuNumeric.zeros(Float64, dims[1], dims[1])
-                @test_throws "Implicit promotion" LinearAlgebra.mul!(C_wider, A, B)
-            end
-        end
-
-        # Integer output with FP input
-        if !(T <: Integer)
-            bad = cuNumeric.zeros(Int, dims[1], dims[1])
-            @test_throws ArgumentError mul!(bad, A, B)
-        end
-    end
-end
diff --git a/test/tests/scoping-advanced.jl b/test/tests/scoping-advanced.jl
deleted file mode 100644
index 87a494195..000000000
--- a/test/tests/scoping-advanced.jl
+++ /dev/null
@@ -1,191 +0,0 @@
-#= 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>
-=#
-
-struct ParamsGS{T<:AbstractFloat}
-    dx::T
-    dt::T
-    c_u::T
-    c_v::T
-    f::T
-    k::T
-
-    # Constructor with default values
-    function ParamsGS{T}(
-        dx::T=one(T), c_u::T=one(T), c_v::T=T(0.3), f::T=T(0.03), k::T=T(0.06)
-    ) where {T<:AbstractFloat}
-        new(dx, dx/5, c_u, c_v, f, k)
-    end
-end
-
-function step(u, v, u_new, v_new, args::ParamsGS)
-    @analyze_lifetimes begin
-        # calculate F_u and F_v functions
-        # currently we don't have NDArray^x working yet.
-        F_u = (
-            (
-                -u[2:(end - 1), 2:(end - 1)] .*
-                (v[2:(end - 1), 2:(end - 1)] .* v[2:(end - 1), 2:(end - 1)])
-            ) + args.f*(1 .- u[2:(end - 1), 2:(end - 1)])
-        )
-        F_v = (
-            (
-                u[2:(end - 1), 2:(end - 1)] .*
-                (v[2:(end - 1), 2:(end - 1)] .* v[2:(end - 1), 2:(end - 1)])
-            ) - (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)]
-            ) ./ args.dx^2 +
-            (
-                u[2:(end - 1), 3:end] - 2*u[2:(end - 1), 2:(end - 1)] +
-                u[2:(end - 1), 1:(end - 2)]
-            ) ./ args.dx^2
-        )
-        v_lap = (
-            (
-                v[3:end, 2:(end - 1)] - 2*v[2:(end - 1), 2:(end - 1)] +
-                v[1:(end - 2), 2:(end - 1)]
-            ) ./ args.dx^2 +
-            (
-                v[2:(end - 1), 3:end] - 2*v[2:(end - 1), 2:(end - 1)] +
-                v[2:(end - 1), 1:(end - 2)]
-            ) ./ args.dx^2
-        )
-
-        # # 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
-        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
-end
-
-# same as above but without @analyze_lifetimes macro
-function step_base(u, v, u_new, v_new, args::ParamsGS)
-    # calculate F_u and F_v functions
-    # currently we don't have NDArray^x working yet.
-    F_u = (
-        (
-            -u[2:(end - 1), 2:(end - 1)] .*
-            (v[2:(end - 1), 2:(end - 1)] .* v[2:(end - 1), 2:(end - 1)])
-        ) + args.f*(1 .- u[2:(end - 1), 2:(end - 1)])
-    )
-    F_v = (
-        (
-            u[2:(end - 1), 2:(end - 1)] .*
-            (v[2:(end - 1), 2:(end - 1)] .* v[2:(end - 1), 2:(end - 1)])
-        ) - (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)]
-        ) ./ args.dx^2 +
-        (
-            u[2:(end - 1), 3:end] - 2*u[2:(end - 1), 2:(end - 1)] +
-            u[2:(end - 1), 1:(end - 2)]
-        ) ./ args.dx^2
-    )
-    v_lap = (
-        (
-            v[3:end, 2:(end - 1)] - 2*v[2:(end - 1), 2:(end - 1)] +
-            v[1:(end - 2), 2:(end - 1)]
-        ) ./ args.dx^2 +
-        (
-            v[2:(end - 1), 3:end] - 2*v[2:(end - 1), 2:(end - 1)] +
-            v[2:(end - 1), 1:(end - 2)]
-        ) ./ args.dx^2
-    )
-
-    # # 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
-    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 gray_scott(FT, n_steps, u_rand, v_rand)
-    N = 100
-    dims = (N, N)
-    args = ParamsGS{FT}()
-    u = cuNumeric.ones(FT, dims)
-    v = cuNumeric.zeros(FT, dims)
-    u_new = cuNumeric.zeros(FT, dims)
-    v_new = cuNumeric.zeros(FT, dims)
-
-    u[1:15, 1:15] = u_rand
-    v[1:15, 1:15] = v_rand
-
-    cuNumeric.disable_gc!(; verbose=false)
-    for n in 1:n_steps
-        step(u, v, u_new, v_new, args)
-        u, u_new = u_new, u
-        v, v_new = v_new, v
-    end
-
-    return u, v
-end
-
-function gray_scott_base(FT, n_steps, u_rand, v_rand)
-    N = 100
-    dims = (N, N)
-    args = ParamsGS{FT}()
-    u = cuNumeric.ones(FT, dims)
-    v = cuNumeric.zeros(FT, dims)
-    u_new = cuNumeric.zeros(FT, dims)
-    v_new = cuNumeric.zeros(FT, dims)
-
-    u[1:15, 1:15] = u_rand
-    v[1:15, 1:15] = v_rand
-
-    cuNumeric.init_gc!()
-    for n in 1:n_steps
-        step_base(u, v, u_new, v_new, args)
-        u, u_new = u_new, u
-        v, v_new = v_new, v
-    end
-
-    return u, v
-end
diff --git a/test/tests/unary_tests.jl b/test/tests/unary_tests.jl
deleted file mode 100644
index 6bc87a582..000000000
--- a/test/tests/unary_tests.jl
+++ /dev/null
@@ -1,113 +0,0 @@
-#= 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>
-=#
-
-# Map functions to their required domains
-const SPECIAL_DOMAINS = Dict(
-    Base.acosh => :greater_than_one,
-    Base.log => :positive,
-    Base.log10 => :positive,
-    Base.log2 => :positive,
-    Base.log1p => :positive, # technicaly anything > -1
-    Base.sqrt => :positive,
-)
-
-function test_unary_operation(func, julia_arr, cunumeric_arr, T)
-    T_OUT = Base.promote_op(func, T)
-
-    # Pre-allocate output arrays
-    cunumeric_in_place = cuNumeric.zeros(T_OUT, size(julia_arr)...)
-
-    # Compute results using different methods
-    julia_res = func.(julia_arr)
-
-    cunumeric_res = func.(cunumeric_arr)
-    cunumeric_in_place .= func.(cunumeric_arr)
-    cunumeric_res2 = map(func, cunumeric_arr)
-
-    allowscalar() do
-        @test cuNumeric.compare(julia_res, cunumeric_in_place, atol(T_OUT), rtol(T_OUT))
-        @test cuNumeric.compare(julia_res, cunumeric_res, atol(T_OUT), rtol(T_OUT))
-        @test cuNumeric.compare(julia_res, cunumeric_res2, atol(T_OUT), rtol(T_OUT))
-    end
-end
-
-skip_on_integer = (Base.acosh, Base.atanh, Base.atan, Base.acos, Base.asin)
-skip_on_bool = (Base.:(-), skip_on_integer...)
-skip_on_complex = (
-    Base.tanh,
-    Base.deg2rad, Base.rad2deg, Base.sign, Base.cbrt,
-    Base.exp2, Base.expm1, Base.log10, Base.log1p, Base.log2,
-    Base.acos, Base.asin, Base.atan, Base.acosh, Base.asinh, Base.atanh,
-)
-
-function test_unary_function_set(func_dict, T, N)
-    default_generator = (T == Bool) ? :uniform : :unit_interval
-
-    @testset "$func" for func in keys(func_dict)
-        if func in skip_on_complex && (T <: Complex)
-            continue
-        end
-
-        # The are only defined for like 3 integers (-1, 0, 1) so just skip them
-        if func in skip_on_integer && (T <: Integer)
-            continue
-        end
-
-        if func in skip_on_bool && (T == Bool)
-            continue
-        end
-
-        domain_type = get(SPECIAL_DOMAINS, func, default_generator)
-
-        # :uniform is the only generator capable of generating bits
-        skip = (T == Bool && domain_type != :uniform)
-        skip && continue
-
-        julia_arr_1D, julia_arr_2D = make_julia_arrays(T, N, domain_type)
-        cunumeric_arr_1D, cunumeric_arr_2D = make_cunumeric_arrays(
-            [julia_arr_1D], [julia_arr_2D], T, N
-        )
-
-        test_unary_operation(func, julia_arr_1D, cunumeric_arr_1D, T)
-        test_unary_operation(func, julia_arr_2D, cunumeric_arr_2D, T)
-    end
-end
-
-function test_unary_reduction_dims(
-    func, julia_arr::AbstractArray{T,N}, cunumeric_arr::NDArray{T,N}
-) where {T,N}
-    allowpromotion(true) do
-        for d in 1:N
-            julia_res = func(julia_arr; dims=d)
-            cunumeric_res = func(cunumeric_arr; dims=d)
-            n = size(julia_arr, d)
-            allowscalar() do
-                @test cuNumeric.compare(
-                    julia_res, cunumeric_res, reduction_atol(T, n), reduction_rtol(T, n)
-                )
-            end
-        end
-
-        # we are testing a multi axis reduction. This will throw a runtime error.
-        # https://github.com/nv-legate/cupynumeric/blob/main/src/cupynumeric/ndarray.cc#L1132
-        if N >= 2
-            @test_throws Exception func(cunumeric_arr, dims=(1, 2))
-        end
-    end
-end
diff --git a/test/tests/util.jl b/test/util.jl
similarity index 61%
rename from test/tests/util.jl
rename to test/util.jl
index fcf21cc0d..b9c3e596c 100644
--- a/test/tests/util.jl
+++ b/test/util.jl
@@ -38,9 +38,19 @@ atol(::Type{I}) where {I<:Integer} = atol(float(I))
 rtol(::Type{Complex{T}}) where {T} = rtol(T)
 atol(::Type{Complex{T}}) where {T} = atol(T)
 
-# Reduction rounding error grows with the number of elements reduced (n).
+# Reduction compares: association/order differences between Julia and cuNumeric.
+# Forward error for summing n terms is O(n ε Σ|x_i|) (Higham), NOT O(n ε |Σ x_i|).
+# So when partial sums cancel, |result| is small and rtol alone is insufficient —
+# the absolute floor must track the *input* magnitude (`scale`, e.g. maximum(abs, x)).
+#
+# CI example Float32 n=10, |sum|≈0.49, absdiff≈6e-5: relative to the result that
+# looks large (~1e-4), but with |x_i|~1e3 the Higham-scale bound is ~n*eps*scale
+# ≈ 1e-3, and 6e-5 sits comfortably under it (different reduction trees, not a bug).
 reduction_rtol(::Type{T}, n) where {T} = rtol(T) * n
-reduction_atol(::Type{T}, n) where {T} = atol(T) * n
+function reduction_atol(::Type{T}, n, scale=1) where {T}
+    FT = float(real(T))
+    return max(atol(T) * n, n * eps(FT) * abs(scale))
+end
 
 is_same(arr1::NDArray, arr2::NDArray) = @allowscalar (arr1 == arr2)[1]
 is_same(arr1::NDArray, arr2::Array) = @allowscalar (arr1 == arr2)[1]
@@ -92,7 +102,7 @@ function make_cunumeric_arrays(julia_arrs_1D, julia_arrs_2D, T, N; count::Int=1)
     return cunumeric_arrs..., cunumeric_arrs_2D...
 end
 
-function safe_isapprox(x, y, rtol, atol)
+function safe_isapprox(x, y, atol, rtol)
     # Handle NaN
     if isnan(x) && isnan(y)
         return true
@@ -108,22 +118,62 @@ function safe_isapprox(x, y, rtol, atol)
         return false
     end
 
-    return isapprox(x, y; rtol=rtol, atol=atol)
+    return isapprox(x, y; atol=atol, rtol=rtol)
+end
+
+function _safe_compare_fail(idx, left, right, atol, rtol)
+    absdiff = abs(left - right)
+    scale = max(abs(left), abs(right))
+    tol = max(atol, rtol * scale)
+    @error "safe_compare mismatch" index=idx left=left right=right absdiff=absdiff atol=atol rtol=rtol tol=tol
+    return false
 end
 
-function safe_compare(x::AbstractArray{T}, y::NDArray{T}, rtol, atol) where {T}
-    for CI in CartesianIndices(x)
-        if !safe_isapprox(x[CI], y[Tuple(CI)...], rtol, atol)
-            println("Failed at index $(Tuple(CI))")
-            println("x[$(Tuple(CI))] = $(x[CI])")
-            println("y[$(Tuple(CI))] = $(y[Tuple(CI)...])")
-            return false
+# Argument order matches cuNumeric.compare(..., atol, rtol) and existing call sites.
+function safe_compare(
+    julia_array::AbstractArray{T1,N}, arr::NDArray{T2,N}, atol, rtol
+) where {T1,T2,N}
+    if cuNumeric.shape(arr) != size(julia_array)
+        @error "safe_compare shape mismatch" left_size=size(julia_array) right_shape=cuNumeric.shape(
+            arr
+        )
+        return false
+    end
+
+    for CI in CartesianIndices(julia_array)
+        x = julia_array[CI]
+        y = arr[Tuple(CI)...]
+        if !safe_isapprox(x, y, atol, rtol)
+            return _safe_compare_fail(Tuple(CI), x, y, atol, rtol)
         end
     end
 
     return true
 end
 
-function safe_compare(x::NDArray{T}, y::AbstractArray{T}, rtol, atol) where {T}
-    return safe_compare(y, x)
+function safe_compare(
+    arr::NDArray{T2,N}, julia_array::AbstractArray{T1,N}, atol, rtol
+) where {T1,T2,N}
+    return safe_compare(julia_array, arr, atol, rtol)
+end
+
+function safe_compare(
+    arr::NDArray{T1,N}, arr2::NDArray{T2,N}, atol, rtol
+) where {T1,T2,N}
+    if cuNumeric.shape(arr) != cuNumeric.shape(arr2)
+        @error "safe_compare shape mismatch" left_shape=cuNumeric.shape(arr) right_shape=cuNumeric.shape(
+            arr2
+        )
+        return false
+    end
+
+    for CI in CartesianIndices(cuNumeric.shape(arr))
+        x = arr[Tuple(CI)...]
+        y = arr2[Tuple(CI)...]
+        if !safe_isapprox(x, y, atol, rtol)
+            return _safe_compare_fail(Tuple(CI), x, y, atol, rtol)
+        end
+    end
+
+    return true
 end
diff --git a/test/tests/scoping.jl b/test/workflows/grayscott.jl
similarity index 66%
rename from test/tests/scoping.jl
rename to test/workflows/grayscott.jl
index 1d6a6ae5b..e81c364f6 100644
--- a/test/tests/scoping.jl
+++ b/test/workflows/grayscott.jl
@@ -17,6 +17,198 @@
  *            Ethan Meitz <emeitz@andrew.cmu.edu>
 =#
 
+struct ParamsGS{T<:AbstractFloat}
+    dx::T
+    dt::T
+    c_u::T
+    c_v::T
+    f::T
+    k::T
+
+    # Constructor with default values
+    function ParamsGS{T}(
+        dx::T=one(T), c_u::T=one(T), c_v::T=T(0.3), f::T=T(0.03), k::T=T(0.06)
+    ) where {T<:AbstractFloat}
+        return new(dx, dx/5, c_u, c_v, f, k)
+    end
+end
+
+function step(u, v, u_new, v_new, args::ParamsGS)
+    @analyze_lifetimes begin
+        # calculate F_u and F_v functions
+        # currently we don't have NDArray^x working yet.
+        F_u = (
+            (
+                -u[2:(end - 1), 2:(end - 1)] .*
+                (v[2:(end - 1), 2:(end - 1)] .* v[2:(end - 1), 2:(end - 1)])
+            ) + args.f*(1 .- u[2:(end - 1), 2:(end - 1)])
+        )
+        F_v = (
+            (
+                u[2:(end - 1), 2:(end - 1)] .*
+                (v[2:(end - 1), 2:(end - 1)] .* v[2:(end - 1), 2:(end - 1)])
+            ) - (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)]
+            ) ./ args.dx^2 +
+            (
+                u[2:(end - 1), 3:end] - 2*u[2:(end - 1), 2:(end - 1)] +
+                u[2:(end - 1), 1:(end - 2)]
+            ) ./ args.dx^2
+        )
+        v_lap = (
+            (
+                v[3:end, 2:(end - 1)] - 2*v[2:(end - 1), 2:(end - 1)] +
+                v[1:(end - 2), 2:(end - 1)]
+            ) ./ args.dx^2 +
+            (
+                v[2:(end - 1), 3:end] - 2*v[2:(end - 1), 2:(end - 1)] +
+                v[2:(end - 1), 1:(end - 2)]
+            ) ./ args.dx^2
+        )
+
+        # # 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
+        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
+end
+
+# same as above but without @analyze_lifetimes macro
+function step_base(u, v, u_new, v_new, args::ParamsGS)
+    # calculate F_u and F_v functions
+    # currently we don't have NDArray^x working yet.
+    F_u = (
+        (
+            -u[2:(end - 1), 2:(end - 1)] .*
+            (v[2:(end - 1), 2:(end - 1)] .* v[2:(end - 1), 2:(end - 1)])
+        ) + args.f*(1 .- u[2:(end - 1), 2:(end - 1)])
+    )
+    F_v = (
+        (
+            u[2:(end - 1), 2:(end - 1)] .*
+            (v[2:(end - 1), 2:(end - 1)] .* v[2:(end - 1), 2:(end - 1)])
+        ) - (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)]
+        ) ./ args.dx^2 +
+        (
+            u[2:(end - 1), 3:end] - 2*u[2:(end - 1), 2:(end - 1)] +
+            u[2:(end - 1), 1:(end - 2)]
+        ) ./ args.dx^2
+    )
+    v_lap = (
+        (
+            v[3:end, 2:(end - 1)] - 2*v[2:(end - 1), 2:(end - 1)] +
+            v[1:(end - 2), 2:(end - 1)]
+        ) ./ args.dx^2 +
+        (
+            v[2:(end - 1), 3:end] - 2*v[2:(end - 1), 2:(end - 1)] +
+            v[2:(end - 1), 1:(end - 2)]
+        ) ./ args.dx^2
+    )
+
+    # # 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
+    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, :]
+    return v_new[end, :] = v[2, :]
+end
+
+function gray_scott(FT, n_steps, u_rand, v_rand)
+    N = 100
+    dims = (N, N)
+    args = ParamsGS{FT}()
+    u = cuNumeric.ones(FT, dims)
+    v = cuNumeric.zeros(FT, dims)
+    u_new = cuNumeric.zeros(FT, dims)
+    v_new = cuNumeric.zeros(FT, dims)
+
+    u[1:15, 1:15] = u_rand
+    v[1:15, 1:15] = v_rand
+
+    cuNumeric.disable_gc!(; verbose=false)
+    for n in 1:n_steps
+        step(u, v, u_new, v_new, args)
+        u, u_new = u_new, u
+        v, v_new = v_new, v
+    end
+
+    return u, v
+end
+
+function gray_scott_base(FT, n_steps, u_rand, v_rand)
+    N = 100
+    dims = (N, N)
+    args = ParamsGS{FT}()
+    u = cuNumeric.ones(FT, dims)
+    v = cuNumeric.zeros(FT, dims)
+    u_new = cuNumeric.zeros(FT, dims)
+    v_new = cuNumeric.zeros(FT, dims)
+
+    u[1:15, 1:15] = u_rand
+    v[1:15, 1:15] = v_rand
+
+    cuNumeric.init_gc!()
+    for n in 1:n_steps
+        step_base(u, v, u_new, v_new, args)
+        u, u_new = u_new, u
+        v, v_new = v_new, v
+    end
+
+    return u, v
+end
+
+#= 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>
+=#
+
 function run_test(op, op_scope, FT, N)
     a = cuNumeric.rand(FT, (N, N))
     b = cuNumeric.rand(FT, (N, N))
@@ -419,3 +611,27 @@ function run_all_ops(FT, N)
 
     return results
 end
+
+@testset "Gray-Scott 2D" begin
+    N = 100
+    @testset verbose = true for T in Base.uniontypes(cuNumeric.SUPPORTED_FLOAT_TYPES)
+        allowscalar() do
+            results = run_all_ops(T, N)
+            for (name, (c_base, c_scoped)) in results
+                @test safe_compare(c_base, c_scoped, atol(T), rtol(T))
+            end
+
+            u_rand = cuNumeric.rand(T, (15, 15))
+            v_rand = cuNumeric.rand(T, (15, 15))
+
+            u, v = gray_scott_base(T, N, u_rand, v_rand)
+            u_scoped, v_scoped = gray_scott(T, N, u_rand, v_rand)
+
+            @test safe_compare(u, u_scoped, atol(T) * N, rtol(T) * 10)
+        end
+    end
+end
+
+@testset "Rewrite pipeline" begin
+    test_scoping_rewrite_pipeline()
+end