diff --git a/Project.toml b/Project.toml
index 9271df1f..281cfb00 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "ReverseDiff"
 uuid = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
-version = "1.15.3"
+version = "1.15.4"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -36,7 +36,8 @@ julia = "1"
 DiffTests = "de460e47-3fe3-5279-bb4a-814414816d5d"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["DiffTests", "FillArrays", "IrrationalConstants", "Test"]
+test = ["DiffTests", "FillArrays", "IrrationalConstants", "StaticArrays", "Test"]
diff --git a/src/api/gradients.jl b/src/api/gradients.jl
index ed92cd4d..fe616651 100644
--- a/src/api/gradients.jl
+++ b/src/api/gradients.jl
@@ -21,7 +21,7 @@ call.
 function gradient(f, input, cfg::GradientConfig = GradientConfig(input))
     tape = GradientTape(f, input, cfg)
     result = construct_result(input_hook(tape))
-    seeded_reverse_pass!(result, tape)
+    result = seeded_reverse_pass!(result, tape)
     empty!(cfg.tape)
     return result
 end
@@ -39,7 +39,7 @@ in it as well.
 """
 function gradient!(result, f, input, cfg::GradientConfig = GradientConfig(input))
     tape = GradientTape(f, input, cfg)
-    seeded_reverse_pass!(result, tape)
+    result = seeded_reverse_pass!(result, tape)
     empty!(cfg.tape)
     return result
 end
@@ -60,7 +60,7 @@ of `f` w.r.t. `input[i].`
 """
 function gradient!(tape::Union{GradientTape,CompiledGradient}, input)
     result = construct_result(input_hook(tape))
-    gradient!(result, tape, input)
+    result = gradient!(result, tape, input)
     return result
 end
 
@@ -77,6 +77,6 @@ in it as well.
 """
 function gradient!(result, tape::Union{GradientTape,CompiledGradient}, input)
     seeded_forward_pass!(tape, input)
-    seeded_reverse_pass!(result, tape)
+    result = seeded_reverse_pass!(result, tape)
     return result
 end
diff --git a/src/api/hessians.jl b/src/api/hessians.jl
index e0107cb9..6fbdf356 100644
--- a/src/api/hessians.jl
+++ b/src/api/hessians.jl
@@ -38,12 +38,13 @@ If `result` is a `DiffResults.DiffResult`, the primal value `f(input)` and the g
 """
 function hessian!(result, f, input::AbstractArray, cfg::HessianConfig = HessianConfig(input))
     ∇f = x -> gradient(f, x, cfg.gradient_config)
-    jacobian!(result, ∇f, input, cfg.jacobian_config)
+    result = jacobian!(result, ∇f, input, cfg.jacobian_config)
     return result
 end
 
 function hessian!(result::DiffResult, f, input::AbstractArray,
                   cfg::HessianConfig = HessianConfig(result, input))
+    # TODO: how to realias `result` here?
     ∇f! = (y, x) -> begin
         gradient_result = DiffResult(zero(eltype(y)), y)
         gradient!(gradient_result, f, x, cfg.gradient_config)
@@ -68,7 +69,7 @@ return the Hessian `H(f)(input)`.
 """
 function hessian!(tape::Union{HessianTape,CompiledHessian}, input::AbstractArray)
     result = construct_result(output_hook(tape), input_hook(tape))
-    hessian!(result, tape, input)
+    result = hessian!(result, tape, input)
     return result
 end
 
@@ -85,11 +86,12 @@ If `result` is a `DiffResults.DiffResult`, the primal value `f(input)` and the g
 """
 function hessian!(result::AbstractArray, tape::Union{HessianTape,CompiledHessian}, input::AbstractArray)
     seeded_forward_pass!(tape, input)
-    seeded_reverse_pass!(result, tape)
+    result = seeded_reverse_pass!(result, tape)
     return result
 end
 
 function hessian!(result::DiffResult, tape::Union{HessianTape,CompiledHessian}, input::AbstractArray)
+    # TODO: how to realias `result` here?
     seeded_forward_pass!(tape, input)
     seeded_reverse_pass!(DiffResult(DiffResults.gradient(result), DiffResults.hessian(result)), tape)
     result = DiffResults.value!(result, func_hook(tape)(input))
diff --git a/src/api/jacobians.jl b/src/api/jacobians.jl
index 132957af..800727a7 100644
--- a/src/api/jacobians.jl
+++ b/src/api/jacobians.jl
@@ -41,7 +41,7 @@ in it as well.
 function jacobian!(result, f, input, cfg::JacobianConfig = JacobianConfig(input))
     tape = JacobianTape(f, input, cfg)
     isa(input, TrackedArray) && empty!(input.tape)
-    jacobian!(result, tape, input)
+    result = jacobian!(result, tape, input)
     empty!(tape.tape)
     return result
 end
@@ -74,7 +74,7 @@ form `f!(output::AbstractArray{<:Real}, input::AbstractArray{<:Real}...)`.
 function jacobian!(result, f!, output, input, cfg::JacobianConfig = JacobianConfig(output, input))
     tape = JacobianTape(f!, output, input, cfg)
     isa(input, TrackedArray) && empty!(input.tape)
-    jacobian!(result, tape, input)
+    result = jacobian!(result, tape, input)
     extract_result_value!(output, output_hook(tape))
     empty!(tape.tape)
     return result
@@ -103,7 +103,7 @@ new `output` values into the tape.
 """
 function jacobian!(tape::Union{JacobianTape,CompiledJacobian}, input)
     result = construct_result(output_hook(tape), input_hook(tape))
-    jacobian!(result, tape, input)
+    result = jacobian!(result, tape, input)
     return result
 end
 
@@ -119,7 +119,7 @@ case the primal value of the target function will be stored in it as well.
 """
 function jacobian!(result, tape::Union{JacobianTape,CompiledJacobian}, input)
     seeded_forward_pass!(tape, input)
-    seeded_reverse_pass!(result, tape)
+    result = seeded_reverse_pass!(result, tape)
     return result
 end
 
diff --git a/src/api/utils.jl b/src/api/utils.jl
index a2469b0e..ef98c123 100644
--- a/src/api/utils.jl
+++ b/src/api/utils.jl
@@ -29,12 +29,12 @@ function seeded_reverse_pass!(result, output::TrackedReal, input, tape)
     unseed!(input)
     seed!(output)
     reverse_pass!(tape)
-    extract_result!(result, output, input)
+    result = extract_result!(result, output, input)
     return result
 end
 
 function seeded_reverse_pass!(result, output::Number, input, tape)
-    extract_result!(result, output)
+    result = extract_result!(result, output)
     return result
 end
 
@@ -58,13 +58,13 @@ function seeded_reverse_pass!(result::AbstractArray, output::AbstractArray, inpu
 end
 
 function seeded_reverse_pass!(result::DiffResult, output::AbstractArray, input::TrackedArray, tape)
-    seeded_reverse_pass!(DiffResults.jacobian(result), output, input, tape)
-    extract_result_value!(result, output)
+    result = seeded_reverse_pass!(DiffResults.jacobian(result), output, input, tape)
+    result = extract_result_value!(result, output)
     return result
 end
 
 function seeded_reverse_pass!(result::Tuple, output::AbstractArray, input::Tuple, tape)
-    for i in eachindex(result)
+    result = map(eachindex(result, input)) do i
         seeded_reverse_pass!(result[i], output, input[i], tape)
     end
     return result
@@ -75,14 +75,14 @@ end
 #####################
 
 function extract_result!(result::Tuple, output, input::Tuple)
-    for i in eachindex(result)
+    result = map(eachindex(result, input)) do i
         extract_result!(result[i], output, input[i])
     end
     return result
 end
 
 function extract_result!(result::Tuple, output)
-    for i in eachindex(result)
+    result = map(eachindex(result)) do i
         extract_result!(result[i], output)
     end
     return result
@@ -111,7 +111,7 @@ function extract_result!(result::DiffResult, output::Number)
 end
 
 function extract_result_value!(result::Tuple, output)
-    for i in eachindex(result)
+    result = map(eachindex(result)) do i
         extract_result_value!(result[i], output)
     end
     return result
diff --git a/test/api/GradientTests.jl b/test/api/GradientTests.jl
index 741bd049..17daece5 100644
--- a/test/api/GradientTests.jl
+++ b/test/api/GradientTests.jl
@@ -5,170 +5,194 @@ using DiffTests, ForwardDiff, ReverseDiff, Test, LinearAlgebra
 include(joinpath(dirname(@__FILE__), "../utils.jl"))
 
 function test_unary_gradient(f, x)
-    test = ForwardDiff.gradient!(DiffResults.GradientResult(x), f, x)
+    @testset "Test unary gradient: f=$f - x::$(typeof(x))" begin
 
-    # without GradientConfig
+        test = ForwardDiff.gradient!(DiffResults.GradientResult(x), f, x)
 
-    test_approx(ReverseDiff.gradient(f, x), DiffResults.gradient(test))
+        # without GradientConfig
 
-    out = similar(x)
-    ReverseDiff.gradient!(out, f, x)
-    test_approx(out, DiffResults.gradient(test))
+        @testset "Without GradientConfig" begin
+            test_approx(ReverseDiff.gradient(f, x), DiffResults.gradient(test))
 
-    result = DiffResults.GradientResult(x)
-    ReverseDiff.gradient!(result, f, x)
-    test_approx(DiffResults.value(result), DiffResults.value(test))
-    test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
+            out = similar(x)
+            ReverseDiff.gradient!(out, f, x)
+            test_approx(out, DiffResults.gradient(test))
 
-    # with GradientConfig
+            result = DiffResults.GradientResult(x)
+            result = ReverseDiff.gradient!(result, f, x)
+            test_approx(DiffResults.value(result), DiffResults.value(test))
+            test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
+        end
 
-    cfg = ReverseDiff.GradientConfig(x)
+        # with GradientConfig
 
-    test_approx(ReverseDiff.gradient(f, x, cfg), DiffResults.gradient(test))
+        @testset "With GradientConfig" begin
+            cfg = ReverseDiff.GradientConfig(x)
 
-    out = similar(x)
-    ReverseDiff.gradient!(out, f, x, cfg)
-    test_approx(out, DiffResults.gradient(test))
+            test_approx(ReverseDiff.gradient(f, x, cfg), DiffResults.gradient(test))
 
-    result = DiffResults.GradientResult(x)
-    ReverseDiff.gradient!(result, f, x, cfg)
-    test_approx(DiffResults.value(result), DiffResults.value(test))
-    test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
+            out = similar(x)
+            ReverseDiff.gradient!(out, f, x, cfg)
+            test_approx(out, DiffResults.gradient(test))
 
-    # with GradientTape
+            result = DiffResults.GradientResult(x)
+            result = ReverseDiff.gradient!(result, f, x, cfg)
+            test_approx(DiffResults.value(result), DiffResults.value(test))
+            test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
+        end
 
-    seedx = rand(eltype(x), size(x))
-    tp = ReverseDiff.GradientTape(f, seedx)
+        # with GradientTape
 
-    test_approx(ReverseDiff.gradient!(tp, x), DiffResults.gradient(test))
+        @testset "With GradientTape" begin
+            seedx = rand(eltype(x), size(x))
+            tp = ReverseDiff.GradientTape(f, seedx)
 
-    out = similar(x)
-    ReverseDiff.gradient!(out, tp, x)
-    test_approx(out, DiffResults.gradient(test))
+            test_approx(ReverseDiff.gradient!(tp, x), DiffResults.gradient(test))
 
-    result = DiffResults.GradientResult(x)
-    ReverseDiff.gradient!(result, tp, x)
-    test_approx(DiffResults.value(result), DiffResults.value(test))
-    test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
+            out = similar(x)
+            ReverseDiff.gradient!(out, tp, x)
+            test_approx(out, DiffResults.gradient(test))
 
-    # with compiled GradientTape
+            result = DiffResults.GradientResult(x)
+            result = ReverseDiff.gradient!(result, tp, x)
+            test_approx(DiffResults.value(result), DiffResults.value(test))
+            test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
+        end
 
-    if length(tp.tape) <= COMPILED_TAPE_LIMIT # otherwise compile time can be crazy
-        ctp = ReverseDiff.compile(tp)
+        # with compiled GradientTape
 
-        test_approx(ReverseDiff.gradient!(ctp, x), DiffResults.gradient(test))
+        @testset "With compiled GradientTape" begin
+            if length(tp.tape) <= COMPILED_TAPE_LIMIT # otherwise compile time can be crazy
+                ctp = ReverseDiff.compile(tp)
 
-        out = similar(x)
-        ReverseDiff.gradient!(out, ctp, x)
-        test_approx(out, DiffResults.gradient(test))
+                test_approx(ReverseDiff.gradient!(ctp, x), DiffResults.gradient(test))
+
+                out = similar(x)
+                ReverseDiff.gradient!(out, ctp, x)
+                test_approx(out, DiffResults.gradient(test))
+
+                result = DiffResults.GradientResult(x)
+                result = ReverseDiff.gradient!(result, ctp, x)
+                test_approx(DiffResults.value(result), DiffResults.value(test))
+                test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
+            end
+        end
 
-        result = DiffResults.GradientResult(x)
-        ReverseDiff.gradient!(result, ctp, x)
-        test_approx(DiffResults.value(result), DiffResults.value(test))
-        test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
     end
 end
 
 function test_ternary_gradient(f, a, b, c)
-    test_val = f(a, b, c)
-    test_a = ForwardDiff.gradient(x -> f(x, b, c), a)
-    test_b = ForwardDiff.gradient(x -> f(a, x, c), b)
-    test_c = ForwardDiff.gradient(x -> f(a, b, x), c)
-
-    # without GradientConfig
-
-    ∇a, ∇b, ∇c = ReverseDiff.gradient(f, (a, b, c))
-    test_approx(∇a, test_a)
-    test_approx(∇b, test_b)
-    test_approx(∇c, test_c)
-
-    ∇a, ∇b, ∇c = map(similar, (a, b, c))
-    ReverseDiff.gradient!((∇a, ∇b, ∇c), f, (a, b, c))
-    test_approx(∇a, test_a)
-    test_approx(∇b, test_b)
-    test_approx(∇c, test_c)
-
-    ∇a, ∇b, ∇c = map(DiffResults.GradientResult, (a, b, c))
-    ReverseDiff.gradient!((∇a, ∇b, ∇c), f, (a, b, c))
-    test_approx(DiffResults.value(∇a), test_val)
-    test_approx(DiffResults.value(∇b), test_val)
-    test_approx(DiffResults.value(∇c), test_val)
-    test_approx(DiffResults.gradient(∇a), test_a)
-    test_approx(DiffResults.gradient(∇b), test_b)
-    test_approx(DiffResults.gradient(∇c), test_c)
-
-    # with GradientConfig
-
-    cfg = ReverseDiff.GradientConfig((a, b, c))
-
-    ∇a, ∇b, ∇c = ReverseDiff.gradient(f, (a, b, c), cfg)
-    test_approx(∇a, test_a)
-    test_approx(∇b, test_b)
-    test_approx(∇c, test_c)
-
-    ∇a, ∇b, ∇c = map(similar, (a, b, c))
-    ReverseDiff.gradient!((∇a, ∇b, ∇c), f, (a, b, c), cfg)
-    test_approx(∇a, test_a)
-    test_approx(∇b, test_b)
-    test_approx(∇c, test_c)
-
-    ∇a, ∇b, ∇c = map(DiffResults.GradientResult, (a, b, c))
-    ReverseDiff.gradient!((∇a, ∇b, ∇c), f, (a, b, c), cfg)
-    test_approx(DiffResults.value(∇a), test_val)
-    test_approx(DiffResults.value(∇b), test_val)
-    test_approx(DiffResults.value(∇c), test_val)
-    test_approx(DiffResults.gradient(∇a), test_a)
-    test_approx(DiffResults.gradient(∇b), test_b)
-    test_approx(DiffResults.gradient(∇c), test_c)
-
-    # with GradientTape
-
-    tp = ReverseDiff.GradientTape(f, (rand(eltype(a), size(a)), rand(eltype(b), size(b)), rand(eltype(c), size(c))))
-
-    ∇a, ∇b, ∇c = ReverseDiff.gradient!(tp, (a, b, c))
-    test_approx(∇a, test_a)
-    test_approx(∇b, test_b)
-    test_approx(∇c, test_c)
-
-    ∇a, ∇b, ∇c = map(similar, (a, b, c))
-    ReverseDiff.gradient!((∇a, ∇b, ∇c), tp, (a, b, c))
-    test_approx(∇a, test_a)
-    test_approx(∇b, test_b)
-    test_approx(∇c, test_c)
-
-    ∇a, ∇b, ∇c = map(DiffResults.GradientResult, (a, b, c))
-    ReverseDiff.gradient!((∇a, ∇b, ∇c), tp, (a, b, c))
-    test_approx(DiffResults.value(∇a), test_val)
-    test_approx(DiffResults.value(∇b), test_val)
-    test_approx(DiffResults.value(∇c), test_val)
-    test_approx(DiffResults.gradient(∇a), test_a)
-    test_approx(DiffResults.gradient(∇b), test_b)
-    test_approx(DiffResults.gradient(∇c), test_c)
-
-    # with compiled GradientTape
-
-    if length(tp.tape) <= COMPILED_TAPE_LIMIT # otherwise compile time can be crazy
-        ctp = ReverseDiff.compile(tp)
-
-        ∇a, ∇b, ∇c = ReverseDiff.gradient!(ctp, (a, b, c))
-        test_approx(∇a, test_a)
-        test_approx(∇b, test_b)
-        test_approx(∇c, test_c)
-
-        ∇a, ∇b, ∇c = map(similar, (a, b, c))
-        ReverseDiff.gradient!((∇a, ∇b, ∇c), ctp, (a, b, c))
-        test_approx(∇a, test_a)
-        test_approx(∇b, test_b)
-        test_approx(∇c, test_c)
-
-        ∇a, ∇b, ∇c = map(DiffResults.GradientResult, (a, b, c))
-        ReverseDiff.gradient!((∇a, ∇b, ∇c), ctp, (a, b, c))
-        test_approx(DiffResults.value(∇a), test_val)
-        test_approx(DiffResults.value(∇b), test_val)
-        test_approx(DiffResults.value(∇c), test_val)
-        test_approx(DiffResults.gradient(∇a), test_a)
-        test_approx(DiffResults.gradient(∇b), test_b)
-        test_approx(DiffResults.gradient(∇c), test_c)
+    @testset "Test ternary gradient: f=$f - a::$(typeof(a)), b::$(typeof(b)), c::$(typeof(c))" begin
+
+        test_val = f(a, b, c)
+        test_a = ForwardDiff.gradient(x -> f(x, b, c), a)
+        test_b = ForwardDiff.gradient(x -> f(a, x, c), b)
+        test_c = ForwardDiff.gradient(x -> f(a, b, x), c)
+
+        # without GradientConfig
+
+        @testset "Without GradientConfig" begin
+            ∇a, ∇b, ∇c = ReverseDiff.gradient(f, (a, b, c))
+            test_approx(∇a, test_a)
+            test_approx(∇b, test_b)
+            test_approx(∇c, test_c)
+
+            ∇a, ∇b, ∇c = map(similar, (a, b, c))
+            ReverseDiff.gradient!((∇a, ∇b, ∇c), f, (a, b, c))
+            test_approx(∇a, test_a)
+            test_approx(∇b, test_b)
+            test_approx(∇c, test_c)
+
+            ∇a, ∇b, ∇c = map(DiffResults.GradientResult, (a, b, c))
+            (∇a, ∇b, ∇c) = ReverseDiff.gradient!((∇a, ∇b, ∇c), f, (a, b, c))
+            test_approx(DiffResults.value(∇a), test_val)
+            test_approx(DiffResults.value(∇b), test_val)
+            test_approx(DiffResults.value(∇c), test_val)
+            test_approx(DiffResults.gradient(∇a), test_a)
+            test_approx(DiffResults.gradient(∇b), test_b)
+            test_approx(DiffResults.gradient(∇c), test_c)
+        end
+
+        # with GradientConfig
+
+        @testset "With GradientConfig" begin
+            cfg = ReverseDiff.GradientConfig((a, b, c))
+
+            ∇a, ∇b, ∇c = ReverseDiff.gradient(f, (a, b, c), cfg)
+            test_approx(∇a, test_a)
+            test_approx(∇b, test_b)
+            test_approx(∇c, test_c)
+
+            ∇a, ∇b, ∇c = map(similar, (a, b, c))
+            ReverseDiff.gradient!((∇a, ∇b, ∇c), f, (a, b, c), cfg)
+            test_approx(∇a, test_a)
+            test_approx(∇b, test_b)
+            test_approx(∇c, test_c)
+
+            ∇a, ∇b, ∇c = map(DiffResults.GradientResult, (a, b, c))
+            (∇a, ∇b, ∇c) = ReverseDiff.gradient!((∇a, ∇b, ∇c), f, (a, b, c), cfg)
+            test_approx(DiffResults.value(∇a), test_val)
+            test_approx(DiffResults.value(∇b), test_val)
+            test_approx(DiffResults.value(∇c), test_val)
+            test_approx(DiffResults.gradient(∇a), test_a)
+            test_approx(DiffResults.gradient(∇b), test_b)
+            test_approx(DiffResults.gradient(∇c), test_c)
+        end
+
+        # with GradientTape
+
+        @testset "With GradientTape" begin
+            tp = ReverseDiff.GradientTape(f, (rand(eltype(a), size(a)), rand(eltype(b), size(b)), rand(eltype(c), size(c))))
+
+            ∇a, ∇b, ∇c = ReverseDiff.gradient!(tp, (a, b, c))
+            test_approx(∇a, test_a)
+            test_approx(∇b, test_b)
+            test_approx(∇c, test_c)
+
+            ∇a, ∇b, ∇c = map(similar, (a, b, c))
+            ReverseDiff.gradient!((∇a, ∇b, ∇c), tp, (a, b, c))
+            test_approx(∇a, test_a)
+            test_approx(∇b, test_b)
+            test_approx(∇c, test_c)
+
+            ∇a, ∇b, ∇c = map(DiffResults.GradientResult, (a, b, c))
+            (∇a, ∇b, ∇c) = ReverseDiff.gradient!((∇a, ∇b, ∇c), tp, (a, b, c))
+            test_approx(DiffResults.value(∇a), test_val)
+            test_approx(DiffResults.value(∇b), test_val)
+            test_approx(DiffResults.value(∇c), test_val)
+            test_approx(DiffResults.gradient(∇a), test_a)
+            test_approx(DiffResults.gradient(∇b), test_b)
+            test_approx(DiffResults.gradient(∇c), test_c)
+        end
+
+        # with compiled GradientTape
+
+        @testset "With compiled GradientTape" begin
+            if length(tp.tape) <= COMPILED_TAPE_LIMIT # otherwise compile time can be crazy
+                ctp = ReverseDiff.compile(tp)
+
+                ∇a, ∇b, ∇c = ReverseDiff.gradient!(ctp, (a, b, c))
+                test_approx(∇a, test_a)
+                test_approx(∇b, test_b)
+                test_approx(∇c, test_c)
+
+                ∇a, ∇b, ∇c = map(similar, (a, b, c))
+                ReverseDiff.gradient!((∇a, ∇b, ∇c), ctp, (a, b, c))
+                test_approx(∇a, test_a)
+                test_approx(∇b, test_b)
+                test_approx(∇c, test_c)
+
+                ∇a, ∇b, ∇c = map(DiffResults.GradientResult, (a, b, c))
+                (∇a, ∇b, ∇c) = ReverseDiff.gradient!((∇a, ∇b, ∇c), ctp, (a, b, c))
+                test_approx(DiffResults.value(∇a), test_val)
+                test_approx(DiffResults.value(∇b), test_val)
+                test_approx(DiffResults.value(∇c), test_val)
+                test_approx(DiffResults.gradient(∇a), test_a)
+                test_approx(DiffResults.gradient(∇b), test_b)
+                test_approx(DiffResults.gradient(∇c), test_c)
+            end
+        end
+
     end
 end
 
@@ -193,10 +217,10 @@ norm_hermitian2(v) = (A = I - 2 * v * transpose(v); norm(transpose(A) * A))
 norm_hermitian3(v) = (A = I - 2 * v * collect(v'); norm(collect(A') * A))
 norm_hermitian4(v) = (A = I - 2 * v * v'; norm(transpose(A) * A))
 norm_hermitian5(v) = (A = I - 2 * v * transpose(v); norm(A' * A))
-norm_hermitian6(v) = (A = (v'v)*I - 2 * v * v'; norm(A' * A))
+norm_hermitian6(v) = (A = (v'v) * I - 2 * v * v'; norm(A' * A))
 
 for f in (norm_hermitian1, norm_hermitian2, norm_hermitian3,
-            norm_hermitian4, norm_hermitian5, norm_hermitian6)
+    norm_hermitian4, norm_hermitian5, norm_hermitian6)
     test_println("VECTOR_TO_NUMBER_FUNCS", f)
     test_unary_gradient(f, rand(5))
 end
diff --git a/test/api/HessianTests.jl b/test/api/HessianTests.jl
index 6dae5018..4667f1a3 100644
--- a/test/api/HessianTests.jl
+++ b/test/api/HessianTests.jl
@@ -12,73 +12,83 @@ ReverseDiff.hessian(DiffTests.mat2num_1, rand(3, 3))
 hess_test_approx(a, b) = test_approx(a, b, 1e-4)
 
 function test_unary_hessian(f, x)
-    test = DiffResults.HessianResult(x)
-    ForwardDiff.hessian!(test, f, x, ForwardDiff.HessianConfig(f, test, x, ForwardDiff.Chunk{1}()))
+    @testset "Test unary Hessian: f=$f - x::$(typeof(x))" begin
+        test = DiffResults.HessianResult(x)
+        ForwardDiff.hessian!(test, f, x, ForwardDiff.HessianConfig(f, test, x, ForwardDiff.Chunk{1}()))
 
-    # without HessianConfig
+        # without HessianConfig
 
-    hess_test_approx(ReverseDiff.hessian(f, x), DiffResults.hessian(test))
+        @testst "Without HessianConfig" begin
+            hess_test_approx(ReverseDiff.hessian(f, x), DiffResults.hessian(test))
 
-    out = similar(DiffResults.hessian(test))
-    ReverseDiff.hessian!(out, f, x)
-    hess_test_approx(out, DiffResults.hessian(test))
+            out = similar(DiffResults.hessian(test))
+            ReverseDiff.hessian!(out, f, x)
+            hess_test_approx(out, DiffResults.hessian(test))
 
-    result = DiffResults.HessianResult(x)
-    ReverseDiff.hessian!(result, f, x)
-    hess_test_approx(DiffResults.value(result), DiffResults.value(test))
-    hess_test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
-    hess_test_approx(DiffResults.hessian(result), DiffResults.hessian(test))
+            result = DiffResults.HessianResult(x)
+            result = ReverseDiff.hessian!(result, f, x)
+            hess_test_approx(DiffResults.value(result), DiffResults.value(test))
+            hess_test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
+            hess_test_approx(DiffResults.hessian(result), DiffResults.hessian(test))
+        end
 
-    # with HessianConfig
+        # with HessianConfig
 
-    cfg = ReverseDiff.HessianConfig(x)
+        @testset "With HessianConfig" begin
+            cfg = ReverseDiff.HessianConfig(x)
 
-    hess_test_approx(ReverseDiff.hessian(f, x, cfg), DiffResults.hessian(test))
+            hess_test_approx(ReverseDiff.hessian(f, x, cfg), DiffResults.hessian(test))
 
-    out = similar(DiffResults.hessian(test))
-    ReverseDiff.hessian!(out, f, x, cfg)
-    hess_test_approx(out, DiffResults.hessian(test))
+            out = similar(DiffResults.hessian(test))
+            ReverseDiff.hessian!(out, f, x, cfg)
+            hess_test_approx(out, DiffResults.hessian(test))
 
-    result = DiffResults.HessianResult(x)
-    cfg = ReverseDiff.HessianConfig(result, x)
-    ReverseDiff.hessian!(result, f, x, cfg)
-    hess_test_approx(DiffResults.value(result), DiffResults.value(test))
-    hess_test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
-    hess_test_approx(DiffResults.hessian(result), DiffResults.hessian(test))
+            result = DiffResults.HessianResult(x)
+            cfg = ReverseDiff.HessianConfig(result, x)
+            result = ReverseDiff.hessian!(result, f, x, cfg)
+            hess_test_approx(DiffResults.value(result), DiffResults.value(test))
+            hess_test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
+            hess_test_approx(DiffResults.hessian(result), DiffResults.hessian(test))
+        end
 
-    # with HessianTape
+        # with HessianTape
 
-    seedx = rand(eltype(x), size(x))
-    tp = ReverseDiff.HessianTape(f, seedx)
+        @testset "With HessianTape" begin
+            seedx = rand(eltype(x), size(x))
+            tp = ReverseDiff.HessianTape(f, seedx)
 
-    hess_test_approx(ReverseDiff.hessian!(tp, x), DiffResults.hessian(test))
+            hess_test_approx(ReverseDiff.hessian!(tp, x), DiffResults.hessian(test))
 
-    out = similar(DiffResults.hessian(test))
-    ReverseDiff.hessian!(out, tp, x)
-    hess_test_approx(out, DiffResults.hessian(test))
+            out = similar(DiffResults.hessian(test))
+            ReverseDiff.hessian!(out, tp, x)
+            hess_test_approx(out, DiffResults.hessian(test))
 
-    result = DiffResults.HessianResult(x)
-    ReverseDiff.hessian!(result, tp, x)
-    hess_test_approx(DiffResults.value(result), DiffResults.value(test))
-    hess_test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
-    hess_test_approx(DiffResults.hessian(result), DiffResults.hessian(test))
+            result = DiffResults.HessianResult(x)
+            result = ReverseDiff.hessian!(result, tp, x)
+            hess_test_approx(DiffResults.value(result), DiffResults.value(test))
+            hess_test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
+            hess_test_approx(DiffResults.hessian(result), DiffResults.hessian(test))
+        end
 
-    # with compiled HessianTape
+        # with compiled HessianTape
 
-    if length(tp.tape) <= 10000 # otherwise compile time can be crazy
-        ctp = ReverseDiff.compile(tp)
+        @testset "With compiled HessianTape" begin
+            if length(tp.tape) <= 10000 # otherwise compile time can be crazy
+                ctp = ReverseDiff.compile(tp)
 
-        hess_test_approx(ReverseDiff.hessian!(ctp, x), DiffResults.hessian(test))
+                hess_test_approx(ReverseDiff.hessian!(ctp, x), DiffResults.hessian(test))
 
-        out = similar(DiffResults.hessian(test))
-        ReverseDiff.hessian!(out, ctp, x)
-        hess_test_approx(out, DiffResults.hessian(test))
+                out = similar(DiffResults.hessian(test))
+                ReverseDiff.hessian!(out, ctp, x)
+                hess_test_approx(out, DiffResults.hessian(test))
 
-        result = DiffResults.HessianResult(x)
-        ReverseDiff.hessian!(result, ctp, x)
-        hess_test_approx(DiffResults.value(result), DiffResults.value(test))
-        hess_test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
-        hess_test_approx(DiffResults.hessian(result), DiffResults.hessian(test))
+                result = DiffResults.HessianResult(x)
+                result = ReverseDiff.hessian!(result, ctp, x)
+                hess_test_approx(DiffResults.value(result), DiffResults.value(test))
+                hess_test_approx(DiffResults.gradient(result), DiffResults.gradient(test))
+                hess_test_approx(DiffResults.hessian(result), DiffResults.hessian(test))
+            end
+        end
     end
 end
 
diff --git a/test/api/JacobianTests.jl b/test/api/JacobianTests.jl
index 507edc7b..c02254bf 100644
--- a/test/api/JacobianTests.jl
+++ b/test/api/JacobianTests.jl
@@ -5,244 +5,274 @@ using DiffTests, ForwardDiff, ReverseDiff, Test
 include(joinpath(dirname(@__FILE__), "../utils.jl"))
 
 function test_unary_jacobian(f, x)
-    test_val = f(x)
-    test = ForwardDiff.jacobian!(DiffResults.JacobianResult(test_val, x), f, x, ForwardDiff.JacobianConfig(f, x))
+    @testset "Test unary Jacobian: f=$f - x::$(typeof(x))" begin
+        test_val = f(x)
+        test = ForwardDiff.jacobian!(DiffResults.JacobianResult(test_val, x), f, x, ForwardDiff.JacobianConfig(f, x))
 
-    # without JacobianConfig
+        # without JacobianConfig
 
-    test_approx(ReverseDiff.jacobian(f, x), DiffResults.jacobian(test))
+        @testset "Without JacobianConfig" begin
+            test_approx(ReverseDiff.jacobian(f, x), DiffResults.jacobian(test))
 
-    out = similar(DiffResults.jacobian(test))
-    ReverseDiff.jacobian!(out, f, x)
-    test_approx(out, DiffResults.jacobian(test))
+            out = similar(DiffResults.jacobian(test))
+            ReverseDiff.jacobian!(out, f, x)
+            test_approx(out, DiffResults.jacobian(test))
 
-    result = DiffResults.JacobianResult(test_val, x)
-    ReverseDiff.jacobian!(result, f, x)
-    test_approx(DiffResults.value(result), DiffResults.value(test))
-    test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
+            result = DiffResults.JacobianResult(test_val, x)
+            result = ReverseDiff.jacobian!(result, f, x)
+            test_approx(DiffResults.value(result), DiffResults.value(test))
+            test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
+        end
 
-    # with JacobianConfig
+        # with JacobianConfig
 
-    cfg = ReverseDiff.JacobianConfig(x)
+        @testset "With JacobianConfig" begin
+            cfg = ReverseDiff.JacobianConfig(x)
 
-    test_approx(ReverseDiff.jacobian(f, x, cfg), DiffResults.jacobian(test))
+            test_approx(ReverseDiff.jacobian(f, x, cfg), DiffResults.jacobian(test))
 
-    out = similar(DiffResults.jacobian(test))
-    ReverseDiff.jacobian!(out, f, x, cfg)
-    test_approx(out, DiffResults.jacobian(test))
+            out = similar(DiffResults.jacobian(test))
+            ReverseDiff.jacobian!(out, f, x, cfg)
+            test_approx(out, DiffResults.jacobian(test))
 
-    result = DiffResults.JacobianResult(test_val, x)
-    ReverseDiff.jacobian!(result, f, x, cfg)
-    test_approx(DiffResults.value(result), DiffResults.value(test))
-    test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
+            result = DiffResults.JacobianResult(test_val, x)
+            result = ReverseDiff.jacobian!(result, f, x, cfg)
+            test_approx(DiffResults.value(result), DiffResults.value(test))
+            test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
+        end
 
-    # with JacobianTape
+        # with JacobianTape
 
-    tp = ReverseDiff.JacobianTape(f, rand(eltype(x), size(x)))
+        @testset "With JacobianTape" begin
+            tp = ReverseDiff.JacobianTape(f, rand(eltype(x), size(x)))
 
-    test_approx(ReverseDiff.jacobian!(tp, x), DiffResults.jacobian(test))
+            test_approx(ReverseDiff.jacobian!(tp, x), DiffResults.jacobian(test))
 
-    out = similar(DiffResults.jacobian(test))
-    ReverseDiff.jacobian!(out, tp, x)
-    test_approx(out, DiffResults.jacobian(test))
+            out = similar(DiffResults.jacobian(test))
+            ReverseDiff.jacobian!(out, tp, x)
+            test_approx(out, DiffResults.jacobian(test))
 
-    result = DiffResults.JacobianResult(test_val, x)
-    ReverseDiff.jacobian!(result, tp, x)
-    test_approx(DiffResults.value(result), DiffResults.value(test))
-    test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
+            result = DiffResults.JacobianResult(test_val, x)
+            result = ReverseDiff.jacobian!(result, tp, x)
+            test_approx(DiffResults.value(result), DiffResults.value(test))
+            test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
+        end
 
-    # with compiled JacobianTape
+        # with compiled JacobianTape
 
-    if length(tp.tape) <= COMPILED_TAPE_LIMIT # otherwise compile time can be crazy
-        ctp = ReverseDiff.compile(tp)
+        @testset "With compiled JacobianTape" begin
+            if length(tp.tape) <= COMPILED_TAPE_LIMIT # otherwise compile time can be crazy
+                ctp = ReverseDiff.compile(tp)
 
-        test_approx(ReverseDiff.jacobian!(ctp, x), DiffResults.jacobian(test))
+                test_approx(ReverseDiff.jacobian!(ctp, x), DiffResults.jacobian(test))
 
-        out = similar(DiffResults.jacobian(test))
-        ReverseDiff.jacobian!(out, ctp, x)
-        test_approx(out, DiffResults.jacobian(test))
+                out = similar(DiffResults.jacobian(test))
+                ReverseDiff.jacobian!(out, ctp, x)
+                test_approx(out, DiffResults.jacobian(test))
 
-        result = DiffResults.JacobianResult(test_val, x)
-        ReverseDiff.jacobian!(result, ctp, x)
-        test_approx(DiffResults.value(result), DiffResults.value(test))
-        test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
+                result = DiffResults.JacobianResult(test_val, x)
+                result = ReverseDiff.jacobian!(result, ctp, x)
+                test_approx(DiffResults.value(result), DiffResults.value(test))
+                test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
+            end
+        end
     end
 end
 
 function test_unary_jacobian(f!, y, x)
-    y_original = copy(y)
-    y_copy = copy(y)
-    test = ForwardDiff.jacobian!(DiffResults.JacobianResult(y_copy, x), f!, y_copy, x)
-
-    # without JacobianConfig
-
-    out = ReverseDiff.jacobian(f!, y, x)
-    test_approx(y, DiffResults.value(test))
-    test_approx(out, DiffResults.jacobian(test))
-    copyto!(y, y_original)
-
-    out = similar(DiffResults.jacobian(test))
-    ReverseDiff.jacobian!(out, f!, y, x)
-    test_approx(y,   DiffResults.value(test))
-    test_approx(out, DiffResults.jacobian(test))
-    copyto!(y, y_original)
-
-    result = DiffResults.JacobianResult(y, x)
-    ReverseDiff.jacobian!(result, f!, y, x)
-    @test DiffResults.value(result) == y
-    test_approx(y, DiffResults.value(test))
-    test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
-    copyto!(y, y_original)
-
-    # with JacobianConfig
-
-    cfg = ReverseDiff.JacobianConfig(y, x)
-
-    out = ReverseDiff.jacobian(f!, y, x, cfg)
-    test_approx(y,   DiffResults.value(test))
-    test_approx(out, DiffResults.jacobian(test))
-    copyto!(y, y_original)
-
-    out = similar(DiffResults.jacobian(test))
-    ReverseDiff.jacobian!(out, f!, y, x, cfg)
-    test_approx(y,   DiffResults.value(test))
-    test_approx(out, DiffResults.jacobian(test))
-    copyto!(y, y_original)
-
-    result = DiffResults.JacobianResult(y, x)
-    ReverseDiff.jacobian!(result, f!, y, x, cfg)
-    @test DiffResults.value(result) == y
-    test_approx(y, DiffResults.value(test))
-    test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
-    copyto!(y, y_original)
-
-    # with JacobianTape
-
-    tp = ReverseDiff.JacobianTape(f!, y, rand(eltype(x), size(x)))
-
-    out = ReverseDiff.jacobian!(tp, x)
-    test_approx(out, DiffResults.jacobian(test))
-
-    out = similar(DiffResults.jacobian(test))
-    ReverseDiff.jacobian!(out, tp, x)
-    test_approx(out, DiffResults.jacobian(test))
-
-    result = DiffResults.JacobianResult(y, x)
-    ReverseDiff.jacobian!(result, tp, x)
-    test_approx(DiffResults.value(result), DiffResults.value(test))
-    test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
-
-    # with compiled JacobianTape
-
-    if length(tp.tape) <= COMPILED_TAPE_LIMIT # otherwise compile time can be crazy
-        ctp = ReverseDiff.compile(tp)
-
-        out = ReverseDiff.jacobian!(ctp, x)
-
-        test_approx(out, DiffResults.jacobian(test))
-
-        out = similar(DiffResults.jacobian(test))
-        ReverseDiff.jacobian!(out, ctp, x)
-        test_approx(out, DiffResults.jacobian(test))
-
-        result = DiffResults.JacobianResult(y, x)
-        ReverseDiff.jacobian!(result, ctp, x)
-        test_approx(DiffResults.value(result), DiffResults.value(test))
-        test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
+    @testset "Test unary Jacobian: f!=$f! - y::$(typeof(y)), x::$(typeof(x))" begin
+        y_original = copy(y)
+        y_copy = copy(y)
+        test = ForwardDiff.jacobian!(DiffResults.JacobianResult(y_copy, x), f!, y_copy, x)
+
+        # without JacobianConfig
+
+        @testset "Without JacobianConfig" begin
+            out = ReverseDiff.jacobian(f!, y, x)
+            test_approx(y, DiffResults.value(test))
+            test_approx(out, DiffResults.jacobian(test))
+            copyto!(y, y_original)
+
+            out = similar(DiffResults.jacobian(test))
+            ReverseDiff.jacobian!(out, f!, y, x)
+            test_approx(y, DiffResults.value(test))
+            test_approx(out, DiffResults.jacobian(test))
+            copyto!(y, y_original)
+
+            result = DiffResults.JacobianResult(y, x)
+            result = ReverseDiff.jacobian!(result, f!, y, x)
+            @test DiffResults.value(result) == y
+            test_approx(y, DiffResults.value(test))
+            test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
+            copyto!(y, y_original)
+        end
+
+        # with JacobianConfig
+
+        @testset "With JacobianConfig" begin
+            cfg = ReverseDiff.JacobianConfig(y, x)
+
+            out = ReverseDiff.jacobian(f!, y, x, cfg)
+            test_approx(y, DiffResults.value(test))
+            test_approx(out, DiffResults.jacobian(test))
+            copyto!(y, y_original)
+
+            out = similar(DiffResults.jacobian(test))
+            ReverseDiff.jacobian!(out, f!, y, x, cfg)
+            test_approx(y, DiffResults.value(test))
+            test_approx(out, DiffResults.jacobian(test))
+            copyto!(y, y_original)
+
+            result = DiffResults.JacobianResult(y, x)
+            result = ReverseDiff.jacobian!(result, f!, y, x, cfg)
+            @test DiffResults.value(result) == y
+            test_approx(y, DiffResults.value(test))
+            test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
+            copyto!(y, y_original)
+        end
+
+        # with JacobianTape
+
+        @testset "With JacobianTape" begin
+            tp = ReverseDiff.JacobianTape(f!, y, rand(eltype(x), size(x)))
+
+            out = ReverseDiff.jacobian!(tp, x)
+            test_approx(out, DiffResults.jacobian(test))
+
+            out = similar(DiffResults.jacobian(test))
+            ReverseDiff.jacobian!(out, tp, x)
+            test_approx(out, DiffResults.jacobian(test))
+
+            result = DiffResults.JacobianResult(y, x)
+            result = ReverseDiff.jacobian!(result, tp, x)
+            test_approx(DiffResults.value(result), DiffResults.value(test))
+            test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
+        end
+
+        # with compiled JacobianTape
+
+        @testset "With compiled JacobianTape" begin
+            if length(tp.tape) <= COMPILED_TAPE_LIMIT # otherwise compile time can be crazy
+                ctp = ReverseDiff.compile(tp)
+
+                out = ReverseDiff.jacobian!(ctp, x)
+
+                test_approx(out, DiffResults.jacobian(test))
+
+                out = similar(DiffResults.jacobian(test))
+                ReverseDiff.jacobian!(out, ctp, x)
+                test_approx(out, DiffResults.jacobian(test))
+
+                result = DiffResults.JacobianResult(y, x)
+                result = ReverseDiff.jacobian!(result, ctp, x)
+                test_approx(DiffResults.value(result), DiffResults.value(test))
+                test_approx(DiffResults.jacobian(result), DiffResults.jacobian(test))
+            end
+        end
     end
 end
 
 function test_binary_jacobian(f, a, b)
-    test_val = f(a, b)
-    test_a = ForwardDiff.jacobian(x -> f(x, b), a)
-    test_b = ForwardDiff.jacobian(x -> f(a, x), b)
-
-    # without JacobianConfig
-
-    Ja, Jb = ReverseDiff.jacobian(f, (a, b))
-    test_approx(Ja, test_a)
-    test_approx(Jb, test_b)
-
-    Ja = similar(a, length(test_val), length(a))
-    Jb = similar(b, length(test_val), length(b))
-    ReverseDiff.jacobian!((Ja, Jb), f, (a, b))
-    test_approx(Ja, test_a)
-    test_approx(Jb, test_b)
-
-    Ja = DiffResults.JacobianResult(test_val, a)
-    Jb = DiffResults.JacobianResult(test_val, b)
-    ReverseDiff.jacobian!((Ja, Jb), f, (a, b))
-    test_approx(DiffResults.value(Ja), test_val)
-    test_approx(DiffResults.value(Jb), test_val)
-    test_approx(DiffResults.gradient(Ja), test_a)
-    test_approx(DiffResults.gradient(Jb), test_b)
-
-    # with JacobianConfig
-
-    cfg = ReverseDiff.JacobianConfig((a, b))
-
-    Ja, Jb = ReverseDiff.jacobian(f, (a, b), cfg)
-    test_approx(Ja, test_a)
-    test_approx(Jb, test_b)
-
-    Ja = similar(a, length(test_val), length(a))
-    Jb = similar(b, length(test_val), length(b))
-    ReverseDiff.jacobian!((Ja, Jb), f, (a, b), cfg)
-    test_approx(Ja, test_a)
-    test_approx(Jb, test_b)
-
-    Ja = DiffResults.JacobianResult(test_val, a)
-    Jb = DiffResults.JacobianResult(test_val, b)
-    ReverseDiff.jacobian!((Ja, Jb), f, (a, b), cfg)
-    test_approx(DiffResults.value(Ja), test_val)
-    test_approx(DiffResults.value(Jb), test_val)
-    test_approx(DiffResults.jacobian(Ja), test_a)
-    test_approx(DiffResults.jacobian(Jb), test_b)
-
-    # with JacobianTape
-
-    tp = ReverseDiff.JacobianTape(f, (rand(eltype(a), size(a)), rand(eltype(b), size(b))))
-
-    Ja, Jb = ReverseDiff.jacobian!(tp, (a, b))
-    test_approx(Ja, test_a)
-    test_approx(Jb, test_b)
-
-    Ja = similar(a, length(test_val), length(a))
-    Jb = similar(b, length(test_val), length(b))
-    ReverseDiff.jacobian!((Ja, Jb), tp, (a, b))
-    test_approx(Ja, test_a)
-    test_approx(Jb, test_b)
-
-    Ja = DiffResults.JacobianResult(test_val, a)
-    Jb = DiffResults.JacobianResult(test_val, b)
-    ReverseDiff.jacobian!((Ja, Jb), tp, (a, b))
-    test_approx(DiffResults.value(Ja), test_val)
-    test_approx(DiffResults.value(Jb), test_val)
-    test_approx(DiffResults.gradient(Ja), test_a)
-    test_approx(DiffResults.gradient(Jb), test_b)
-
-    # with compiled JacobianTape
-
-    if length(tp.tape) <= COMPILED_TAPE_LIMIT # otherwise compile time can be crazy
-        ctp = ReverseDiff.compile(tp)
-
-        Ja, Jb = ReverseDiff.jacobian!(ctp, (a, b))
-        test_approx(Ja, test_a)
-        test_approx(Jb, test_b)
-
-        Ja = similar(a, length(test_val), length(a))
-        Jb = similar(b, length(test_val), length(b))
-        ReverseDiff.jacobian!((Ja, Jb), ctp, (a, b))
-        test_approx(Ja, test_a)
-        test_approx(Jb, test_b)
-
-        Ja = DiffResults.JacobianResult(test_val, a)
-        Jb = DiffResults.JacobianResult(test_val, b)
-        ReverseDiff.jacobian!((Ja, Jb), ctp, (a, b))
-        test_approx(DiffResults.value(Ja), test_val)
-        test_approx(DiffResults.value(Jb), test_val)
-        test_approx(DiffResults.gradient(Ja), test_a)
-        test_approx(DiffResults.gradient(Jb), test_b)
+    @testset "Test binary Jacobian: f=$f - a::$(typeof(a)), b::$(typeof(b))" begin
+        test_val = f(a, b)
+        test_a = ForwardDiff.jacobian(x -> f(x, b), a)
+        test_b = ForwardDiff.jacobian(x -> f(a, x), b)
+
+        # without JacobianConfig
+
+        @testset "Without JacobianConfig" begin
+            Ja, Jb = ReverseDiff.jacobian(f, (a, b))
+            test_approx(Ja, test_a)
+            test_approx(Jb, test_b)
+
+            Ja = similar(a, length(test_val), length(a))
+            Jb = similar(b, length(test_val), length(b))
+            ReverseDiff.jacobian!((Ja, Jb), f, (a, b))
+            test_approx(Ja, test_a)
+            test_approx(Jb, test_b)
+
+            Ja = DiffResults.JacobianResult(test_val, a)
+            Jb = DiffResults.JacobianResult(test_val, b)
+            (Ja, Jb) = ReverseDiff.jacobian!((Ja, Jb), f, (a, b))
+            test_approx(DiffResults.value(Ja), test_val)
+            test_approx(DiffResults.value(Jb), test_val)
+            test_approx(DiffResults.gradient(Ja), test_a)
+            test_approx(DiffResults.gradient(Jb), test_b)
+        end
+
+        # with JacobianConfig
+
+        @testset "With JacobianConfig" begin
+            cfg = ReverseDiff.JacobianConfig((a, b))
+
+            Ja, Jb = ReverseDiff.jacobian(f, (a, b), cfg)
+            test_approx(Ja, test_a)
+            test_approx(Jb, test_b)
+
+            Ja = similar(a, length(test_val), length(a))
+            Jb = similar(b, length(test_val), length(b))
+            ReverseDiff.jacobian!((Ja, Jb), f, (a, b), cfg)
+            test_approx(Ja, test_a)
+            test_approx(Jb, test_b)
+
+            Ja = DiffResults.JacobianResult(test_val, a)
+            Jb = DiffResults.JacobianResult(test_val, b)
+            (Ja, Jb) = ReverseDiff.jacobian!((Ja, Jb), f, (a, b), cfg)
+            test_approx(DiffResults.value(Ja), test_val)
+            test_approx(DiffResults.value(Jb), test_val)
+            test_approx(DiffResults.jacobian(Ja), test_a)
+            test_approx(DiffResults.jacobian(Jb), test_b)
+        end
+
+        # with JacobianTape
+
+        @testset "With JacobianTape" begin
+            tp = ReverseDiff.JacobianTape(f, (rand(eltype(a), size(a)), rand(eltype(b), size(b))))
+
+            Ja, Jb = ReverseDiff.jacobian!(tp, (a, b))
+            test_approx(Ja, test_a)
+            test_approx(Jb, test_b)
+
+            Ja = similar(a, length(test_val), length(a))
+            Jb = similar(b, length(test_val), length(b))
+            ReverseDiff.jacobian!((Ja, Jb), tp, (a, b))
+            test_approx(Ja, test_a)
+            test_approx(Jb, test_b)
+
+            Ja = DiffResults.JacobianResult(test_val, a)
+            Jb = DiffResults.JacobianResult(test_val, b)
+            (Ja, Jb) = ReverseDiff.jacobian!((Ja, Jb), tp, (a, b))
+            test_approx(DiffResults.value(Ja), test_val)
+            test_approx(DiffResults.value(Jb), test_val)
+            test_approx(DiffResults.gradient(Ja), test_a)
+            test_approx(DiffResults.gradient(Jb), test_b)
+        end
+
+        # with compiled JacobianTape
+
+        @testset "With compiled JacobianTape" begin
+            if length(tp.tape) <= COMPILED_TAPE_LIMIT # otherwise compile time can be crazy
+                ctp = ReverseDiff.compile(tp)
+
+                Ja, Jb = ReverseDiff.jacobian!(ctp, (a, b))
+                test_approx(Ja, test_a)
+                test_approx(Jb, test_b)
+
+                Ja = similar(a, length(test_val), length(a))
+                Jb = similar(b, length(test_val), length(b))
+                ReverseDiff.jacobian!((Ja, Jb), ctp, (a, b))
+                test_approx(Ja, test_a)
+                test_approx(Jb, test_b)
+
+                Ja = DiffResults.JacobianResult(test_val, a)
+                Jb = DiffResults.JacobianResult(test_val, b)
+                (Ja, Jb) = ReverseDiff.jacobian!((Ja, Jb), ctp, (a, b))
+                test_approx(DiffResults.value(Ja), test_val)
+                test_approx(DiffResults.value(Jb), test_val)
+                test_approx(DiffResults.gradient(Ja), test_a)
+                test_approx(DiffResults.gradient(Jb), test_b)
+            end
+        end
     end
 end
 
diff --git a/test/compat/CompatTests.jl b/test/compat/CompatTests.jl
index 9dc08ed8..f509fcdd 100644
--- a/test/compat/CompatTests.jl
+++ b/test/compat/CompatTests.jl
@@ -1,13 +1,33 @@
 module CompatTests
 
-using FillArrays, ReverseDiff, Test
+using DiffResults, FillArrays, StaticArrays, ReverseDiff, Test
 
-@test ReverseDiff.gradient(fill(2.0, 3)) do x
-  sum(abs2.(x .- Zeros(3)))
-end == fill(4.0, 3)
+@testset "FillArrays" begin
+    @test ReverseDiff.gradient(fill(2.0, 3)) do x
+        sum(abs2.(x .- Zeros(3)))
+    end == fill(4.0, 3)
 
-@test ReverseDiff.gradient(fill(2.0, 3)) do x
-  sum(abs2.(x .- (1:3)))
-end == [2, 0, -2]
+    @test ReverseDiff.gradient(fill(2.0, 3)) do x
+        sum(abs2.(x .- (1:3)))
+    end == [2, 0, -2]
+end
 
-end
\ No newline at end of file
+sumabs2(x) = sum(abs2, x)
+
+@testset "StaticArrays" begin
+    @testset "Gradient" begin
+        x = MVector{2}(3.0, 4.0)
+        result = ReverseDiff.gradient!(DiffResults.GradientResult(x), sumabs2, x)
+        @test_broken x == [3.0, 4.0]
+        @test_broken DiffResults.value(result) == 25.0
+    end
+
+    @testset "Hessian" begin
+        x = MVector{2}(3.0, 4.0)
+        result = ReverseDiff.hessian!(DiffResults.HessianResult(x), sumabs2, x)
+        @test_broken x == [3.0, 4.0]
+        @test_broken DiffResults.value(result) == 25.0
+    end
+end
+
+end