diff --git a/lib/SimpleNonlinearSolve/.github/workflows/CI.yml b/lib/SimpleNonlinearSolve/.github/workflows/CI.yml
index 8855ab61a..3bc260f64 100644
--- a/lib/SimpleNonlinearSolve/.github/workflows/CI.yml
+++ b/lib/SimpleNonlinearSolve/.github/workflows/CI.yml
@@ -32,9 +32,11 @@ jobs:
             ${{ runner.os }}-
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
+        with:
+          annotate: true
         env:
           GROUP: ${{ matrix.group }}
-          JULIA_NUM_THREADS: 11
+          JULIA_NUM_THREADS: "auto"
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v3
         with:
diff --git a/lib/SimpleNonlinearSolve/test/23_test_problems.jl b/lib/SimpleNonlinearSolve/test/23_test_problems.jl
index 18ac5e9bd..73bb304ec 100644
--- a/lib/SimpleNonlinearSolve/test/23_test_problems.jl
+++ b/lib/SimpleNonlinearSolve/test/23_test_problems.jl
@@ -1,4 +1,4 @@
-using SimpleNonlinearSolve, LinearAlgebra, NonlinearProblemLibrary, DiffEqBase, Test
+using SimpleNonlinearSolve, LinearAlgebra, NonlinearProblemLibrary, DiffEqBase, XUnit
 
 problems = NonlinearProblemLibrary.problems
 dicts = NonlinearProblemLibrary.dicts
@@ -37,49 +37,51 @@ function test_on_library(problems, dicts, alg_ops, broken_tests, ϵ = 1e-4;
     end
 end
 
-@testset "SimpleNewtonRaphson 23 Test Problems" begin
-    alg_ops = (SimpleNewtonRaphson(),)
+@testset "23 Test Problems" begin
+    @testcase "SimpleNewtonRaphson 23 Test Problems" begin
+        alg_ops = (SimpleNewtonRaphson(),)
 
-    # dictionary with indices of test problems where method does not converge to small residual
-    broken_tests = Dict(alg => Int[] for alg in alg_ops)
-    broken_tests[alg_ops[1]] = []
+        # dictionary with indices of test problems where method does not converge to small residual
+        broken_tests = Dict(alg => Int[] for alg in alg_ops)
+        broken_tests[alg_ops[1]] = []
 
-    test_on_library(problems, dicts, alg_ops, broken_tests)
-end
+        test_on_library(problems, dicts, alg_ops, broken_tests)
+    end
 
-@testset "SimpleTrustRegion 23 Test Problems" begin
-    alg_ops = (SimpleTrustRegion(),)
+    @testcase "SimpleTrustRegion 23 Test Problems" begin
+        alg_ops = (SimpleTrustRegion(),)
 
-    # dictionary with indices of test problems where method does not converge to small residual
-    broken_tests = Dict(alg => Int[] for alg in alg_ops)
-    broken_tests[alg_ops[1]] = [3, 6, 15, 16, 21]
+        # dictionary with indices of test problems where method does not converge to small residual
+        broken_tests = Dict(alg => Int[] for alg in alg_ops)
+        broken_tests[alg_ops[1]] = [3, 6, 15, 16, 21]
 
-    test_on_library(problems, dicts, alg_ops, broken_tests)
-end
+        test_on_library(problems, dicts, alg_ops, broken_tests)
+    end
 
-@testset "SimpleDFSane 23 Test Problems" begin
-    alg_ops = (SimpleDFSane(),)
+    @testcase "SimpleDFSane 23 Test Problems" begin
+        alg_ops = (SimpleDFSane(),)
 
-    broken_tests = Dict(alg => Int[] for alg in alg_ops)
-    broken_tests[alg_ops[1]] = [1, 2, 3, 4, 5, 6, 11, 21]
+        broken_tests = Dict(alg => Int[] for alg in alg_ops)
+        broken_tests[alg_ops[1]] = [1, 2, 3, 4, 5, 6, 11, 21]
 
-    test_on_library(problems, dicts, alg_ops, broken_tests)
-end
+        test_on_library(problems, dicts, alg_ops, broken_tests)
+    end
 
-@testset "SimpleBroyden 23 Test Problems" begin
-    alg_ops = (SimpleBroyden(),)
+    @testcase "SimpleBroyden 23 Test Problems" begin
+        alg_ops = (SimpleBroyden(),)
 
-    broken_tests = Dict(alg => Int[] for alg in alg_ops)
-    broken_tests[alg_ops[1]] = [1, 5, 11]
+        broken_tests = Dict(alg => Int[] for alg in alg_ops)
+        broken_tests[alg_ops[1]] = [1, 5, 11]
 
-    test_on_library(problems, dicts, alg_ops, broken_tests)
-end
+        test_on_library(problems, dicts, alg_ops, broken_tests)
+    end
 
-@testset "SimpleKlement 23 Test Problems" begin
-    alg_ops = (SimpleKlement(),)
+    @testcase "SimpleKlement 23 Test Problems" begin
+        alg_ops = (SimpleKlement(),)
 
-    broken_tests = Dict(alg => Int[] for alg in alg_ops)
-    broken_tests[alg_ops[1]] = [1, 2, 4, 5, 11, 12, 22]
+        broken_tests = Dict(alg => Int[] for alg in alg_ops)
+        broken_tests[alg_ops[1]] = [1, 2, 4, 5, 11, 12, 22]
 
-    test_on_library(problems, dicts, alg_ops, broken_tests)
+        test_on_library(problems, dicts, alg_ops, broken_tests)
+    end
 end
diff --git a/lib/SimpleNonlinearSolve/test/Project.toml b/lib/SimpleNonlinearSolve/test/Project.toml
index f690826d0..9077234a6 100644
--- a/lib/SimpleNonlinearSolve/test/Project.toml
+++ b/lib/SimpleNonlinearSolve/test/Project.toml
@@ -12,6 +12,7 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+XUnit = "3e3c03f2-1a94-11e9-2981-050a4ca824ab"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [compat]
diff --git a/lib/SimpleNonlinearSolve/test/adjoint.jl b/lib/SimpleNonlinearSolve/test/adjoint.jl
index 27c746981..cf2ef38d7 100644
--- a/lib/SimpleNonlinearSolve/test/adjoint.jl
+++ b/lib/SimpleNonlinearSolve/test/adjoint.jl
@@ -1,6 +1,6 @@
-using ForwardDiff, SciMLSensitivity, SimpleNonlinearSolve, Test, Zygote
+using ForwardDiff, SciMLSensitivity, SimpleNonlinearSolve, XUnit, Zygote
 
-@testset "Simple Adjoint Test" begin
+@testcase "Simple Adjoint Test" begin
     ff(u, p) = u .^ 2 .- p
 
     function solve_nlprob(p)
diff --git a/lib/SimpleNonlinearSolve/test/basictests.jl b/lib/SimpleNonlinearSolve/test/basictests.jl
index 5938197c7..c46230476 100644
--- a/lib/SimpleNonlinearSolve/test/basictests.jl
+++ b/lib/SimpleNonlinearSolve/test/basictests.jl
@@ -1,5 +1,5 @@
 using AllocCheck, LinearSolve, SimpleNonlinearSolve, StaticArrays, Random,
-    LinearAlgebra, Test, ForwardDiff, DiffEqBase
+    LinearAlgebra, XUnit, ForwardDiff, DiffEqBase
 import PolyesterForwardDiff
 
 _nameof(x) = applicable(nameof, x) ? nameof(x) : _nameof(typeof(x))
@@ -25,33 +25,29 @@ const TERMINATION_CONDITIONS = [
     AbsSafeTerminationMode(), RelSafeBestTerminationMode(), AbsSafeBestTerminationMode(),
 ]
 
-# --- SimpleNewtonRaphson tests ---
-
-@testset "$(alg)" for alg in (SimpleNewtonRaphson, SimpleTrustRegion)
-    # Eval else the alg is type unstable
-    @eval begin
-        function benchmark_nlsolve_oop(f, u0, p = 2.0; autodiff = nothing)
-            prob = NonlinearProblem{false}(f, u0, p)
-            return solve(prob, $(alg)(; autodiff), abstol = 1e-9)
-        end
+function benchmark_nlsolve_oop(f::F, u0, p = 2.0; solver) where {F}
+    prob = NonlinearProblem{false}(f, u0, p)
+    return solve(prob, solver; abstol = 1e-9)
+end
+function benchmark_nlsolve_iip(f!::F, u0, p = 2.0; solver) where {F}
+    prob = NonlinearProblem{true}(f!, u0, p)
+    return solve(prob, solver; abstol = 1e-9)
+end
 
-        function benchmark_nlsolve_iip(f, u0, p = 2.0; autodiff = nothing)
-            prob = NonlinearProblem{true}(f, u0, p)
-            return solve(prob, $(alg)(; autodiff), abstol = 1e-9)
-        end
-    end
+# --- SimpleNewtonRaphson tests ---
 
+@testcase "$(alg)" for alg in (SimpleNewtonRaphson, SimpleTrustRegion)
     @testset "AutoDiff: $(_nameof(autodiff))" for autodiff in (AutoFiniteDiff(),
         AutoForwardDiff(), AutoPolyesterForwardDiff())
         @testset "[OOP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
             u0 isa SVector && autodiff isa AutoPolyesterForwardDiff && continue
-            sol = benchmark_nlsolve_oop(quadratic_f, u0; autodiff)
+            sol = benchmark_nlsolve_oop(quadratic_f, u0; solver = alg(; autodiff))
             @test SciMLBase.successful_retcode(sol)
             @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
         end
 
         @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
-            sol = benchmark_nlsolve_iip(quadratic_f!, u0; autodiff)
+            sol = benchmark_nlsolve_iip(quadratic_f!, u0; solver = alg(; autodiff))
             @test SciMLBase.successful_retcode(sol)
             @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
         end
@@ -67,16 +63,11 @@ end
 
 # --- SimpleHalley tests ---
 
-@testset "SimpleHalley" begin
-    function benchmark_nlsolve_oop(f, u0, p = 2.0; autodiff = nothing)
-        prob = NonlinearProblem{false}(f, u0, p)
-        return solve(prob, SimpleHalley(; autodiff), abstol = 1e-9)
-    end
-
+@testcase "SimpleHalley" begin
     @testset "AutoDiff: $(_nameof(autodiff))" for autodiff in (AutoFiniteDiff(),
         AutoForwardDiff())
         @testset "[OOP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
-            sol = benchmark_nlsolve_oop(quadratic_f, u0; autodiff)
+            sol = benchmark_nlsolve_oop(quadratic_f, u0; solver = SimpleHalley(; autodiff))
             @test SciMLBase.successful_retcode(sol)
             @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
         end
@@ -92,27 +83,17 @@ end
 
 # --- SimpleBroyden / SimpleKlement / SimpleLimitedMemoryBroyden tests ---
 
-@testset "$(_nameof(alg))" for alg in [SimpleBroyden(), SimpleKlement(), SimpleDFSane(),
+@testcase "$(_nameof(alg))" for alg in [SimpleBroyden(), SimpleKlement(), SimpleDFSane(),
     SimpleLimitedMemoryBroyden(), SimpleBroyden(; linesearch = Val(true)),
     SimpleLimitedMemoryBroyden(; linesearch = Val(true))]
-    function benchmark_nlsolve_oop(f, u0, p = 2.0)
-        prob = NonlinearProblem{false}(f, u0, p)
-        return solve(prob, alg, abstol = 1e-9)
-    end
-
-    function benchmark_nlsolve_iip(f, u0, p = 2.0)
-        prob = NonlinearProblem{true}(f, u0, p)
-        return solve(prob, alg, abstol = 1e-9)
-    end
-
     @testset "[OOP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
-        sol = benchmark_nlsolve_oop(quadratic_f, u0)
+        sol = benchmark_nlsolve_oop(quadratic_f, u0; solver = alg)
         @test SciMLBase.successful_retcode(sol)
         @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
     end
 
     @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
-        sol = benchmark_nlsolve_iip(quadratic_f!, u0)
+        sol = benchmark_nlsolve_iip(quadratic_f!, u0; solver = alg)
         @test SciMLBase.successful_retcode(sol)
         @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
     end
@@ -126,16 +107,11 @@ end
 end
 
 @testset "Newton Fails" begin
-    function benchmark_nlsolve_oop(f, u0, p, alg)
-        prob = NonlinearProblem{false}(f, u0, p)
-        return solve(prob, alg; abstol = 1e-9)
-    end
-
     u0 = [-10.0, -1.0, 1.0, 2.0, 3.0, 4.0, 10.0]
     p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 
-    for alg in (SimpleDFSane(), SimpleTrustRegion(), SimpleHalley())
-        sol = benchmark_nlsolve_oop(newton_fails, u0, p, alg)
+    @testcase "$(alg)" for alg in (SimpleDFSane(), SimpleTrustRegion(), SimpleHalley())
+        sol = benchmark_nlsolve_oop(newton_fails, u0, p; solver = alg)
         @test SciMLBase.successful_retcode(sol)
         @test all(abs.(newton_fails(sol.u, p)) .< 1e-9)
     end
@@ -144,7 +120,7 @@ end
 # --- Allocation Checks ---
 
 ## SimpleDFSane needs to allocate a history vector
-@testset "Allocation Checks: $(_nameof(alg))" for alg in (SimpleNewtonRaphson(),
+@testcase "Allocation Checks: $(_nameof(alg))" for alg in (SimpleNewtonRaphson(),
     SimpleHalley(), SimpleBroyden(), SimpleKlement(), SimpleLimitedMemoryBroyden(),
     SimpleTrustRegion(), SimpleDFSane(), SimpleBroyden(; linesearch = Val(true)),
     SimpleLimitedMemoryBroyden(; linesearch = Val(true)))
@@ -173,7 +149,7 @@ end
 
 # --- Interval Nonlinear Problems ---
 
-@testset "Interval Nonlinear Problem: $(alg)" for alg in (Bisection(), Falsi(), Ridder(),
+@testcase "Interval Nonlinear Problem: $(alg)" for alg in (Bisection(), Falsi(), Ridder(),
     Brent(), ITP(), Alefeld())
     tspan = (1.0, 20.0)
 
@@ -215,7 +191,8 @@ end
     end
 end
 
-@testset "Tolerance Tests Interval Methods: $(alg)" for alg in (Bisection(), Falsi(), ITP())
+@testcase "Tolerance Tests Interval Methods: $(alg)" for alg in (Bisection(), Falsi(),
+    ITP())
     tspan = (1.0, 20.0)
     probB = IntervalNonlinearProblem(quadratic_f, tspan, 2.0)
     tols = [0.1, 0.01, 0.001, 0.0001, 1e-5, 1e-6, 1e-7]
@@ -228,7 +205,7 @@ end
     end
 end
 
-@testset "Tolerance Tests Interval Methods: $(alg)" for alg in (Ridder(), Brent())
+@testcase "Tolerance Tests Interval Methods: $(alg)" for alg in (Ridder(), Brent())
     tspan = (1.0, 20.0)
     probB = IntervalNonlinearProblem(quadratic_f, tspan, 2.0)
     tols = [0.1] # Ridder and Brent converge rapidly so as we lower tolerance below 0.01, it converges with max precision to the solution
@@ -241,7 +218,7 @@ end
     end
 end
 
-@testset "Flipped Signs and Reversed Tspan: $(alg)" for alg in (Alefeld(), Bisection(),
+@testcase "Flipped Signs and Reversed Tspan: $(alg)" for alg in (Alefeld(), Bisection(),
     Falsi(), Brent(), ITP(), Ridder())
     f1(u, p) = u * u - p
     f2(u, p) = p - u * u
diff --git a/lib/SimpleNonlinearSolve/test/cuda.jl b/lib/SimpleNonlinearSolve/test/cuda.jl
index fddf751de..7708ae5c4 100644
--- a/lib/SimpleNonlinearSolve/test/cuda.jl
+++ b/lib/SimpleNonlinearSolve/test/cuda.jl
@@ -1,4 +1,4 @@
-using SimpleNonlinearSolve, StaticArrays, CUDA, Test
+using SimpleNonlinearSolve, StaticArrays, CUDA, XUnit
 
 CUDA.allowscalar(false)
 
@@ -6,12 +6,10 @@ f(u, p) = u .* u .- 2
 f!(du, u, p) = du .= u .* u .- 2
 
 @testset "Solving on GPUs" begin
-    for alg in (SimpleNewtonRaphson(), SimpleDFSane(), SimpleTrustRegion(), SimpleBroyden(),
-        SimpleLimitedMemoryBroyden(), SimpleKlement(), SimpleHalley(),
-        SimpleBroyden(; linesearch = Val(true)),
+    @testcase "$(alg)" for alg in (SimpleNewtonRaphson(), SimpleDFSane(),
+        SimpleTrustRegion(), SimpleBroyden(), SimpleLimitedMemoryBroyden(), SimpleKlement(),
+        SimpleHalley(), SimpleBroyden(; linesearch = Val(true)),
         SimpleLimitedMemoryBroyden(; linesearch = Val(true)))
-        @info "Testing $alg on CUDA"
-
         # Static Arrays
         u0 = @SVector[1.0f0, 1.0f0]
         probN = NonlinearProblem{false}(f, u0)
@@ -27,12 +25,13 @@ f!(du, u, p) = du .= u .* u .- 2
         @test maximum(abs, sol.resid) ≤ 1.0f-6
 
         # Regular Arrays Inplace
-        alg isa SimpleHalley && continue
-        u0 = [1.0, 1.0]
-        probN = NonlinearProblem{true}(f!, u0)
-        sol = solve(probN, alg; abstol = 1.0f-6)
-        @test SciMLBase.successful_retcode(sol)
-        @test maximum(abs, sol.resid) ≤ 1.0f-6
+        if !(alg isa SimpleHalley)
+            u0 = [1.0, 1.0]
+            probN = NonlinearProblem{true}(f!, u0)
+            sol = solve(probN, alg; abstol = 1.0f-6)
+            @test SciMLBase.successful_retcode(sol)
+            @test maximum(abs, sol.resid) ≤ 1.0f-6
+        end
     end
 end
 
@@ -44,7 +43,8 @@ end
 @testset "CUDA Kernel Launch Test" begin
     prob = NonlinearProblem{false}(f, @SVector[1.0f0, 1.0f0])
 
-    for alg in (SimpleNewtonRaphson(), SimpleDFSane(), SimpleTrustRegion(), SimpleBroyden(),
+    @testcase "$(alg)" for alg in (SimpleNewtonRaphson(), SimpleDFSane(),
+        SimpleTrustRegion(), SimpleBroyden(),
         SimpleLimitedMemoryBroyden(), SimpleKlement(), SimpleHalley(),
         SimpleBroyden(; linesearch = Val(true)),
         SimpleLimitedMemoryBroyden(; linesearch = Val(true)))
diff --git a/lib/SimpleNonlinearSolve/test/cuda/Project.toml b/lib/SimpleNonlinearSolve/test/cuda/Project.toml
index 36c63059b..4fdee40a7 100644
--- a/lib/SimpleNonlinearSolve/test/cuda/Project.toml
+++ b/lib/SimpleNonlinearSolve/test/cuda/Project.toml
@@ -3,3 +3,4 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+XUnit = "3e3c03f2-1a94-11e9-2981-050a4ca824ab"
diff --git a/lib/SimpleNonlinearSolve/test/forward_ad.jl b/lib/SimpleNonlinearSolve/test/forward_ad.jl
index f545ccb0c..78059c2f0 100644
--- a/lib/SimpleNonlinearSolve/test/forward_ad.jl
+++ b/lib/SimpleNonlinearSolve/test/forward_ad.jl
@@ -1,4 +1,4 @@
-using ForwardDiff, SimpleNonlinearSolve, StaticArrays, Test, LinearAlgebra
+using ForwardDiff, SimpleNonlinearSolve, StaticArrays, XUnit, LinearAlgebra
 import SimpleNonlinearSolve: AbstractSimpleNonlinearSolveAlgorithm
 
 test_f!(du, u, p) = (@. du = u^2 - p)
@@ -35,7 +35,7 @@ __compatible(::AbstractSimpleNonlinearSolveAlgorithm, ::Val{:iip}) = true
 __compatible(::AbstractSimpleNonlinearSolveAlgorithm, ::Val{:oop}) = true
 __compatible(::SimpleHalley, ::Val{:iip}) = false
 
-@testset "ForwardDiff.jl Integration: $(alg)" for alg in (SimpleNewtonRaphson(),
+@testcase "ForwardDiff.jl Integration: $(alg)" for alg in (SimpleNewtonRaphson(),
     SimpleTrustRegion(), SimpleHalley(), SimpleBroyden(), SimpleKlement(), SimpleDFSane())
     us = (2.0, @SVector[1.0, 1.0], [1.0, 1.0], ones(2, 2), @SArray ones(2, 2))
 
diff --git a/lib/SimpleNonlinearSolve/test/least_squares.jl b/lib/SimpleNonlinearSolve/test/least_squares.jl
index bc801421f..6c48a9807 100644
--- a/lib/SimpleNonlinearSolve/test/least_squares.jl
+++ b/lib/SimpleNonlinearSolve/test/least_squares.jl
@@ -1,4 +1,4 @@
-using SimpleNonlinearSolve, LinearAlgebra, Test
+using SimpleNonlinearSolve, LinearAlgebra, XUnit
 
 true_function(x, θ) = @. θ[1] * exp(θ[2] * x) * cos(θ[3] * x + θ[4])
 
@@ -14,7 +14,8 @@ end
 θ_init = θ_true .+ 0.1
 prob_oop = NonlinearLeastSquaresProblem{false}(loss_function, θ_init, x)
 
-for solver in [SimpleNewtonRaphson(AutoForwardDiff()), SimpleGaussNewton(AutoForwardDiff()),
+@testcase "Solver: $(solver)" for solver in [SimpleNewtonRaphson(AutoForwardDiff()),
+    SimpleGaussNewton(AutoForwardDiff()),
     SimpleNewtonRaphson(AutoFiniteDiff()), SimpleGaussNewton(AutoFiniteDiff())]
     sol = solve(prob_oop, solver)
     @test norm(sol.resid) < 1e-12
diff --git a/lib/SimpleNonlinearSolve/test/matrix_resizing_tests.jl b/lib/SimpleNonlinearSolve/test/matrix_resizing_tests.jl
index 455aac91a..63a0a1735 100644
--- a/lib/SimpleNonlinearSolve/test/matrix_resizing_tests.jl
+++ b/lib/SimpleNonlinearSolve/test/matrix_resizing_tests.jl
@@ -1,4 +1,4 @@
-using SimpleNonlinearSolve
+using SimpleNonlinearSolve, XUnit
 
 ff(u, p) = u .* u .- p
 u0 = ones(2, 3)
@@ -6,7 +6,7 @@ p = 2.0
 vecprob = NonlinearProblem(ff, vec(u0), p)
 prob = NonlinearProblem(ff, u0, p)
 
-@testset "$(alg)" for alg in (SimpleKlement(), SimpleBroyden(), SimpleNewtonRaphson(),
+@testcase "$(alg)" for alg in (SimpleKlement(), SimpleBroyden(), SimpleNewtonRaphson(),
     SimpleDFSane(), SimpleLimitedMemoryBroyden(; threshold = Val(2)), SimpleTrustRegion())
     @test vec(solve(prob, alg).u) ≈ solve(vecprob, alg).u
 end
diff --git a/lib/SimpleNonlinearSolve/test/runtests.jl b/lib/SimpleNonlinearSolve/test/runtests.jl
index c9ff996dc..5d9216ae3 100644
--- a/lib/SimpleNonlinearSolve/test/runtests.jl
+++ b/lib/SimpleNonlinearSolve/test/runtests.jl
@@ -1,4 +1,4 @@
-using Pkg, SafeTestsets, Test
+using Pkg, SafeTestsets, XUnit
 
 const GROUP = get(ENV, "GROUP", "All")
 
@@ -8,18 +8,18 @@ function activate_env(env)
     Pkg.instantiate()
 end
 
-@time @testset "SimpleNonlinearSolve.jl" begin
+@testset runner=ParallelTestRunner() "SimpleNonlinearSolve.jl" begin
     if GROUP == "All" || GROUP == "Core"
-        @time @safetestset "Basic Tests" include("basictests.jl")
-        @time @safetestset "Forward AD" include("forward_ad.jl")
-        @time @safetestset "Matrix Resizing Tests" include("matrix_resizing_tests.jl")
-        @time @safetestset "Least Squares Tests" include("least_squares.jl")
-        @time @safetestset "23 Test Problems" include("23_test_problems.jl")
-        @time @safetestset "Simple Adjoint Tests" include("adjoint.jl")
+        @safetestset "Basic Tests" include("basictests.jl")
+        @safetestset "Forward AD" include("forward_ad.jl")
+        @safetestset "Matrix Resizing Tests" include("matrix_resizing_tests.jl")
+        @safetestset "Least Squares Tests" include("least_squares.jl")
+        @safetestset "23 Test Problems" include("23_test_problems.jl")
+        @safetestset "Simple Adjoint Tests" include("adjoint.jl")
     end
 
     if GROUP == "CUDA"
         activate_env("cuda")
-        @time @safetestset "CUDA Tests" include("cuda.jl")
+        @safetestset "CUDA Tests" include("cuda.jl")
     end
 end